【数据融合】基于AIS和雷达的多传感器航迹融合matlab代码

1 简介

AIS船舶自动识别系统能提供精确的船舶位置与属性信息,通过信息融合技术能有效弥补雷达测向精度不高的影响。为了对更大范围内的海域进行监视,本文对升空AIS与雷达信息融合技术进行了研究,以提高平台对海的跟踪探测能力。本文主要涉及的工作:1、对升空平台AIS系统出现的问题进行研究。本文介绍了AIS系统的工作原理与SOTDMA协议过程,分析了升空AIS监视范围与出现信号冲突的问题,并且根据信号区别接收原理得到了时隙冲突、检测容量与发现概率之间的关系。分析分区域接收提高发现概率的可行性,并在matlab下仿真验证。2、通过对AIS数据报的编解码原理进行学习,分析了AIS主要的数据报文类型,并且通过软件实现AIS数据包的解析,并在电子海图上显示,为融合中心利用AIS信息提供支撑。3、本文基于信息融合系统的理论知识,分析了数据融合的模型与基本过程。利用拟合的航向信息对航迹平滑,提升海面目标航迹质量。在不同的测量噪声环境下,从插值精度与算法耗时方面仿真对比分析了最小二乘与拉格朗日插值方法的优缺点。推导存在测量噪声的基于AIS的三维空间配准算法,利用kalman滤波对系统误差进行估计,仿真分析了该算法配准精度与收敛速度。提出了基于DS理论的灰色航迹关联方法,该方法通过利用冗余的历史数据的组合证据进行判决,提高了正确关联概率。通过将其与统计双门限、模糊双门限关联算法在正确关联率、运行效率、适用环境等方面进行对比,为操作员选择算法提供参考。

2 部分代码

function [gmm,variance,likelihood,isactive] = PHD_update(gmm,z,isactive,cst)

% PHD_UPDATE executes the update step of a GM-PHD filter on a given
%           Gaussian Mixture.
%
% IN       gmm       input mixture
%           cst       the constant structure
%
% OUT       gmm       updated mixture
%
% AUTHOR   Isabel Schlangen, (c) 2016

Nmeas = size(z,1);                % number of measurements
active_pred = find(isactive);     % collect all components

gmmtable = repmat(gmm(active_pred),1,Nmeas);

%% update all active components
normtable = cst.pFA*ones(Nmeas,1);
% ASSOCIATIONS
for jj=1:Nmeas
   for ii=1:length(active_pred)
       % innovation:
       y = z(jj,:)'-cst.H*gmm(active_pred(ii)).m;
       S = cst.H * gmm(active_pred(ii)).C * cst.H' + cst.R;

       % Kalman gain:
       K = (gmm(active_pred(ii)).C * cst.H')/S;
       
       gate = y'/S*y;
       
       % discard the combination if we want gating and m is outside of the
       % threshold around z:
       if ~(cst.gating && gate>cst.tgate) 
           if cst.Nz == 2
               detS = S(1,1)*S(2,2)-S(2,1)*S(1,2);
           else
               detS = det(S);
           end        
           q = exp(-0.5 * gate)/sqrt((2*pi)^(cst.Nz)*abs(detS));

           %update the Gaussian:
           gmmtable(ii,jj).w = q * cst.pD * gmm(active_pred(ii)).w;
           gmmtable(ii,jj).m = gmm(active_pred(ii)).m + K*y;
           gmmtable(ii,jj).C = (eye(cst.Nx)-K*cst.H) * gmm(active_pred(ii)).C;
           gmmtable(ii,jj).i = 1;
       else
           gmmtable(ii,jj).w = 0;
           gmmtable(ii,jj).i = 0;
       end
   end
   % normalise per measurement
   normtable(jj) = cst.pFA + sum([gmmtable(:,jj).w]);
   for ii=1:length(active_pred)
       gmmtable(ii,jj).w = gmmtable(ii,jj).w/normtable(jj);
       if gmmtable(ii,jj).w > 1
           keyboard;
       end
   end
end
if sum(sum([gmmtable(:).w]))>Nmeas
   keyboard;
end

% MISSED DETECTIONS
mol_mdterm = 0;
for ii=active_pred
   mol_mdterm = mol_mdterm + cst.pD*gmm(ii).w; %for the likelihood
   gmm(ii).w = (1-cst.pD)*gmm(ii).w;           %for weight update
end

likelihood = exp(-mol_mdterm)*prod(normtable);

%% variance
allweights =  zeros(size(gmmtable));
for ii=1:size(gmmtable,1)
  allweights(ii,:) = [gmmtable(ii,:).w]; 
end
variance = sum([gmm(ii).w]) + sum(sum(allweights.*(1-allweights)));

%% save all components that are still active for the next iteration
inactive = find(~isactive);  % collect indices of free space in gmm
activateind = 0;
for ii=1:length(active_pred)
   for jj=1:Nmeas
       if gmmtable(ii,jj).w > cst.pruning*cst.tprune
           activateind = activateind+1;
           if activateind>length(inactive)
               error('Maximum number of Gauss components reached.\n');            
           end
           gmm(inactive(activateind)) = gmmtable(ii,jj);
           isactive(inactive(activateind)) = 1;
       end
   end
end

%% Pruning
active_up = find(isactive);
prunedweights = 0;
for ii=active_up
   if gmm(ii).w<cst.pruning*cst.tprune
       prunedweights = prunedweights + gmm(ii).w;
       gmm(ii).w = 0;
       gmm(ii).i = 0;
       isactive(ii) = 0;
   end
end
active_pruned = find(isactive);
ncomp = length(active_pruned);
for ii=active_pruned
   gmm(ii).w = gmm(ii).w + prunedweights/ncomp;
end

% fprintf('MD: %g AS: %g, pruned:%g, total:%g\n',sum([gmm(active_pred).w]),sum([gmm(setdiff(active_up,active_pred)).w]),prunedweights,sum([gmm([gmm.i]==1).w]));

%% Merging
[gmm,isactive] = merging(gmm,isactive,cst);

% fprintf('total weight after merging: %g\n', sum([gmm([gmm.i]==1).w]));

3 仿真结果

4 参考文献

[1]陈芳淮. 升空AIS与雷达航迹信息融合技术研究. Diss. 电子科技大学.