简介:本资源面向船舶航迹研究与Matlab聚类学习者,复现论文《基于轨迹聚类的船舶异常行为识别研究》的核心流程,解决航迹数据提取、聚类与偏离预测的完整实现问题。包内共20个文件,以14个m脚本为主,辅以mat数据、zip子包与png图示,整体约4.32MB,涵盖DBSCAN聚类、改进Hausdorff距离计算、航迹提取、阈值分类与阈值寻优等模块,可直接运行。已有766人学习下载。读者可借此掌握DBSCAN距离算法与H距离计算方法,理解航迹分类大致流程,并基于聚类结果完成航迹偏离预测,经检测准确率较高。船舶航迹方向的研究者可替换数据与相关模块,将其拓展为自己的模型;需要学习DBSCAN的读者也可将其作为应用案例,快速上手聚类与预测的完整链路。
1. 从一堆乱麻般的AIS轨迹说起:这个项目到底在干什么
如果你手头有一批船舶AIS轨迹数据,每条轨迹几百到几千个点,几十上百条船混在一起,想按航行模式自动分成几类——比如近岸穿梭的、直线过境的、锚泊打转的——直接上DBSCAN大概率会翻车。原因很简单:DBSCAN默认用欧氏距离算两点相似度,而两条轨迹的采样点数量、时间戳、航向角往往对不齐,直接算距离等于拿苹果比橘子。这个Matlab项目就是冲着这个痛点来的:它把轨迹间的相似度度量换成了改进的Hausdorff距离,再喂给DBSCAN做聚类,让“形状相似但采样点不对齐”的轨迹也能被正确归到一类。整套代码在Matlab环境下运行,适合做船舶行为分析、航迹预测预处理、海上交通流挖掘的开发者直接拿来改。你不需要从零推导距离公式,但得知道每个参数动了之后聚类结果会怎么变——这正是下面要拆清楚的事。
2. 核心模块拆解:改进Hausdorff距离怎么算、DBSCAN怎么接
2.1 为什么原始Hausdorff距离在航迹上不好用
原始Hausdorff距离的定义是:集合A中每个点到集合B中最近点的距离的最大值,再取双向的最大。放到船舶轨迹上,它有两个致命问题。第一,它对离群点极其敏感——一条轨迹里只要有一个GPS跳变点,整个距离就被拉爆。第二,它只取“最远最近点”,忽略了轨迹的整体走向,两条形状完全不同的轨迹可能因为局部某一段靠近而得到一个很小的距离值。
改进思路通常从两个方向走:一是用平均Hausdorff距离替代最大距离,把“最远”换成“平均”,降低离群点影响;二是引入航向、航速等运动特征做加权,让距离度量不只依赖空间坐标。这个项目采用的是加权平均Hausdorff距离,权重向量在代码里是一个可配置参数。
function d = improvedHausdorff(trajA, trajB, w) % trajA, trajB: N×3 矩阵,列分别为 x, y, 航向角(弧度) % w: 1×3 权重向量,对应位置x、位置y、航向的权重 % 返回改进后的双向平均Hausdorff距离 % 计算A中每个点到B中每个点的加权距离矩阵 distMat = zeros(size(trajA,1), size(trajB,1)); for i = 1:size(trajA,1) for j = 1:size(trajB,1) dx = trajA(i,1) - trajB(j,1); dy = trajA(i,2) - trajB(j,2); dtheta = abs(trajA(i,3) - trajB(j,3)); % 航向角环形修正,避免359°和1°被算成358° dtheta = min(dtheta, 2*pi - dtheta); distMat(i,j) = w(1)*abs(dx) + w(2)*abs(dy) + w(3)*dtheta; end end % 单向平均Hausdorff:A中每个点到B的最近距离,再取平均 hAB = mean(min(distMat, [], 2)); hBA = mean(min(distMat, [], 1)); d = max(hAB, hBA); % 双向取大,保持对称性 end这段代码的关键在三个地方。第一,距离矩阵是逐点算的,时间复杂度O(N×M),轨迹点多了会慢,实际用的时候一般先做抽稀。第二,航向角做了环形修正,min(dtheta, 2*pi - dtheta)这行如果漏掉,359°和1°会被算成358°的差距,聚类结果直接乱掉。第三,最后取max(hAB, hBA)而不是平均,是为了保持距离度量的对称性——如果取平均,d(A,B)和d(B,A)虽然相等,但会丢失“单向最大偏离”的信息。
权重向量w一般怎么设?如果只关心空间位置,设成[1, 1, 0];如果航向很重要(比如分进出港航道),航向权重可以给到0.5到1之间,具体看数据里航向的数值范围。位置坐标如果是经纬度,记得先转成平面坐标(比如UTM),否则纬度和经度的实际距离尺度不一样。
2.2 DBSCAN的输入不是点,是距离矩阵
标准DBSCAN的实现通常接收一个N×2的点坐标矩阵,内部自己算欧氏距离。但这个项目里,每个“点”是一条完整轨迹,轨迹之间的相似度由上面的改进Hausdorff距离给出。所以需要把DBSCAN改造成接收预计算的距离矩阵。
function [labels, coreIdx] = dbscanOnDistMat(D, eps, minPts) % D: N×N 对称距离矩阵,D(i,j)为轨迹i和j之间的改进Hausdorff距离 % eps: 邻域半径阈值 % minPts: 核心点的最小邻居数 % labels: N×1 聚类标签,-1表示噪声 N = size(D, 1); labels = zeros(N, 1); % 0表示未访问 clusterId = 0; coreIdx = []; for i = 1:N if labels(i) ~= 0 continue; % 已处理过 end % 找i的eps邻域 neighbors = find(D(i, :) <= eps & (1:N) ~= i); if length(neighbors) < minPts labels(i) = -1; % 暂标噪声,后续可能被归入某簇 continue; end % i是核心点,新建簇 clusterId = clusterId + 1; labels(i) = clusterId; coreIdx = [coreIdx; i]; % 队列扩展 seedSet = neighbors(:)'; while ~isempty(seedSet) q = seedSet(1); seedSet(1) = []; if labels(q) == -1 labels(q) = clusterId; % 边界点归入当前簇 end if labels(q) ~= 0 continue; end labels(q) = clusterId; qNeighbors = find(D(q, :) <= eps & (1:N) ~= q); if length(qNeighbors) >= minPts seedSet = [seedSet, qNeighbors(:)']; end end end end这个实现和标准DBSCAN在逻辑上完全一致,区别只在于邻域查询从“算坐标距离”变成了“查距离矩阵”。有两个参数需要你根据数据调:eps和minPts。eps的选取有个经验做法——把距离矩阵中所有非对角元素拿出来排序,画一条k-距离曲线(k取minPts-1),拐点位置就是比较合适的eps。minPts一般取轨迹总数的1%到5%,太少会导致大量小簇,太多会把不同类合并。
注意:距离矩阵是对称的,但如果你在改进Hausdorff里用了非对称的权重或者单向计算,记得手动对称化
D = (D + D')/2,否则DBSCAN的邻域查询会出问题。
2.3 完整调用链路:从原始AIS数据到聚类标签
把上面两个模块串起来,主流程大概长这样:
% 1. 读取原始轨迹数据,假设每条轨迹是一个结构体 % trajData: 1×M cell,每个cell是Ni×3矩阵 [x, y, heading] load('ais_trajectories.mat'); % 2. 轨迹抽稀,降低计算量(Douglas-Peucker或等间隔采样) trajData = cellfun(@(t) downsampleTraj(t, 50), trajData, 'UniformOutput', false); % 3. 计算两两之间的改进Hausdorff距离 M = length(trajData); D = zeros(M, M); w = [1, 1, 0.5]; % 位置权重1,航向权重0.5 for i = 1:M for j = i+1:M D(i,j) = improvedHausdorff(trajData{i}, trajData{j}, w); D(j,i) = D(i,j); end end % 4. 用k-距离曲线辅助选eps k = 4; % minPts = 5 sortedD = sort(D(D>0)); kdist = sortedD(1:k:end); figure; plot(kdist); xlabel('轨迹序号'); ylabel('k-距离'); % 观察拐点,手动设定eps % 5. 运行DBSCAN eps = 0.8; % 根据上一步拐点调整 minPts = 5; [labels, coreIdx] = dbscanOnDistMat(D, eps, minPts); % 6. 可视化:不同颜色画不同簇的轨迹 plotClusters(trajData, labels);抽稀那一步很关键。原始AIS数据一分钟一个点,一条轨迹几千个点,两两算距离是O(N²)的复杂度,不抽稀跑不动。抽稀到50到100个点通常能保留轨迹形状,同时把计算量降一个数量级。w的权重需要根据你的数据特点调——如果航向数据噪声大,把航向权重降到0.1甚至0;如果轨迹空间分布密集但航向差异明显,航向权重可以提到1。
3. 跑起来之前:环境、数据格式和参数调试
3.1 Matlab版本和依赖
这个项目对Matlab版本没有特别苛刻的要求,R2018a及以上都能跑。用到的函数都是基础函数:矩阵运算、cellfun、sort、plot,不需要额外工具箱。如果你用的是Octave,大部分代码也能跑,但cellfun的UniformOutput参数在旧版Octave里行为略有差异,建议在Matlab里跑。
数据输入格式是cell数组,每个cell是一个Ni×3的矩阵,三列分别是x坐标、y坐标、航向角(弧度)。如果你的原始数据是经纬度,先转平面坐标;如果航向是角度制,先转弧度。这些预处理代码不在项目主文件里,需要你自己写一小段。
3.2 参数调试的实操顺序
不要一上来就调DBSCAN的eps和minPts,先把距离矩阵的质量确认好。具体顺序是:
第一步,固定权重w,跑一遍距离矩阵,看距离值的分布。如果大部分距离值都集中在很小的范围,说明权重设置让轨迹之间区分度不够,需要调整权重或者检查数据是否已经归一化。
第二步,画k-距离曲线。minPts先设一个大概值(比如5),看曲线拐点是否明显。拐点不明显说明距离度量本身区分度不够,回到第一步。
第三步,固定eps,调minPts。从3开始往上加,观察聚类数量的变化。通常存在一个区间,minPts在这个区间内变化时聚类数量稳定,取区间中值。
第四步,固定minPts,微调eps。eps增大,簇会合并;eps减小,簇会分裂。根据你对“一类”的定义来定。
% 参数扫描示例:观察不同eps下的聚类数量 epsList = 0.4:0.1:1.2; numClusters = zeros(size(epsList)); for k = 1:length(epsList) labels = dbscanOnDistMat(D, epsList(k), 5); numClusters(k) = length(unique(labels(labels>0))); end plot(epsList, numClusters, '-o'); xlabel('eps'); ylabel('聚类数量');这段扫描代码能帮你快速看到eps对结果的影响趋势。如果聚类数量在某个eps区间内保持稳定,那个区间就是比较鲁棒的选择。
3.3 结果验证:怎么知道聚对了
聚类没有绝对正确的答案,但有几个检查手段。一是可视化,把不同簇的轨迹用不同颜色画在地图上,看空间上是否合理分离。二是算轮廓系数,虽然它是为欧氏距离设计的,但在预计算距离矩阵上也能用——把距离矩阵传给silhouette函数,看平均轮廓系数是否大于0.5。三是抽样检查,从每个簇里随机抽几条轨迹,看它们的航向、航速变化模式是否一致。
提示:如果某个簇的轨迹数量特别大,占了总数的一半以上,大概率是
eps设大了,把不同类合并了。反过来,如果噪声点(label=-1)超过30%,说明eps太小或者minPts太大。
4. 避坑与排查:跑这个项目最容易翻车的五个地方
4.1 航向角没做环形修正,聚类结果完全乱套
现象:两条航向分别是359°和1°的轨迹,被算成距离很大,分到了不同簇。原因:直接用角度差绝对值,359-1=358,但实际只差2°。解决:在距离计算里加环形修正,dtheta = min(dtheta, 2*pi - dtheta)。如果航向是角度制,先转弧度再算。
4.2 经纬度直接当平面坐标用,距离尺度失真
现象:聚类结果在地图上看起来明显不对,东西方向和南北方向的距离被扭曲。原因:经纬度是球面坐标,1度经度和1度纬度的实际距离不同(除非在赤道)。解决:先用utmgeoid或手动公式把经纬度转成UTM平面坐标,再算距离。如果数据范围小,也可以用简单的等距圆柱投影近似。
4.3 轨迹点数量差异太大,平均Hausdorff被长轨迹主导
现象:一条1000个点的轨迹和一条50个点的轨迹算距离,结果偏向长轨迹。原因:平均Hausdorff对每个点取平均,长轨迹的每个点都参与,短轨迹的贡献被稀释。解决:先对所有轨迹做等间隔抽稀,统一到相近的点数(比如都抽到50到100个点),再算距离。
4.4 距离矩阵太大导致内存溢出
现象:跑1000条轨迹时Matlab报内存不足。原因:1000×1000的double矩阵是8MB,不算大,但如果轨迹本身很大(每条几千个点),在算距离矩阵时的中间变量会爆。解决:分块计算距离矩阵,每次算100×100的块,存到稀疏矩阵或直接写入文件。或者先对轨迹做聚类预筛选,把明显不相似的轨迹对跳过。
4.5 eps和minPts用默认值,结果没法看
现象:直接拿网上教程里的eps=0.5, minPts=10跑,要么全是一个簇,要么全是噪声。原因:不同数据的距离尺度完全不同,没有通用参数。解决:老老实实画k-距离曲线,根据拐点定eps;minPts从3开始试,观察聚类数量稳定性。这一步没有捷径。
5. 二次开发:换距离度量、加运动特征、接预测模型
这个项目的骨架是“距离矩阵 + DBSCAN”,所以二次开发最自然的切入点就是换距离度量。改进Hausdorff只是其中一种选择,你完全可以换成DTW(动态时间规整)、LCSS(最长公共子序列)或者Fréchet距离。换的时候只需要改improvedHausdorff函数,保持输入输出格式一致(输入两条轨迹矩阵,输出一个标量距离),DBSCAN部分不用动。
function d = dtwDistance(trajA, trajB) % 用DTW替代Hausdorff,适合采样率不一致的轨迹 Na = size(trajA,1); Nb = size(trajB,1); Dmat = inf(Na+1, Nb+1); Dmat(1,1) = 0; for i = 1:Na for j = 1:Nb cost = norm(trajA(i,1:2) - trajB(j,1:2)); Dmat(i+1,j+1) = cost + min([Dmat(i,j+1), Dmat(i+1,j), Dmat(i,j)]); end end d = Dmat(Na+1, Nb+1); endDTW的好处是允许时间轴上的伸缩,两条轨迹即使一条快一条慢也能对齐。代价是计算复杂度O(N×M),比Hausdorff的O(N×M)在常数项上更大,而且DTW距离不满足三角不等式,DBSCAN的邻域查询逻辑虽然不受影响,但k-距离曲线的解释会略有不同。
另一个改造方向是往距离度量里加更多运动特征。现在的权重向量只有三维(x, y, 航向),你可以扩展到五维甚至更多:航速、加速度、转向率。每加一维,权重向量的长度加一,距离计算里多一项加权。但要注意,维度越高,权重调参的难度越大,而且高维空间里距离的区分度会下降(维度灾难)。我一般会先做特征相关性分析,把高度相关的特征去掉,保留三到四个独立特征就够了。
如果你做这个聚类的最终目的是航迹预测,那聚类标签可以作为预测模型的输入特征。比如用LSTM做轨迹预测时,把每条轨迹的聚类标签做one-hot编码,拼接到输入序列后面,能让模型区分不同航行模式。这一步的接口很简单:labels是N×1的向量,labels(i)对应第i条轨迹的簇编号,直接拼到你的特征矩阵里就行。
还有一个实用技巧:增量聚类。船舶轨迹数据是持续到来的,每次新来一批轨迹就重新跑全量DBSCAN太慢。可以把已有簇的核心点保存下来,新轨迹先和核心点算距离,如果落在某个核心点的eps邻域内,直接归入该簇;否则再跑局部DBSCAN。这样能把计算量从O(N²)降到O(N×K),K是核心点数量。
注意:增量聚类会引入“概念漂移”问题——随着时间推移,新的轨迹模式可能出现,已有簇的核心点不能代表新数据。我一般会每隔一段时间(比如一周)做一次全量重聚类,用新结果替换旧的核心点集合。
从那以后我每次换距离度量或者加新特征,都强制走一遍“距离分布检查 → k-距离曲线 → 参数扫描”这三步,不再凭感觉设eps。希望帮到你。
本文还有配套的精品资源,点击获取