简介：这是一套基于改进的豪斯多夫距离与DBSCAN的船舶航迹聚类Matlab实现，复现了《基于轨迹聚类的船舶异常行为识别研究》论文核心内容，适合学习聚类算法、从事船舶航迹分析或异常行为识别的研究者参考。项目共20个文件、压缩包约4.32MB，主体为14个可直接运行的.m脚本，另含.mat航迹数据、说明文档md、结果图片与子压缩包，覆盖航迹提取、DBSCAN聚类、H距离计算、阈值寻优及航迹偏离预测等完整流程。目前已有766人学习。通过该案例，既能掌握豪斯多夫距离与DBSCAN结合的改进聚类思路，也能理解航迹分类、聚类中心提取和预测阈值寻优的实际操作方法。代码模块划分清晰，替换数据或相关模块后，还可拓展为自定义船舶航迹模型，适合作为论文复现和课程设计的实用参考。
1. 船舶航迹聚类：真正决定成败的不是DBSCAN，而是提前一步的距离定义
DBSCAN作为密度聚类算法的原理几句话就能讲完，可真到船舶航迹聚类这个任务上，最先翻车的往往不是Eps和MinPts，而是最不起眼的“轨迹间距离”。AIS船舶轨迹是典型的不等长序列：有的船十分钟回传一批密集点，有的半小时才补一条稀疏航迹；定位信号还带着几十米到上百米的漂移。把这些数据直接喂进DBSCAN，欧氏距离根本没法定长对齐，聚类结果自然糊成一团。
问题收敛到一个被许多实现验证过的答案上：用Hausdorff距离衡量两条轨迹的几何接近程度，再把距离矩阵交给DBSCAN。它不要求轨迹等长，也不要求时间对齐，只问“两条轨迹占据的航迹形状有多像”。搞懂这一步，才算真正理解航迹聚类为什么能work——聚类算法只是最后装配的外壳，真正的机制藏在距离函数里。
2. 轨迹相似度的底层逻辑：为什么逐点对齐在AIS数据上行不通
2.1 不等长轨迹进入距离函数之前的第一个问题
一条AIS轨迹在程序里长什么样？存在Matlab的结构体或cell数组里，每个元素是一个 m×2 的坐标矩阵，m 因航次长短、上报策略而完全不一致。有的轨迹只有二十个点，是一条从锚地到泊位的短距机动；有的轨迹有上千个点，是一条完整的跨区航行记录。它们之间没有天然的“第i个点对应第i个点”关系，因为采样时刻、航速、转向点位置全部是错开的。
很多第一次做航迹聚类的人会尝试把轨迹规整成定长向量：取最长的轨迹做基准，短的补零或复制末尾点；或者统一重采样到100个点再做逐点欧氏距离。这两条路都有硬伤。补零等于人为制造一段“静止在某个坐标上”的假轨迹，会在聚类时把有靠泊行为的航迹和没有靠泊行为的航迹强行拉到一起；重采样到等长虽然解决了长度问题，但等长之后的逐点比较仍然要求“两条船在同一语义位置上对齐”——实际情况是轨迹A第20个点是转向点，轨迹B第20个点可能还在直线航行段，逐点相减算出来的距离包含了大量语义错位。
所以在航迹聚类里，距离函数的第一条准则不是“怎么让长度相等”，而是“怎么在两个不需要对齐的集合之间定义远近”。这个出发点直接把问题推向集合级距离度量，而不是序列级距离度量。
2.2 逐点对齐与DTW的失效场景
DTW（动态时间规整）常常被当作变长序列聚类的万能解，但在船舶AIS数据上它有两个具体问题。第一个问题是计算量。DTW要构造两条轨迹的点对点累积代价矩阵，复杂度 O(mn)，五千条轨迹两两计算就是两千五百万次DTW，即使轨迹平均只有一百个点，这个开销在普通机器上也已经到分钟甚至小时级。更大的问题是它的时间单调约束：DTW要求匹配路径沿时间轴单调推进，可AIS轨迹的时间戳并不干净。船舶在港内停泊时，AIS会反复上报几乎重叠的位置点；在卫星覆盖不好的海域，回传间隔可能拉长到几十分钟。这些时间异常会让DTW把“停留段”疯狂膨胀，把本应相近的两条轨迹匹配出诡异的大距离。
另一种逐点方案是离散Fréchet距离，它比DTW更贴合空间形状比较，但核心问题依旧：它本质上是“找一条同时穿过两条轨迹的绳子”，对离群点和采样密度极敏感，一条轨迹上多一个被噪声拉偏的GPS点，Fréchet距离就会被这个点锁死。相比之下，Hausdorff距离不做任何时间或次序假设，它只关心几何上的最大近邻偏差，计算上也能完全矩阵化（pdist2 一步出结果），在工程实现上是更务实的选择。
2.3 Hausdorff距离的几何定义：双向最大最小距离
Hausdorff距离定义在两个点集之间。设有轨迹 A={a1,...,am}，轨迹 B={b1,...,bn}，先定义单向距离：
h(A,B) = max_{a∈A} min_{b∈B} ||a-b||
含义是：对A里的每一个点，找它在B里的最近邻居，取这些最近距离里的最大值。直觉化地说，就是“A当中离B最远的那一个点，离B有多远”。这个值刻画的是A被B覆盖得有多充分。反过来 h(B,A) 刻画B被A覆盖的程度，两者通常不相等，所以完整的Hausdorff距离取两者的最大值：
H(A,B) = max( h(A,B), h(B,A) )
在航迹语境下，H(A,B)小，说明两条轨迹虽然点数和稀疏程度不同，但它们的几何形状在空间上互相覆盖得很好；H(A,B)大，则说明至少有一条轨迹有相当一部分区域落在另一条的覆盖范围之外。
用这个度量，轨迹A是一条从某锚地A到B泊位的直线，轨迹B是一条绕了远路的折线，它们的Hausdorff距离会由折线偏离直线的最大幅度决定；而两条平行但一前一后航行的轨迹，只要方向一致、间距不大，Hausdorff距离就会很小。这正好符合船舶航迹聚类的需求：聚在一起的是走同一走廊、空间形态接近的航迹，而不是时间上恰好同步的航迹。
3. DBSCAN的密度机制：从距离矩阵到航迹簇
3.1 Eps邻域与MinPts：密度可达是怎么传递的
拿到轨迹两两之间的Hausdorff距离矩阵后，DBSCAN的工作就变成了纯粹的密度问题。算法只有两个参数：邻域半径Eps和最小点数MinPts。对任意一条轨迹，统计在Eps范围内存在多少条其他轨迹；如果数量达到MinPts，这条轨迹就是核心轨迹。
核心轨迹的密度可达关系沿着链传递：轨迹X的Eps邻域里有核心轨迹Y，Y的邻域里又有核心轨迹Z，那么即使Z和X之间的距离超过了Eps，X和Z也会通过Y这条链被划进同一个簇。这是DBSCAN区别于K-means的底层机制——簇的形状不是凸的，可以沿着航道自然弯折。航道聚类场景里这非常关键：一条近海航道可能有一个弧形的拐弯，航迹在拐弯前后首尾距离很远，但中间被大量连续航迹衔接，只有基于传递性的密度聚类才能把这段弧形完整地圈出来。
Eps的取值单位在这里直接跟着距离函数的单位走。Hausdorff距离由经纬度坐标算出来时是海里或公里，Eps设为0.5就意味着只把空间上互相覆盖在0.5海里范围内的轨迹当邻居。可以先跑一个直方图看轨迹间距离分布的拐点，再定Eps，避免纯玄学调参。
3.2 核心轨迹与边界轨迹：离群航迹如何被分离
在一个簇里，轨迹分两类角色。核心轨迹是内部支柱，特点是邻居多；边界轨迹刚好达标或略低于MinPts，但它落在某个核心轨迹的Eps邻域内，因此被吸收进簇。真正有意思的是噪声点：既不是核心，也不在任何核心的邻域内。
放到船舶航迹数据里，噪声轨迹通常对应三类异常行为：偏离主航道的渔船作业轨迹、部分只上报了零星几个点的残缺轨迹，以及因为AIS信号丢失而出现大段跳变的轨迹。这些轨迹的Hausdorff距离到任何簇都很远，DBSCAN给它们贴上噪声标签而不是强行归入最近的簇，这种“宁缺毋滥”的行为是无监督聚类里难得的硬边界。MinPts设得越大，孤立的轨迹对越容易被判定为噪声；设得太小，两三条偶然接近的轨迹就会串成一个没有实际意义的微簇，所以MinPts一般至少取5到10，用来过滤掉“偶遇型”轨迹对。
3.3 MATLAB里用预计算距离矩阵跑DBSCAN的实现结构
使用统计与机器学习工具箱时，代码非常短：
% trajectories: 1×N cell 数组，每个元素是 m×2 的经纬度坐标矩阵 numTraj = numel(trajectories); D = zeros(numTraj); for i = 1:numTraj for j = i+1:numTraj D(i,j) = improvedHausdorffTraj(trajectories{i}, trajectories{j}, params); D(j,i) = D(i,j); end end % 预计算距离矩阵直接喂给 dbscan idx = dbscan(D, eps, minPts, 'Distance', 'precomputed'); % idx 为 N×1 标签向量，-1 表示噪声（离群航迹）核心点在于'Distance','precomputed'：dbscan 不再自行计算向量间的距离，而是直接查表。这样做的工程价值是把“轨迹距离计算”和“密度聚类”彻底解耦，距离函数怎么改都不影响聚类逻辑。
没有工具箱时，DBSCAN可以用图连通分量近似实现。先构造邻接矩阵 A = (D <= eps)，再找出所有邻居数不少于 MinPts 的核心轨迹；对每个核心轨迹做广度优先搜索，沿着核心到核心的边扩张簇边界。这个思路很多开源实现都在用，本质上是把密度可达翻译成图上的连通性，效果与标准DBSCAN基本一致。
4. 改进的Hausdorff距离：针对AIS噪声、时间错位与方向语义的三处修改
4.1 原始Hausdorff被单点噪声绑架：分位数截断
原始Hausdorff距离有一个众所周知的脾气：它由最大值决定。A里只要有一个点离B特别远，不管其他点匹配得多好，整个距离都被这一个点拉爆。AIS轨迹恰恰是离群点高发数据——船舶定位受多径效应和传感器漂移影响，一条正常航迹里偶尔混入一两个偏离航线几百米的野点，几乎无法避免。用原始Hausdorff计算时，这两条实际上属于同一簇的轨迹，距离可能比它们到另一条完全不同航线的轨迹还大。
常见的改进思路是把单向距离里的 max 换成高分位点。把 min(D, [], 2) 得到的是“每个A点到B的最近距离”序列，对这个序列取95%分位数，等效于扔掉最不匹配的那5%的A点，让距离由绝大多数点的覆盖质量决定。这其中有个隐含假设：AIS噪声在轨迹点中的占比远低于5%。实际数据里这个假设通常成立，所以分位数截断是性价比最高的鲁棒化手段。
4.2 时间轴错位导致几何失真：重采样到统一时间网格
Hausdorff距离虽然不要求轨迹长度一致，但原始轨迹的采样密度对距离结果有隐性影响。一条轨迹在弯道段密集拐点，另一条在直航段稀疏分布，单向距离里“每个点找最近邻”对密集段会赋予过高的权重——密集段的点数量多，等于变相提高了弯道在距离里的话语权。而且AIS的采样间隔不固定，不能直接拿原始点集比形状。
处理办法是先重采样到固定点数。把每条轨迹的累计弧长归一化到0到1之间，再用线性插值或样条插值取出N个等弧长位置上的坐标点（N一般取16到32）。等弧长重采样而不是等时间重采样，是为了保证轨迹的几何形态不变：船舶在弯道减速、直道加速，原始点密度与几何凸凹无关，按时间采样会把密集的慢速段重复计入，按弧长采样则让每个轨迹点均匀地代表一段走过的路。重采样之后的Hausdorff距离不再受轨迹点数量差异干扰，成为纯粹的形状距离。
4.3 对向航行的语义差异：方向感知惩罚项的引入
两条轨迹几何上几乎重合，但一条从某近海锚地开往港区，另一条反向从港区驶出，它们的Hausdorff距离可能非常小。同一个聚簇算法如果只看空间形状，会把对向航行的两条轨迹划进同一组。这在很多场景下是错误语义：港区调度里，进港航迹和出港航迹对应完全不同的行为模式，混在一起会直接污染后续的航迹预测。
改进方案是给距离加上方向惩罚项。取每条轨迹首点到末点的位移向量，计算夹角余弦cosθ；同向时cosθ≈1，惩罚趋近于0；垂直时惩罚为w；对向时cosθ≈-1，惩罚为2w。修正后的距离写成：
d = d_geo + w * (1 - cosθ)
这里还有一个可配置细节：如果任务关心的是“这条走廊而不是那条走廊”，对向航行应该共享同一走廊，就需要把 cosθ 取绝对值，让同向和对向的惩罚都为0，只惩罚方向差异较大的轨迹。两种做法分别对应“区分航向”和“区分航路走廊”两种业务语义，需要作为参数暴露出来。
4.4 改进距离的MATLAB实现与参数实测
把三处修改合起来，完整实现如下：
function d = improvedHausdorffTraj(trajA, trajB, params) % trajA, trajB: m×2 和 n×2 的 [x, y] 轨迹点 % params.resampleN : 重采样点数，一般 16~32 % params.q : 单向距离的分位数，一般 0.90~0.98 % params.w : 方向惩罚权重，取典型轨迹间距的 0.1~0.3 % params.mode : 'directed' 区分对向航行 / 'symmetric' 合并对向 % 1. 等弧长重采样 trajA = resampleByArcLength(trajA, params.resampleN); trajB = resampleByArcLength(trajB, params.resampleN); % 2. 点间欧氏距离矩阵 D = pdist2(trajA, trajB); % resampleN×resampleN % 3. 分位数截断的单向 Hausdorff hAB = quantile(min(D, [], 2), params.q); % 每个A点找最近B点，再取高分位 hBA = quantile(min(D, [], 1), params.q); % 每个B点找最近A点，再取高分位 dGeo = max(hAB, hBA); % 双向取最大，保持对称 % 4. 方向惩罚项 vA = trajA(end, :) - trajA(1, :); vB = trajB(end, :) - trajB(1, :); cosTheta = dot(vA, vB) / (norm(vA) * norm(vB) + eps); if strcmp(params.mode, 'symmetric') cosTheta = abs(cosTheta); end d = dGeo + params.w * (1 - cosTheta); endresampleByArcLength 的内部逻辑是对累计弧长做插值：
function traj = resampleByArcLength(traj, N) segLen = sqrt(sum(diff(traj).^2, 2)); cumLen = [0; cumsum(segLen)]; total = cumLen(end); s = linspace(0, total, N)'; traj = [interp1(cumLen, traj(:,1), s), ... interp1(cumLen, traj(:,2), s)]; end关键的参数含义：q=0.95 意味着允许5%的轨迹点是孤立野点，对这5%完全视而不见；q=0.90 更激进，保留的形状特征更少，适合噪声严重的AIS数据；q=1.0 则是原始Hausdorff。w 的取值要看 dGeo 的典型量级，如果两条相近航迹的 dGeo 在0.1海里量级，而反向航迹的 cosθ 是-1，那么 w 取0.05就已经足够把反向距离拉开一倍，不需要设太大，否则会把方向差异放大到覆盖几何差异。
5. 船舶航迹聚类的常见问题与避坑：从参数到数据清洗的排查记录
5.1 Eps选不好：全簇合一或全盘噪声
现象：聚类结果要么所有轨迹聚成一个巨大的簇，要么几乎全是噪声点，中间过渡状态完全不存在。
原因：Hausdorff距离分布往往呈长尾形态，大量轨迹间距离集中在一个很窄的区间，Eps在某个临界值前后，邻居数量剧烈跳变。Eps取得比临界值大一点点，密度图就连成片；小一点点，核心轨迹全部不够MinPts，所有轨迹都被判为噪声。
解决：先算完距离矩阵后，把每条的 k 近邻距离（k 取 MinPts-1）排序画曲线，找曲线由缓变陡的拐点作为Eps基准，再在这个值附近做小范围网格搜索，用轮廓系数圈定最优区间。不要凭对“海里”的直觉直接填一个值，因为Hausdorff距离是“最大覆盖偏差”，和普通点间距不是一个量纲。
5.2 单个GPS漂移点让整条轨迹失去归属
现象：肉眼明显属于同一条航道的十几条轨迹，聚类后被劈成两簇，仔细检查发现簇边界上恰好有一条轨迹带了一个离航线几百米的野点。
原因：q 设成了1.0，退化成了原始Hausdorff。最大值对离群点完全没有抵抗力，一个野点就把这条轨迹推向所有其他轨迹的覆盖范围之外。
解决：把 q 降到0.92到0.95。注意 q 不能低于0.85，否则两条航向垂直但几何上偶有接近的轨迹也会被算成近邻，簇的区分度下降。折中办法是同时对原始轨迹做一次中值滤波，把位置序列里偏离局部邻域超过阈值（比如50米）的坐标点直接剔除。
5.3 重采样方式不当导致轨迹形态畸变
现象：重采样后的轨迹在图上明显“缩水”，原本船在锚地附近绕的小圈被抹成一个点；或者弯道被拉直，聚类结果与航图判读的航路分布不符。
原因：重采样按时间轴插值，而AIS时间戳在慢速航行段非常密集、在快速段稀疏。时间等间隔插值会把慢速段的点按时间均匀摊开，绕泊位的圆周运动被压缩成一小段原地抖动的点，几何形状失真。
解决：改成按累计弧长归一化做插值（4.4节的resampleByArcLength），确保每个插值点对应等长的航迹段。弧长插值也有边界情况：轨迹里有静默或长时间停泊时，连续重复的GPS点会产生零长度弧段，需要在插值前先剔除重复点。
5.4 航向角跨环绕边界引发的角度差错误
现象：方向惩罚项计算出异常大的角度差，方向完全相同的两条轨迹被误判为“几乎垂直”。检查数据发现一条轨迹的航向角是358度，另一条是2度，直接相减得到356度，余弦值接近1，但逻辑不对。
原因：如果方向项用航向角序列或者首尾方位角相减，角度的周期性问题没有处理。直接在角度值上做减法，358和2的实际差是4度，但数值差356度，惩罚被错误放大。
解决：不要用角度相减，始终用单位向量的点积算余弦——点积自带周期性，358度和2度的向量几乎完全同向，cosθ≈1。如果确实需要角度值，就用wrapToPi或mod(angDiff+180, 360)-180把差值折叠到[-180, 180]区间，再做后续计算。
5.5 距离矩阵里的NaN、零点和不对称
现象：距离矩阵对角线是0，但非对角线出现NaN；或者 D(i,j) 和 D(j,i) 数值不一致；部分轨迹和所有轨迹的距离都是0。
原因：三处来源。第一，经纬度数据里有NAN或Inf坐标，pdist2算出来的距离带NaN，quantile 对 NaN 敏感，一传十十传百；第二，AIS缺失值被填成(0,0)，这在经纬度坐标里恰好落在某处海面上，会让轨迹与无数轨迹“距离接近”；第三，单向Hausdorff的hAB和hBA本身天然不对称，如果最后只做了单向距离就放进矩阵，聚类结果会依赖遍历顺序。
解决：清洗阶段把坐标不在合理经纬度范围（如经纬度绝对值超过阈值或有重复固定值）的点整段剔除；填充的(0,0)直接丢弃。不对称问题在计算时强制合并：H = max(hAB, hBA)，同时在下三角复制时保证 D(j,i)=D(i,j)，构建一个对称且非负的矩阵再交给dbscan。
6. 改进距离的设计取舍：鲁棒性、粒度与可解释性的平衡点
分位数q、方向权重w、重采样点数N，这三个参数构成了整个改进距离的取舍面。q决定你有多信任AIS数据：q越接近1，距离越忠实于轨迹细节，但也被噪声绑架的风险越高；q越小，轨迹的尖角和转向特征被抹掉，两簇容易过度合并。N决定几何描述的精细度：N=16时轨迹只是粗略的折线，两条贴近但在某处有小幅度绕行的轨迹会合并；N=64则开始区分这些小细节，但也开始把定位噪声当成真实几何。没有绝对正确的组合，只有匹配数据质量的选择。
一个务实的验证顺序是先固定MinPts=5，用k近邻拐点找到Eps范围，然后在{0.90, 0.95, 0.98}×{0.0, 0.05, 0.1}的网格上用轮廓系数选优。轮廓系数在航迹聚类里尤其好用，因为Hausdorff距离本身就是度量，簇内距离和簇间距离的比值可以直接解释。我一般在选定参数后还会抽查簇中心的轨迹：把簇内所有轨迹重采样后取逐点中位数，画出一条“代表航迹”，直接和电子海图对照语义。这一步能暴露纯数字指标看不见的问题——比如某个簇其实是两条反向航道被过分鲁棒的q值强行焊在一起。
聚类结果进入船舶航迹预测的下游任务时，我发现方向感知的距离对预测质量的提升比预想中大。簇中心航迹本身就能作为航迹预测的先验模板，新来一条轨迹先按改进Hausdorff距离归属到最近的簇，再用该簇历史航迹推演后续走向。从那以后，我每次做航迹聚类都强制走一遍完整流程：先清洗坐标和缺失值，再等弧长重采样，然后分位数鲁棒化，最后检查方向惩罚的对向语义。这套习惯帮我少踩了很多坑，也希望帮到你。
本文还有配套的精品资源，点击获取
