简介针对基于改进Hausdorff距离的DBSCAN船舶航迹聚类与异常行为识别问题这份Matlab源码项目提供了完整可运行实现适合从事船舶航迹分析、聚类算法应用或异常检测研究的学生和工程师。项目完整复现论文《基于轨迹聚类的船舶异常行为识别研究》的核心流程覆盖航迹数据提取、DBSCAN聚类、聚类中心提取、Hausdorff距离计算、预测阈值寻优及偏离预警等环节。资源共20个文件以14个.m脚本和1个.mat数据为核心另含2个结果/画图zip压缩包、2个png说明/联系方式图片及1份md说明文档整体仅4.32MB可直接运行并根据自身数据替换模块。已有765人学习下载适合希望以真实案例掌握DBSCAN变体距离计算和航迹聚类流程的读者。借助该案例不仅可理清航迹分类与阈值寻优的关键思路还能基于聚类结果做航迹偏离预测扩展为个人异常行为识别模型实用性强。1. 为什么普通DBSCAN聚类船舶航迹总是一团糟Harsdorf距离的引入逻辑做船舶航迹聚类时最让人头疼的不是算法不会写而是“看起来在一条航线上的船算出来居然不是一类”。AIS数据里每条轨迹是长短不一的点序列有的船一小时报一次位置有的五秒报一次同一条航线上船可能因为避让、风流影响偏离航道几百米。直接把坐标点喂给DBSCAN要么因为轨迹长度不同导致距离计算失败要么因为局部抖动太大把本该聚成一类的轨迹拆成好几块。这就是标题里“Matlab基于改进的Harsdorf距离的DBSCAN船舶航迹聚类”要解决的核心问题先用Harsdorf距离规范术语是Hausdorff距离度量两条轨迹的整体形状差异再把距离矩阵交给DBSCAN做密度聚类。这个方案适合手里有AIS或船舶轨迹数据、想做航线挖掘、港口交通流分析或异常行为检测的从业者也适合刚接触轨迹聚类、想避开“点级聚类”那些坑的初学者。它的价值在于把“轨迹相似度怎么算”和“聚类怎么跑”两件事彻底拆开每一层都能独立调优。2. 航迹相似度怎么算从标准Hausdorff到改进Hausdorff距离2.1 标准Hausdorff距离的定义和它的病态问题Hausdorff距离度量的是两个点集之间的最大不匹配程度。给定两条轨迹A和B标准双向Hausdorff距离定义为H(A, B) max(h(A, B), h(B, A))其中h(A, B) max(min(||a - b||))表示A中每个点到B中最近点的距离的最大值。直观理解就是两条轨迹里最离谱的那个点决定了它们之间的距离。这个性质在轨迹聚类里既是优点也是灾难。优点是它对轨迹形状的整体差异非常敏感两条轨迹只要主体走向不同距离就会很大缺点是它对单个离群点毫无抵抗力。实际AIS数据里一条轨迹上只要有一个定位跳变点或者靠港绕行点哪怕其余几百个点和另一条轨迹完全重合计算出来的Hausdorff距离也会被这个坏点拉到极大值聚类时这条船就会被误判成不属于任何航线。我在处理长江口AIS数据时踩过这个坑。两条轨迹从吴淞口到绿华山锚地前半段几乎平行其中一条因为避让渔船绕了一个半径约1.5海里的弯结果标准Hausdorff距离从正常的0.3海里飙到1.8海里DBSCAN直接把它们分成两簇。后来我把绕行段之外的点全部参与计算发现98%的点的最近距离都小于0.1海里。这说明标准定义把“最大误差”当成“整体差异”对轨迹聚类这种本身就含噪的数据来说过于苛刻。2.2 改进Hausdorff距离的三种常见形式要解决标准定义的病态问题最常见的做法是把“最大值”换成“平均值”或“分位数”。我在实际项目中主要用以下三种改进形式它们各有适用场景。第一种是平均Hausdorff距离Modified Hausdorff Distance, MHD。把h(A, B)从max改成对A中所有点取平均即h_mhd(A, B) mean(min(||a - b||))。双向同样取两个方向的均值。这样单个离群点的影响被稀释到全部点里轨迹整体越接近距离越小。MHD的缺点是它对局部密度的差异不敏感一条轨迹在某个区域点特别密、另一条点特别稀时平均值会被密集段的重复计算带偏。第二种是截断Hausdorff距离Quantile/Truncated Hausdorff Distance。先把A中每个点到B的最近距离排序然后取第90%或95%分位数而不是最大值。这样做既保留了“尾部惩罚”的思想又不会被极端离群点完全控制。分位数P的选取要看数据噪声比例我一般先用箱线图看一下最近距离分布的异常点比例再决定P取0.9还是0.95。如果异常点比例超过5%取0.9更稳。第三种是带时间权重的改进Hausdorff。船舶轨迹不只是空间路径还有时间属性。两条船在同一片海域但时间差几个小时的轨迹业务上往往不算同一批航次。做法是把每个点的最近距离乘一个时间衰减因子比如exp(-|Δt|/T)T是时间尺度参数单位小时。这样空间接近但时间错开的轨迹会被拉开距离。这个改进在港口交通流分析里特别有用能区分“同向但不同时段的班轮”和“同一时段的编队”。2.3 Matlab里实现改进Hausdorff距离重采样是前提写改进Hausdorff距离之前必须先处理两条轨迹长度不一致的问题。虽然Hausdorff距离定义上不需要轨迹等长但为了后续能复用pdist2这类向量化计算我习惯先把每条轨迹重采样到固定点数N。N的选择很关键太小会丢失航道弯曲特征太大则增加计算量。我一般按轨迹长度和航速来定N取64或128能保留大部分转向特征。% 轨迹重采样到N个点保持原始形状 function traj_resampled resampleTrajectory(traj, N) % traj: n*2矩阵[经度 纬度]已转平面坐标 % 计算累计弧长 dx diff(traj(:,1)); dy diff(traj(:,2)); seg_len sqrt(dx.^2 dy.^2); cum_len [0; cumsum(seg_len)]; total_len cum_len(end); % 按等弧长生成新采样位置 query linspace(0, total_len, N); traj_resampled(:,1) interp1(cum_len, traj(:,1), query, linear); traj_resampled(:,2) interp1(cum_len, traj(:,2), query, linear); % 处理重复点导致cum_len非严格递增的情况 if any(diff(cum_len) 0) traj_resampled traj_resampled(~isnan(traj_resampled(:,1)), :); end end这段代码的思路是先算轨迹累计弧长再按弧长等间隔采样。用interp1插值时如果原始轨迹里有原地停留的重复点AIS数据里很常见船速为0但坐标不变cum_len会出现相等值interp1会报错或者插出NaN。所以我在最后加了一行防呆处理把NaN行过滤掉。实际使用中我还会加一个判断如果traj点数小于N直接用原始点加末尾点重复补全避免插值外推。重采样之后改进Hausdorff距离函数就简单了。我写了一个支持三种模式的函数方便切换对比效果。% 计算两条轨迹的改进Hausdorff距离 function d mhd_distance(A, B, method, pct) % A, B: N*2矩阵重采样后的轨迹 % method: std标准双向, mhd平均, quantile截断 % pct: quantile模式的分位数如0.9 D pdist2(A, B); % N*N欧氏距离矩阵 h_AB min(D, [], 2); % A中每个点到B的最近距离 h_BA min(D, [], 1); % B中每个点到A的最近距离 switch method case std d max(max(h_AB), max(h_BA)); case mhd d (mean(h_AB) mean(h_BA)) / 2; case quantile d (quantile(h_AB, pct) quantile(h_BA, pct)) / 2; otherwise error(未知方法); end endpdist2计算两条轨迹所有点对之间的欧氏距离然后分别沿行和列取最小值得到两个方向的最近距离序列。标准模式用max取最大值平均模式用mean取均值截断模式用quantile取分位数。这里要特别注意pdist2默认用欧氏距离但船舶轨迹坐标如果直接用电离度经纬度一个纬度差对应的海里数和经度差完全不同算出来的距离在南北方向和东西方向不可比。所以坐标必须先转成平面坐标我一般用墨卡托投影的本地近似即把经度乘以cos(平均纬度)做缩放单位统一成海里。这一步不在函数里做而是数据预处理时完成。2.4 四种距离定义的效果对比和选型建议距离定义对离群点敏感度对局部密度差异敏感度计算开销适用场景标准Hausdorff极高低低形状差异极大的轨迹分类平均MHD低中低AIS轨迹聚类推荐首选截断Quantile中可控低低噪声点比例较高的数据时间加权MHD低中中需要区分时间维度的交通流分析我在实际项目里的选择逻辑是如果只是做航线挖掘不考虑时间维度直接用平均MHD如果数据里有比较多的异常轨迹比如定位漂移、靠港绕行用截断Quantile配合P0.9如果业务上要区分不同批次的编队航行用时间加权MHD。选型时还要考虑计算量平均MHD和截断Quantile的耗时代价几乎一样时间加权MHD多一个时间差计算轨迹数量超过一万条时会明显变慢。另外提醒一点无论选哪种改进形式重采样的点数N要固定否则两条轨迹点数不同pdist2算出来的距离矩阵不是方阵虽然不影响min操作但后续如果要对距离做缓存或并行计算统一N会省很多麻烦。3. 把Hausdorff距离矩阵喂给DBSCANMatlab最小可跑通实现3.1 为什么不能直接调用dbscan函数Matlab从R2019b开始提供内置的dbscan函数但它接受的是特征矩阵内部用欧氏距离找邻域不适合直接接自定义距离。虽然dbscan函数支持Distance参数可以传precomputed并输入距离矩阵但它的距离矩阵要求是方阵且对角线为0这对标准Hausdorff没问题对改进Hausdorff也没问题——只要我们自己保证对称性。不过实际用下来我发现内置版本在距离矩阵很大超过几千条轨迹时占内存比较凶而且它内部用KDTree加速的逻辑对自定义距离不生效。所以我的习惯是轨迹数量在5000条以内用自己写的主循环超过5000条再考虑内置版本或者分块计算。自己写DBSCAN主循环的另一个好处是可以方便地输出每个点的类型核心点、边界点、噪声点这在排查聚类效果时非常有用。内置dbscan返回的聚类标签里噪声点标为0但核心点和边界点不区分排查参数时少了一层信息。% 基于距离矩阵的DBSCAN主循环 function [labels, core_flags] dbscan_dist(D, eps, minPts) % D: n*n距离矩阵需满足D(i,j) D(j,i)对角线为0 % eps: 邻域半径 % minPts: 邻域内最少点数 n size(D, 1); labels zeros(n, 1); % 0表示未分类-1表示噪声 core_flags false(n, 1); % 是否核心点 cluster_id 0; % 先统计每个点的邻域点数判断是否核心点 for i 1:n neighbor_count sum(D(i,:) eps) - 1; % 减1去掉自身 if neighbor_count minPts core_flags(i) true; end end % 核心点扩张聚类 for i 1:n if labels(i) ~ 0 || ~core_flags(i) continue; end cluster_id cluster_id 1; labels(i) cluster_id; seed find(core_flags D(i,:) eps labels 0); seed seed(:); idx 1; while idx length(seed) cur seed(idx); if labels(cur) 0 labels(cur) cluster_id; end if core_flags(cur) new_neighbors find(D(cur,:) eps labels 0); seed [seed; new_neighbors(:)]; %#okAGROW end idx idx 1; end end % 剩余未访问点标记为噪声 labels(labels 0) -1; end这段代码的关键在于先用一次全量扫描确定核心点再用核心点做BFS扩张。注意第二个循环里我只从核心点开始聚类边界点邻域点数不够但不是噪声会在扩张过程中被吸收。seed队列的写法用了一个动态增长的向量Matlab里这会被效率警告但如果轨迹数量在几千条级别实际耗时可接受。如果想提速可以预分配seed数组或者用队列结构体但代价是代码变复杂。我在项目中还发现一个细节统计邻域点数时用sum(D(i,:) eps) - 1减1是为了去掉自身D(i,i)0必然满足条件。如果距离矩阵没有严格保证对角线为0这里排查半天都找不出问题所以构造距离矩阵时务必确认对角线。3.2 完整的聚类流程从AIS原始轨迹到聚类结果整个流程分四步坐标转换和轨迹清洗、重采样、计算距离矩阵、DBSCAN聚类。我通常把前两步合成一个预处理脚本后两步合成一个聚类脚本中间用mat文件传递数据这样调试时不用重复跑耗时的距离计算。% 主脚本跑通完整聚类流程 % 输入tracks是元胞数组每个元素是n*2矩阵 [lng, lat] % 输出labels, 每个元胞对应的聚类标签 % Step1: 坐标转换为平面坐标本地墨卡托近似 lat0 mean(cellfun((t) mean(t(:,2)), tracks)); lng_scale cosd(lat0); for k 1:length(tracks) tracks{k}(:,1) tracks{k}(:,1) * lng_scale * 60; % 经度转海里 tracks{k}(:,2) tracks{k}(:,2) * 60; % 纬度转海里 end % Step2: 重采样到固定长度N64 N 64; tracks_rs cell(size(tracks)); for k 1:length(tracks) tracks_rs{k} resampleTrajectory(tracks{k}, N); end % Step3: 计算改进Hausdorff距离矩阵 n_tracks length(tracks_rs); D zeros(n_tracks, n_tracks); for i 1:n_tracks for j i1:n_tracks d mhd_distance(tracks_rs{i}, tracks_rs{j}, mhd, 0.9); D(i,j) d; D(j,i) d; end end % Step4: 用距离矩阵跑DBSCAN eps 2.0; % 单位海里根据距离分布设定 minPts 5; [labels, core_flags] dbscan_dist(D, eps, minPts);坐标转换里的lng_scale cosd(lat0)把经度差乘以纬度的余弦再统一乘以60把度转成海里。这个近似在纬度跨度小于200海里的区域内误差可以忽略如果研究区域跨好几个纬度带比如中国沿海从北到南建议分区域分别处理或者用更严格的墨卡托投影函数。重采样和距离计算分别循环复杂度是O(n_tracks^2 * N^2)三千条轨迹、N64时pdist2每调用一次要算4096个点对距离整体耗时约几分钟属于可以接受的范围。如果把N减到32耗时能缩短到四分之一代价是曲线细节丢失拐弯急的航段可能被拉直。3.3 这里有几个必须同步调整的参数第一个是重采样点数N。我见过有人直接取所有轨迹点数的中位数结果长短轨迹差异大的数据里长轨迹被严重压缩短轨迹被过度插值。建议N取轨迹点数分布的50%分位数附近然后统一。第二个是距离矩阵的存储精度我默认用double但如果轨迹数量超过一万条D矩阵占内存超过800MB可以考虑临时转single存储计算时再转回double精度损失对聚类结果影响很小。第三个是eps的初始值不要拍脑袋。我一般是先跑一次距离矩阵把D上三角不包括对角线拉平画出直方图和累计分布看距离分布的自然分层点在哪里那个点就是eps的首选值。这个方法比按比例取中位数可靠因为距离分布通常是多峰的不同峰对应不同尺度的轨迹差异强行取中位数会把多尺度结构压平。% 距离分布分析确定eps的参考范围 triu_idx find(triu(ones(n_tracks), 1)); dist_vals D(triu_idx); % 画累计分布曲线 [f, x] ecdf(dist_vals); plot(x, f); grid on; xlabel(Hausdorff距离海里); ylabel(累计概率); % 找斜率突变点即曲线从陡变缓的位置 % 视觉检查后取累计概率0.1到0.3之间的平缓段起始值作为eps候选ecdf画出来的曲线如果有一段明显变缓说明轨迹距离存在一个“类内差异”和“类间差异”的分界带。eps取在变缓段的起始值附近聚类结果最稳定。如果曲线全程平滑没有明显拐点说明轨迹差异是连续谱这时候改距离定义比如从MHD换成截断Quantile比调eps更有效。3.4 pdist2和暴力循环之间怎么选上面用的双层循环调mhd_distance每个距离都要算一次pdist2重复计算了每条轨迹和其他轨迹的距离过程中间结果。另一个做法是预先算好所有重采样后轨迹的点坐标拼接成一个大矩阵用pdist2一次算出所有轨迹点之间的欧氏距离再通过reshape和min操作提取最近距离。% 向量化方案一次pdist2批量计算所有轨迹间MHD距离 % 把所有轨迹点拼成 n_tracks*N 行的大矩阵 all_points zeros(n_tracks * N, 2); for k 1:n_tracks all_points((k-1)*N1 : k*N, :) tracks_rs{k}; end D_all pdist2(all_points, all_points); % 提取每条轨迹点块之间的最近距离 block_idx reshape(1:n_tracks*N, N, n_tracks); % 初始化距离矩阵 D_vec zeros(n_tracks, n_tracks); for i 1:n_tracks for j i1:n_tracks blockD D_all((i-1)*N1:i*N, (j-1)*N1:j*N); h_ij min(blockD, [], 2); h_ji min(blockD, [], 1); D_vec(i,j) (mean(h_ij) mean(h_ji)) / 2; D_vec(j,i) D_vec(i,j); end end这个方案的优点是一次pdist2把整个大矩阵的欧氏距离算完之后的操作都是矩阵切片在轨迹数量超过两千条时速度明显优于双层循环调函数。缺点是大矩阵D_all是(n_tracks*N)的方阵三千条轨迹、N64时就是19.2万行D_all要占约590GB内存根本放不下。所以实际中这个向量化方案只适合小规模数据几百条轨迹以内或者分块计算把轨迹分成若干批每批之间算距离再合并。我在项目里一般用分块的思路每批包含500条轨迹内部分块pdist2两个批次之间单独算交叉距离这样内存和速度都能兼顾。分块还有一个额外好处是方便用parfor并行每个worker只需要加载自己负责的距离块。4. 让DBSCAN参数不再靠猜eps和minPts的标定方法与自适应流程4.1 为什么船舶航迹聚类不能用默认参数很多教程里默认minPts 2 * 维度对二维数据就是minPts4eps取k-dist曲线的拐点。这在点云聚类里问题不大但轨迹聚类完全不一样。我们聚类的基本单元是“整条轨迹”不是“单个点”。一条轨迹代表一个样本样本之间距离的尺度和坐标系下的物理距离直接相关但密度分布却受航路分布影响极大。繁忙航道里同向轨迹的Hausdorff距离可能集中在0.5到1.5海里开阔海域里散货船的随意航行轨迹距离可能分布在5到20海里。用一个全局eps处理整个研究海域结果往往是繁忙区域聚成一坨稀疏区域全部变成噪声点。我在处理某沿海港口进出港轨迹时遇到过同一个eps1.5海里在进港航道的分道通航带里能聚出漂亮的3条交通流但到了港外锚地散乱的等待锚泊轨迹全部被判成噪声锚地明明是有聚集特征的。后来我把研究区域划分成港内、近海、锚地三个子区域分别标定eps才得到合理的簇结构。这说明轨迹聚类的eps不是算法参数而是业务参数必须结合具体海域的交通密度来标定。4.2 k-dist曲线在轨迹距离矩阵上的适用性k-dist曲线的做法是对每个样本找它到第k近邻的距离k minPts - 1按降序排列后画曲线。曲线拐点对应的距离值就是eps候选。这个方法对点云有效对距离矩阵同样有效但有个前提距离矩阵必须是真实反映样本密度的。如果用的距离定义本身对离群点敏感比如标准Hausdorffk-dist曲线会非常平滑拐点模糊根本找不到稳定的eps。所以我都是先确定距离定义用改进版再画k-dist图。% 画k-dist图来标定eps function [kdist, sorted_idx] plot_kdist(D, minPts) k minPts - 1; n size(D, 1); kdist_all zeros(n, 1); for i 1:n d_row D(i,:); d_row(i) inf; % 排除自身 sorted_d sort(d_row); kdist_all(i) sorted_d(k); end [kdist, sorted_idx] sort(kdist_all, descend); plot(1:n, kdist, .-); xlabel(样本序号按k-dist降序); ylabel(第k近邻距离海里); end画完曲线后理想情况是看到一个明显的“膝盖”形状左侧一段陡降然后突然变平缓。拐点对应的y值就是eps。这里要提醒k minPts - 1所以minPts取5时看的是第4近邻距离。minPts越大曲线越平滑拐点越不明显minPts越小曲线越毛糙拐点越锐利。我一般先用minPts5画出线目测拐点后再用minPts3和minPts7各画一次如果三次标出的eps差别不大在20%以内说明这个eps稳定可信。如果差别很大说明数据本身的簇结构不明显此时要回头检查距离定义而不是硬调参数。4.3 分位点法比目测拐点更可复现的eps设定k-dist曲线的目测拐点有一个问题不同人看到的位置不一样代码评审时说不清。更可复现的做法是直接取k-dist值的某个分位数。比如把k-dist序列画完后取它的90%分位数、95%分位数作为eps候选然后对比聚类结果。% 用分位点设定eps并快速评估聚类规模 eps_candidates [quantile(kdist, 0.80), quantile(kdist, 0.85), quantile(kdist, 0.90)]; for eps_test eps_candidates labels_tmp dbscan_dist(D, eps_test, 5); n_clusters length(unique(labels_tmp(labels_tmp 0))); n_noise sum(labels_tmp -1); fprintf(eps%.2f 簇数量%d 噪声点%d\n, eps_test, n_clusters, n_noise); end这个做法的逻辑是取90%分位数意味着允许10%的样本在邻域外这部分通常对应轨迹差异极大的极端数据取之后作为eps能让大多数正常轨迹先聚起来。实际经验中eps取k-dist的85%到90%分位数时簇结构最稳定。如果取的太高比如95%以上容易把几个本来独立的簇合并取太低比如70%以下会把一个簇从中间切开。用分位点法的好处是每一步都有明确计算依据复现实验时只要约定“eps kdist的p90”换数据、换人都能拿到一样的参数。4.4 分层聚类和自适应处理多密度区域的实用技巧如果研究区域确实包含密度差异明显的多个子区域全局eps怎么调都不理想那就别硬调了直接分层处理。我的做法是先做一次粗聚类把eps设成全局的80%分位数minPts设6跑一遍得到几个大簇和大量噪声点。然后对噪声点单独做一次精聚类把eps降为原来的40%minPts降到3只看噪声点之间的距离矩阵。这样两层下来繁忙区域和稀疏区域都能得到合理的簇。% 分层聚类先粗后细 eps_coarse quantile(kdist, 0.8); labels_coarse dbscan_dist(D, eps_coarse, 6); noise_idx find(labels_coarse -1); if ~isempty(noise_idx) D_noise D(noise_idx, noise_idx); kdist_noise sort(min(D_noise diag(inf), [], 2), descend); eps_fine quantile(kdist_noise, 0.8); labels_fine dbscan_dist(D_noise, eps_fine, 3); % 把精聚类标签映射回原标签空间 labels_coarse(noise_idx(labels_fine 0)) max(labels_coarse) labels_fine(labels_fine 0); end注意精聚类这一步minPts从6降到3是有讲究的。粗聚类时minPts大能把零星轨迹视为噪声保证簇的纯度精聚类时minPts小能在稀疏区域捕捉到“三五条船同路”的弱聚集模式。如果精聚类时还保持minPts6大多数小簇都达不到核心点要求精聚类等于白跑。还有一点精聚类用的eps_fine只基于噪声点的距离分布计算不要沿用全局的eps因为噪声点的距离分布和核心区域的分布本来就是两回事。分层聚类跑完后记得检查粗聚类被精聚类重新划分的簇这些簇在业务上往往是“同一航向但时间分散”的模式如果恰好是你要找的交通流特征就保留如果只是数据噪声就再把它们标回噪声点。4.5 minPts和轨迹数量、噪声比例的联动关系minPts的设定不能只看数据维度还要看轨迹总数和预期噪声比例。一个经验公式是minPts ≈ ln(n_tracks)但我觉得这太粗略。我更倾向于用业务语义来定你想让几艘船组成一个“可识别的航次模式”如果3艘船同向且距离接近就算一条交通流minPts3如果要求至少5艘才算minPts5。这和点云聚类里“核心点最少要有k个邻居”的含义完全不同轨迹聚类里minPts代表的是“成为一类航迹模式的最少样本数”。另外要留意minPts对噪声点比例的影响。minPts越小噪声点越少但簇的数量会碎片化minPts越大噪声点越多大簇更稳定。我在做港口航行安全分析时通常把噪声点比例控制在10%到20%之间低于10%说明聚类过松高于20%说明聚类过紧。这个比例可以直接作为参数调优的量化目标跑一组(minPts, eps)的组合记录每个组合的噪声点比例找最接近15%的那组参数。如果所有组合的噪声点比例都远超20%说明距离定义里混入了太多异常轨迹应该先清洗数据而不是一味加大eps。5. DBSCANHausdorff轨迹聚类的5个常见翻车点和排查记录5.1 翻车点一聚出来的“航线”是一条贯穿全图的之字折线现象聚类结果里某个簇的轨迹首尾相连在地图上形成一条从起点到终点横跨整个研究区域的之字形折线看起来像是把完全不同的航段硬接在一起。原因这类轨迹通常是同一艘船在一段时间内的完整航程中间经过了多次转向和靠泊。用Hausdorff距离衡量时它和任何一条子航段的距离都不算太大于是一个长航程轨迹作为“桥梁”把好几条短航段轨迹串进同一个簇。DBSCAN的核心点扩张机制在这种情况下会放大这个效应长航程轨迹是核心点它的邻居包括A航段和B航段于是A和B被分到一个簇。解决在距离定义里加入航向变化惩罚。具体做法是重采样时保留每个点的瞬时航向角然后在Hausdorff距离里增加一项航向差异项让转向次数、角度差大的轨迹距离拉大。我一般用d_final d_mhd λ * mean(|θ_A - θ_B|)其中θ_A是轨迹A每个点的航向θ_B是轨迹B对应点的航向用最近点匹配λ取0.5到1之间。加了这个惩罚项后长航程的之字轨迹和短航段的距离明显增大不会被错误桥接。5.2 翻车点二同一批AIS数据跑两遍每次聚类结果都不一样现象代码没有改任何参数只是重新运行了一遍脚本聚类结果图看起来差不多但簇的编号变了某些边缘轨迹一会儿归A簇一会儿归噪声点。原因Matlab内置dbscan或者自己写的循环里如果用了parfor或者在循环里依赖了随机数种子比如某些版本的pdist2在多线程模式下会有浮点累加顺序不同距离矩阵会有微小数值差异。当某条轨迹正好落在eps边界附近时微小的数值抖动就会改变它的归属。这不是算法逻辑错误而是数值稳定性问题。解决固定随机种子只是治标更关键的是让结果对边界不敏感。检查距离矩阵的对角线是否严格为0检查eps的精度。我习惯把eps的精度控制在两位小数比如1.52海里然后用eps ± 0.1做敏感性测试看有多少轨迹的标签发生变化。如果变化超过5%说明eps落在了一个不稳定区间需要微调eps到更平稳的位置。另外跑聚类时把距离矩阵保存下来用同一个D重复聚类确认DBSCAN本身是否稳定如果同一个D的结果都不一致那说明代码有随机性优先修代码。5.3 翻车点三内置dbscan报错“Distance matrix must be symmetric”现象调用内置dbscan(D, eps, minPts, Distance, precomputed)时Matlab报错说距离矩阵不对称。原因内置dbscan对距离矩阵的对称性检查非常严格要求norm(D - D)精确为0。但我们的距离矩阵在计算时如果用了浮点运算D(i,j)和D(j,i)可能存在1e-12级别的差异虽然不是业务问题但不满足内置函数的校验。另一个常见原因是重采样时某条轨迹有NaN导致pdist2返回NaNNaN在比较和判断时会破坏对称性。解决要么自己写主循环我上面的dbscan_dist就没有这个限制要么在调用内置函数前强制对称化D (D D) / 2同时用isnan检查距离矩阵里是否有NaN有NaN就把对应轨迹从集合里剔除。我的习惯是统一用自己写的dbscan_dist它对距离矩阵的要求只有对角线非负、其余位置为正数宽松得多。5.4 翻车点四聚类结果全部合并成一个簇航线区分度为零现象无论eps怎么调小聚类结果始终是一个大簇加上极少数噪声点完全看不出多条航线的结构。原因距离定义有问题。如果你用的是标准Hausdorff距离而数据里存在一条特别长的轨迹它和所有轨迹的距离都很大但剩余轨迹之间的标准Hausdorff距离都很小于是大部分轨迹被聚在一起。还有一种情况是坐标没转平面坐标经纬度直接参与计算在低纬度海区经度方向的物理距离被严重低估导致所有轨迹看起来都很近。解决先检查距离矩阵的数值范围。把距离矩阵的最大值、中位数、90%分位数打印出来如果在同一片研究区域里最大距离和90%分位数相差不超过3倍说明距离定义区分度不够改用MHD或截断Quantile。如果是坐标问题确认lng_scale cosd(lat0)这一步有没有正确执行纬度方向乘以60后经度方向也应该乘以60*lng_scale而不是只乘lng_scale。我把这个错误犯过一次排查了两天才发现经度方向少乘了60导致东西向距离被压缩了60倍聚类自然全糊在一起。5.5 翻车点五轨迹数量过万时距离矩阵算到一半内存爆了现象代码在计算距离矩阵的双层循环里跑了几分钟后Matlab直接报错“Out of memory”或者系统变得极卡风扇狂转。原因双层循环里每次都调用pdist2临时变量会被反复创建和释放Matlab的内存碎片化严重。另外如果用了向量化方案把D_all一次性算出来内存占用是(n_tracks*N)^2 * 8字节一万条轨迹、N64时D_all要占3.3TB显然不可能。问题出在单机内存上限而不是算法逻辑。解决分块计算是唯一出路。我的做法是把轨迹分成块每块500条两个块之间计算交叉距离后立即存入D矩阵的对应区域释放临时变量。同时把pdist2的调用限制在块内不要一次性算整个大矩阵。另外可以合理利用Matlab的single类型距离矩阵用single存储内存直接减半只在计算最小值和均值时临时转double。如果一万条轨迹确实要处理建议先降低重采样点数N到32或者用tall数组配合datastore但那个复杂度太高不如分块实在。我的经验是单机Matlab用分块方案处理一万条轨迹、N32距离矩阵计算时间大约20到30分钟属于可以过夜的批量任务。6. 聚类结果怎么验证才算数内部指标、外部对照和可视化三个习惯聚类跑完不是终点验证才是让结果能拿去写报告的那一步。我习惯同时做三个层面的验证缺一不可。第一个层面是内部指标。DBSCAN不像KMeans那样有明确的SSE曲线我常用的是轮廓系数Silhouette Coefficient但直接调用Matlab的silhouette函数时会发现它要求传入原始特征矩阵不认距离矩阵。这时要自己写对每个样本计算它与自己簇内所有样本的平均距离a再计算它与最近邻簇的平均距离b轮廓系数s (b - a) / max(a, b)。噪声点不参与轮廓系数计算因为它们的归属是-1没有簇内距离。整体轮廓系数超过0.4说明簇结构清晰低于0.2说明eps或距离定义还有问题。注意轮廓系数对团状簇有效对任意形状的轨迹簇可能偏低所以还要结合下面的外部对照。% 基于距离矩阵计算整体轮廓系数噪声点除外 function sc silhouette_from_dist(D, labels) valid labels 0; Dv D(valid, valid); lv labels(valid); n size(Dv, 1); s_all zeros(n, 1); for i 1:n same_cluster lv lv(i); % 簇内平均距离 a mean(Dv(i, same_cluster (1:n) ~ i)); % 找最近邻簇 other_clusters unique(lv(lv ~ lv(i))); b_min inf; for k 1:length(other_clusters) idx_other lv other_clusters(k); b_k mean(Dv(i, idx_other)); b_min min(b_min, b_k); end s_all(i) (b_min - a) / max(a, b_min); end sc mean(s_all); end第二个层面是外部对照。聚类结果最好能和真实业务标签做对比比如已知某条航线上的固定班轮或者已知的锚地分区。如果没有现成标签我常用的办法是拿一周的AIS数据聚类再用下一周的数据做时间外验证同一片海域、同一批起终点聚类出的航线结构应该基本一致。如果跨周验证时航线中心线偏移超过20%说明聚类结果里混入了临时性干扰比如临时交通管制、恶劣天气绕航需要进一步清洗。第三个层面是可视化验证。不要只看散点图要把每条簇的轨迹中心线画出来。中心线求法很简单对簇内所有重采样轨迹的对应点做中位数得到一条代表航线再画一个以中心线为中心、宽度为轨迹距离标准差的范围带。如果两条中心线交叉或者范围带大面积重叠说明这两个簇实际上是一条航线被切开了。这个可视化检查比任何指标都直观我每次聚类完必做能快速发现用轮廓系数看不出来的问题。做完这三个层面的验证后聚类结果才算真正可用。最后说一个我的习惯每次跑完聚类不管结果好坏我都会把eps、minPts、距离定义、N、噪声点比例五个参数连同聚类结果的截图存成一个版本记录。这个看起来繁琐但当你换了一批数据、调了三个月参数后回头找当年能用的参数组合时这份记录就是后悔药。没有它所有调参经验都要重来一遍。希望这些记录方式对你有用也祝你跑出来的聚类结果一次比一次干净。本文还有配套的精品资源点击获取