红外图像弱小目标

红外图像弱小目标PF-TBD算法源程序

% 将粒子滤波算法应用于红外弱小目标TBD问题，验证其检测、跟踪目标的有效性

clear;clc

% 粒子数目

N = 2000;

% 采样时间

T = 1;

% 仿真结束时间(采样总帧数)

T_end = 30;

% 假定目标从某一特定帧开始出现，然后在另一特定帧消失

T_ap = 6;

T_dp = 24;

% 采样时间序列

SimTime = zeros(floor(T_end/T),1);

% 分辨单元数目

N_x = 32; % 横向分布单元数目

M_y = 32; % 纵向分布单元数目

% 分辨单元的宽度

Delta_X = 1;

Delta_Y = 1;

% 传感器的模糊参数值

SIGMA = 0.7;

% 目标初始出现概率

mu = 0.05;

% 目标速度区间

vmin = 0.2;

vmax = 1;

% 目标强度(灰度值)区间

Imin = 10;

Imax = 30;

% 抽样阈值(在大于r_th的区域内均匀分布)

r_th = 2.5;

% 扩散因子(目标影响相邻分辨单元的程度)

p = 3;

% 目标Markov 过程转移概率相关参数

Pb = 0.05;

Pd = 0.05;

% 转移矩阵的表达式

PI_T = [1-Pb,Pb

Pd,1-Pd];

% 转移矩阵PI 的行数(列数)

PI_s = size(PI_T,1);

% 系统状态转移矩阵

Phi = [1,T,0,0,0

0,1,0,0,0

0,0,1,T,0

0,0,0,1,0

0,0,0,0,1];

% 系统噪声协方差矩阵中的目标状态和灰度幅值噪声强度

q1 = 0.001; % 目标状态变化强度

q2 = 0.01; % 目标灰度值变化强度

% 系统噪声协方差矩阵

Q = [q1*T^3/3,q1*T^2/2,0,0,0

q1*T^2/2,q1*T, 0,0,0

0,0,q1*T^3/3,q1*T^2/2,0

0,0,q1*T^2/2, q1*T, 0

0,0,0,0,q2*T];

% 系统观测噪声

R = 1.5^2;

%*************************************

% 变量取值初始化过程

%*************************************

% 定义滤波初值(假定目标出现时的初值)

X = [4.2,0.45,7.2,0.25,20]';

% 整个粒子的集合，其中第6位代表当前时刻E判决标识位，而第7位为上一时刻值

X_PF = zeros(7,N);

% 单个粒子状态值

X_PF_i = zeros(7,1);

% 初始时假定每个粒子的权值为均匀的

w_i = 1/N*zeros(1,N);

% 预测粒子的均值及其协方差

X_PF_mean = zeros(5,1);

Px_PF = Q;

Xtru = zeros(5,floor(T_end/T));

Xest_PF = zeros(5,floor(T_end/T));

% 经PF估计出的目标出现概率

Prob_PF = zeros(1,floor(T_end/T));

% 每帧红外图像的观测值(以每个像素点为基准)

Measure_i = zeros(N_x,M_y);

H_i = zeros(N_x,M_y);

% 为了计算各粒子的权值，生成初始时刻的红外观测图像

for i=1:N_x

for j = 1:M_y

Measure_i(i,j) = sqrt(R)*randn(1);

end

end

% 当考虑不设置阈值的情况时，此时假定目标会在所有的观测单元中随机出现；

% 当考虑阈值时，在大于门限的区域内均匀分布(如何实现)

% 由于假定目标初始出现概率为0.05，则可认为只有100个粒子的E值为1，其余均为0

for i=1:20:2000

% 假定每20个粒子出现一个E0=1的粒子

X_PF_i = zeros(7,1);

% 在测量值超过门限的区域中均匀地抽取粒子，首先找出超过门限的区域

[PosRow,PosCol] = find(Measure_i>r_th);

NPoints = length(PosRow);

TarPos = ceil(rand(1)*NPoints);

% 抽取粒子

X_PF_i(1) = PosRow(TarPos);

X_PF_i(2) = vmin + (vmax-vmin)*rand(1);

X_PF_i(3) = PosCol(TarPos);

X_PF_i(4) = vmin + (vmax-vmin)*rand(1);

X_PF_i(5) = Imin + (Imax-Imin)*rand(1);

X_PF_i(6) = 1;

X_PF_i(7) = 0;% 保存的是上一时刻的E判决值

X_PF(:,i) = X_PF_i;

end