简介这套Matlab项目复现了论文《基于轨迹聚类的船舶异常行为识别研究》的核心流程围绕改进的豪斯多夫距离与DBSCAN聚类实现航迹数据提取、聚类分析、聚类中心提取、基于豪斯多夫距离的航迹预测及预测阈值寻优等环节适合船舶航迹研究人员、航运异常行为分析者和DBSCAN算法学习者作为工程案例参考。压缩包共20个文件以14个可直接运行的m脚本为主体覆盖DBSCAN、豪斯多夫距离计算、航迹提取、阈值分类与寻优模块另有1个mat航迹数据、2个辅助zip压缩包、2个png效果图和1个md说明文档包体仅4.32MB。已有766人学习/下载。读者可在运行中理解DBSCAN距离度量和豪斯多夫距离的计算逻辑掌握航迹聚类整体流程替换自身船舶数据或相应模块即可拓展为自定义聚类与偏离预测模型代码完整、结构清晰便于二次开发。1. 为什么船舶航迹聚类要用“改进的 Harsdorf 距离 DBSCAN”手上有一批 AIS 轨迹数据少则几十条、多则上千条航线都在同一片海域里来回跑。人工一条条看能看出几条主走廊却说不清哪些轨迹真正属于同一类航行模式。航迹聚类的思路是把“一条轨迹”整体当成一个对象先算轨迹两两之间的距离再用聚类算法归组。标题里的“改进的 Harsdorf 距离”就是负责衡量两条航迹有多像的尺子。需要先说清一点学术资料里这个距离更常见的拼写是 Hausdorff工程笔记里写成 Harsdorf 的变体也不少本文按标题口径统一写“Harsdorf 距离”。而 DBSCAN 不需要预先指定聚成几类还能把离群的轨迹当噪声直接丢出去这两点对船舶航迹数据特别实用。适合谁手上有 AIS / 轨迹数据想知道船都在走哪几条固定走廊、哪些轨迹算异常、想按航路把船分组的人。2. 先搞清楚原始 Hausdorff 距离为什么不适合航迹聚类改进改在哪儿2.1 一条轨迹就是一个点集两个距离定义差的不是一点点把轨迹的时间维度先放一边一条航迹就是平面上的一个点集。要算轨迹 A 和轨迹 B 之间的距离常见做法是先定义一个“点到轨迹的距离”A 中的某个点 a 到 B 的距离是 a 与 B 中所有点的欧氏距离的最小值。原始 Hausdorff 距离在此基础上取的是“最坏情况”。单向 Hausdorff 距离是 A 中所有点到 B 的距离的最大值也就是 A 里离 B 最远的那个点它离 B 有多远。双向 Hausdorff 再把 B 到 A 的单向距离也算一遍取两个单向距离中的较大值。这个定义看上去对称、精细但它对离群点极其敏感。一条十公里长的正常航迹只要中间有一个点因为 AIS 丢星漂出去两公里max 这个操作就会把这个两公里直接算成整条轨迹的距离。于是这条轨迹和谁比都显得很远聚类时它要么被单独拎出来要么被算成噪声。船舶轨迹数据里这种野点几乎避免不了进出港时船速骤降、GPS 多径、设备重启每一条轨迹都能给你找出几个离群点。原始 Hausdorff 距离在这种数据上会让距离矩阵整体“虚高”聚类层次全被抹平完全看不出密度差异。2.2 改进方案分位数截断 双向平均 航向惩罚改进的 Harsdorf 距离要做三件事都是冲着原始定义的缺陷去的。第一把 max 换成 q 分位数将 A 中每个点到 B 的最近距离排个序取第 q 百分位的那个值作为单向距离。q 通常取 95意思是“允许最差的 5% 野点不参与计算”。这一下就把丢星、锚泊漂移这类问题压住了。第二把单向距离改成双向平均A 到 B 的 q 分位距离和 B 到 A 的 q 分位距离各算各的再取平均保持对称性。第三加一个可选的航向惩罚项两条几何上重叠的轨迹如果方向完全相反空间距离会非常小但语义上一条进港、一条出港肯定不该归成一类。航向惩罚的具体做法是算完双向平均距离后加上一项 lambda 乘以两条轨迹平均航向差的归一化值。航向差的绝对值取 0 到 180 度再除以 180压到 0 到 1 区间。lambda 通常从 0.3 起步设到 0.5 以上时反方向轨迹的距离会被显著拉开。要注意的是这个惩罚项应该只在两条轨迹空间上有重叠的区域才有意义一条在北、一条在南航向差再大也不能说明任何问题。所以工程上更稳妥的做法是先判断两条轨迹的包络矩形是否相交不相交就跳过航向惩罚只保留空间距离。2.3 为什么不直接用 DTW 或等间隔重采样有人会问轨迹对齐之后用欧氏距离或者用动态时间规整 DTW不是更常规吗这里对比一下就知道改进 Harsdorf 的优势在哪。等间隔重采样加欧氏距离实现最简单但要求两条轨迹点数一致、时间对齐而且只要有一个野点欧氏距离会被直接拉偏。DTW 对时间偏移非常鲁棒能处理船速不一致导致的轨迹伸缩但计算复杂度高轨迹一多距离矩阵开销吃不消。改进 Harsdorf 距离既不要求两条轨迹等长也不要求同采样率算的是“整体形状的接近程度”复杂度比 DTW 低一个量级还有 q 分位数兜底抗野点。方法优点不适合的场景等间隔重采样 欧氏距离实现快、容易理解轨迹长短不一、野点多、时间不同步DTW 动态时间规整对时间伸缩鲁棒轨迹数量大时计算量不可控改进 Harsdorf 距离不要求等长、抗野点、计算量适中完全忽略时间维度信息时间维度被丢掉确实是改进 Harsdorf 的代价。如果聚类结果里船速明显不同的轨迹被分到一起说明光靠空间形状不够需要回到距离函数里加 SOG 差异项或者按航次时长先做一轮规则过滤。这也是为什么很多人试完这个方案后又回去在距离函数里加参数而不是换聚类算法。3. 把 AIS 轨迹变成可计算的距离矩阵Matlab 数据流水线3.1 原始 AIS 记录切分成独立航次AIS 数据在数据库里通常是按报文一条条存的一个 MMSI 对应一条船但一条船会跨多个航次。如果不做切分船在港口停一天甚至跨一周的轨迹会被连成一条超长轨迹聚类结果会出现“跨区域长航线”这种毫无语义的簇。切分的原则有两个时间间隔和船速。同一艘船两条报文时间差超过阈值常见取 5 分钟说明中间要么关机要么停靠应该切一刀速度持续低于某个值如 SOG 小于 0.5 节超过 10 分钟也判定为停泊前后拆开。function trajCell splitAIS(rec) % rec: table列为 MMSI, time, lon, lat, sog, cog rec sortrows(rec, {MMSI, time}); t rec.time; dt [inf; seconds(diff(t))]; mmsiChange [true; diff(rec.MMSI) ~ 0]; timeGap dt 300; % 间断超过 300 秒切一刀 cutIdx find(mmsiChange | timeGap); cutIdx [cutIdx; height(rec) 1]; trajCell cell(length(cutIdx) - 1, 1); for i 1:length(cutIdx) - 1 seg rec(cutIdx(i):cutIdx(i 1) - 1, :); if height(seg) 10 % 少于 10 个点的轨迹不要 trajCell{i} seg; end end trajCell trajCell(~cellfun(isempty, trajCell)); end这里dt 300是时间切分的核心阈值300 秒对航速 10 节以上的船来说已经能开出 1.5 公里足以说明断开。代码里另一个容易被忽略的细节是height(seg) 10这个过滤条件少于 10 个报文的片段大概率是设备开机测试或者轨迹太短对聚类没有任何贡献只会增加距离矩阵的计算量。如果后续要用航向惩罚项切分时还要把sog、cog两列保留下来不要只留经纬度和时间。3.2 经纬度投影到本地平面距离不能直接用度数算经纬度是角度不是平面长度。纬度 1 度大约 111 公里经度 1 度在赤道也是约 111 公里但在 60 度纬度处只有约 55 公里。如果直接把经纬度差值当欧氏距离用高纬度区域的轨迹间距会被系统性放大聚类结果跟着变形。所以距离计算之前必须把经纬度投影到以米为单位的局部平面坐标。研究区域在几百公里以内时用一个简单的等距圆柱投影就够了以研究海域中心为原点经度方向乘以中心纬度的余弦修正。function [x, y] ll2local(lon, lat, lon0, lat0) % 经纬度转局部米制坐标区域尺度 500 km 时误差可接受 R 6371000; % 地球平均半径米 x R * deg2rad(lon - lon0) .* cos(deg2rad(lat0)); y R * deg2rad(lat - lat0); end投影中心lon0, lat0的选取会影响全部轨迹的坐标值但对距离计算的影响是整体的平移和缩放聚类结果不太会因此改变。需要留意的是cos(deg2rad(lat0))里的纬度必须用中心纬度而不是每条点的纬度否则投影不是平面变换距离就不再具有平移不变性。研究区域跨度超过 5 个经度时建议改用墨卡托投影Matlab 里可以直接调用 Mapping Toolbox 的projfwd没有这个工具箱就保持等距圆柱投影但聚类结果的边界簇要重新核对一遍。3.3 实现改进的 Harsdorf 距离矩阵距离矩阵是整个流水线的核心产出矩阵里的每个元素 D(i,j) 表示第 i 条轨迹和第 j 条轨迹之间的距离。改进 Harsdorf 距离的实现分三步求两条轨迹点集两两之间的欧氏距离矩阵对每一行取最小值得到“每个点离对方轨迹最近的距离”再对这个最近距离向量取 q 分位数。两个方向的 q 分位数都算出来取平均加上可选的航向惩罚。function d improvedHarsdorf(P, Q, q, lambda) % P, Q: n x 3 矩阵列分别为 x, y, cog航向角度 if isempty(P) || isempty(Q) d inf; return; end Dpq pdist2(P(:, 1:2), Q(:, 1:2)); % 轨迹点位两两距离 d1q quantile(min(Dpq, [], 2), q / 100); d2q quantile(min(Dpq, [], 1), q / 100); d 0.5 * (d1q d2q); if lambda 0 c1 mean(P(:, 3)); c2 mean(Q(:, 3)); dc abs(mod(c1 - c2 180, 360) - 180) / 180; d d lambda * dc; end endq取值 95 表示“最差的 5% 点被忽略”这个值在绝大多数航迹数据上都稳定除非你的轨迹里野点占比超过 5%那应该先回头清理数据而不是把 q 调低到 80。lambda是航向惩罚权重0 表示关闭需要区分对向航线时从 0.3 开始调。lambda 超过 0.5 后空间上完全重叠但方向相反的轨迹会被拉到 0.5 公里以上的距离差聚类结果会明显碎掉不建议一上来就给高分。两两距离矩阵的计算用pdist2在轨迹点数几百时没问题但整体复杂度是 O(K^2 * N * M)K 是轨迹数量N、M 是轨迹点数。所以下一步必须考虑抽稀。抽稀是控制计算量的关键。常见做法是对每条轨迹按弧长等间隔取 50 到 100 个点而不是按时间等间隔。按时间等间隔会保留大量低速段进出港的点而高速直线段只有零星几个点形状表达不均匀。实现也不复杂先对轨迹点求累计弧长再用interp1在等距弧长位置上插值。4. 跑 DBSCAN 并交互调参eps 和 minPts 怎么定4.1 距离矩阵喂进 DBSCAN工具箱版和手写版DBSCAN 的核心只有两个参数eps 是邻域半径minPts 是判定核心点所需的最少邻居数。Matlab 在 R2019a 之后的 Statistics Toolbox 里直接内置了dbscan函数可以传入预计算距离矩阵。如果没有带工具箱手写实现也不算复杂逻辑就是遍历未访问点找邻域内的点够 minPts 就扩张簇。function labels dbscanPrecomputed(D, eps, minPts) n size(D, 1); labels zeros(n, 1); % 0 表示噪声 visited false(n, 1); clusterId 0; for i 1:n if visited(i), continue; end visited(i) true; neigh find(D(i, :) eps); if numel(neigh) minPts labels(i) 0; % 邻居不够先记为噪声 else clusterId clusterId 1; labels(i) clusterId; queue neigh; while ~isempty(queue) j queue(1); queue(1) []; if ~visited(j) visited(j) true; labels(j) clusterId; jneigh find(D(j, :) eps); if numel(jneigh) minPts queue union(queue, jneigh); end end end end end end这段代码的关键在visited标记。访问过的点即使后来被别的簇扩展到了也不再改变归属这保证了每个点只可能属于一个簇。噪声点的处理是 DBSCAN 的特点初始标记为噪声的点如果后来被某个核心点扩展到状态会被改成该簇的点。所以不要在算法中途看到 label 为 0 就急着清理要等全部跑完。D(i, :) eps这条判断要求距离矩阵里不能有 NaN否则find会把 NaN 当不满足条件跳过邻域莫名其妙变小聚类结果完全不可复现。4.2 k-距离图选 eps别再靠猜eps 是 DBSCAN 里最容易翻车的参数很多人直接把 0.5 或者 1 填进去结果聚类出来 80% 是噪声。船舶轨迹距离的量级可能是几百米到几百公里eps 必须跟距离矩阵的量纲匹配。经典做法是画 k-距离图取 minPts 作为 k对每条轨迹找到它到所有其他轨迹的第 k 小距离不含自身按降序排列画折线曲线拐点就是合理的 eps。minPts 3; kDist sort(D, 2); kDist kDist(:, minPts 1); % 第一列是自身距离 0跳过 plot(sort(kDist, descend), o-); xlabel(轨迹编号按距离降序); ylabel(sprintf(第 %d 近邻距离, minPts)); grid on;拐点的判断有主观成分经验是找曲线从陡峭变平缓的“膝盖”位置对应横坐标即轨迹编号纵坐标即 eps。如果画出来是一条平滑斜线没有明显拐点说明数据里根本没有清晰的密度层次这时候不要硬找 eps应该回头检查距离函数是不是把轨迹长度差异抹平了或者轨迹数据本身质量太差。另一种常见误用是把 k 直接设成 minPts 的 2 倍这会让你看到的是“到第 6 近邻的距离”密度估计被过度平滑拐点更不明显。k 就取 minPts这是标准做法。4.3 用轮廓系数和噪声比例判断聚类结果好不好调完 eps 和 minPts不能只看图好看。量化指标最常用的两个轮廓系数和噪声比例。Matlab 的silhouette函数需要特征矩阵但我们的特征是轨迹间的距离矩阵可以先做多维缩放把距离矩阵嵌入二维平面再算轮廓系数。Y mdscale(D, 2); % 距离矩阵嵌入二维坐标 s silhouette(Y, labels); mean(s(labels 0)) % 噪声样本不参与评价 sum(labels 0) / numel(labels) % 噪声比例轮廓系数范围在 -1 到 1大于 0.25 说明簇内距离明显小于簇间距离结果可以接受低于 0.25 先怀疑距离函数里的 q 或 lambda 参数而不是急着调 eps。噪声比例参考区间是 5% 到 15%太低说明 eps 可能偏大簇之间边界模糊太高说明 eps 偏小真实密度结构被拆碎。这两个指标一配合比单看聚类图可靠得多。另外提醒一句mdscale在轨迹数量小于 10 条时意义不大数据量不够就别强行跑聚类先把轨迹预处理做好。5. 船舶轨迹聚类避坑与排查5 条血泪经验5.1 高纬度距离变形经纬度直接当平面坐标用现象低纬度海域聚类结果正常一换到高纬度海域同一批船型的轨迹聚类结果整个“糊”在一起簇边界完全对不上地理常识。原因经纬度差值是角度高纬度地区相同的经度差对应更短的实际距离。直接用经纬度算欧氏距离相当于在高纬把实际距离缩小了DBSCAN 按实际距离设定的 eps 在角度空间里对应了更大的范围聚类自然过松。解决所有距离计算前统一走ll2local投影聚类结果的坐标展示也基于投影后的点。研究区域跨 10 个纬度以上考虑分区投影不要试图用一个中心点覆盖整个大区域。5.2 eps 用默认值 0.5结果全是噪声现象聚类结果跑出来地图上只有几个孤零零的点其余轨迹全部标成噪声或者反过来所有轨迹并成一个簇。原因eps 的量纲必须匹配距离矩阵。轨迹距离动辄几十公里eps 填 0.5 等于要求轨迹之间距离小于 500 米大多数真实轨迹都不可能满足于是全是噪声。而 eps 填太大又会把整个海域的所有轨迹连到一起。解决先画 k-距离图看拐点再定 eps。如果两个看起来合理的候选值一个全噪声、一个全一簇说明数据里没有明显的密度层次回去检查轨迹预处理和距离函数。5.3 对向航线被并到同一簇现象进港和出港的船走同一条水道轨迹在空间上完全重合聚类结果把它们归成了同一类下游统计航次流量时方向直接抵消。原因改进 Harsdorf 距离在 lambda 0 时只看空间形状对“重叠但反向”的两条轨迹计算出的距离很小DBSCAN 自然认为它们是邻居。解决在距离函数里打开航向惩罚lambda 从 0.3 起步。如果业务上方向敏感也可以对轨迹做方向归一化统一把轨迹整理成从起点到终点的方向反向轨迹会被自动分开但这种方法会破坏后续的代表轨迹提取语义不如直接加航向惩罚干净。5.4 一条船多个航次被当成一条轨迹现象聚类结果里出现一条超长轨迹跨了港口到港口甚至横跨整个研究区域两侧和周围的簇完全对不上。原因航次切分失效。时间阈值设得太大船在港口停几个小时被当成正常航行或者 SOG 持续低值判停泊的阈值写错导致停泊前后没有被切开。解决切分代码里同时启用时间间断和速度双条件不要只依赖一个阈值。切完以后按轨迹总时长和总里程做分布统计把明显偏长的轨迹单独挑出来重看切分参数。5.5 距离矩阵算到天荒地老算法复杂度失控现象500 条轨迹每条轨迹 3000 个点跑了两小时还没出距离矩阵Matlab 占内存 16GB 以上。原因两两距离复杂度 O(K^2 * N * M)K500、NM3000 时单次计算量是 10^12 量级不管用什么矩阵加速都救不回来。解决轨迹先抽稀到 50 到 100 个点再算距离。抽稀后精度损失对聚类影响很小因为改进 Harsdorf 本身用的是分位数去掉的点本来就有一定冗余。计算端用parfor按行并行避免一次性构造超大矩阵。D inf(K, K); parfor i 1:K row zeros(1, K); for j i 1:K row(j) improvedHarsdorf(trajCell{i}, trajCell{j}, q, lambda); end D(i, :) row; end D D D;parfor里每轮只写一行D(i, :)这是 Matlab 并行切片允许的写入方式不需要用parpool手动开池。D提前用inf填充对角线保持 inf后面跑 DBSCAN 时邻域判断 eps不会把自身算进去省一次显式的去零处理。如果你发现parfor在远程服务器上起不来先检查是不是用了动态数据结构cell 数组最好先预分配再填。6. 让聚类结果直接用于航路分析三个小习惯6.1 每个簇用 medoid 轨迹当代表航路聚类完成后散点图能看轮廓但报告里需要的是一条可读的“代表航路”。对每个簇找出簇内到其他所有轨迹的 Harsdorf 距离之和最小的那条轨迹作为 medoid它是最能代表整簇几何特征的轨迹。idx find(labels 3); % 以簇 3 为例 subD D(idx, idx); [~, pos] min(sum(subD, 2)); centerTraj trajCell{idx(pos)}; % 代表航路用 medoid 而不是平均值是因为轨迹平均需要逐点对齐两条轨迹点数和坐标都不同平均出来的“船迹”很可能不再是任何一条真实航路甚至可能穿过陆地。medoid 是实实在在存在的轨迹查 MMSI、查时间、查航行意图都有据可依。这个习惯帮我少解释了无数次“这条线是哪里来的”。6.2 按平均 SOG、COG 和起终点给簇贴语义标签聚类结果最终要落到业务词汇上。对每个簇统计平均航速、平均航向、起点区域和终点区域自动生成“进港 1 号走廊”“过境直线航路”这类标签下游做航路流量统计和异常检测时直接按标签过滤。数据量少时每个簇画一张轨迹叠加图和一个统计表就够数据量大时对簇内轨迹点做二维核密度估计提取高密度骨架线比直接堆散点直观得多。6.3 一个自己的实验管理习惯我每次改完距离函数或 DBSCAN 参数都会把三个数记下来k-距离图的拐点位置、eps 的最终取值、噪声比例和平均轮廓系数。改的是 Harsdorf 距离里的 q 还是 lambda还是改的 eps 和 minPts分开记录。这个习惯帮我避免了最典型的翻车——换了一批海域数据后同样参数跑出来一团糟却完全想不起上一次是怎么调出来的。聚类本来就是“距离定义 参数选择”的双层问题只记参数不记距离函数等于没有实验记录。希望这些细节能帮你在自己的航迹数据上少走几趟弯路也希望你下一次调参时第一件事是画 k-距离图而不是猜 eps。本文还有配套的精品资源点击获取