一、模型简介

1.1问题背景

作者复现了《基于轨迹聚类的船舶异常行为识别研究_胡智辉》的文章内容，对模型感兴趣的朋友可以下载文章深入研究。本文主要是分享复现过程中的算法实现。

1.2具体内容

本模型通过处理AIS数据得到船舶轨迹，通过计算改进的 Hausdorff 距离，对轨迹进行DNSCAN聚类，得到不同类型的航迹簇。主要包含以下内容

AIS数据的预处理

船舶轨迹分割

船舶轨迹相似度度量

船舶轨迹表达方式

船舶轨迹相似度量方法

改进的 Hausdorff 距离

船舶轨迹聚类及轨迹提取

基于改进DBSCAN算法轨迹聚类

船舶典型轨迹的提取

船舶轨迹聚类效果对比

二、代码分享

具体的代码原理比较繁琐，笔者时间精力有限，难以详尽描述，有需要的可以私信咨询。这里就直接分享部分源码

2.1主函数

% clc
% close all
% clear all
% % %%
% % global noisenum
% % % %% 数据路径
% % % path='AIS1246-1.CSV';
% % % %% 计算H距离
% % % [data1,px,tempMapDis,tempMapV,tempMapHeading]=myHuasdoff(path);
% % % Hdist = (tempMapDis + tempMapV + tempMapHeading) / 3;
% % % %% 取距离中的较小值
% % % for i = 1 : size(Hdist,1)
% % %     for j = i + 1 : size(Hdist,2)
% % %         max1 = Hdist(i,j);
% % %         if max1 < Hdist(j,i)
% % %             max1 = Hdist(j,i);
% % %         end
% % %         Hdist(i,j) = max1;
% % %         Hdist(j,i) = max1;
% % %     end
% % % end
% load data.mat
% %% 设置聚类参数
% eps=0.056;
% MinNum=18;
% %% DBscan 聚类
% cluster_label=DBSCAN(Hdist,eps,MinNum);
% %% 计算聚类数
% cluster_num=max(cluster_label);
% %% 计算轨迹平均位置
% meanpos=zeros(length(cluster_label),2);
% for i=1:length(cluster_label)
%     meanpos(i,1)=mean(data1(data1(:,1)==i,2));
%     meanpos(i,2)=mean(data1(data1(:,1)==i,3));
% end
% %% 画图
% % %% plot trace
% figure
% Plottrace(data1,cluster_label)
% title(sprintf(['DBSCAN聚类结果（航迹示意图） \n邻域距离 = ' num2str(eps) ', 最小样本数 = ' num2str(MinNum) '\n航迹类别:' num2str(cluster_num) '类,噪声航迹:' num2str(noisenum) '条']));
% xlabel('经度')
% ylabel('纬度')
% %% plot meanpos
% % figure
% % Plotmeanpos(meanpos, cluster_label);
% % title(sprintf(['DBSCAN聚类结果（航迹平均坐标示意图） \n邻域距离 = ' num2str(eps) ', 最小样本数 = ' num2str(MinNum) '\n航迹类别:' num2str(cluster_num) '类,噪声航迹:' num2str(noisenum) '条']));
% % xlabel('经度')
% % ylabel('维度')
% % 逐类计数典型轨迹
% filename='典型轨迹数据.xlsx';
% DXcolor={'black','magenta'};
% DXindexlist=[15,36];
% Legendslist={'主航道出口船舶典型轨迹','主航道进口船舶典型轨迹'};
% 
% for k=1:cluster_num
%     data_ind=[];
%     temp_trace_ind=find(cluster_label==k);
%     for t=1:length(temp_trace_ind)
%         data_ind1=find(data1(:,1)==temp_trace_ind(t));
%         data_ind=[data_ind;data_ind1];
%     end
%     %     Legendsname = ['第' num2str(k) '类航迹的典型航迹'];
%     Legendsname=Legendslist{k};
%     [DXfinalData{k},DXdatanew{k},DXpx{k},DXfData{k},DXMapDis{k},DXMapV{k},DXMapHeading{k}]=DXGJ(data1(data_ind,:),DXcolor{k},Legendsname,DXindexlist(k));
%     datawriter(filename,DXfData{k},Legendsname)% end
% end
% 
% % %%
% % addpath '识别数据集'
% % filename_train_in1='进口正常行为数据.xlsx';
% % filename_train_in0='进口异常行为数据.xlsx';
% addpath '正常与异常数据'
% filename_train_in1='进口正常轨迹数据.xlsx';
% filename_train_in0='进口异常轨迹数据.xlsx';
% [data_in1,finalData_in1]=makeData(filename_train_in1);
% [data_in0,finalData_in0]=makeData(filename_train_in0);
% stepNum=20;
% DistThre=linspace(1,250,stepNum);
% SogThre=linspace(1,8,stepNum);
% CogThre=linspace(1,8,stepNum);
% %
% % DistThre=1:0.01:250;
% DXTrace=DXfData{2};
% % ACC_P=[];
% % ACC_S=[];
% % ACC_C=[];
% for i=1:stepNum
%     tDistThre=DistThre(i);
%     tSogThre=SogThre(i);
%     tCogThre=CogThre(i);
%     [ACC_P(i),ACC_S(i),ACC_C(i)]=drawErrerRate(finalData_in1,finalData_in0,DXTrace,tDistThre,tSogThre,tCogThre);
% end
% %% 
filename='阈值.xlsx';
data=xlsread(filename);
figure
plot(DistThre,[ACC_P.r],'-gp','DisPlayName','准确率')
hold on
plot(DistThre,1-[ACC_P.r],'-rd','DisPlayName','误警率')
xlabel('位置偏差')
legend
figure
plot(SogThre,[ACC_S.r],'-gp','DisPlayName','准确率')
hold on
plot(SogThre,1-[ACC_S.r],'-rd','DisPlayName','误警率')
xlabel('航速偏差')
legend
figure
plot(CogThre,[ACC_C.r],'-gp','DisPlayName','准确率')
hold on
plot(CogThre,1-[ACC_C.r],'-rd','DisPlayName','误警率')
xlabel('航向偏差')
legend
%% 2类：进口 1类出口
global DistThre SogThre CogThre
DistThre=60;
SogThre=2.8;
CogThre=3.8;
addpath '正常与异常数据'
filename_in='验证进口道船舶异常行为.xlsx';
tracetype=input('请选择要验证的船舶类型：1（出口）/ 2（进口）');
try
DXTrace=DXfData{tracetype};
catch
    disp('输入错误')
    return
end
checkTrace(filename_in,DXTrace);```

2.1 计算H距离

function [data1,px,tempMapDis,tempMapV,tempMapHeading]=myHuasdoff(datapath)
data1 = xlsread(datapath); 
[sy,sx] = size(data1);
tempx = 1;tempy = 1;
global finalData;
finalData{1,1} = data1(1,:);
px(1,1) = 1;
px(1,2) = data1(1,1);
for i = 2 : sy
    if data1(i,1) ~= data1(i - 1,1)
        tempx = tempx + 1;
        px(tempx,1) = tempx;px(tempx,2) = data1(i,1);
        finalData{tempx,1} = data1(i,:);
        tempy = 1;
    else
        tempy = tempy + 1;
        finalData{tempx,tempy} = data1(i,:);
    end
end
global n;
global yy;
xx = size(finalData,1);
yy = size(finalData,2);
n = xx;
tempMapDis = zeros(n,n);
tempMapV = zeros(n,n);
tempMapHeading = zeros(n,n);
%%计算最小值
for i = 1 : xx
    for j = 1 : xx
       ......
    end
end
%%归一化
for i = 1 : n
    for j = 1 : n
... ...
        end
    end
end
disGy = disGy - disMin;
vGy = vGy - vMin;
headingGy = headingGy - headingMin;
for i = 1 : n
    for j = 1 : n
 ......
    end
end

end

2.2 DBSCAN聚类

function [labels, isnoise]=DBSCAN(D,epsilon,MinPts)
    C=0;
    n=size(D,1);
    labels=zeros(n,1);
%     D=pdist2(X,X);
%     D=X;
    visited=false(n,1);
    isnoise=false(n,1);
    for i=1:n
        if ~visited(i)
            visited(i)=true;
            
            Neighbors=RegionQuery(i);
            if numel(Neighbors)<MinPts
                % X(i,:) is NOISE
                isnoise(i)=true;
            else
                C=C+1;
                ExpandCluster(i,Neighbors,C);
            end
            
        end
    
    end
    
    function ExpandCluster(i,Neighbors,C)
        labels(i)=C;
        
        k = 1;
        while true
            j = Neighbors(k);
            
            if ~visited(j)
                visited(j)=true;
                Neighbors2=RegionQuery(j);
                if numel(Neighbors2)>=MinPts
                    Neighbors=[Neighbors Neighbors2];   %#ok
                end
            end
            if labels(j)==0
                labels(j)=C;
            end
            
            k = k + 1;
            if k > numel(Neighbors)
                break;
            end
        end
    end
    
    function Neighbors=RegionQuery(i)
    end

2.3提取典型轨迹

%求同一簇典型轨迹
function [finalData1,data1,px,fData,tempMapDis,tempMapV,tempMapHeading]=DXGJ(data1,DXcolor,DXlegend,DXindex)
% function main()
%%数据预处理
% data1 = xlsread(path); %xls xlsx
%    data1 = csvread('1.csv'); %csv
[sy,sx] = size(data1);
tempx = 1;tempy = 1;
global finalData;
finalData=[];
finalData{1,1} = data1(1,:);
px(1,1) = 1;
px(1,2) = data1(1,1);
for i = 2 : sy
    if data1(i,1) ~= data1(i - 1,1)
        tempx = tempx + 1;
        px(tempx,1) = tempx;px(tempx,2) = data1(i, 1);
        finalData{tempx,1} = data1(i,:);
        tempy = 1;
    else
        tempy = tempy + 1;
        finalData{tempx,tempy} = data1(i,:);
    end
end

finalData1 = finalData;
%画出船舶轨迹线
% for i = 1 : length(data1(:,1)) - 1
%     if data1(i,1) == data1(i + 1,1)
%         line([data1(i,2) data1(i + 1,2)],[data1(i,3) data1(i + 1,3)]);
%         hold on
%     end
% end
% hold on

global a;
a=DXindex;
global xx;
global yy;
xx = size(finalData,1);
yy = size(finalData,2);
tempX = [];
tempY = [];
tempH = [];
tempV = [];

%%计算交点
for k = 1 : yy
    if isempty(finalData{a,k})
        break;
    end
    [t1,t2,t3,t4] = toIntersection(k);
    if t1 == -1
        continue;
    end
    tempX = [tempX,t1];
    tempY = [tempY,t2];
    tempH = [tempH,t3];
    tempV = [tempV,t4];
end
hold on;
plot(tempX,tempY,'Color',DXcolor,'LineWidth',1,'DisPlayName',DXlegend) %典型轨迹的划线
hold on;
fData ={};
for i = 1 : size(tempX,2)
    temp = [tempX(i),tempY(i),tempH(i),tempV(i)];
    fData{1,i} = temp;
    temp = [1,temp];
    finalData{xx + 1,i} = temp;
end
xx = size(finalData,1);
yy = size(finalData,2);
n = xx;
tempMapDis = zeros(1,n - 1);
tempMapV = zeros(1,n - 1);
tempMapHeading = zeros(1,n - 1);

%%计算最小值
for j = 1 : xx - 1
    totalDis = 0;
    totalV = 0;
    totalHeading = 0;
    len = 0;
    for k = 1 : yy
        if isempty(finalData{xx,k})
            break;
        end
        min1 = 999999999;
        minl = -1;
        for l = 1 : yy
            if isempty(finalData{j,l})
                break;
            end
            if min1 > osDis(xx,k,j,l)
                min1 = osDis(xx,k,j,l);
                minl = l;
            end
        end
        [t1,t2,t3] = toMinDis(xx,k,j,minl);
        t2 = abs(t2 - finalData{xx,k}(4));
        t3 = abs(t3 - finalData{xx,k}(5));
        totalDis = totalDis + t1;
        totalHeading = totalHeading + t2;
        totalV = totalV + t3;
        len = len + 1;
    end
    tempMapDis(1,j) = totalDis / len;
    tempMapV(1,j) = totalV / len;
    tempMapHeading(1,j) = totalHeading / len;
end
end

2.4 根据聚类结果 进行航迹预测

function [ACC_P,ACC_S,ACC_C]=drawErrerRate(finalData_in1,finalData_in0,DXTrace,DistThre,SogThre,CogThre)
% DistThre=60;
% SogThre=2.8;
% CogThre=3.8;
% filename_train_in1='进口正常轨迹数据.xlsx';
% filename_train_in0='进口异常轨迹数据.xlsx';
% [data_in1,finalData_in1]=makeData(filename_train_in1);
% [data_in0,finalData_in0]=makeData(filename_train_in0);
% pred_label_in1_P=zeros(length(finalData_in1),1);
% pred_label_in0_P=zeros(length(finalData_in0),1);
% pred_label_in1_S=zeros(length(finalData_in1),1);
% pred_label_in0_S=zeros(length(finalData_in0),1);
% pred_label_in1_C=zeros(length(finalData_in1),1);
% pred_label_in0_C=zeros(length(finalData_in0),1);
TP_P=0;
FP_P=0;
TN_P=0;
FN_P=0;

TP_S=0;
FP_S=0;
TN_S=0;
FN_S=0;

TP_C=0;
FP_C=0;
TN_C=0;
FN_C=0;
for ind=1:size(finalData_in1,1)
    currentTrace=finalData_in1(ind,:);

    % DXTrace=DXfData{t};
    [DevPosition,DevSog,DevCog]=getDev(currentTrace,DXTrace);
        pred_label_in_P=zeros(length(DevPosition),1);
    pred_label_in_S=zeros(length(DevSog),1);
    pred_label_in_C=zeros(length(DevCog),1);
    for node=1:length(DevPosition)
        if DevPosition(node)<DistThre
            pred_label_in_P(node)=1;
        end
        if DevSog(node)<SogThre
            pred_label_in_S(node)=1;
        end
        if DevCog(node)<CogThre
            pred_label_in_C(node)=1;
        end
    end
    %% 找到偏离点
    errerind_P=find(pred_label_in_P==0);
    diff1_P=diff(diff(errerind_P));
    if isempty(find(diff1_P==0))
        TP_P=TP_P+1;
%         max(DevPosition)
    else
        FN_P=FN_P+1;
    end
    %%
    %% 找到偏离点
    errerind_S=find(pred_label_in_S==0);
    diff1_S=diff(diff(errerind_S));
    if isempty(find(diff1_S==0))
        TP_S=TP_S+1;
    else
        FN_S=FN_S+1;
    end
    %%
    %% 找到偏离点
    errerind_C=find(pred_label_in_C==0);
    diff1_C=diff(diff(errerind_C));
    if isempty(find(diff1_C==0))
        TP_C=TP_C+1;
    else
        FN_C=FN_C+1;
    end
end

%%
for ind=1:size(finalData_in0,1)
    currentTrace=finalData_in0(ind,:);

    % DXTrace=DXfData{t};
    [DevPosition,DevSog,DevCog]=getDev(currentTrace,DXTrace);
        pred_label_in_P=zeros(length(DevPosition),1);
    pred_label_in_S=zeros(length(DevSog),1);
    pred_label_in_C=zeros(length(DevCog),1);
    for node=1:length(DevPosition)
        if DevPosition(node)<DistThre
            pred_label_in_P(node)=1;
        end
        if DevSog(node)<SogThre
            pred_label_in_S(node)=1;
        end
        if DevCog(node)<CogThre
            pred_label_in_C(node)=1;
        end
    end
    %% 找到偏离点
    errerind_P=find(pred_label_in_P==0);
    diff1_P=diff(diff(errerind_P));
    if isempty(find(diff1_P==0))
        FP_P=FP_P+1;
    else
        TN_P=TN_P+1;
    end
    %%
    %% 找到偏离点
    errerind_S=find(pred_label_in_S==0);
    diff1_S=diff(diff(errerind_S));
    if isempty(find(diff1_S==0))
        FP_S=FP_S+1;
    else
        TN_S=TN_S+1;
    end
    %%
    %% 找到偏离点
    errerind_C=find(pred_label_in_C==0);
    diff1_C=diff(diff(errerind_C));
    if isempty(find(diff1_C==0))
        FP_C=FP_C+1;
    else
        TN_C=TN_C+1;
    end
end

ACC_P.r=(TP_P+TN_P)/(TP_P+TN_P+FP_P+FN_P);
ACC_S.r=(TP_S+TN_S)/(TP_S+TN_S+FP_S+FN_S);
ACC_C.r=(TP_C+TN_C)/(TP_C+TN_C+FP_C+FN_C);
ACC_P.TP=TP_P;
ACC_P.TN=TN_P;
ACC_P.FP=FP_P;
ACC_P.FN=FN_P;

ACC_S.TP=TP_S;
ACC_S.TN=TN_S;
ACC_S.FP=FP_S;
ACC_S.FN=FN_S;

ACC_C.TP=TP_C;
ACC_C.TN=TN_C;
ACC_C.FP=FP_C;
ACC_C.FN=FN_C;

三、结果展示

本文的过程相对来说比较复杂，结果也比较多，主要将后期可视化的结果进行展示，完整的数据处理过程和相关结果见第四部分代码分享。

3.1DBSCAN聚类结果

3.2 典型航迹提取

3.2 航迹阈值分类

3.2 最佳分类阈值

四、代码分享

源码链接：C币下载
需要免费分享的关注+点赞+收藏的 私聊领取
完整的项目内容
包含：
模型文件和航迹数据
DBSCAN代码
H距离计算代码
航迹提取代码
基于H距离的阈值分类代码
阈值寻优代码
均可直接运行
可以结合这个案例，学习DBSCAN距离算法和H距离的计算方法，掌握了航迹分类的大致流程，取得了较好的聚类效果，并基于航迹聚类进行航迹偏离预测，经检测准确率较高。
需要学习DBSCAN的朋友可以当作一个应用案例来学习，船舶航迹研究方向的朋友可以通过替换数据和相关模块实现拓展成自己的模型实现

参考文献
[1]:基于轨迹聚类的船舶异常行为识别研究
胡智辉
[2]:Typical trajectory extraction method for ships based on AIS data
and trajectory clustering
Huaipeng Li∗
Beijing Industry University, Beijing, China