基于HDBSCAN的3D航空点云聚类问题：层状聚类优化诉求

问题背景 #

我有15000个欧洲上空飞机的点数据，包含纬度、经度、海拔信息，目前在做法国区域内的点聚类与可视化，目标是得到类球形聚类，但现在的HDBSCAN聚类结果完全按海拔分层——水平距离极远但海拔接近的点被归为一类，呈现层状分布。调整HDBSCAN参数后仍无改善，求可行的参数方案或替代算法，让聚类同时兼顾垂直（海拔）和水平（经纬度）维度。

现有代码 #

import easygui
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from matplotlib.path import Path
from matplotlib.patches import PathPatch
from shapely.geometry import Point, Polygon
import matplotlib.path as mplPath
import plotly.graph_objs as go
from plotly.offline import plot
from plotly.subplots import make_subplots
from sklearn.cluster import DBSCAN
from sklearn.cluster import KMeans
import hdbscan
from sklearn.mixture import GaussianMixture

# 读取Excel数据
df = pd.read_excel("1finalnidata.xlsx", sheet_name='50vzorek', usecols=["time", "icao24", "lat", "lon", "geoaltitude"], na_values=["NULL"])

# 提取数据列
k2 = df["time"].values
m2 = df["icao24"].values
n2 = df["lat"].values
b2 = df["lon"].values
z2 = df["geoaltitude"].values

# 读取空域数据
df = pd.read_excel("FinalniProstoryPYT.xlsx", sheet_name='Sheet1', usecols=["aid", "lat", "lon", "seznam"])

id = df["aid"].values
lat = df["lat"].values
lon = df["lon"].values
sez = df["seznam"].values
sez = sez[:38]

# 初始化空域坐标列表
airspaceLON = []
airspaceLAT = []

check = None
poly = []
planedatID = []
planedatAS = []
soux = []
souy = []
souz = []
print("3")

# 筛选空域内的飞机点
for i in range(0, len(sez)):
    check = sez[i]
    for j in range(0, len(id)):
        if id[j] == check:
            airspaceLON.append(lon[j])
            airspaceLAT.append(lat[j])

        if j == 8861:
            poly = list(zip(airspaceLAT, airspaceLON))
            poly_path = mplPath.Path(np.array(poly))

            for x in range(0, len(n2)):
                point = (n2[x], b2[x])
                if poly_path.contains_point(point):
                    planedatID.append(m2[x])
                    planedatAS.append(check)
                    soux.append(n2[x])
                    souy.append(b2[x])
                    souz.append(z2[x])

                poly.clear()
                airspaceLON.clear()
                airspaceLAT.clear()

matice5 = np.column_stack((planedatID, planedatAS))
souPuvBod = list(zip(souy, soux))

# 提取法国区域的多边形
polyg1 = []
polyg2 = []

for j in range(0, len(id)):
    if id[j] == "FRA":
        polyg1.append(lon[j])
        polyg2.append(lat[j])

poly1 = list(zip(polyg1, polyg2))

polygon = Polygon(np.array(poly1))
poly_path1 = mplPath.Path(np.array(poly1))

# 筛选法国区域内的飞机点
bod1=[]
bod2=[]
bod3=[]

for x in range(0, len(souPuvBod)):
    point = (souy[x], soux[x])
    if poly_path1.contains_point(point):
        bod1.append(souy[x])
        bod2.append(soux[x])
        bod3.append(souz[x])

data = np.column_stack((bod1, bod2, bod3))

# 运行HDBSCAN聚类
clusterer = hdbscan.HDBSCAN(min_cluster_size=40, min_samples=5, metric='euclidean')
labels = clusterer.fit_predict(np.column_stack((bod1, bod2, bod3)))

# 统计聚类数量
num_clusters = len(set(labels)) - (1 if -1 in labels else 0)

# 生成3D可视化
fig = make_subplots(rows=1, cols=1, specs=[[{'type':'scatter3d'}]])

# 生成聚类颜色
colors = ['rgb({},{},{})'.format(int(255/num_clusters*i), 0, int(255/num_clusters*(num_clusters-i-1))) for i in range(num_clusters)]

# 添加每个聚类的散点
for i in range(num_clusters):
    x = []
    y = []
    z = []
    for j in range(len(bod1)):
        if labels[j] == i:
            x.append(bod1[j])
            y.append(bod2[j])
            z.append(bod3[j])
    fig.add_trace(go.Scatter3d(x=x, y=y, z=z, mode='markers', marker=dict(color=colors[i], size=3)),row=1, col=1)

# 设置布局
fig.update_layout(title='Data for Specific Polygon in 3D', scene=dict(xaxis=dict(title='Longitude'), yaxis=dict(title='Latitude'), zaxis=dict(title='Altitude')))
plot(fig, filename='3d_scatter_with_polygon.html')  # 保存为HTML文件

核心问题分析 #

经纬度的数值范围（法国经度约2°W-8°E，纬度42°N-51°N）与海拔（0-12000米）的量级差异极大，欧几里得距离会被海拔的大数值主导，导致聚类仅关注海拔维度，忽略水平距离。

解决方案 #

1. 数据标准化/归一化 #

统一三个特征的尺度，让每个维度对距离计算的贡献权重一致，这是最基础且有效的修复手段：

from sklearn.preprocessing import StandardScaler

# 标准化数据：经度、纬度、海拔
scaler = StandardScaler()
scaled_data = scaler.fit_transform(np.column_stack((bod1, bod2, bod3)))

# 使用标准化后的数据重新运行HDBSCAN
clusterer = hdbscan.HDBSCAN(min_cluster_size=40, min_samples=5, metric='euclidean')
labels = clusterer.fit_predict(scaled_data)

2. 自定义地理空间距离度量 #

使用更贴合航空场景的三维距离计算方式：先通过Haversine公式计算地面水平距离，再结合海拔差得到真实三维距离，避免经纬度与海拔的尺度冲突：

import math

def haversine_3d(point1, point2):
    # point1: [lon, lat, alt], point2: [lon, lat, alt]
    lon1, lat1, alt1 = point1
    lon2, lat2, alt2 = point2
    
    # 计算地面Haversine距离（单位：米）
    R = 6371000  # 地球半径（米）
    phi1 = math.radians(lat1)
    phi2 = math.radians(lat2)
    delta_phi = math.radians(lat2 - lat1)
    delta_lambda = math.radians(lon2 - lon1)
    
    a = math.sin(delta_phi/2)**2 + math.cos(phi1)*math.cos(phi2)*math.sin(delta_lambda/2)**2
    c = 2*math.atan2(math.sqrt(a), math.sqrt(1-a))
    ground_dist = R * c
    
    # 计算海拔差
    alt_dist = abs(alt1 - alt2)
    
    # 返回三维欧几里得距离
    return math.sqrt(ground_dist**2 + alt_dist**2)

# 传递自定义距离给HDBSCAN，需设置algorithm='brute'
clusterer = hdbscan.HDBSCAN(min_cluster_size=40, min_samples=5, 
                            metric=haversine_3d,
                            algorithm='brute')
labels = clusterer.fit_predict(np.column_stack((bod1, bod2, bod3)))

3. 尝试替代聚类算法 #

如果HDBSCAN调整后仍达不到预期，可以尝试以下算法：

• DBSCAN：同样需要先标准化数据，重点调整eps（邻域半径）和min_samples参数：

from sklearn.cluster import DBSCAN

scaled_data = StandardScaler().fit_transform(np.column_stack((bod1, bod2, bod3)))
dbscan = DBSCAN(eps=0.5, min_samples=10)  # 根据数据分布调整eps值
labels = dbscan.fit_predict(scaled_data)

• 高斯混合模型（GMM）：基于概率的聚类算法，天然适合类球形分布的数据集，需要提前指定聚类数量：

from sklearn.mixture import GaussianMixture

scaled_data = StandardScaler().fit_transform(np.column_stack((bod1, bod2, bod3)))
gmm = GaussianMixture(n_components=10, random_state=42)  # n_components根据实际情况调整
labels = gmm.fit_predict(scaled_data)

额外调优建议 #

• 先可视化标准化后的三维数据，确认三个维度的尺度已统一；
• 对于HDBSCAN，可尝试缩小min_cluster_size（比如从40降至20），或设置cluster_selection_epsilon限制聚类的最大距离；
• 若使用自定义距离度量，可给海拔差添加权重系数（如alt_dist * 0.1），让水平距离的优先级更高，更符合航空场景的聚类逻辑。

内容的提问来源于stack exchange，提问作者Martin Kavka