简介：本资源面向船舶航迹研究与Matlab聚类学习者，复现论文《基于轨迹聚类的船舶异常行为识别研究》的核心流程，解决航迹提取、聚类与偏离预测的实现问题。包内共20个文件，以14个m脚本为主，辅以mat数据、zip子包及png图示，整体约4.32MB，涵盖DBSCAN聚类、改进Hausdorff距离计算、航迹提取、阈值分类与寻优等模块，代码完整可直接运行。已有766人学习下载，适合作为DBSCAN与H距离计算的应用案例。读者可借此掌握航迹分类的完整流程，理解聚类中心提取与基于豪斯多夫距离的航迹预测方法，并通过阈值寻优获得较高检测准确率；船舶航迹方向的研究者可替换数据与相关模块，将案例拓展为自有模型，也可用于异常行为识别与航迹偏离预警的入门实践。
1. 从一次航迹交叉误判说起：Hausdorff 距离为什么在船舶聚类里会翻车
AIS 数据里最让人头疼的场景之一，是两条船在狭窄水道里对遇、交叉、再分开。你拿 DBSCAN 去跑聚类，本意是把同一条航道、同一种航行模式的轨迹归到一起，结果两条明明属于不同航路的船，因为中间有一段空间上贴得很近，被硬生生塞进同一个簇。这不是 DBSCAN 的错，是距离度量选错了。
经典 DBSCAN 默认用欧氏距离逐点比较，但船舶航迹是变长序列，两条轨迹点数不一样、采样间隔不一样、航速不一样，逐点欧氏距离根本对不齐。工程上更常见的做法是先把轨迹当点集，用 Hausdorff 距离衡量两条轨迹的"最大不匹配程度"。问题在于，标准 Hausdorff 距离取的是单向最大最小距离的最大值，它对局部离群点极其敏感——一条轨迹上只要有一个 AIS 跳点，整个距离就被这个点绑架，聚类结果直接崩。
这就是"改进的 Hausdorff 距离"要解决的问题：把最大值改成排序后的分位数，或者做平均，把个别异常点的影响压下去，同时保留 Hausdorff 对形状差异的敏感度。配合 DBSCAN 的密度可达逻辑，就能把"形状相似、空间邻近"的航迹聚成簇，而不是被一两个跳点带偏。
这篇笔记面向三类人：手上有一批 AIS 航迹点、想用 Matlab 做无监督聚类的航运数据分析从业者；做过 DBSCAN 但发现欧氏距离效果差、想换度量的算法工程师；以及需要把聚类结果接到航路提取、异常检测下游任务的人。下面从距离定义、DBSCAN 改造、参数调优到避坑，一步步给出可复现的 Matlab 实现路径。
2. 改进 Hausdorff 距离：从最大值到分位数，到底改了什么
2.1 标准 Hausdorff 距离的定义与它的致命短板
给定两条轨迹点集 $A={a_1,...,a_m}$ 和 $B={b_1,...,b_n}$，标准 Hausdorff 距离定义为：
$$H(A,B)=\max\left(\sup_{a\in A}\inf_{b\in B}|a-b|,\ \sup_{b\in B}\inf_{a\in A}|a-b|\right)$$
翻译成人话：对 A 里每个点，找它到 B 的最近距离，取这些最近距离里的最大值；B 对 A 同理；两个方向取大的那个。它衡量的是"两条轨迹最不匹配的那一处有多远"。
这个定义在图像模板匹配里很好用，因为图像是稠密点集，个别噪声点会被邻域平滑掉。但船舶航迹是稀疏、带噪、非均匀采样的序列，AIS 跳点、信号丢失补点、锚泊时的点堆积都会制造离群点。只要 A 中有一个点离 B 很远，$\sup$ 就把这个距离顶到天上，两条本该相似的轨迹距离被拉到极大，DBSCAN 里就永远达不到 eps 邻域，聚不到一起。
我早期直接用pdist2算双向最大最小再取 max，跑出来的簇数量对不上，人工核对发现同一航路的船被拆成七八个簇，就是被跳点坑的。
2.2 改进思路：分位数 Hausdorff 与平均 Hausdorff
改进的核心是把"取最大值"这个操作换掉。常见两种做法：
分位数 Hausdorff（Quantile Hausdorff）：把单向的最近距离集合排序，取第 $k$ 分位而不是最大值。设 $d_A(i)=\inf_{b\in B}|a_i-b|$，把它排序后取 $Q_p$ 分位：
$$H_q(A,B)=\max\left(Q_p{d_A(i)},\ Q_p{d_B(j)}\right)$$
$p$ 取 0.9 或 0.95 时，能容忍约 5%~10% 的离群点，同时保留对整体形状差异的敏感。$p=1$ 就退化成标准 Hausdorff。
平均 Hausdorff（Modified/Average Hausdorff）：直接对最近距离集合取均值：
$$H_{avg}(A,B)=\max\left(\frac{1}{m}\sum_{i=1}^{m}d_A(i),\ \frac{1}{n}\sum_{j=1}^{n}d_B(j)\right)$$
平均版本对离群点最鲁棒，但会弱化形状差异——两条整体偏移但形状一致的轨迹，平均距离可能也不小，这反而是我们想要的（偏移就该分开）。缺点是它对"局部形状突变"不敏感，比如一条直线航迹和一条带急转弯的航迹，平均下来可能差不多。
工程上我一般先用分位数版本，$p$ 从 0.95 往下调，观察簇数量和轮廓系数；如果数据跳点特别多（比如内河 AIS 丢包严重），再切到平均版本。两者在 Matlab 里都是几行的事，不用纠结理论最优，用聚类指标说话。
2.3 在 Matlab 里实现两种改进距离
下面给出可直接复用的函数。输入是两条轨迹的坐标矩阵（N×2 或 N×3），输出标量距离。
function d = hausdorff_quantile(A, B, p) % A, B: N×2 或 N×3 轨迹点坐标 % p: 分位数，0<p<=1，p=1 退化为标准 Hausdorff if nargin < 3 p = 0.95; end % 计算 A 中每个点到 B 的最近距离 dA = min(pdist2(A, B), [], 2); % m×1 % 计算 B 中每个点到 A 的最近距离 dB = min(pdist2(B, A), [], 2); % n×1 % 排序后取分位数，替代 max qA = quantile(dA, p); qB = quantile(dB, p); d = max(qA, qB); endfunction d = hausdorff_average(A, B) % 平均 Hausdorff：对最近距离取均值再取双向最大 dA = min(pdist2(A, B), [], 2); dB = min(pdist2(B, A), [], 2); d = max(mean(dA), mean(dB)); end逻辑说明：pdist2(A,B)返回 m×n 的距离矩阵，min(...,[],2)沿第二维取最小，得到每个点到另一条轨迹的最近距离。quantile对向量排序后插值取分位，比手写排序更稳。参数p是唯一需要调的：0.9 偏鲁棒，0.98 偏严格，建议从 0.95 起步。
提示：如果轨迹点带时间戳或航向，可以把距离改成加权形式，比如空间距离乘一个航向差惩罚项，但那是另一个话题，先把纯空间版本跑通。
2.4 距离矩阵的批量计算与内存控制
DBSCAN 需要两两距离，N 条轨迹就是 N×N 矩阵。N=500 时约 12.5 万次距离计算，每次内部还有 pdist2，直接双重循环会慢到怀疑人生。做法是先算好上三角，再对称填充：
N = numel(traj); % traj 是 cell 数组，每个元素是一条轨迹 D = zeros(N, N); for i = 1:N-1 for j = i+1:N D(i,j) = hausdorff_quantile(traj{i}, traj{j}, 0.95); D(j,i) = D(i,j); end end参数说明：traj用 cell 存变长轨迹，避免用 NaN 填充的矩阵。N 超过 2000 时内存和耗时都会吃紧，这时要么先做轨迹降采样（Douglas-Peucker 抽稀到 50~100 点），要么改用 KD-tree 加速最近邻查询。我一般先抽稀再算，精度损失很小，速度提升一个数量级。
3. 把距离塞进 DBSCAN：Matlab 实现与参数怎么定
3.1 为什么不用 Matlab 自带 dbscan
Matlab 的 Statistics and Machine Learning Toolbox 里有dbscan(X, epsilon, minpts)，但它只接受数值特征矩阵，内部用欧氏距离。你要用自定义的 Hausdorff 距离，有两条路：一是把距离矩阵转成pdist格式喂给linkage那套层次聚类（但那不是 DBSCAN），二是自己写 DBSCAN 主循环。我选后者，因为逻辑透明、好调参、好加缓存。
DBSCAN 的核心就两个概念：eps 邻域内点数 ≥ minpts 的点是核心点；核心点密度可达的点归入同一簇；不属于任何簇的点是噪声。用预计算的距离矩阵 D，判断邻域就是查 D 的第 i 行哪些列 ≤ eps。
3.2 基于距离矩阵的 DBSCAN 主循环
function [labels, coreIdx] = dbscan_custom(D, eps, minpts) % D: N×N 对称距离矩阵，对角线为 0 % eps: 邻域半径 % minpts: 核心点最小邻居数（含自身） N = size(D, 1); labels = zeros(N, 1); % 0 表示未访问/噪声 visited = false(N, 1); clusterId = 0; coreIdx = false(N, 1); for i = 1:N if visited(i), continue; end visited(i) = true; neighbors = find(D(i, :) <= eps); % 含自身 if numel(neighbors) < minpts labels(i) = -1; % 暂标噪声，后续可能被核心点吸收 continue; end % i 是核心点，开新簇 clusterId = clusterId + 1; coreIdx(i) = true; labels(i) = clusterId; seedSet = neighbors(neighbors ~= i); while ~isempty(seedSet) q = seedSet(1); seedSet(1) = []; if ~visited(q) visited(q) = true; qNeighbors = find(D(q, :) <= eps); if numel(qNeighbors) >= minpts coreIdx(q) = true; % 把新发现的核心点邻居加入种子集 newSeeds = qNeighbors(labels(qNeighbors) <= 0); seedSet = unique([seedSet; newSeeds(:)]); end end if labels(q) <= 0 labels(q) = clusterId; % 边界点或未分类点归入当前簇 end end end end逻辑说明：外层遍历每个未访问点，判断是否核心点。核心点开新簇，用种子集做广度扩展，扩展时对每个未访问邻居再判断核心性，是核心点就把它的邻居也拉进种子集。labels(q) <= 0这个条件保证噪声点（-1）和未分类点（0）都能被正确归入簇，而已经属于其他簇的点不会被抢走。
参数说明：eps是距离阈值，量纲和你的 Hausdorff 距离一致（通常是米或度）；minpts一般取 4~10，航迹聚类里我常用 5。labels输出中 -1 是噪声，正整数是簇编号。
3.3 eps 和 minpts 的确定：k-距离曲线与经验区间
eps 定不好，DBSCAN 要么全是一个簇，要么全是噪声。标准做法是画k-距离曲线：对每个点，取它到第 k 近邻的距离（k = minpts），排序后画曲线，拐点处就是 eps 的合理值。
k = 5; kdist = zeros(N, 1); for i = 1:N sortedD = sort(D(i, :)); kdist(i) = sortedD(k + 1); % 第1个是自身(0)，第k+1个是第k近邻 end kdist = sort(kdist, 'descend'); plot(kdist); grid on; xlabel('轨迹点序号(按k距离降序)'); ylabel(['第' num2str(k) '近邻距离']); title('k-距离曲线，拐点即eps候选');拐点位置用肉眼找，或者用最大曲率法自动定位。船舶航迹里，如果轨迹已经抽稀到 50~100 点、坐标用米制投影，eps 常见落在 200~800 米区间；用经纬度的话要先把度换算成米（1 度纬度约 111km，经度乘 cos(纬度)）。
注意：eps 和距离度量强耦合。你换了分位数 p，距离整体尺度会变，eps 必须重新标定，不能沿用上一组参数。
3.4 聚类结果怎么验证：轮廓系数与人工抽查
无监督聚类没有标签，验证靠内部指标加人工。轮廓系数（Silhouette）在自定义距离下也能算，只要你有距离矩阵：
function s = silhouette_custom(D, labels) % 基于预计算距离矩阵的轮廓系数 valid = labels > 0; D = D(valid, valid); labels = labels(valid); N = numel(labels); s = zeros(N, 1); for i = 1:N same = (labels == labels(i)); same(i) = false; if ~any(same), s(i) = 0; continue; end a = mean(D(i, same)); otherClusters = unique(labels(~same)); b = inf; for c = otherClusters' idx = (labels == c); b = min(b, mean(D(i, idx))); end s(i) = (b - a) / max(a, b); end s = mean(s); end轮廓系数越接近 1 越好，0 附近说明簇间重叠。但别迷信这个数——航迹聚类里，两条平行航路的船形状相似、距离近，轮廓系数天然偏低，这时候要靠人工把簇画到地图上核对。我一般把每个簇的轨迹用不同颜色画出来，叠加航道底图，一眼就能看出聚得对不对。
4. 避坑与排查：航迹聚类里最容易翻车的五件事
4.1 现象：所有轨迹聚成一个簇，eps 怎么调都没用
原因：距离没有归一化，或者坐标单位混用。AIS 原始经纬度下，1 度约 111km，两条相距 500 米的船距离只有 0.0045 度，而 Hausdorff 距离算出来也是这个量级，eps 设 0.01 就把所有船都圈进去了。另一个常见原因是轨迹点数差异太大，长轨迹的最近距离集合被大量小值拉低。
解决：统一投影到米制平面（UTM 或高斯克吕格），或者用 haversine 算球面距离。同时先对轨迹做等间隔重采样，让每条轨迹点数接近，再算 Hausdorff。我习惯在预处理阶段就把经纬度转成局部平面坐标，后面所有距离都是米，eps 直接按米设，心里有数。
4.2 现象：同一航路的船被拆成很多小簇
原因：分位数 p 取太高（接近 1），退化成标准 Hausdorff，个别跳点把距离顶大；或者 minpts 设太大，密度不够的点全成噪声。
解决：把 p 降到 0.9 甚至 0.85，观察簇合并情况；minpts 从 5 降到 3 试试。同时检查数据里有没有明显的 AIS 跳点——速度突变超过 50 节、位置瞬间跳几海里的点，预处理阶段就该剔除。跳点不除，再鲁棒的距离也救不回来。
4.3 现象：程序跑得极慢，N=1000 要跑十几分钟
原因：双重循环里每次调用pdist2都重新分配内存，且 Hausdorff 内部对每条轨迹重复计算最近邻。Matlab 的循环本身不慢，慢在内存分配和重复计算。
解决：先把所有轨迹抽稀到统一长度（比如 64 点），用矩阵化方式批量算。或者对每条轨迹预计算一个 KD-tree（createns），最近邻查询用knnsearch，比pdist2快数倍。另一个技巧是距离矩阵用single精度存，内存减半，精度对聚类足够。
4.4 现象：换了分位数 p，聚类结果完全变样，无法复现
原因：p 改变了距离尺度，eps 没跟着调。距离整体变小，原来的 eps 就显得过大，簇被合并；距离变大则相反。
解决：把 p 和 eps 绑定成一组参数，每次调 p 都重新跑 k-距离曲线定 eps。更稳的做法是先把距离矩阵归一化到 [0,1]（除以最大距离），eps 用归一化后的值（比如 0.1~0.3），这样换 p 时 eps 的相对含义不变。我在项目里会把 (p, eps, minpts) 三个参数一起记录，方便回溯。
4.5 现象：噪声点里混着真正有价值的异常航迹
原因：DBSCAN 把密度不够的点一律标 -1，但有些异常航迹（比如偏离主航道的船）本身就是分析目标，不该被当垃圾扔掉。
解决：不要把 -1 直接丢弃，单独存一份，后续用孤立森林或人工核查。另外可以调低 minpts 让更多点进入簇，把真正的稀疏异常留给下游。聚类是手段不是目的，噪声集合往往藏着最有意思的样本。
5. 进阶技巧：用层次化 eps 和轨迹重采样把聚类质量再抬一档
前面跑通的是单尺度 DBSCAN，但真实航迹数据往往多尺度共存：主干航道上船密集，支流航道船稀疏，用一个全局 eps 必然顾此失彼。我后来常用的做法是层次化 eps：先粗聚类拿到大簇，再对每个大簇内部用更小的 eps 做二次聚类，把同一航道里不同航行模式（上行/下行、快船/慢船）分开。
具体操作：第一轮 eps 取 k-距离曲线拐点的 1.5 倍，minpts=5，拿到粗簇；对每个簇内的轨迹子集重新算距离矩阵，eps 取该子集 k-距离拐点的 0.8 倍，再跑一次dbscan_custom。这样主干道不会被过度拆分，支流也不会被主干道的大 eps 吞掉。代价是参数变多，但每一层的 eps 都有本层数据支撑，比拍脑袋强。
另一个提升点是把轨迹重采样做扎实。原始 AIS 采样间隔从几秒到几分钟不等，直接算 Hausdorff 会让采样密的轨迹主导距离。我一般用弧长等间隔重采样到固定点数：
function P = resample_traj(P, n) % 按弧长等间隔重采样到 n 个点 d = [0; cumsum(sqrt(sum(diff(P).^2, 2)))]; if d(end) == 0, P = repmat(P(1,:), n, 1); return; end qi = linspace(0, d(end), n)'; P = [interp1(d, P(:,1), qi), interp1(d, P(:,2), qi)]; end参数说明：n取 64 或 128，兼顾形状保真和计算量。重采样后所有轨迹等长，Hausdorff 的最近距离集合长度一致，分位数 p 的含义更稳定，聚类结果的可复现性明显变好。
验证改进效果，我习惯做两件事：一是把聚类前后的簇内平均距离和簇间最小距离打出来，看分离度有没有提升；二是固定随机种子（DBSCAN 本身确定性，但重采样和抽稀可能引入随机），跑三次看簇数量波动。波动超过 10% 就说明参数太敏感，得回去调 p 或 eps。
最后说个血泪教训：别一上来就追求复杂距离。我见过太多人（包括早期的自己）在 Hausdorff 上叠了一堆加权、DTW、形状上下文，结果预处理没做好，跳点没除、坐标没投影，再花哨的距离也是白搭。先把数据洗干净、坐标统一、轨迹重采样，用最朴素的分位数 Hausdorff 跑一版基线，再谈改进。这套流程我在几个内河和近海航迹项目里反复用过，稳。希望帮到你。
本文还有配套的精品资源，点击获取
