简介这套Matlab项目基于改进的豪斯多夫距离实现DBSCAN船舶航迹聚类复现了《基于轨迹聚类的船舶异常行为识别研究》论文中的关键流程。面向船舶航迹研究、异常行为识别以及希望将密度聚类用于实际数据的开发者可直接运行或替换数据做扩展也可将航迹提取与聚类模块迁移到其他轨迹分析场景。压缩包共20个文件以14个M脚本为主辅以MAT数据文件、结果图、说明文档和示例压缩包整体仅4.32MB结构清晰便于定位各算法模块并配有使用说明帮助快速上手。已有765人学习下载代码完整覆盖航迹数据提取、DBSCAN聚类、聚类中心提炼、基于豪斯多夫距离的阈值分类与阈值寻优并在此基础上实现航迹偏离预测经测试准确率较高。既能作为DBSCAN和豪斯多夫距离的落地案例学习也能为船舶异常行为识别提供可复用的实验框架。1. 换一把量航迹的尺子Harsdorf 距离改进 DBSCAN 的 Matlab 实现值在哪AIS 数据越攒越多真正想回答的问题却始终绕不开一个哪几条船走的是同一条航线传统 DBSCAN 直接把经纬度点丢进去分出来的是点团不是航迹。这份 Matlab 项目做了个替换——把每艘船一整条轨迹当成一个对象用 Harsdorf 距离去量轨迹之间的“不相似程度”再喂给改进后的 DBSCAN 做航迹聚类。跑出来的结果是整条整条的航迹簇港口交通流、异常绕行、固定班轮一眼就能看出来。适合手里有 AIS 数据、想按轨迹而不是按点做分析的从业者和研究生。2. 为什么是 Harsdorf 距离轨迹相似度度量的选型逻辑2.1 欧氏距离度量轨迹的三个死穴做航迹聚类之前最容易被惯性带偏的点是直接用 DBSCAN 聚类坐标点。DBSCAN 的标准输入是 N 个样本点每个样本点是一组坐标邻域查询用欧氏距离。但船迹不是点是变长的点序列直接套点聚类至少有三个问题。第一轨迹长度天然不对齐。一艘远洋船从宁波到新加坡AIS 报文能积累几千个点另一艘港内拖轮作业半小时可能只有几十个点。把两个长度差两个数量级的序列硬塞进欧氏距离要先对齐长度对齐本身就是个灾难。第二时间不对齐。两艘船走同一条航路一艘凌晨过点、一艘中午过点时间戳完全不同。按时间索引对齐后匹配上的是完全不同空间位置的点算出来的欧氏距离大得离谱真实的空间相似性被时间维度搅浑。第三空间上小偏移会被放大。两艘船前后相差 200 米平行进港形状几乎一样但点对点欧氏距离对平移极其敏感算出来两条轨迹“差距很大”。航迹聚类要的是整体形态相似不是逐点坐标吻合。所以这个项目把度量对象从“点”换成“整条轨迹”用 Harsdorf 距离来刻画两条轨迹之间的最大不相似程度。Harsdorf 距离不要求两条轨迹等长、不要求时间对齐只需要两串空间点之间存在一个距离定义。这个特性让它成为轨迹聚类的第一选择。2.2 Harsdorf 距离定义与 Matlab 函数实现Harsdorf 距离的思想很朴素对轨迹 A 上的每个点找到它在轨迹 B 上最近的那个点记下这段距离A 上所有点都算完之后取最大值得到 A 到 B 的单向距离。反过来再算一遍 B 到 A 的单向距离两者取大就是双向 Harsdorf 距离。function d hausdorff_dist(traj1, traj2) % 双向 Harsdorf 距离 % traj1, traj2: N×2 或 N×3 矩阵前两列为平面坐标x, y % 返回标量 d两条轨迹之间的最大不相似度单位与坐标一致 % A - B 单向A 上每点到 B 的最小距离再取最大 d1 0; for i 1:size(traj1, 1) diff traj2(:, 1:2) - traj1(i, 1:2); dists sqrt(sum(diff.^2, 2)); d1 max(d1, min(dists)); end % B - A 单向B 上每点到 A 的最小距离再取最大 d2 0; for j 1:size(traj2, 1) diff traj1(:, 1:2) - traj2(j, 1:2); dists sqrt(sum(diff.^2, 2)); d2 max(d2, min(dists)); end d max(d1, d2); end两个单向循环的逻辑完全对称先找一个方向上的“最大缝隙”再看另一个方向。d 越小说明两条轨迹整体形态越接近d 越大说明至少有一个方向的某个点在另一条轨迹上找不到近邻也就是存在明显偏移或分叉。这段代码的时间复杂度是 O(N×M)N 和 M 是两条轨迹的点数。几百条轨迹两两互算还能接受上千条就要做优化。常见做法是先用 pdist2 把逐点距离矩阵一次性算出来再对每行取最小、对列取最大可以省掉内部的平方开方循环。另外这里默认输入是平面坐标如果手上是经纬度先把经纬度投影成墨卡托平面坐标再来算否则直接用经纬度算欧氏距离在高纬度地区误差会越来越大。2.3 轨迹预处理重采样、清洗与坐标转换Harsdorf 距离不怕轨迹长短不一但怕脏数据。AIS 原始报文里常见三类问题坐标越界、速度为零的长期抛锚点、以及信号丢失导致的时间断层。这些噪点会直接污染距离矩阵所以预处理阶段要做的事不少。第一步是清洗把纬度小于 -90 或大于 90、经度超出合理海域范围、以及对地航速为 0 但持续超过阈值的点先摘出来判断是抛锚还是设备异常。第二步是重采样AIS 报文间隔本身不固定一般用等时间间隔线性插值把每条轨迹统一到相同采样密度。function trajOut resample_traj(lon, lat, t, dt) % 按固定时间间隔 dt 对轨迹点做线性插值重采样 % lon, lat: 原始经度/纬度序列t: 时间戳秒 % 返回 N×3 矩阵x, y, t t0 t(1); t1 t(end); tNew (t0:dt:t1); lonNew interp1(t, lon, tNew, linear); latNew interp1(t, lat, tNew, linear); % 经纬度转平面坐标等距圆柱投影适合中小范围海域 % 1 度纬度约 111320 米经度按平均纬度缩放 lat0 mean(latNew); x lonNew * 111320 * cosd(lat0); y latNew * 111320; trajOut [x, y, tNew]; enddt 的取值决定了轨迹的细节保留程度。60 秒一个点适合大多数沿海航迹分析能保留航向变化的主要形态同时把噪声点的影响压下去。dt 太小会放大 AIS 抖动dt 太大又会让小尺度绕行被插值抹平。我一般先画几条典型轨迹看 60 秒插值后有没有明显失真再决定要不要收紧到 30 秒。坐标转换这一步经常被忽略。AIS 给的是经纬度Harsdorf 距离要的是平面坐标。上面代码用的是等距圆柱投影在港区这类几十公里范围的场景精度够用如果做跨洋航迹聚类建议换墨卡托投影避免高纬度纬线间距失真。3. 改进 DBSCAN 的核心实现从距离矩阵到聚类结果3.1 全轨迹距离矩阵的批量计算预处理完的轨迹集合是一组 cell 数组每个元素是 N×3 的矩阵。改进后的 DBSCAN 不再接收坐标点而是接收轨迹之间的距离矩阵。矩阵的每个元素 D(i,j) 表示第 i 条轨迹和第 j 条轨迹之间的 Harsdorf 距离。function D build_traj_dist_matrix(trajs) % 批量计算所有轨迹两两之间的 Harsdorf 距离 % trajs: 1×K cell每个元素是 resample 后的 N×3 轨迹矩阵 % 返回 K×K 距离矩阵single 类型节省内存 K numel(trajs); D zeros(K, K, single); for i 1:K for j i1:K d hausdorff_dist(trajs{i}, trajs{j}); D(i, j) d; D(j, i) d; end D(i, i) 0; end end这段有两个细节值得留意。第一只算上三角再对称赋值避免每条轨迹对算两遍直接省一半计算时间。第二矩阵用 single 存储。K1500 条轨迹时double 矩阵是 1500×1500×8 字节等于 18MBsingle 直接减半。真跑到上万条轨迹时这个矩阵会是几百 MB 级别分块计算或写 parfor 并行就是刚需了。% 多核并行版本的核心改动parfor 取代内层 for parpool(local, 6); for i 1:K parfor j i1:K ... end end用 parfor 时注意Harsdorf 函数本身计算量不小在 8 核机器上开 parfor 能明显提速但如果轨迹数量只有几十条每次计算的耗时不足以抵消并行调度的开销反而更慢。我一般阈值设成 K 大于 500 才开并行。3.2 基于距离矩阵的 DBSCAN 主循环标准 DBSCAN 的核心是三件事找邻域、判定核心点、扩展簇。改进版把邻域查询从“以某点为圆心画半径 eps 的圆”换成“从距离矩阵取一行找出所有小于 eps 的轨迹下标”。算法骨架完全不变变的只有邻域的定义方式。function labels dbscan_trajectory(D, eps, minPts) % 基于 Harsdorf 距离矩阵的改进 DBSCAN % D: K×K 距离矩阵 % eps: 距离阈值两条轨迹距离小于 eps 视为邻域轨迹 % minPts: 核心轨迹的最少邻域轨迹数 % labels: K×1-1 表示噪声1 表示簇编号 K size(D, 1); labels zeros(K, 1); clusterId 0; for i 1:K if labels(i) ~ 0 continue; end % 邻居查找距离矩阵第 i 行中 eps 的所有轨迹 neighbors find(D(i, :) eps); if numel(neighbors) minPts labels(i) -1; % 噪声轨迹 continue; end clusterId clusterId 1; labels(i) clusterId; seeds neighbors; while ~isempty(seeds) idx seeds(1); seeds(1) []; if labels(idx) -1 % 边界轨迹归入当前簇 labels(idx) clusterId; elseif labels(idx) 0 labels(idx) clusterId; nbrs find(D(idx, :) eps); if numel(nbrs) minPts % 核心点扩展种子集合 seeds [seeds; nbrs]; end end end end end主循环里有个小坑find(D(i, :) eps) 会把轨迹自身也包含进去因为 D(i, i) 是 0必然小于 eps。这在判定 minPts 时相当于白送一个计数。如果 minPts 取 2实际要求是“自己加另一个邻居”语义上说得通如果 minPts 取 1任何一条非孤立轨迹都变成核心轨迹几乎不会产生噪声。所以我对 minPts 的建议是至少取 2最好取 3 到 5把自邻域的影响稀释掉。3.3 聚类评估与可视化聚类跑完要回答一个现实问题这次分簇分得好不好对航迹聚类来说轮廓系数是性价比最高的评估指标Matlab 自带 silhouette 函数但这里不能直接用原始坐标必须用距离矩阵喂进去。% 距离矩阵版本silhouette 允许传距离矩阵而不用原始数据 s silhouette(labels, D); meanSil mean(s); fprintf(平均轮廓系数: %.3f\n, meanSil);轮廓系数在 -1 到 1 之间越接近 1 说明簇内轨迹之间 Harsdorf 距离显著小于簇间距离。平均轮廓系数 0.5 以上算是可接受0.7 以上是相当清晰的聚类。可视化也和平常画散点不一样轨迹是折线按簇分组上色画线figure; hold on; colors lines(max(labels)); for i 1:numel(trajs) if labels(i) -1 plot(trajs{i}(:, 1), trajs{i}(:, 2), Color, [0.7 0.7 0.7], LineWidth, 0.5); else plot(trajs{i}(:, 1), trajs{i}(:, 2), Color, colors(labels(i), :), LineWidth, 1); end end灰色轨迹就是 DBSCAN 判定为噪声的船彩色轨迹按簇聚合。这一步非常直观画完之后港口哪些航路是主交通流哪些是孤立绕行比任何统计数字都清楚。我习惯在画图代码后面加一段打印每个簇的轨迹条数、平均轨迹长度、簇内最大 Harsdorf 距离。这三个数字配合轮廓系数能快速暴露参数问题。4. Eps 和 MinPts 不调好聚类全是噪点参数调试三板斧4.1 K 距离图定 Eps拐点法实操DBSCAN 最让新手头疼的就是 eps 和 minPts 没有标准答案。但航迹聚类有更具体的问题eps 的物理含义是“两条轨迹之间的最大允许偏移量”这个量纲是米。因此要比通用 DBSCAN 更直观但选值仍然需要数据辅助。K 距离图是选 eps 最经典的手段核心思路是把每条轨迹到其最近第 k 条轨迹的距离从小到大排列画成曲线曲线上拐点位置的纵坐标就是建议的 eps。function kDist k_distance_curve(D, k) % 计算 K 距离图每条轨迹到第 k 近轨迹的距离 % D: 距离矩阵k: 近邻序号一般用 minPts-1 或 minPts % 返回按降序排列的 K 距离序列 K size(D, 1); kDist zeros(K, 1); for i 1:K row D(i, :); row(i) inf; % 排除自身 srow sort(row); kDist(i) srow(min(k, K-1)); end kDist sort(kDist, descend); plot(kDist, LineWidth, 1.5); xlabel(轨迹序号); ylabel([第 , num2str(k), 近 Harsdorf 距离]); end注意排除自身距离 D(i,i)0否则每个轨迹的最近邻都是自己整条曲线全被 0 带偏。k 的取值跟随 minPtsminPts3 时 k 也取 3曲线拐点会比较明显。调用后观察曲线形态如果曲线像一条平滑下降的斜线拐点不清晰说明轨迹之间的距离分布很均匀没有天然簇结构如果曲线前段急速下降、中段出现明显平台或拐弯拐点的纵坐标就是 eps 的直接候选值。我习惯取拐点纵坐标乘以 0.8 到 1.2 的范围再下去细调。4.2 网格搜索配轮廓系数K 距离图只能给出候选区间真正定值还是得靠网格搜索。航迹聚类的搜索空间不大eps 在候选值附近取 5 到 8 个值minPts 取 2、3、5 三个值跑 20 组左右每组记录簇数、噪声占比、平均轮廓系数。function bestParams grid_search_eps_minpts(D, epsRange, minPtsRange) % 简单网格搜索遍历 eps 和 minPts 组合返回轮廓系数最优的一组 % epsRange: 候选 eps 向量minPtsRange: 候选 minPts 向量 bestScore -1; bestParams []; for e epsRange for m minPtsRange labels dbscan_trajectory(D, e, m); valid labels 0; nClusters max(labels); nNoise sum(labels -1); if nClusters 2 continue; % 全聚成一个簇淘汰 end s silhouette(labels, D); score mean(s); fprintf(eps%.3f minPts%d - 簇数%d 噪点%d 轮廓%.3f\n, ... e, m, nClusters, nNoise, score); if score bestScore bestScore score; bestParams [e, m]; end end end end每次跑完打印那几个关键数字不是为了刷屏而是为了看趋势轮廓系数最优的那组参数附近如果簇数和噪点占比变化剧烈说明聚类结构不稳定如果轮廓系数最优时噪点占比超过 30%大概率是 eps 偏小或 minPts 偏大边界轨迹全被误杀。关于网格搜索有一段血泪经验别追求轮廓系数的唯一最大值要取一个“平坦区间”的中间值。轮廓系数在 2 到 3 个相邻 eps 取值上都很接近时取中间那个 eps这种参数对数据噪声的鲁棒性最好。4.3 批量跑次验证稳定性网格搜索跑完一轮不等于参数就能直接用。真实 AIS 数据集往往有时间切片比如某港口一周的数据可以按天切、按班次切。同一个参数在不同切片上能不能跑出结构一致的簇是判断参数真正有效的标准。% 对每天的数据分别聚类统计簇数分布 days unique(t_raw(:, 4)); % 假设第 4 列是日期编号 nClustersPerDay zeros(numel(days), 1); for i 1:numel(days) idx t_raw(:, 4) days(i); D_day build_traj_dist_matrix(trajs(idx)); labels_day dbscan_trajectory(D_day, bestEps, bestMinPts); nClustersPerDay(i) max(labels_day); end disp(nClustersPerDay);簇数在一天内波动超过正负 30%说明参数过拟合了某一个特定数据分布就得把搜索范围放宽重来。这一步没有玄学本质是检验参数对数据密度差异的容忍度。船舶轨迹聚类有个特性工作日与周末、白天与夜间的船舶密度差异非常大一套参数能同时适应这两种密度才算合格。另外要特别提一下 minPts 的取值边界。轨迹聚类场景下 minPts 很少需要超过 5因为同一条航路在某一时间段内真正并行的船就那么多。minPts 设太大小规模船队直接全部变噪声设太小两三条离群轨迹拉到一起就成簇。港口场景我从 minPts3 起步远洋航线场景从 minPts2 起步然后按 K 距离图和轮廓系数往回调。参数推荐范围实用经验epsK 距离图拐点 ±20%偏大一条龙大簇偏小全是噪点minPts2 ~ 5港口取 3远洋取 2密度不均时取 5重采样间隔30s ~ 120s60s 起步看航向变化是否被抹平距离矩阵精度single2000 条轨迹时比 double 省一半内存5. 航迹聚类避坑手册五条实测踩坑记录5.1 一条轨迹被拆成两簇AIS 丢包陷阱现象同一条船连续 10 小时的完整航迹聚类结果里首尾竟然分属两个簇中间还夹着噪点标记。原因AIS 信号在近岸基站覆盖边缘经常丢失轨迹中间会出现一两个小时的空窗。重采样时线性插值会在这个空窗期拉出一条直线这条虚拟线段和任何真实轨迹都不相似Harsdorf 距离对离群点又是取最大值一条轨迹直接被自己的虚假线段推离原簇。解决预处理的清洗阶段先做断点检测相邻 AIS 时刻间隔超过阈值我一般用 30 分钟就沿断点把轨迹切成两段再重采样。这既保住了原始轨迹的大部分有效信息又避免插值虚空段污染距离矩阵。% 断点检测时间间隔超过阈值则切分 tDiff diff(t); cutIdx find(tDiff 30 * 60);5.2 长轨迹主导相似度Harsdorf 的 max 放大效应现象聚类结果里最大的一簇几乎全部是远洋长航迹近海短航迹七零八落全是噪点轮廓系数倒是很高但一看簇内轨迹形态两条船根本不在同一片海域。原因Harsdorf 距离的定义取的是最大值长轨迹航程几千公里在航路末端的任何一个离群点都会产生一个极大距离值。短轨迹与长轨迹的 Harsdorf 距离被这个极值主导短轨迹之间即使形态一致整体相似度反而不如长轨迹内部的相互“包容”。解决换成截断式 Harsdorf 距离。把两个方向的点距离序列先排序只取前 95% 的距离值再取最大砍掉长尾离群点的影响。这在 Matlab 里改动很小在 hausdorff_dist 内部把 max 换成对排序后 95% 分位的 max但效果显著。5.3 Out of Memory距离矩阵吃到内存爆掉现象轨迹数量攒到 2000 条时Matlab 直接报 Out of Memory窗口卡死。原因距离矩阵本身是 2000×2000single 类型只有 16MB不算大。但 build_traj_dist_matrix 里每次调用 hausdorff_dist 都会生成临时矩阵2000×2000 次调用的临时分配累计起来非常可观内存碎片化严重。解决先预分配 D 矩阵并固定为 single内部函数里用局部变量减少重复分配批量计算时每 100 条轨迹存一次中间结果到磁盘跑完之后合并。再不够就上 parfor把内层循环分到多个 worker每个 worker 只持有部分轨迹数据。5.4 中文注释乱码编码问题导致的代码不可读现象打开项目源码中文注释全部变成“锟斤拷”一类乱码部分脚本直接无法运行。原因Matlab 2023 中文版默认编辑器编码是 GBK而项目文件保存为 UTF-8文件头没有被正确识别中文字符逐个字节被错误解析。解决先在 Matlab 主页预设里把“语言和位置”调成 “English (United States)” 或指定 UTF-8 编码再重新打开文件。更稳妥的做法是运行前执行一次feature(DefaultCharacterSet, UTF-8)或者在文件开头用英文注释。从那以后我再也不在有中文注释的 Matlab 工程里直接双击打开了。5.5 NaN 传染轨迹里一个坏点毁掉全部距离现象聚类结果里所有轨迹标签全是 0部分输出直接 NaN找半天没发现代码逻辑有问题。原因AIS 原始数据里有 NaN 坐标没被清洗掉hausdorff_dist 内部的 sum 和 min 遇到 NaN 会把结果一并带成 NaNNaN 再顺着距离矩阵传播给 DBSCANfind(D(i,:)eps) 永远返回空集。解决预处理阶段强制加一行trajOut(~isfinite(trajOut(:, 1:2))) NaN; trajOut(any(isnan(trajOut(:, 1:2)), 2), :) [];宁缺毋滥坏的坐标点直接删行绝对不能让 NaN 进入距离计算。validIdx isfinite(traj(:, 1)) isfinite(traj(:, 2)); traj traj(validIdx, :);6. 验证与进阶合成轨迹先跑通真实数据再上船6.1 合成数据验证参数调完别急着上真实 AIS 数据。先用合成轨迹把整套流程验证一遍这一步能省掉大量排查时间。造两组已知答案的轨迹一组走直线带微小幅度的随机弯曲另一组走 S 形路线把这两组数据混在一起跑完聚类后看能否完整分回两组。% 合成轨迹A 组直线抖动B 组 S 形 for i 1:20 t (0:0.1:10); noiseA randn(size(t)) * 0.02; trajs{i} [t, 0.2*t noiseA, t]; end for i 21:40 t (0:0.1:10); trajs{i} [t, sin(t) randn(size(t)) * 0.02, t]; end跑完如果 A 组和 B 组被干净分开说明距离度量、聚类参数都站得住。合成都分不开真实数据只会更糟。这个方法能帮你把整套流程的项目逻辑从黑匣子变成可控实验。6.2 两个可行的进阶方向第一个进阶方向是截断 Harsdorf 距离。前面避坑里提到长轨迹主导问题把 max 换成排序后按 95% 分位截断对离群点和 AIS 漂移的鲁棒性明显更好代价是多一次排序计算时间上升约 20%换来的聚类质量提升值这个价。第二个进阶方向是轨迹分段聚类。船舶行为并不总是“一条船 一条完整行为”航行、锚泊、靠泊三个阶段形态差异极大。先把轨迹按速度或转向率切成航段再对航段做 Harsdorf 距离聚类出来的簇更贴近“行为”而不是“航次”。这两步做完基本就能把这份工程扩展到港口交通流分析之外的异常绕行检测、渔船偷捕识别等高价值场景。从那以后我每次换数据集都强制先用合成样本验证一遍参数再放开跑。这个习惯帮我躲过了好几次“看起来像模像样、实际全在乱聚”的翻车现场。希望帮到你。本文还有配套的精品资源点击获取
