1. 程式人生 > 其它 >【運動學】基於matlab Singer模型演算法機動目標跟蹤【含Matlab原始碼 1157期】

【運動學】基於matlab Singer模型演算法機動目標跟蹤【含Matlab原始碼 1157期】

一、簡介

目標跟蹤是利用各種型別的感測器獲得的關於目標資訊,對目標的真實狀態以及未來的狀態進行估計和預測的一門技術。目標跟蹤技術在軍事以及民用等諸多領域裡都有著廣泛的應用。隨著電子技術以及計算機技術的進步,各種新的技術與理論被應用到目標跟蹤領域中來,目標跟蹤技術也逐漸的發展成為一門跨學科,跨行業,多層面的技術。目標跟蹤技術研究的主要內容是從感測器中獲得的關於目標的不精確的資訊中,準確的估計和預測到目標的真實資訊,因此需要各種濾波器對收集到的資料進行濾波處理。卡爾曼(Kalman)濾波作為一種效能優良的濾波演算法,在目標跟蹤領域裡得到廣泛的應用,但是濾波器是以一定的目標跟蹤的模型為基礎的。因此目標跟蹤的研究物件也主要包括跟蹤模型和濾波演算法這兩個方面。在這兩個方面,國內外的眾多學者對此進行了深入的研究,並取得了較為豐碩的成果。作為一種新型的資料融合演算法,互動式多模型(IMM)演算法由於其優良的跟蹤效果,較寬的跟蹤頻帶,在最近幾年得到了足夠的重視

二、原始碼

%**********利用Singer模型演算法對機動目標進行跟蹤*************
function [xx5,yy5,ex5,exv5]=singer(T,r,N)
clc
clear
close all
% %***************模擬條件********************
T=2;      %雷達掃描週期
r=10000;  %量測誤差方差
N=50;%Monte Carlo模擬次數
%alpha=1/60;%機動時間常數的倒數,即機動頻率
F=[1 T (1/2)*T^2 0 0 0;
   0 1 T  0 0 0;
   0 0 1  0 0 0;
   0 0 0  1 T (1/2)*T^2; 
   0 0 0  0 1 T;
   0 0 0  0 0 1];%狀態轉移矩陣
H=[1 0 0 0 0 0;
   0 0 0 1 0 0];%量測矩陣
sigmax=r;%X方向量測噪聲方差
sigmay=r;%Y方向量測噪聲方差
R=[sigmax 0;
   0 sigmay];%量測噪聲協方差
%sigmaax=0.01;%X方向目標加速度方差
%sigmaay=0.01;%Y方向目標加速度方差
qq11=T^5/20;
qq12=T^4/8;
qq13=T^3/6;
qq22=T^3/3;
qq23=T^2/2;
qq33=T;
qq44=T^5/20;
qq45=T^4/8;
qq46=T^3/6;
qq55=T^3/3;
qq56=T^2/2;
qq66=T;
Q=[qq11 qq12 qq13 0    0    0;
   qq12 qq22 qq23 0    0    0;
   qq13 qq23 qq33 0    0    0;
   0    0    0    qq44 qq45 qq46;
   0    0    0    qq45 qq55 qq56;
   0    0    0    qq46 qq56 qq66];%過程噪聲協方差
for j=1:N
    [x,y,zx,zy,NN]=target_movement;
    load target_movement_out
    z=[zx';zy'];
    X=[z(1,3) (z(1,3)-z(1,2))/T (z(1,3)-2*z(1,2)+z(1,1))/T^2 z(2,3) (z(2,3)-z(2,2))/T (z(2,3)-2*z(2,2)+z(2,1))/T^2]';%狀態向量初始化
    %濾波協方差初始化
    P11=R(1,1);
    P12=R(1,1)/T;
    P13=R(1,1)/T^2;
    P22=2*R(1,1)/T^2;
    P23=3*R(1,1)/T^3;
    P33=6*R(1,1)/T^4;
    P44=R(2,2);
    P45=R(2,2)/T;
    P46=R(2,2)/T^2;
    P55=2*R(2,2)/T^2;
    P56=3*R(2,2)/T^3;
    P66=6*R(2,2)/T^4;
    P=[P11 P12 P13 0   0   0;
        P12 P22 P23 0   0   0;
        P13 P23 P33 0   0   0;
        0   0   0   P44 P45 P46;
        0   0   0   P45 P55 P56;
        0   0   0   P46 P56 P66];
    MX(:,3)=X; 
    EX(j,3)=(X(1)-x(3)).^2;%x方向位置初始方差
    EXv(j,3)=(X(2)-vvx(3)).^2;%x方向速度初始方差
    EY(j,3)=(X(4)-y(3)).^2;%y方向位置初始方差
    EYv(j,3)=(X(5)-vvy(3)).^2;%y方向速度初始方差
    for i=4:NN
        x1=F*X;
        z1=H*x1;
        P1=F*P*F'+Q;
        S=H*P1*H'+R;
        v=z(:,i)-z1;
        W=P1*H'*inv(S);
        X=x1+W*v;
        P=P1-W*S*W';
        Mv=v'*inv(S)*v;
        MX(:,i)=X;
        MEX(:,i,j)=MX(:,i);
        EX(j,i)=(X(1)-x(i)).^2;%x方向位置初始方差
        EXv(j,i)=(X(2)-vvx(i)).^2;%x方向速度初始方差
        EY(j,i)=(X(4)-y(i)).^2;%y方向位置初始方差
        EYv(j,i)=(X(5)-vvy(i)).^2;%x方向速度初始方差        
    end
end
function [x,y,zx,zy,NN]=target_movement
%函式定義:生成目標運動的真實值和測量值
% %***************模擬條件*******************************************************
T=2;      %雷達掃描週期
r=10000;  %量測誤差方差
x0=2000;%目標在X軸方向的起始位置
y0=10000;%目標在Y軸方向的起始位置
xv0=0;%目標在X軸方向的起始速度
yv0=-15;%目標在Y軸方向的起始速度
NN=500;%取樣點數
x=zeros(NN,1);%X軸位置初始化
y=zeros(NN,1);%Y軸位置初始化
x(1)=x0;%X軸初始位置
y(1)=y0;%Y軸初始位置
vx(1)=xv0;%X軸初始速度
vy(1)=yv0;%Y軸初始速度
for i=1:NN-1
    if i<200
        ax=0;
        ay=0;
        vx(i+1)=vx(i)+ax*T;
        vy(i+1)=vy(i)+ay*T;
    elseif (i>=200)&(i<=300)
        ax=15/200;
        ay=15/200;
        vx(i+1)=vx(i)+ax*T;
        vy(i+1)=vy(i)+ay*T;
    elseif (i>300)&(i<=500)
        ax=0;
        ay=0;
        vx(i+1)=vx(i)+ax*T;
        vy(i+1)=vy(i)+ay*T;     
    end
    x(i+1)=x(i)+vx(i)*T+0.5*ax*T^2+0.5*0*T^2*randn;%X軸的動態方程
    y(i+1)=y(i)+vy(i)*T+0.5*ay*T^2+0.5*0*T^2*randn;%Y軸的動態方程
end
%***************產生量測噪聲********************
nx=100*randn(NN,1);
ny=100*randn(NN,1);
%***************量測值**************************
zx=x+nx;
zy=y+ny;
vvx=vx;
vvy=vy;
save target_movement_out vvx vvy

%i=1:NN;
%k=4:1:NN;
%l=4:1:NN;
%figure(1)
%plot(x,y,'-dm');
%title('目標運動軌跡')
%xlabel('x方向')
%ylabel('y方向')
%legend('目標運動軌跡')

三、執行結果





四、備註

版本:2014a

完整程式碼或代寫加QQ912100926