Matlab实现AIS航迹聚类：改进Hausdorff距离与DBSCAN参数调优

简介：面向船舶航迹聚类与异常行为识别研究的Matlab完整工程包，复现了基于改进Hausdorff距离的DBSCAN聚类方法。资源以论文《基于轨迹聚类的船舶异常行为识别研究》为蓝本，完整覆盖航迹数据提取、聚类中心获取、基于H距离的航迹预测及阈值寻优等关键环节，并包含轨迹滤波、交叉点提取、误差率计算等子功能，代码经整理可一键运行。包内共20个文件，以14个M脚本为核心，涵盖DBSCAN聚类、H距离计算、阈值分类、航迹提取绘制等核心函数，另含mat数据文件、结果压缩包、说明文档与预览图片，整体约4.32MB，模块边界清晰，便于替换数据或修改算法拓展至自己的模型。目前已有766人学习下载。借助该案例可系统掌握DBSCAN与Hausdorff距离的实现细节，理解航迹分类的完整流程，并基于聚类结果进行航迹偏离预测，对船舶航迹方向的研究者和需要DBSCAN应用案例的Matlab学习者都很有参考价值。

1. 这个项目到底解决了什么：让AIS航迹聚类不再是一堆散点拼图

拿到一整个港区或某片海域的AIS报文，最常见的尴尬是：点画出来很密，但看不出几条船在走什么固定路线。聚类停在散点级别，结果永远是“这里船多、那里船少”，对航道规划、船舶行为分析、航迹预测都没什么用。这个Matlab项目解决的是把“点聚类”升级成“轨迹聚类”：先用改进的Hausdorff距离度量两条航迹的相似程度，再交给DBSCAN去分簇，最后得到的是若干条典型航路。适合正在做船舶航迹聚类、航路提取、异常行为识别的人，尤其是手里有AIS数据但不知道从哪一步开始处理的读者。标题里拼写的Harsdorf其实是Hausdorff（豪斯多夫），下文中统一用正确拼写。这类项目最值得啃的不是聚类调用本身，而是距离怎么定义、噪声怎么压、两个聚类参数怎么定。

2. 数据与预处理链路：从AIS原始报文到干净的多段航迹

2.1 AIS原始数据长什么样：每一行是一条带状态的点

AIS报文最常见的落地方式是CSV或TXT，每行记录一个时刻的船舶状态。核心字段一般是：MMSI（船号）、经度、纬度、SOG（对地航速）、COG（对地航向）、时间戳。有时候还有船型、船长、艏向，但对聚类来说，MMSI、经纬度、时间戳三样就够起步了。这个项目里读取数据的典型代码块如下：

function [trajDB] = loadAISData(filepath) % 读取AIS CSV，按MMSI和时间排序，返回结构体数组 opts = detectImportOptions(filepath); opts.VariableNamingRule = 'preserve'; T = readtable(filepath, opts); % 列名做一次容错映射，常见命名如 lon/lat 或 longitude/latitude if ismember('longitude', T.Properties.VariableNames) T.Properties.VariableNames{'longitude'} = 'Lon'; T.Properties.VariableNames{'latitude'} = 'Lat'; end % 剔除明显异常：经纬度越界、SOG为负、时间缺失 valid = T.Lon > -180 & T.Lon < 180 & ... T.Lat > -90 & T.Lat < 90 & ... T.SOG >= 0 & ~ismissing(T.Time); T = T(valid, :); % 按船、按时间排序，为后续轨迹切分打基础 T = sortrows(T, {'MMSI', 'Time'}); trajDB = T; end

这段代码做的事是先容忍列名差异，再剔除越界脏点，最后按MMSI加时间排序。参数上需要注意：如果数据里时间字段是字符串，需要在外面先统一成datetime类型，否则sortrows排序出来会是字典序，导致轨迹乱序。detectImportOptions对带表头的CSV很好用，但对几百MB的大文件建议手动指定列号读取，速度差一个量级。

2.2 航迹切分与重采样：断点判断比想象中更讲究

一条船不可能一天24小时都处于“连续运动”状态。定位会丢星、船会靠泊、AIS设备会间歇性关机，这些都会让轨迹断裂。切分逻辑上常用的两个阈值是时间间隔和距离间隔：如果两条相邻记录的时间差超过设定值（比如10分钟），或位移超过设定值（比如5海里），就判定为两条航迹。这个项目里切分后还要做时间重采样，把不规则报文变成等间隔轨迹，否则后面算Hausdorff距离时两条轨迹点的对应关系会非常混乱。

function [segments] = splitTrajectories(T, maxTimeGap, maxDistGap) % 输入排序后的AIS表，输出航迹元胞数组 mmsiList = unique(T.MMSI); segments = cell(0, 1); for i = 1:length(mmsiList) idx = T.MMSI == mmsiList(i); sub = T(idx, :); if height(sub) < 5 continue; % 点数太少，不足以构成航迹，直接跳过 end % 按时间间隔与球面距离间隔切分 dT = hours(diff(sub.Time)); dDist = distance(sub.Lat(1:end-1), sub.Lon(1:end-1), ... sub.Lat(2:end), sub.Lon(2:end)); gapIdx = find(dT > maxTimeGap | dDist > maxDistGap); % 切分点位置 + 首尾索引，生成子轨迹 breaks = [1; gapIdx+1; height(sub)+1]; for j = 1:length(breaks)-1 seg = sub(breaks(j):breaks(j+1)-1, :); if height(seg) >= 5 segments{end+1, 1} = seg; %#ok<AGROW> end end end end

切分阈值不是拍脑袋出来的。港口内作业船舶的AIS中断可能只有两分钟，而远洋航行中断十几分钟也很正常。我第一次跑这套代码时用了一个固定的5分钟阈值，结果港内拖轮被切得稀碎。后来改成按船速动态调整：慢速船用小时间窗，快速船用大时间窗。这是个值得注意的经验点。

2.3 等距重采样与坐标归一化：把轨迹变成同长度的向量

切分后的轨迹长度参差不齐，有的50个点，有的200个点。直接算Hausdorff距离虽然不要求两条轨迹点数严格相等，但点数差异过大会被离群点带偏。这个项目的做法是把每条轨迹重采样到固定点数N，这里N=64。做法是对轨迹的总长度做归一化，然后按弧长等比例取点，而不是按时间等间隔取点——因为船在转弯时会减速，时间等间隔会让转弯处的点更密，弧长等比例则能保证形状特征被均匀保留。

function [traj] = resampleTrajectory(seg, N) % 按弧长等比例重采样，输出 N x 2 的经纬度序列 lat = seg.Lat; lon = seg.Lon; % 球面距离累加求累计弧长 d = distance(lat(1:end-1), lon(1:end-1), lat(2:end), lon(2:end)); cumLen = [0; cumsum(d)]; totalLen = cumLen(end); if totalLen == 0 traj = [lat(1), lon(1)] * ones(1, N); % 静止轨迹退化处理 return; end % 目标弧长位置，差值得到对应坐标 sTarget = linspace(0, totalLen, N)'; traj = [interp1(cumLen, lat, sTarget, 'linear'), ... interp1(cumLen, lon, sTarget, 'linear')]; end

重采样之后，每一条轨迹都是 N x 2 的矩阵。这里还有一个细节：纬度差1度对应的地面距离约111公里，经度差1度对应的距离则随着纬度变化，直接拿经纬度差当作空间距离会造成高纬度地区轨迹被“横向拉伸”。这套代码后面算距离时用的是球面距离函数，避免了这个坑。同时，为了后续可视化时坐标轴不变形，也可以用米制单位（UTM投影）做坐标变换，但球面距离已经够用。

3. 核心算法拆解：改进Hausdorff距离如何让DBSCAN在航迹上落地

3.1 标准Hausdorff距离为什么不够用：最大最小值太脆弱

两条轨迹的Hausdorff距离定义为：先求A中每个点到B中最近点的距离，取所有“最近距离”里的最大值；再从B向A算一遍，最终取两者中更大的那个。公式上看很直观，描述的是两条曲线整体偏离的最大程度。它的致命弱点是单个离群点或一条轨迹上个别异常定位点，就能把距离值拉得很大。AIS数据里野值多，经常出现一个跳点，比如突然偏移两三海里，标准Hausdorff距离立即失真。DBSCAN又对距离值极度敏感，距离矩阵稍微偏大，聚类就会退化成一堆噪声点。

改进方向常见的有两种。一种是改用单向Hausdorff距离（Directed Hausdorff），只取一个方向，减少部分对称性带来的冗余计算；另一种是使用MHD（Modified Hausdorff Distance），核心思想是把最大最小值替换成平均最近距离：A中每个点到B的最近距离求平均，再与B到A的平均值取平均。这个做法对单点跳变的容忍度明显提升。本项目的改进点就是在MHD基础上，叠加航向夹角加权。

3.2 改进距离实现：平均最近距离加航向加权

这个项目的距离计算函数拆开看分三步：先算点集间的球面最近距离矩阵；再在点对距离基础上叠加航向差异惩罚；最后把两个方向的结果对称平均，得到距离值。航向信息来自AIS原始报文中的COG字段。航向加权的意义在于：两条空间上接近但方向相反的船，比如进出港航道上的对向航行，在航路聚类里应该被分到不同簇。如果只用空间位置，进出港双向往返航线会混成一条，无法支撑后续的航迹预测。

function [d] = trajDistanceMHD(P, Q, w_cog) % 改进Hausdorff距离：平均最近距离 + 航向夹角加权 % P, Q: N x 4 矩阵 [Lon, Lat, COG, SOG] % w_cog: 航向权重系数，0表示不启用航向加权 Np = size(P, 1); Nq = size(Q, 1); % 第一步：球面距离矩阵，这里用haversine公式 lon1 = P(:,1) * pi/180; lat1 = P(:,2) * pi/180; lon2 = Q(:,1) * pi/180; lat2 = Q(:,2) * pi/180; dlon = lon2' - lon1; % Np x Nq dlat = lat2' - lat1; a = sin(dlat/2).^2 + cos(lat1) * cos(lat2') .* sin(dlon/2).^2; distMat = 6371 * 2 * atan2(sqrt(a), sqrt(1-a)); % 单位km % 第二步：航向夹角差异矩阵，归一化到[0,1] cog1 = P(:,3) * pi/180; cog2 = Q(:,3) * pi/180; % 角度差取绝对值并折叠到[0, pi] angDiff = abs(cog2' - cog1); angDiff = min(angDiff, 2*pi - angDiff); angDiff = angDiff / pi; weightMat = 1 + w_cog * angDiff; % 夹角越大权重越高 % 第三步：距离矩阵叠加权重，再算平均最近距离 M = distMat .* weightMat; % A->B 方向平均最近距离 minAB = min(M, [], 2); % B->A 方向平均最近距离 minBA = min(M, [], 1); d = 0.5 * (mean(minAB) + mean(minBA)); end

参数说明里最值得关注的是w_cog。设成0时退化成空间MHD距离；设成1时航向差异最大可以把距离放大到原来的2倍；设成0.5时方向相反的轨迹距离增加约50%。我用0.3左右跑港口进出港数据效果比较好——进出双向能分开，又不会把交叉相遇的航道拆得过碎。计算距离矩阵时有一个明显可以加速的地方：提前计算好每条轨迹的重采样点坐标，把P和Q都转成米制坐标，然后用pdist2计算欧氏距离矩阵，最后再转成公里单位，能省掉大量三角函数的重复调用。

3.3 DBSCAN聚类主流程：距离矩阵决定一切

DBSCAN聚类本身不直接消费轨迹数据，它消费的是上面算出来的两两距离矩阵。核心流程是：指定一个邻域半径Eps和一个最小点数MinPts，然后遍历所有轨迹段：如果某条轨迹的Eps邻域内轨迹数大于等于MinPts，就把它标为核心轨迹；核心轨迹密度可达的轨迹归入同一簇；不满足核心条件也不可达的标为噪声。工程实现上有两种路径：一种是调用Matlab自带的dbscan函数，但需要注意它接受的是点坐标矩阵而不是距离矩阵，需要用squareform把距离矩阵转成向量再调用；另一种是自己写循环版DBSCAN，便于二次开发和加日志。

function [labels, coreFlags] = dbscanTrajectories(D, eps, minPts) % D: N x N 距离矩阵 % eps: 邻域半径 % minPts: 最小邻域轨迹数 N = size(D, 1); labels = zeros(N, 1); coreFlags = false(N, 1); clusterId = 0; % 用稀疏矩阵找邻域，避免大规模全矩阵多次索引 neighborIdx = cell(N, 1); for i = 1:N neighborIdx{i} = find(D(i, :) <= eps); end % 先标记核心点 for i = 1:N if numel(neighborIdx{i}) >= minPts coreFlags(i) = true; end end % BFS遍历密度可达簇 for i = 1:N if labels(i) ~= 0 || ~coreFlags(i) continue; end clusterId = clusterId + 1; labels(i) = clusterId; queue = neighborIdx{i}; while ~isempty(queue) j = queue(1); queue(1) = []; if labels(j) ~= 0 continue; end labels(j) = clusterId; if coreFlags(j) queue = [queue, neighborIdx{j}]; %#ok<AGROW> end end end % 边界点归属于最早访问它的核心点簇 labels(labels == 0 & coreFlags) = clusterId + 1; labels(labels == 0) = -1; % 噪声点 end

这段代码是教学版BFS实现，逻辑清晰但有三处实际使用时需要注意。第一，省内存但费时间的写法是上面这样；如果矩阵规模在500条以上，建议直接用Matlab内置的dbscan函数加squareform，速度提高数倍。第二，边界点归属顺序会影响最终簇编号但不影响分簇结果，因为这条轨迹最终落在哪个簇是由第一个访问它的核心点决定的。第三，minPts通常设置为3到5，设置成1会导致几乎所有互相碰得到的轨迹都连成一团，聚类结果退化成连通域剖分。

4. Matlab工程实现：模块划分、运行步骤与参数调优

4.1 目录结构与模块职责：一个脚本对应一个环节

这个项目的代码组织比较规整，适合直接照着改成自己的数据管线。根目录下分四类文件：数据读取脚本、预处理函数、距离计算和聚类主函数、可视化脚本。具体目录结构如下。

project_root/ ├── data/ │ ├── ais_raw.csv # 原始AIS报文 │ └── trajSet.mat # 预处理输出的航迹集合 ├── src/ │ ├── loadAISData.m # 读取原始AIS数据 │ ├── splitTrajectories.m # 按时间和距离阈值切分航迹 │ ├── resampleTrajectory.m # 等弧长重采样 │ ├── trajDistanceMHD.m # 改进Hausdorff距离 │ ├── dbscanTrajectories.m # DBSCAN聚类 │ └── plotClusterResult.m # 聚类结果可视化 ├── runCluster.m # 主流程脚本，一键运行 └── config.m # 全局参数配置文件

主流程脚本是runCluster.m，它的执行顺序就是上面目录里从上到下的调用顺序。config.m里放的是全局参数，包括前面提到的切分阈值、重采样点数、航向权重系数、Eps和MinPts。所有参数集中在一个文件里的好处是调参时不用改函数内部代码，缺点是参数多了以后容易出现“改了config但函数里没有引用到”的情况。我一般会在config.m最后加一段断言，检查每个参数是否在合理范围，避免跑到一半才发现Eps是负数。

4.2 主流程脚本拆解：一键跑通的骨架

runCluster.m是理解整个项目运行的入口。它做的事情可以用五步概括：加载数据、切分航迹、重采样、算距离矩阵、聚类并画图。下面给出去掉绘图细节后的核心骨架。

% runCluster.m % 一键运行脚本：AIS原始数据 -> 航迹集合 -> 距离矩阵 -> DBSCAN簇 clear; clc; close all; run('config.m'); % 第一步：读取原始AIS数据 fprintf('[1/5] 读取AIS数据...\n'); T = loadAISData('data/ais_raw.csv'); % 第二步：切分航迹，参数来自config fprintf('[2/5] 切分航迹...\n'); segments = splitTrajectories(T, cfg.maxTimeGap, cfg.maxDistGap); fprintf('共切分出 %d 条航迹段\n', length(segments)); % 第三步：重采样到统一点数 fprintf('[3/5] 航迹重采样...\n'); trajSet = cell(length(segments), 1); for i = 1:length(segments) trajSet{i} = resampleTrajectory(segments{i}, cfg.resampleN); end % 第四步：计算两两距离矩阵（核心耗时步骤） fprintf('[4/5] 计算改进Hausdorff距离矩阵...\n'); numTraj = length(trajSet); distMat = zeros(numTraj, numTraj); for i = 1:numTraj for j = i+1:numTraj d = trajDistanceMHD(trajSet{i}, trajSet{j}, cfg.w_cog); distMat(i, j) = d; distMat(j, i) = d; end end % 第五步：DBSCAN聚类并保存结果 fprintf('[5/5] DBSCAN聚类中...\n'); [labels, coreFlags] = dbscanTrajectories(distMat, cfg.eps, cfg.minPts); save('output/clusterResult.mat', 'labels', 'coreFlags', 'distMat', 'trajSet');

这段脚本的定位是所有后续改动的“主灶台”。替换数据时只需要改第一个loadAISData的路径；调整距离定义时只需要改trajDistanceMHD的内部实现；换成自适应Eps时只需要把第五步前面的定值eps改成计算得到的动态值。需要特别提醒的是第四步，双重循环的计算复杂度是O(N^2)，每条轨迹距离计算内部又是两个点集的双重循环，总计算量是O(N^2 * 64^2)。数据量上到几千条轨迹后，这一步能跑到让人怀疑人生。我后面在避坑章节里会写怎么处理。

4.3 参数表与调优顺序：先粗再细，别一上来就抠Eps

聚类的成败很大程度上不取决于算法实现，而是参数是否合理。这个项目涉及的核心参数按影响面排序如下表。调参顺序建议是：先确定重采样点数和切分阈值，再看距离矩阵的分布，最后才调Eps和MinPts。顺序反了会陷入组合爆炸——Eps调半天发现是切分参数不对导致的距离分布异常。

第1章提到的“翻车”高发区就在eps上。有些初学者直接用距离矩阵中位数的一半当eps，结果簇数奇多。比较稳的做法是先画出所有轨迹两两距离的直方图，找一个“谷底”作为eps初值——直方图双峰形状时选两峰之间的谷底，单峰时选中位数。参数设置完成后，可视化脚本会把每条轨迹按簇标成不同颜色，簇中心轨迹可以用各簇内距离中位数最小的那条轨迹代表，之后的航迹预测就在簇中心轨迹上做外推。

5. 避坑与排查：跑聚类代码时最常见的几个实际问题

5.1 聚类结果全是噪声点，画出来一片灰

这是一个极其常见的情况。现象是运行完成后labels数组里大多数值是-1，地图上所有轨迹都是灰色，只有零星两三条有颜色。原因通常有两层：一是距离矩阵里大部分轨迹间距离偏大，而eps取小了；二是改进Hausdorff距离里航向加权系数w_cog太大，导致空间上接近的轨迹因为航向差异被放大到超过eps。解决方法是先关掉航向加权看基线结果，即把w_cog改成0，看纯空间距离下的聚类是否正常。如果基线正常、加了航向就不行，降低w_cog到0.1附近。如果基线本身就不行，那就看距离直方图，重新选eps位置，选在双峰谷底而不取中位数。

5.2 同一份数据跑出不同簇数，结果不稳定

现象是第二次运行、或换了台电脑运行后，聚类簇数和边界点归属发生变化。原因通常出在DBSCAN的边界点归属和求解顺序上。教学版BFS实现里，邻域搜索顺序取决于轨迹在矩阵里的排列顺序，而轨迹排列顺序又取决于切分后segments元胞数组的顺序。如果切分时用了并行计算、或者数据读取时表排序不稳定，顺序轻微变化就会让边界点被归到相邻簇。解决方法是两步：一是对segments内的轨迹做一次按起点经纬度的稳定排序，保证输入顺序确定；二是用Matlab内置dbscan代替自己写的BFS版本，内置函数对边界点的处理逻辑是确定的。如果排查后发现顺序没问题簇数还不稳定，再看数据量——低于20条轨迹时DBSCAN本身就不稳定，最好多收集一段时间数据。

5.3 计算距离矩阵时内存和CPU被吃满，进度条一动不动

现象是跑到第四步，Matlab主界面左下角显示“正忙”，CPU单核满载，内存占用持续上涨，甚至中途内存溢出。原因很直接：双重循环未向量化且内部反复调用球面距离三角函数，大O复杂度太高。解决方法是三条路并行推进：第一，预计算每条轨迹的经纬度弧度值，把cos(lat)和sin(lat)存下来，不要在距离函数内重复计算；第二，用单精度single存储distMat——距离矩阵对精度要求不高，单精度能让内存减半；第三，对大规模数据使用Matlab Parallel Computing Toolbox把外层循环改成parfor。另外还有一个技巧：如果轨迹点数不多，可以把重采样点全部拼接成一个大矩阵，用pdist2一次性算出所有轨迹间点对距离，再做分块归约，速度通常能提升5到10倍。

5.4 重采样后轨迹长度一致了，但相似轨迹距离还是很大

现象是两条肉眼看起来几乎重合的航迹，距离矩阵里数值却很大，聚类分不到同一簇。原因往往是轨迹形状相同但方向相反，或者重采样起点没对齐。方向相反的场景已经可以通过航向加权解决；但不能依赖权重，更稳妥的做法是在切分阶段就判断轨迹的方向变化范围，把航向反转超过180度的轨迹直接切断成两段。重采样起点没对齐的原因则是对轨迹首尾端点衔接和循环性问题处理不当——航道环形路径比如防波堤绕行时，两条同向轨迹一个从中间切入、一个从起点切入，按弧长重采样后点集合看起来差异巨大。解决方法是算距离前先做一次端点对齐：以轨迹首尾点的连线为参考，尝试头尾翻转，取距离较小的那个结果作为最终距离。这也是改进Hausdorff距离在船舶场景里最常见的一个工程化补丁。

6. 二次开发思路：把聚类结果变成航迹预测和个性化距离函数

这套项目最大的实用价值在聚类之后：得到簇和簇中心轨迹后，可以做船舶典型航线提取、异常轨迹检测、短时航迹预测。我实际做过的扩展是把每一簇内所有轨迹聚合，计算每条轨迹对簇中心的MHD距离，取距离最小的轨迹作为簇中心，然后对簇中心轨迹做外推，实现“下一时刻船大概会走到哪个位置”的预测。具体做法是把簇中心轨迹按时间序列建模，用滑动窗口多项式拟合位置序列，外推未来几分钟的位置。预测误差评估用船位经纬度均方根误差，效果上比直接用原始AIS点做RTS平滑预测更稳定，原因是聚类消除了异常轨迹的干扰。

二次开发还有两个更直接的方向。一是把固定的eps改成自适应值：计算所有轨迹的k-最近邻距离，画出k-dist曲线，取曲线拐点作为eps。这个做法能解决不同时间段AIS疏密程度不同导致的聚类结果不可比问题。二是扩展距离函数：当前距离只用了经纬度、航向和航速，加入船舶类型、船长维度后，同一类船型的轨迹更容易聚合在一起，适合做分船型的航线行为分析。我当时把SOG导入了距离计算，用归一化速度差乘以一个较小权重，做出来之后对慢速作业船和快速过路船的区分明显更好。从那以后我每次做轨迹聚类都强制自己先看一眼距离矩阵直方图，再决定距离函数和eps，不再一上来就跑全文。希望帮到你。

本文还有配套的精品资源，点击获取