简介面向船舶航迹研究与Matlab聚类学习者的实战代码包围绕改进Hausdorff距离的DBSCAN航迹聚类展开复现《基于轨迹聚类的船舶异常行为识别研究》的关键流程。内容覆盖航迹数据提取、DBSCAN聚类、聚类中心提取、H距离计算、基于H距离的阈值分类与阈值寻优并延伸到航迹偏离预测与准确率检验适合作为聚类算法应用案例也可为船舶航迹方向的模型拓展提供基础。压缩包共20个文件约4.32MB以14个m脚本为核心配合mat数据、zip绘图与结果包、png示意图及md说明模块划分清晰可直接运行对照对数据预处理、距离参数影响与结果可视化也有相应脚本支撑。已有768人学习下载。读者可据此掌握DBSCAN与H距离计算、航迹分类的大致流程和预测阈值调优思路并通过替换数据与相关模块迁移到自己的研究场景。1. 一条航迹一条曲线为什么传统 DBSCAN 在船舶轨迹上总翻车某开发者丢给我一份 AIS 导出表两万多条航迹想用 matlab 做聚类再往下接船舶航迹预测。他直接调了内置的dbscan把经纬度当两维特征喂进去结果一次聚出几百个碎片簇一次又全部塞进同一个簇噪声点占了六成。这不是参数没调好是模型的假设错了DBSCAN 的邻域是「点与点」的欧氏距离而一条航迹是长度不一、采样不均的点序列两条航迹之间的「远近」根本不是某一个点能代表的。所以这套工程的第一块骨头是把距离度量从「点对点」换成「曲线对曲线」——也就是标题里说的 Harsdorf一般写作 Hausdorff豪斯多夫距离并且要对它做改进。下面按数据清洗、距离改造、聚类主干、踩坑、二次开发这条线把这个 matlab 船舶航迹聚类工程拆开讲。2. 工程结构与数据契约从原始 AIS 点到可聚类的特征矩阵拿到一份航迹聚类代码先别急着看算法先看它怎么把一张 AIS 表变成「航迹集合」。这一步错了后面距离算得再精致也白搭。这个工程的定位很明确离线批处理输入一段时间的 AIS 轨迹记录输出每条航迹的簇标签、簇中心航迹和可视化图服务于后续的航迹预测与异常航迹筛查。2.1 为什么是 MATLAB矩阵运算、工具箱与可视化闭环选 matlab 不是情怀是这条流水线的每一步都吃矩阵。AIS 数据本质是「航迹 ID 时间 经纬度 运动状态」的长表清洗阶段是分组过滤距离阶段是一个 N×N 的成对矩阵聚类阶段反复切片查邻域最后还要把几千条折线画在同一张图上比对效果。matlab 的groupsummary、pdist2、quantile、table加上geobubble/plot这一套基本不用引第三方库就能闭环对做算法验证的人来说省下来的时间比运行速度更值钱。代价也要说清楚。N×N 距离矩阵在 matlab 里默认是doubleN8000 就是 8000×8000×8 字节约 512MBN 再翻一倍直接吃掉 2GB这是后面避坑章节里第一个要处理的现场。另外 matlab 没有原生的空间索引加速DBSCAN 的邻域查询只能靠预计算矩阵硬扛所以这个工程的设计前提是「中等规模、离线跑、可接受几分钟」。要上几十万条航迹得换语言或者分块这一点在架构上就得认。工程目录大概长这样不追求花哨但每一层职责清晰ShipTrackCluster/ ├── data/ │ ├── raw/ % 原始 AIS 导出csv / xlsx只读 │ └── proc/ % 清洗后的 .mat一条航迹一个 cell ├── src/ │ ├── preprocess/ │ │ ├── loadAIS.m % 读表 字段归一化 │ │ ├── cleanTrack.m % 单条航迹的合法性过滤 │ │ └── resampleTrack.m % 按弧长重采样到固定点数 │ ├── distance/ │ │ ├── latlon2xy.m % 经纬度 - 平面米制 │ │ ├── hausdorffStd.m % 经典 Hausdorff作对照 │ │ └── hausdorffImproved.m % 改进版工程实际使用 │ ├── cluster/ │ │ ├── dbscanByDist.m % 距离矩阵驱动的 DBSCAN 主干 │ │ └── kdistCurve.m % k-距离曲线辅助定 eps │ └── eval/ │ ├── clusterMetrics.m % 噪声率、簇数、轮廓系数 │ └── plotTracks.m % 按簇着色画航迹 ├── config/ │ └── params.m % 所有阈值集中在这里 ├── main_cluster.m % 一键跑通入口 └── README.mdconfig/params.m单独拎出来是我强烈建议保留的习惯。航迹聚类的阈值特别多最小点数、最小航程、重采样点数、分位数、eps、minPts散落在各个函数里改一次参数要翻五个文件血泪经验。2.2 数据契约AIS 轨迹表的字段约定与清洗口径上下游要对接就得先定字段。这个工程假定输入是一张宽表最少包含下面几列列名允许不同但要在loadAIS.m里映射过来字段含义用途常见坑MMSI船舶标识航迹分组键存在非船舶 MMSI需按号段过滤BaseDateTime时间戳排序、切分航迹时区不统一混了本地时间LAT / LON纬度 / 经度几何位置个别记录为 0 或 91是设备默认值SOG对地航速过滤停泊点单位是节不是 km/hCOG对地航向方向项校验0~360359 与 1 相邻不能直接相减Status航行状态剔除锚泊/系泊编码各家不一需按文档映射清洗的核心逻辑是「先按船分组再按时间排序再做几何与业务过滤」。这一层不要写成一个大函数拆成单航迹级别的cleanTrack更好测。下面是分组清洗的主干function tracks loadAIS(csvFile, cfg) % 读入原始 AIS 表按 MMSI 分组输出 cell 数组每个元素为 N×2 的经纬度序列 T readtable(csvFile); T.Properties.VariableNames normalizeNames(T.Properties.VariableNames); % 1) 基础合法性经纬度范围、时间有效性 valid T.LAT -90 T.LAT 90 T.LON -180 T.LON 180 ... ~isnat(T.BaseDateTime) T.SOG 0 T.SOG 60; T T(valid, :); % 2) 按船 时间排序保证航迹的方向是有序的 T sortrows(T, {MMSI, BaseDateTime}); % 3) 时间切分相邻两点间隔超过 cfg.maxGapMin 分钟就断开成新航迹 g findgroups(T.MMSI); tracks {}; for k 1:max(g) sub T(g k, :); if height(sub) cfg.minPoints, continue; end segStart [true; minutes(diff(sub.BaseDateTime)) cfg.maxGapMin]; segId cumsum(segStart); for s unique(segId) seg sub(segId s, :); trk cleanTrack([seg.LON, seg.LAT], seg.SOG, cfg); if ~isempty(trk), tracks{end1} trk; end %#okAGROW end end fprintf(清洗完成保留 %d 条航迹\n, numel(tracks)); end这段的关键在三个参数。cfg.maxGapMin决定「断轨」的粒度AIS 岸基接收有盲区同一条船中间断两小时很常见间隔阈值定小了会把一条远洋航迹切成七八段定大了会把两次不同航次拼成一条假航迹我一般从 30 分钟起调看数据的时间间隔分布再定。cfg.minPoints是航迹最少点数低于这个数的航迹没法算形状距离只能当噪声丢掉常用 20~30。cleanTrack里再做速度突跳剔除和停泊点裁剪把 SOG 长期小于 0.5 节的端点段切掉因为锚泊船的轨迹本质上是一个点团参与曲线聚类只会污染结果。提示切分之前一定要按时间排序。AIS 导出表经常是按接收时间而非船舶时间排列的不排序直接diff算时间间隔得到的断轨结果完全是随机的这个坑非常隐蔽。2.3 重采样把变长序列压成可比较的固定点数清洗完的航迹长度从几十点到上千点不等直接算 Hausdorff 虽然不要求等长但点密度差异会让「最近点距离」失真——采样密的航迹天然更容易靠近对方。所以工程里加了一步按弧长重采样把所有航迹统一到cfg.resampleN个点默认 50。50 个点描述一条几十公里的航迹形状保真度够用同时把后面点对距离矩阵的规模压到可控范围。function P resampleTrack(P, n) % 按累积弧长等距重采样P 为 N×2 平面坐标 segLen sqrt(sum(diff(P).^2, 2)); % 每段长度 s [0; cumsum(segLen)]; % 累积弧长 if s(end) eps, P repmat(P(1,:), n, 1); return; end si linspace(0, s(end), n); % 目标弧长序列 P [interp1(s, P(:,1), si, linear), interp1(s, P(:,2), si, linear)]; end重采样用弧长而不是时间是因为我们要比较的是「形状」不是「速度剖面」。同一条航线上一艘快船和一艘慢船时间采样点分布完全不同但几何形状一致按弧长对齐后距离才会小。这也是后面簇内求平均航迹的基础所有航迹点数一致逐点取中位数就是一条合法的代表航迹。代码里有个边界情况要处理如果一条航迹所有点重合s(end)接近 0插值会报错或者产生 NaN所以先判一下直接复制首点。这种退化航迹在港口附近的数据里并不罕见。3. 距离度量改造改进 Hausdorff 到底改了什么这是整个工程的核心。DBSCAN 只是一个「按密度找簇」的框架它对「密度」的全部理解都来自那个距离函数。距离函数换了聚类结果就是另一套东西。3.1 经典 Hausdorff 的缺陷一个离群点毁掉整对航迹经典定义很简洁两条点集 A、B 之间的 Hausdorff 距离等于双向「最大最小距离」的较大者。H(A,B) max( h(A,B), h(B,A) ) h(A,B) max_{a∈A} min_{b∈B} ||a - b||直白说就是「A 里最孤单的那个点到 B 的最近距离」。问题立刻暴露取的是 max只要有一个点偏出去整条航迹的距离就被这个点绑架。AIS 数据里这种点遍地都是——GPS 跳变、靠泊时的小幅漂移、接收机解算异常产生的「飞点」。一条航迹只要混进一个离岸两公里的坏点它和任何航迹的 Hausdorff 距离都不会小于两公里聚类时必然被踢成噪声。第二个缺陷是不含方向信息。两条平行反向的航迹一去一回点集在空间上高度重叠经典 Hausdorff 距离接近 0会被聚成一类但业务上它们是相反方向的航次聚在一起毫无意义。第三个缺陷是对长度差异不敏感一条 100 公里的长航迹和它开头 5 公里的短航迹Hausdorff 距离可能很小因为短航迹的点全都落在长航迹附近。这三个缺陷在点云配准里影响不大但放到航迹聚类上每一条都会直接毁掉结果。3.2 改进方案分位数截断 平均项 方向惩罚改进思路不复杂就是把「一个最坏点说了算」改成「整体形状说了算但保留对局部差异的敏感」。第一项是分位数截断。把max换成quantile(·, q)q 取 0.9 或者 0.95。这样最远的那 5%~10% 点对不参与决策飞点被自然过滤同时保留了「大部分点都要靠近」的约束不会像平均距离那样被整体拉平。第二项是平均距离项。分位数虽然抗噪但对「大部分点都稍微偏一点」这种情况不够敏感加上双向平均距离可以补充整体贴合度。第三项是方向惩罚。用首尾点连线的主方向夹角余弦把方向差异折算成一个和距离同量纲的惩罚项乘上平均距离作为尺度因子保证不同长度的航迹之间可比。三项加权求和权重放在params.m里。默认给 0.5 / 0.4 / 0.1这个配比的意思是「形状主导、整体贴合次之、方向做微调」。如果业务上方向非常重要比如分上行下行航道把方向权重提到 0.2~0.3同时把 eps 相应放大。function d hausdorffImproved(A, B, cfg) % A, B : N×2 / M×2 的平面米制坐标已重采样通常等长 % cfg.q 分位数截掉最远的一批点对默认 0.9 % cfg.wMax 截断 Hausdorff 项权重 % cfg.wMean 双向平均距离项权重 % cfg.wDir 方向惩罚项权重 if nargin 3, cfg params(); end D pdist2(A, B); % N×M 点对距离只算一次 minAB min(D, [], 2); % 每个 a 到 B 的最近距离 minBA min(D, [], 1); % 每个 b 到 A 的最近距离 % 项 1分位数截断的双向 Hausdorff hMax max(quantile(minAB, cfg.q), quantile(minBA, cfg.q)); % 项 2双向平均最近距离 mAvg (mean(minAB) mean(minBA)) / 2; % 项 3主方向夹角用向量叉积/点积算避免 atan2 的分支问题 vA A(end,:) - A(1,:); vB B(end,:) - B(1,:); cosT dot(vA, vB) / (norm(vA) * norm(vB) eps); dirPen 1 - cosT; % 0 同向2 反向 % 方向项乘上尺度因子保证量纲一致 d cfg.wMax * hMax cfg.wMean * mAvg cfg.wDir * dirPen * max(mAvg, eps); endpdist2只调一次min(D,[],1)直接复用同一个矩阵的转置语义省掉一次 N×M 的重算。对 50 个点的重采样结果来说单次距离计算是 2500 次浮点运算一万条航迹两两比较是五千万次matlab 里跑几分钟量级可以接受。3.3 距离矩阵预计算与对称性校验DBSCAN 需要在任意两点之间查距离不可能每次现算。所以工程里先一次性算出 N×N 的上三角再对称填充。这里有个必须做的自检改进距离因为分位数取值和方向项的存在理论上d(A,B)和d(B,A)应该相等但实现里只要有一处取整或者索引写反就会不对称而 DBSCAN 是默认「距离对称」的不对称的矩阵会让聚类结果随点的遍历顺序剧烈变化表现为「跑两次结果差很多」很难查。function D buildDistMatrix(tracks, cfg) n numel(tracks); D zeros(n, n, single); % 单精度省一半内存 for i 1:n-1 for j i1:n d hausdorffImproved(tracks{i}, tracks{j}, cfg); D(i,j) d; D(j,i) d; % 强制对称不留隐患 end end % 对称性自检随机抽 20 对验证实现级别的一致性 idx randperm(n, min(20, n)); maxAsym 0; for i idx for j idx d1 hausdorffImproved(tracks{i}, tracks{j}, cfg); d2 hausdorffImproved(tracks{j}, tracks{i}, cfg); maxAsym max(maxAsym, abs(d1 - d2)); end end fprintf(距离矩阵最大非对称偏差%.6f 米\n, maxAsym); endsingle是有意为之。距离的量级在米单精度有效数字 7 位对几十公里尺度的航迹来说精度到厘米级完全够用内存直接砍半。那个自检循环只跑 400 对几秒钟的事但能挡住后面 90% 的玄学问题。4. DBSCAN 核心改造用距离矩阵替换欧氏邻域距离有了聚类主干反而是最简单的部分——只要把「邻域」的定义从「坐标差小于 eps」换成「查距离矩阵第 i 行小于 eps 的列」其余逻辑和标准 DBSCAN 一模一样。4.1 算法主干核心点判定、簇扩展与边界点归属主干分三段先一次性把每个点的 eps 邻域算出来存成 cell同时标记核心点然后遍历未访问点核心点开新簇并做 BFS 扩展扩展过程中遇到已标记为噪声的点把它「回收」成边界点。function [labels, isCore] dbscanByDist(D, eps, minPts) % D : N×N 对称距离矩阵 % eps : 邻域半径单位与 D 一致米 % labels : N×10 未访问-1 噪声正整数为簇编号 n size(D, 1); labels zeros(n, 1); isCore false(n, 1); neighbors cell(n, 1); % 第一遍邻域预计算 核心点判定 for i 1:n neighbors{i} find(D(i, :) eps); isCore(i) numel(neighbors{i}) minPts; end clusterId 0; for i 1:n if labels(i) ~ 0, continue; end if ~isCore(i) labels(i) -1; % 先当噪声扩展时可能被回收 continue; end clusterId clusterId 1; labels(i) clusterId; queue neighbors{i}; head 1; while head numel(queue) j queue(head); head head 1; if labels(j) -1 labels(j) clusterId; % 边界点回收 continue; end if labels(j) ~ 0, continue; end labels(j) clusterId; if isCore(j) queue [queue; neighbors{j}]; % 核心点才继续向外扩 end end end end参数只有两个但含义和点聚类完全不同。eps是「两条航迹算不算邻居」的阈值单位是米经验值取「同一航道内不同航次之间的典型距离」通常 500~2000 米取决于航道宽度和航迹尺度。minPts是「一个簇至少要有多密」船舶航迹是有明显主航道的同一个 OD 对的航次重复率高minPts 取 3~5 就够取大了会把冷门航线全判噪声。if labels(j) -1那个分支的顺序不能反。噪声点被回收后必须continue否则会继续走「未访问就标记并扩展」的逻辑等于把一个边界点当核心点用簇会长出不该有的尾巴。4.2 参数整定k-距离曲线怎么读eps靠猜是最容易翻车的地方。标准做法是画 k-距离曲线对每个点取它第 k 近的距离k 一般取 minPts-1降序排列画出来曲线从陡变缓的拐点就是建议 eps。原理是拐点左边是「簇内的点」右边是「离群点」落在拐点附近作为半径能把密集区包住、把稀疏区排除。function epsHat kdistCurve(D, minPts, doPlot) Dsort sort(D, 2); % 每行升序第 1 列是自身距离 0 k max(minPts - 1, 1); kdist sort(Dsort(:, k 1), descend); if doPlot figure; plot(kdist, LineWidth, 1.2); grid on; xlabel(航迹序号按第 k 近邻距离降序); ylabel(sprintf(第 %d 近邻距离 / 米, k)); title(k-距离曲线拐点处即建议 eps); end % 用最大曲率点近似拐点二阶差分绝对值最大处 d2 abs(diff(kdist, 2)); [~, pos] max(d2); epsHat kdist(pos); fprintf(建议 eps ≈ %.1f 米minPts%d\n, epsHat, minPts); end自动找拐点只是给个起点不能全信。实际调的时候我会把这个建议值上下浮动 20% 各跑一次看噪声率和簇数的变化曲线eps 稍微放大噪声率就骤降、簇数骤减说明已经在大簇合并得往回收eps 缩小噪声率暴涨说明簇开始碎裂。稳定区间通常就在那个平台段。4.3 结果评估噪声率、簇数与只对非噪声点算的轮廓系数航迹聚类没有标签评估只能靠内部指标加人工看图。工程里跑三个数噪声率-1 的占比、簇数、轮廓系数。轮廓系数要特别注意必须排除噪声点再算因为噪声点本来就不属于任何簇把它们算进去会让系数虚低误判成「聚类效果差」。function m clusterMetrics(D, labels, trackLens) noiseRate mean(labels -1); nCluster numel(unique(labels(labels 0))); idx labels 0; if nCluster 2 m struct(noiseRate, noiseRate, nCluster, nCluster, silhouette, NaN); return; end Dv D(idx, idx); L labels(idx); s silhouetteFromDist(Dv, L); % 距离矩阵版的轮廓系数 % 顺便统计每个簇的平均航程太长/太短的簇要人工复核 meanLen zeros(nCluster, 1); for c 1:nCluster meanLen(c) mean(trackLens(labels c)); end m struct(noiseRate, noiseRate, nCluster, nCluster, ... silhouette, mean(s), meanLen, meanLen); end指标之外一定要画图。按簇编号取一组区分度高的颜色把每条航迹画成折线叠在底图上噪声点单独用灰色。这一步是唯一能看出「簇是不是业务上合理」的手段如果某一个簇把上行和下行的航迹混在一起指标再漂亮也得回去调方向权重。某实验室的同学就是因为只看轮廓系数没看图交上去的结果里两个明显不同的航道被合并成一类返工了整整两天。5. 避坑记录航迹聚类跑不通的五个典型现场这一节全部来自实际调试。现象、原因、解决三件事写清楚照着对号入座基本能覆盖八成问题。5.1 现象所有航迹聚成一个簇或者几乎全是噪声原因几乎总是eps的单位搞错了。前处理阶段如果把经纬度直接丢进pdist2没做投影距离的量级是「度」0.01 度大约是 1.1 公里可参数表里填的 eps 是 800米于是所有点对距离都远大于 800全部判噪声反过来如果 eps 填了 0.5那 0.5 度以内的点全成邻居整片海域合成一个簇。解决是把投影强制放在距离函数入口作为第一行代码而不是靠调用方自觉。function XY latlon2xy(lat, lon, lat0, lon0) % 以 (lat0, lon0) 为原点的局部等距投影输出米制平面坐标 R 6371000; x deg2rad(lon - lon0) .* R .* cosd(lat0); y deg2rad(lat - lat0) .* R; XY [x, y]; end参考点取整片数据的中心经纬度。研究区域跨度在一两百公里以内等距近似的误差在千分之一量级对聚类阈值的影响可以忽略跨省甚至跨洋的数据要换墨卡托或者按纬度分带处理否则高纬度的东西向距离会被压缩得厉害。5.2 现象跑到一半报Out of memory或者机器卡死原因是对称矩阵乘法式的暴力预计算。N15000 的double距离矩阵是 1.8GBmatlab 在算上三角赋值的时候还会临时分配峰值翻倍普通办公机直接跪。解决办法有三条按顺序试第一矩阵用single内存砍半对米级距离精度无损第二分块计算每次处理 2000 条块内算完就写盘或者写回全局矩阵避免同时驻留两份大矩阵第三先做粗筛用低维特征起点、终点、平均航向、航程快速排除明显不可能同簇的组合只对候选对算精细距离。第三条在航迹数量上万时效果最明显能把候选对压掉七成以上。注意zeros(n, n, single)只是声明单精度pdist2的返回值仍然是 double赋值得靠自动转换。如果中间再做D D这类运算会临时升回 double内存还是会爆。5.3 现象同一条数据跑两次聚类结果不一样DBSCAN 本身是确定性的前提是点和邻域的处理顺序固定。结果不稳定通常来自两个地方一是航迹集合的顺序如果前面用unique或者findgroups之后没有排序两次读盘的顺序可能不同取决于文件系统的返回顺序边界点归属就会变二是距离矩阵不对称d(i,j)和d(j,i)差一点点某些点对恰好卡在 eps 边界上这次算邻居下次不算。解决办法是双管齐下航迹入库前按 MMSI 加起始时间做一次确定性排序保证每次跑的顺序完全一致距离矩阵构建时强制D(j,i) D(i,j)并且加 3.3 节那段自检。如果业务上允许还可以加一步「多次扰动投票」对航迹索引做几次随机置换各跑一遍用共现矩阵决定最终簇归属能明显缓解边界点的抖动代价是耗时成倍。5.4 现象簇是聚出来了但明显上行和下行的航迹混在一起方向项没起作用。先去查dirPen的尺度因子如果mAvg很小比如两条航迹高度重合平均距离只有几十米方向惩罚乘以这个因子之后几乎为零等于没加。再查首尾点连线的主方向如果航迹有回头弯首尾点可能几乎重合vA的模接近 0cosT变成噪声。解决是方向项不用首尾点改用整个序列的加权主方向或者用相邻点方向向量的平均。function u meanHeading(A) % 用相邻点方向向量的单位化和避免首尾点重合导致的退化 v diff(A); v v ./ (vecnorm(v, 2, 2) eps); u sum(v, 1); u u / (norm(u) eps); end同时把方向惩罚的尺度因子从mAvg换成max(mAvg, cfg.minScale)给一个下限防止平均距离过小时方向项被压没。minScale取 eps 的十分之一左右比较合适。5.5 现象短航迹全被判成噪声簇的数量比预期少很多原因是距离度量的尺度不匹配。一条 30 个点、5 公里的短航迹和一条 50 个点、100 公里的长航迹算距离即使它们起点重合、方向一致长航迹的中后段离短航迹很远平均距离项被拉得很大直接超过 eps。解决是在前处理里做长度分层按航程把所有航迹分成短、中、长三档每档单独聚类各自用适合的 eps。另一种做法是在距离计算前对两条航迹做局部对齐——只比较重叠部分也就是取长航迹中距短航迹首尾最近的那一段参与计算。分层简单可靠工程上我更推荐先分层效果不好再上对齐。% 按航程分层各层独立跑聚类 edges [0, 20, 100, inf]; % 公里 for b 1:numel(edges)-1 sel trackLenKm edges(b) trackLenKm edges(b1); if nnz(sel) 10, continue; end D buildDistMatrix(tracks(sel), cfg); [lab, ~] dbscanByDist(D, cfg.epsByBand(b), cfg.minPts); % ... 保存该层结果 end分层之后每层的 eps 独立整定短航迹层用小 eps长航迹层用大 eps噪声率会立刻降下来。6. 二次开发把簇结果接到航迹预测与批量跑数聚类本身不是终点业务上真正要用的是「给定一条新航迹判断它属于哪条习惯航线然后照着这条航线的历史走势做预测」。所以二次开发的第一件事是把簇转成可检索的代表航迹。6.1 从簇到预测代表航迹与最近邻匹配因为前面已经统一重采样到 50 个点簇内求代表航迹的直接做法是逐点取中位数——比取平均更抗异常航次。中位数航迹相当于这个簇的「典型航线」后续预测就拿它做模板。function refs buildRefTracks(tracks, labels, cfg) nCluster max(labels); refs struct(id, {}, track, {}, count, {}); for c 1:nCluster members tracks(labels c); if numel(members) cfg.minPts, continue; end % 所有成员已在预处理阶段重采样到 cfg.resampleN 点 P cell2mat(cellfun((x) reshape(x, 1, [], 2), members, UniformOutput, false)); ref squeeze(median(P, 1)); % 逐点中位数 refs(end1) struct(id, c, track, ref, count, numel(members)); %#okAGROW end end匹配新航迹就是算它到所有代表航迹的改进 Hausdorff 距离取最小且小于cfg.eps的那个簇都超过阈值就判为未知航线走异常告警。这一步可以直接复用第 3 章的距离函数不用再写一套。改造时最容易忽略的是「代表航迹也要重采样」。如果新航迹是实时增量数据点数跟 50 对不上pdist2出来的矩阵维度没问题Hausdorff 允许不等长但平均距离项会被采样密度影响判断结果会偏。所以在匹配入口再做一次重采样成本极低。6.2 封装成可复用函数与批量跑数最后一步是把整条链路包成一个函数参数全部走params()返回一个 table方便直接写 csv 交给下游。function result runCluster(csvFile, cfg) tracks loadAIS(csvFile, cfg); tracks cellfun((t) resampleTrack(t, cfg.resampleN), tracks, UniformOutput, false); D buildDistMatrix(tracks, cfg); [labels, isCore] dbscanByDist(D, cfg.eps, cfg.minPts); m clusterMetrics(D, labels, cellfun((t) arclength(t), tracks)); result table((1:numel(tracks)), labels, isCore, ... VariableNames, {TrackId, ClusterId, IsCore}); writetable(result, fullfile(cfg.outDir, cluster_result.csv)); fprintf(簇数 %d噪声率 %.2f%%轮廓系数 %.3f\n, ... m.nCluster, m.noiseRate*100, m.silhouette); end批量跑数的时候建议把cfg存一份到输出目录和结果一起归档。我吃过亏一个月后回看某次实验结果指标很漂亮但完全想不起当时的 eps 和分位数是多少只能整批重跑。从那以后我每次跑聚类都强制把参数快照写盘一个save(fullfile(outDir,cfg.mat),cfg)就能解决成本几乎为零。改造别人的代码时也一样先在入口处把参数存档再去动算法出了问题至少有后悔药。希望帮到你。本文还有配套的精品资源点击获取