matlab_检测前跟踪（TBD），通过多帧回波数据积 累和联合处理，可以显著提高雷达的微弱目标检测跟踪性能

% DP-TBD for single target tracking % %假设过程噪声和量测噪声均为0均值高斯白噪声 % 作者：郭振 % 时间：2014/3/19 %% clear all clear clc close all % % rand('state',sum(100*clock)); % Shuffle the pack! % randn('state',sum(100*clock)); % Shuffle the pack! %% step1 initialization K=10; %帧数 ss = 5; % state size X=[x,dot_x, y, dot_y, z] os = 3; % observation size Y=[x ,y ,I] X=zeros(ss,K); %状态 I=zeros(K); %各帧状态对应的值函数 S=zeros(K); %存储各帧之间的状态转移关系 Nx=100; %径向距离分辨单元 Ny=100; %切向距离分辨单元 DeltaX=200; %设一个分辨单元大小为20 DeltaY=200; %% simulation condition T_step=2.5; F = [1 T_step 0 0 0; 0 1 0 0 0; 0 0 1 T_step 0; 0 0 0 1 0;0 0 0 0 1]; % 目标状态转移矩阵 Q_CV=[T_step^3/3 T_step^2/2 0 0 0 ;T_step^2/2 T_step 0 0 0;... % 过程噪声方差矩阵 0 0 T_step^3/3 T_step^2/2 0;0 0 T_step^2/2 T_step 0; 0 0 0 0 1/T_step]; q_CV=10; % 连续过程噪声的能量谱密度 q_I=10; % 目标回波信号能量变换速度的功率谱密度 Q=Q_CV*diag([q_CV q_CV q_CV q_CV q_I]); Sigma=Q; %状态为5维时的过程噪声 FFourDim=F(1:4,1:4); processNoise =sample_gaussian(zeros(length(Q),1),Q,K);% Note ' here! state_noise_samples is: ss-Total_time Mach=340; initx = [20000 -1.5*Mach 20000 -1.6*Mach 20000]'; x(:,1) =initx; % Initial state. for t=2:K x(:,t)=F*x(:,t-1)+processNoise(:,t); % For CV Model. end % figure(1) %目标运动轨迹 % plot(x(1,:),x(3,:)) %% 产生量测（像素点平面数据） Theta=1; % 背景噪声方差 BGNosie=BackgroundNoise(Nx,Ny,K,Theta);% 背景噪声 Data_target=zeros(Nx,Ny,K); SNR=10; Power_noise_av=1; Power_target_av=10^(SNR/10)*Power_noise_av; % The average of Target Power. s=1; % extent of blurring Lx=1; % x方向上的目标影响力的衰减速度 Ly=1; % y方向上的目标影响力的衰减速度 for t=1:K IndexX=ceil(x(1,t)/DeltaX); IndexY=ceil(x(3,t)/DeltaY); Buf=Power_target_av*exp(1i*rand(1)*2*pi); %Yiwei' definition Data_target(IndexX,IndexY,t)=Data_target(IndexX,IndexY,K)+Buf; % for i=(IndexX-1):(IndexX+1) % for j=(IndexY-1):(IndexY+1) % if(0<=i)&(i<=100)&(0<=j)&(j<=100) %Salmond's definition % Buf=DeltaX^2*Power_target_av/2/pi/(s)^2*... % exp((-(i*DeltaX-x(1,t))^2-(j*DeltaY-x(3,t))^2)/2/(s)^2) % Data_target(i,j,t)=Data_target(i,j,t)+Buf; % Yiwei's definition % Buf=Power_target_av*exp(1i*rand(1)*2*pi)*... % exp((-(i*DeltaX+0.5-x(1,t))^2*Lx/2/DeltaX-(j*DeltaY+0.5-x(3,t))^2*Ly/2/DeltaY)); % Data_target(i,j,t)=Data_target(i,j,t)+Buf; % end % end % end end DataScan=BGNosie+Data_target; % 目标幅度与背景噪声组成量测的幅相信息，量测的位置位于分辨单元的中心位置 Data_Add=abs(DataScan); Data_Add=sum(DataScan,3); % mesh(abs(Data_Add)) % % figure(2) % image(abs(Data_Add)) %% step2 recurrence accumulation % 求使值函数最大的 %% 航迹数据结构 Point.Value=0; Point.LastPointPosition=[0 0]; Point.LastPointIndex=0; Point.State=[0 0 0 0]; Point.Intensity=0; Point.Position=[0 0]; Point.OtherPointSetIndex=0; Track.Step=0; Track.Value=0; %航迹的值函数 % CandidateTrack; % DeclaredTrack; %% 第一帧数据单独处理 % 选取强度超过一定阈值的单元为航迹头单元 DataScanBuf=abs(DataScan(:,:,1)); DetectThreshold=5; %幅度检测阈值 [MaxBufOne RowIndexBuf]=max(DataScanBuf); % 寻找每一列中最大元素 [MaxBufTwo ColumnIndex]=max(MaxBufOne); % 寻找矩阵中最大元素 RowIndex=RowIndexBuf(ColumnIndex); Flag=0; while(MaxBufTwo>DetectThreshold) %将该航迹头数据保存 Point.Value=0; Point.LastPointPosition=[0 0]; Point.State=[(RowIndex-0.5)*DeltaX -1.5*Mach (ColumnIndex-0.5)*DeltaY -1.6*Mach MaxBufTwo]'; Xk=Point.State; IndexX=ceil(Xk(1)/DeltaX); IndexY=ceil(Xk(3)/DeltaY); Point.Position=[IndexX IndexY]; Point.Intensity=DataScan(RowIndex,ColumnIndex,1); Point.OtherPointSetIndex=1; Track.PointSet(1).Point=Point; % 这个变量其实没用 Track.PointSet(1).OtherPointSet(1)=Point; Track.Step=1; if Flag==0 CandidateTrack(1)=Track; Flag=1; else CandidateTrack=[CandidateTrack Track]; end % 下一个航迹头处理 DataScanBuf(RowIndex,ColumnIndex)=0; %将处理过的最大元素值赋0 [MaxBufOne RowIndexBuf]=max(DataScanBuf); % 寻找每一列中最大元素 [MaxBufTwo ColumnIndex]=max(MaxBufOne); % 寻找矩阵中最大元素 RowIndex=RowIndexBuf(ColumnIndex); % MaxBufTwo为最大元素的值 end %% 从第二帧数据开始迭代贝叶斯处理 for k=2:K % 目标状态的预测,将预测得到的16个状态值存储在OtherPoint中 %暂时做单目标处理 if (CandidateTrack(1).Step==1) LastPoint=CandidateTrack(1).PointSet(1).Point; Xk=LastPoint.State; [XX] = StatePrediction(F,Xk,T_step); LastPointValue=LastPoint.Value; MaxValue=-10000000; % LastIndexX=ceil(Xk(1)/DeltaX); LastIndexY=ceil(Xk(3)/DeltaY); for i=1:16 %对16个预测状态进行值函数计算 % 计算该状态的值函数 IndexX=ceil(XX(1,i)/DeltaX); IndexY=ceil(XX(3,i)/DeltaY); Xkk=F*Xk; LogPdfxkk_xk=1/((2*pi)^2*sqrt(det(Sigma)))*exp(-0.5*(XX(:,i)-Xkk)'/Sigma*(XX(:,i)-Xkk)); LogPdfzkk_xkk=1/(sqrt(2*pi))*exp(-(abs(DataScan(IndexX,IndexY,k))-Power_target_av)^2/2); %量测噪声的分布，这里的表示有问题!!!!!!!!!!!! %LogPdfzkk_xkk=LogPdfzkk_xkk1*LogPdfzkk_xkk2; Value=log10(LogPdfxkk_xk)+log10(LogPdfzkk_xkk)+LastPointValue; if(Value>MaxValue) MaxValue=Value; MaxValueIndex=i; MaxIndexX=IndexX; MaxIndexY=IndexY; end Point.Value=Value; Point.LastPointPosition=[LastIndexX LastIndexY]; Point.State=XX(:,i); Point.Intensity=DataScan(IndexX,IndexY,k); Point.Position=[IndexX IndexY]; Point.LastPointIndex=1; Point.OtherPointSetIndex=i; CandidateTrack(1).PointSet(k).OtherPointSet(i)=Point; %将点的信息存储 end % 第k时刻的状态的更新 %i=numel(CandidateTrack.PointSet); %%%????? CandidateTrack(1).PointSet(k).Point= CandidateTrack(1).PointSet(k).OtherPointSet(MaxValueIndex); %将最大值赋予状态，这个变量无用，可以记录当前最大状态 CandidateTrack(1).Step=CandidateTrack(1).Step+1; CandidateTrack(1).Value=MaxValue; %% CandidateTrack.Step>1 else CellFlag=zeros(Nx,Ny); %指示该单元内是否有目标预测状态存在 LastCellStorage=cell(Nx,Ny); % 胞元数组 % PredictionCellStorage=cell(1); PredictionPointNumber=0; OtherPointSetNumber=numel(CandidateTrack(1).PointSet(k-1).OtherPointSet); for m=1:OtherPointSetNumber % 对上一状态OtherPointSet中的点（上限为16个点）进行状态预测 LastPoint=CandidateTrack(1).PointSet(k-1).OtherPointSet(m); %将上一状态点取出 Xk=LastPoint.State; %状态预测 [XX] = StatePrediction(F,Xk,T_step); LastIndexX= LastPoint.Position(1); LastIndexY= LastPoint.Position(2); % LastIndexX= ceil(Xk(1)/DeltaX); % LastIndexY= ceil(Xk(3)/DeltaY); for i=1:16 %对预测的16个状态点进行操作 % XX(2,i)=(XX(1,i)-Xk(1))/T_step; %更新x方向速度 % XX(4,i)=(XX(3,i)-Xk(3))/T_step; %更新y方向速度 % %对预测的状态进行存储，避免重叠存储 IndexX=ceil(XX(1,i)/DeltaX); IndexY=ceil(XX(3,i)/DeltaY); CellFlag(IndexX,IndexY)=Ce