MATLAB改进Hausdorff距离DBSCAN船舶航迹聚类实战

简介：本资源面向船舶航迹研究与Matlab聚类算法学习者，复现论文《基于轨迹聚类的船舶异常行为识别研究》的核心流程，解决航迹数据提取、聚类与偏离预测的建模问题。包内共20个文件，以14个m脚本为主体，辅以mat数据、zip子包及png图示，整体约4.32MB，涵盖DBSCAN聚类、改进Hausdorff距离计算、航迹提取、阈值分类与寻优等模块，可直接运行。已有766人学习下载。读者可借此掌握DBSCAN距离度量与H距离计算方法，理解航迹分类的完整流程，并基于聚类中心开展航迹偏离预测，经检测准确率较高。船舶航迹方向的研究者可替换数据与相关模块，将其拓展为自有模型，适合作为聚类算法的应用案例学习。

1. 从一次航迹交叉误判说起：这套 MATLAB 聚类方案到底在解决什么

去年帮一个做港口监控的朋友看数据，他手里有几百条 AIS 航迹，想按航行模式自动分组。他先用的是最朴素的欧氏距离加 DBSCAN，结果两条走向完全一致、只是横向错开几百米的航道，被硬生生拆成了两类；而两条在交叉口短暂重叠、实际航向差 90 度的船，反倒被塞进同一簇。问题不在 DBSCAN，在距离度量。欧氏距离只看点位远近，不看方向，航迹这种带时序和走向的数据，用它就是拿错尺子量东西。

这套「MATLAB 基于改进 Hausdorff 距离的 DBSCAN 船舶航迹聚类」资源，核心就干一件事：把航迹之间的相似度度量从「点对点距离」换成「考虑方向与形状的改进 Hausdorff 距离」，再喂给 DBSCAN 做密度聚类。它适合三类人：手里有 AIS 或雷达航迹数据、想做航行模式挖掘的工程人员；用 MATLAB 做交通流或海事方向课题、需要一套能跑通并改得动的研究者；以及被欧氏距离坑过、想搞清楚航迹相似度到底该怎么算的人。整套东西是 MATLAB 脚本加函数的形式，不依赖深度学习工具箱，普通配置的机器就能跑，改参数、换数据都直接。

2. 改进 Hausdorff 距离：为什么不能直接用经典版本

2.1 经典 Hausdorff 的致命短板

经典 Hausdorff 距离的定义是：集合 A 中每个点到集合 B 的最近距离，取最大值；反过来 B 到 A 也取最大值，两者再取大。写成公式就是 H(A,B) = max( sup_{a∈A} inf_{b∈B} d(a,b), sup_{b∈B} inf_{a∈A} d(a,b) )。它对「最不像的那一对点」极其敏感——只要有一条航迹里有一个离群点，整个距离就被这个点绑架。航迹数据里 GPS 漂移、个别跳点太常见了，直接套经典 Hausdorff，聚类结果会被几个噪声点带偏。

更麻烦的是，经典 Hausdorff 只关心「点集之间的几何接近程度」，完全不看航迹的走向。两条平行但反向的航迹，点集可能挨得很近，Hausdorff 距离很小，可它们在实际语义上是完全相反的航行模式。这就是我朋友那个案例里交叉口误判的根源。

2.2 改进思路：方向加权 + 分位数截断

这套资源里的改进主要落在两个地方，我拆开讲。

第一是方向加权。对每条航迹的每个点，除了坐标，还算一个局部航向角（用相邻点差分得到）。在计算点对距离时，不只用欧氏距离，而是把航向差也纳入：d_improved = d_euclid + λ · (1 - cos(Δθ))。λ 是方向权重系数，Δθ 是两点航向角之差。当两条航迹走向一致时 cos(Δθ)≈1，惩罚项趋近 0；走向相反时 cos(Δθ)≈-1，惩罚项拉满。这样反向航迹的距离被显著放大，不会再被误聚。

第二是分位数截断。把 sup（取最大）换成取距离集合的某个高分位数，比如 90% 分位。这样单个离群点最多影响那 10% 的尾部，不会主导整个距离值。分位数参数 q 一般取 0.85～0.95，q 越大越接近经典 Hausdorff、越敏感，q 越小越鲁棒但可能过度平滑。

2.3 距离矩阵的 MATLAB 实现

下面这段是距离计算的核心函数，我按资源里的结构重写并加了注释，可以直接放进你的工程。

function D = improvedHausdorff(trajA, trajB, lambda, q) % trajA, trajB: N×3 矩阵，列为 [x, y, heading(弧度)] % lambda: 方向权重系数，建议 0.3~0.8 % q: 分位数，建议 0.85~0.95 % D: 改进 Hausdorff 距离（标量） % 计算 A 中每个点到 B 中每个点的改进距离 dAB = pairDist(trajA, trajB, lambda); dBA = pairDist(trajB, trajA, lambda); % 每个点到对方集合的最近距离 hAB = min(dAB, [], 2); % N_A × 1 hBA = min(dBA, [], 2); % N_B × 1 % 分位数截断代替 max D = max(quantile(hAB, q), quantile(hBA, q)); end function d = pairDist(A, B, lambda) % A: M×3, B: K×3，返回 M×K 改进距离矩阵 M = size(A, 1); K = size(B, 1); d = zeros(M, K); for i = 1:M dx = A(i,1) - B(:,1); dy = A(i,2) - B(:,2); euclid = sqrt(dx.^2 + dy.^2); dtheta = A(i,3) - B(:,3); dirPenalty = 1 - cos(dtheta); % 方向惩罚项 d(i,:) = euclid + lambda * dirPenalty; end end

逻辑说明：pairDist算出两个航迹点集之间的全连接改进距离矩阵，欧氏距离负责位置接近度，1 - cos(Δθ)负责方向一致性惩罚。improvedHausdorff里对每个点取到对方集合的最近距离，再用quantile取分位数，最后两个方向取大。参数上，lambda调大方向权重更强，适合航道方向差异明显的场景；q调小抗噪更强，但太小（低于 0.8）会让距离区分度下降。我一般先用 lambda=0.5、q=0.9 跑一版看聚类轮廓，再微调。

提示：航向角要做归一化到 [-π, π]，否则 cos 计算会因 2π 跳变出错。资源里在预处理阶段用wrapToPi处理过，自己接数据时别忘了这一步。

3. DBSCAN 参数整定：eps 和 MinPts 到底怎么定

3.1 为什么 DBSCAN 配改进距离更稳

DBSCAN 靠密度可达性聚类，不需要预先指定簇数，还能把低密度区域标成噪声——这对航迹数据很合适，因为真实航行模式的数量事先不知道，而且总有零星异常航迹。但 DBSCAN 对 eps（邻域半径）极其敏感，而 eps 的合理取值直接取决于距离度量的尺度。换成改进 Hausdorff 后，距离分布和欧氏距离完全不同，不能沿用老参数，必须重新整定。

3.2 k-距离曲线定 eps

标准做法是画 k-距离曲线：对每个点，算它到第 k 近邻的距离，排序后画出来，拐点位置就是 eps 的参考值。k 一般取 MinPts。下面这段直接可用。

% D: N×N 距离矩阵（用 improvedHausdorff 两两算出来） N = size(D, 1); MinPts = 5; % 先设一个初值 k = MinPts; sortedKDist = zeros(N, 1); for i = 1:N di = D(i, :); di(i) = inf; % 排除自身 di = sort(di); sortedKDist(i) = di(k); % 第 k 近邻距离 end sortedKDist = sort(sortedKDist, 'descend'); plot(sortedKDist); xlabel('样本序号（按k-距离降序）'); ylabel(['第' num2str(k) '近邻距离']); grid on; % 找拐点：曲线从陡降转为平缓的位置

逻辑说明：di(i)=inf排除自身距离，排序后取第 k 个就是 k-距离。降序排列后，曲线拐点（曲率最大处）对应的纵坐标就是 eps 建议值。实际操作时我会把拐点值上下浮动 20% 各跑一次，看哪个簇数稳定、噪声比例合理（一般噪声占比 5%～15% 可接受）。MinPts 的经验规则是 MinPts ≥ 维度+1，航迹特征维度不高，取 4～8 都行，数据量大就往上调。

3.3 聚类主流程与参数落点

% 假设 trajs 是 cell 数组，每个元素是一条 N×3 航迹 n = numel(trajs); D = zeros(n, n); for i = 1:n for j = i+1:n D(i,j) = improvedHausdorff(trajs{i}, trajs{j}, 0.5, 0.9); D(j,i) = D(i,j); end end epsVal = 0.8; % 由 k-距离曲线拐点确定 MinPts = 5; labels = dbscan(D, epsVal, MinPts, 'Distance', 'precomputed'); % labels: 0 表示噪声，正整数为簇编号

逻辑说明：dbscan用'precomputed'模式直接吃距离矩阵，避免它内部再算一遍欧氏距离。epsVal和MinPts就是上面整定出来的值。跑完看labels里 0 的比例和簇的数量，如果噪声超过 20%，说明 eps 偏小；如果几乎所有点挤成一簇，eps 偏大。这一步没有银弹，就是反复试，但有了 k-距离曲线打底，通常三五次就能收敛。

注意：距离矩阵是对称的，循环里只算上三角再镜像赋值，能省一半计算量。航迹条数上百以后，这个优化很实在。

4. 避坑与排查：我踩过的五个真实问题

4.1 现象：聚类结果每次跑都不一样

原因：DBSCAN 本身是确定性的，但如果你的航迹点顺序在预处理时被打乱，或者用了随机采样做加速，距离矩阵就会变。另一个隐蔽原因是航向角计算依赖点顺序，顺序一变航向就变。

解决：预处理阶段固定航迹点顺序，按时间戳排序后再算航向。如果做了降采样，用等间隔抽取而不是随机抽取。跑之前把随机种子固定rng(42)，虽然 DBSCAN 不用随机数，但上游如果有任何随机环节都能锁住。

4.2 现象：两条明显不同的航迹距离却很小

原因：方向权重 lambda 设得太小，或者航向角没做归一化，导致 cos 惩罚项失效。还有一种情况是两条航迹长度差异极大，短航迹的点全被长航迹的某个局部区域「吸」住了。

解决：先把 lambda 提到 0.5 以上试；检查航向角是否用wrapToPi归一化；对长度差异大的航迹，考虑先做重采样统一到相近点数，或者在距离里加长度惩罚项。资源里对长度差异有处理，但如果你换数据，这点要自己盯。

4.3 现象：eps 稍微一动，簇数从 3 跳到 15

原因：距离分布本身没有明显的密度分层，说明改进距离的区分度不够，或者数据里噪声太多把密度结构糊掉了。

解决：先回去调距离参数，把 lambda 和 q 组合扫一遍，看距离分布的直方图有没有双峰或明显间隙。如果数据噪声确实大，先做一道航迹清洗（去掉点数过少、速度异常的航迹）再聚类。别硬调 eps，那是治标。

4.4 现象：计算距离矩阵慢到无法接受

原因：双重循环加逐点对计算，复杂度是 O(n² · L²)，n 是航迹数，L 是平均点数。n 上百、L 上千时，MATLAB 纯循环扛不住。

解决：把pairDist里的内层循环向量化（上面代码已经用矩阵运算替代了逐点循环）；对航迹做等间隔降采样，把 L 压到 100～200；如果 n 超过 500，考虑先用欧氏距离粗筛，只对可能邻近的航迹算改进 Hausdorff。我一般会在距离计算前加一个进度打印，不然跑着跑着不知道卡在哪。

4.5 现象：噪声点里混着真实的小众航行模式

原因：DBSCAN 把低密度区域一律标噪声，但有些真实模式样本就是少，密度天然低。

解决：调小 MinPts 让密度门槛降低，或者对标记为噪声的点做二次聚类（用更小的 eps 再跑一遍 DBSCAN）。另一个思路是聚类后对噪声点做人工抽检，确认是不是真的异常。资源里默认参数偏向保守，小众模式容易被吞，这点在做实际业务时要留意。

5. 进阶技巧：用轮廓系数和可视化把结果验到位

聚类跑完不算完，得验证。我习惯用两个手段：轮廓系数和航迹叠加图。

轮廓系数衡量每个点跟自己簇的紧密度 vs 跟最近邻簇的分离度，取值 [-1, 1]，越接近 1 越好。但注意，轮廓系数默认用欧氏距离算，对改进 Hausdorff 聚类的评估只能做参考，不能当唯一标准。我的做法是用它做相对比较：同一份数据，不同参数组合下谁的轮廓系数高，谁就更可能合理。

% labels 来自 dbscan，D 是改进 Hausdorff 距离矩阵 sil = silhouette((1:numel(labels))', labels, 'Distance', 'precomputed'); % 注意：silhouette 的 precomputed 模式需要传入方阵距离 meanSil = mean(sil(~isnan(sil))); fprintf('平均轮廓系数: %.3f\n', meanSil);

逻辑说明：silhouette在'precomputed'模式下直接吃距离矩阵，但要求距离矩阵是方阵且对角线为 0。如果你的 D 对角线不是 0，先D(1:n+1:end) = 0。meanSil排除 NaN（噪声点标签为 0 时可能产生 NaN）。这个值我一般拿来横向比参数，不追求绝对值多高。

可视化更直观。把每条航迹按簇编号上色，画在同一张经纬度图上，噪声点用灰色。一眼就能看出聚类是不是把走向一致的航迹归到了一起、交叉口有没有误判。

figure; hold on; colors = lines(max(labels)); for i = 1:numel(trajs) if labels(i) == 0 plot(trajs{i}(:,1), trajs{i}(:,2), 'Color', [0.7 0.7 0.7]); else plot(trajs{i}(:,1), trajs{i}(:,2), 'Color', colors(labels(i),:)); end end xlabel('经度'); ylabel('纬度'); title('航迹聚类结果'); grid on;

参数上，lines(n)生成 n 种区分度较好的颜色，噪声统一灰色。如果簇数超过 7，颜色会开始接近，可以换turbo或parula色图。画完重点看三处：簇内航迹走向是否一致、簇间是否有明显方向差异、噪声点是不是真的零散异常。

一个我常用的调参习惯：固定 MinPts=5，把 eps 从 k-距离拐点的 0.8 倍扫到 1.2 倍，步长 0.05，每个值记录簇数、噪声比例、平均轮廓系数，列成表挑最优。这套流程走下来，基本不会出现「跑出来不知道对不对」的情况。从那以后我每次拿到新航迹数据，都强制先跑一遍 k-距离曲线和参数扫描，再动手解释结果——跳过这步的代价，就是被错误聚类带偏结论。希望这套东西能帮到你。

本文还有配套的精品资源，点击获取