雷达无线电系列（二）经典CFAR算法图文解析与实现（matlab）

一，CFAR基础知识介绍

简介

　　恒虚警检测技术是指雷达系统在保持虚警概率恒定条件下对接收机输出的信号与噪声作判别以确定目标信号是否存在的技术。

前提

　　由于接收机输出端中肯定存有噪声(包括大气噪声、人为噪声、内部噪声和杂波等)，而信号一般是叠加在噪声上的。这就需要在接收机输出的噪声或信号加噪声条件下，采用检测技术判别是否有目标信号。

误差概率

　　任何形式的判决必然存在着两种误差概率：发现概率和虚警概率。当接收机输出端存在目标回波信号，而判决时判为有目标的概率为Pd，判为无目标的概率为1-Pad。当接收机输出端只有噪声时，而判为有目标的概率为Pfa。由于噪声是随机变量，其特征可用概率密度函数表示，因此信号加噪声也是一随机变量 [1]

具体过程

　　恒虚警检测器首先对输入的噪声进行处理后确定一个门限，将此门限与输入端信号相比，如输入端信号超过了此门限，则判为有目标，否则，判为无目标。

方法

　　为了使系统的虚警概率保持恒定，此门限是随输入噪声变化而进行快速的自适应调整的，噪声处理方法是随噪声的不同分布而异，因而恒虚警检测技术包括恒虚警处理技术和目标检测技术两大部分。恒虚警处理技术包括快门限、慢门限处理技术、目标检测包括似然比检测、二进制检测、序贯检测和非参量检测等 [1]。

 

二，经典CFAR算法解析与性能图谱展示

 

　　①均值类CFRA

　　核心思想是通过对参考窗内采样数据取平均来估计背景功率。

　　CA-CFAR、GO-CFAR、SO-CFAR算法这三个是最经典的均值类CFAR算法，后续有其变形log-CFAR算法但是原理是相同的。

 

图1 均值类CFAR检测器原理图

 

 

 　　②统计有序CFAR

　　核心思想：通过对参考窗内的数据由小到大排序选取其中第K个数值假设其为杂波背景噪声。

图2 OS-CFAR检测器原理图

 

　　③自适应CFAR

　　针对不同的杂波选用不同的决策方法来进行。

图3 HCE-CFAR检测器原理图

 

图4 VI-CFAR检测器原理图

 

　　④杂波图CFAR

　　经典的空域杂波图CFAR算法，假设每个距离元和方位角上的平方率值遵从独立同分布，通过对时间维上的数据进行遗忘迭代来拟合杂波背景噪声。

　　

图5 杂波图CFAR检测器原理图

 

三，CFRA发展现状

 　　均值类-统计有序类-自适应-杂波图-变换域-压缩感知-神经网络

 

四，经典CFAR算法实现源码（matlab）

 cfar_ac.m　　单元平均恒虚警算法实现

优点：损失率最少的一种算法；

缺点：多目标遮掩，杂波边缘性能也欠佳；

function [ index, XT ] = cfar_ac( xc, N, pro_N, PAD)
%   假设回波服从高斯分布
%   
alpha=N.*(PAD.^(-1./N)-1);
index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
XT=zeros(1,length(index));

for i=index
    cell_left=xc(1,i-N/2-pro_N/2:i-pro_N/2-1);
    cell_right=xc(1,i+pro_N/2+1:i+N/2+pro_N/2);
    Z=(sum(cell_left)+sum(cell_right))./N;
    
    XT(1,i-N/2-pro_N/2)=Z.*alpha;
end

end

test_cfar_ca.m

%% 程序初始化
clc;clear all;close all;

%% 均匀背景噪声（单目标&多目标）
% shape=[200];
% variance=200;
% noise_db=20;
% noise_p=10.^(noise_db./10);
% show_out=0;
% [ xc ] = env_uniform(variance,  shape, noise_db,show_out);

% 单目标
% SNR1=15;    signal1_p=10.^(SNR1./10).*noise_p;
% xc(1,90)=signal1_p;
% 
% 多目标
% SNR1=20;signal1_p=10.^(SNR1./10).*noise_p;
% xc(1,90)=signal1_p;
% SNR2=15;signal2_p=10.^(SNR2./10).*noise_p;
% xc(1,102)=signal2_p;


%% 杂波边缘背景噪声（单目标&多目标）
shape=[100,200];
variance=200;
noise_db=[20,30];
noise_p=10.^(noise_db./10);
show_out=0;
[ xc ] = env_edge(variance,  shape, noise_db,show_out);

% 多目标
SNR1=15;signal1_p=10.^(SNR1./10).*noise_p(1,end);
SNR2=12;signal2_p=10.^(SNR2./10).*noise_p(1,end);
SNR3=8;signal3_p=10.^(SNR3./10).*noise_p(1,end);
SNR4=5;signal4_p=10.^(SNR4./10).*noise_p(1,end);

loc1=randi([43,44],1,1);
xc(1,loc1)=signal1_p;
loc2=randi([46,48],1,1);
xc(1,loc2)=signal3_p;
loc3=randi([50,53],1,1);
xc(1,loc3)=signal2_p;
loc4=randi([55,58],1,1);
xc(1,loc4)=signal1_p;
loc5=randi([90,93],1,1);    %接近杂波区，但是依然在
%     xc(1,loc5)=signal1_p;
xc(1,loc5)=signal2_p;     %这种杂波边缘有目标的情况，还需要优化
%     xc(1,loc5)=signal3_p;

loc6=randi([102,108],1,1);
xc(1,loc6)=signal3_p;

%% 算法结果&图谱显示
N=36;
pro_N=2;
PAD=10^(-4);
[ index, XT ] = cfar_ac( xc, N, pro_N, PAD);

figure;
plot(10.*log(abs(xc))./log(10)),hold on;
plot(index,10.*log(abs(XT))./log(10)),hold on;

图6 CA-CFAR效果图

 

 

cfar_go.m　　最大选择恒虚警算法实现

优点：杂波边缘区域虚警概率降低

缺点：多目标遮掩

function [ index, XT ] = cfar_go( xc, N, pro_N, PAD)
%   假设回波服从高斯分布
%   alpha赋值有些问题，一个比较复杂的高次函数
alpha=N.*(PAD.^(-1./N)-1);

index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
XT=zeros(1,length(index));

for i=index
    cell_left=xc(1,i-N/2-pro_N/2:i-pro_N/2-1);
    cell_right=xc(1,i+pro_N/2+1:i+N/2+pro_N/2);
    Z=max([mean(cell_left),mean(cell_right)]);
    
    XT(1,i-N/2-pro_N/2)=Z.*alpha;
end

end

test_cfar_go.m

%% 程序初始化
clc;clear all;close all;

%% 均匀背景噪声（单目标&多目标）
% shape=[200];
% variance=200;
% noise_db=20;
% noise_p=10.^(noise_db./10);
% show_out=0;
% [ xc ] = env_uniform(variance,  shape, noise_db,show_out);

% 单目标
% SNR1=15;    signal1_p=10.^(SNR1./10).*noise_p;
% xc(1,90)=signal1_p;
% 
% 多目标
% SNR1=20;signal1_p=10.^(SNR1./10).*noise_p;
% xc(1,90)=signal1_p;
% SNR2=15;signal2_p=10.^(SNR2./10).*noise_p;
% xc(1,102)=signal2_p;


%% 杂波边缘背景噪声（单目标&多目标）
shape=[100,200];
variance=200;
noise_db=[20,30];
noise_p=10.^(noise_db./10);
show_out=0;
[ xc ] = env_edge(variance,  shape, noise_db,show_out);

% 多目标
SNR1=15;signal1_p=10.^(SNR1./10).*noise_p(1,end);
SNR2=12;signal2_p=10.^(SNR2./10).*noise_p(1,end);
SNR3=8;signal3_p=10.^(SNR3./10).*noise_p(1,end);
SNR4=5;signal4_p=10.^(SNR4./10).*noise_p(1,end);

loc1=randi([43,44],1,1);
xc(1,loc1)=signal1_p;
loc2=randi([46,48],1,1);
xc(1,loc2)=signal3_p;
loc3=randi([50,53],1,1);
xc(1,loc3)=signal2_p;
loc4=randi([55,58],1,1);
xc(1,loc4)=signal1_p;
% loc5=randi([90,93],1,1);    %接近杂波区，但是依然在
% xc(1,loc5)=signal1_p;
% xc(1,loc5)=signal2_p;     %这种杂波边缘有目标的情况，还需要优化
% xc(1,loc5)=signal3_p;

loc6=randi([102,108],1,1);
xc(1,loc6)=signal3_p;

%% 算法结果&图谱显示
N=36;
pro_N=2;
PAD=10^(-4);
[ index, XT ] = cfar_go( abs(xc), N, pro_N, PAD);

figure;
plot(10.*log(abs(xc))./log(10)),hold on;
plot(index,10.*log(abs(XT))./log(10)),hold on;

 

图7 GO-CFAR效果图

 

 

cfar_so.m　　最小选择恒虚警算法实现

优点：多目标效果有改进；

缺点：杂波边缘区域虚警概率提升；

function [ index, XT ] = cfar_so( xc, N, pro_N, PAD)
%   假设回波服从高斯分布
%   alpha赋值有些问题，一个比较复杂的高次函数
alpha=N.*(PAD.^(-1./N)-1);

index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
XT=zeros(1,length(index));

for i=index
    cell_left=xc(1,i-N/2-pro_N/2:i-pro_N/2-1);
    cell_right=xc(1,i+pro_N/2+1:i+N/2+pro_N/2);
    Z=min([mean(cell_left),mean(cell_right)]);
    
    XT(1,i-N/2-pro_N/2)=Z.*alpha;
end

end

test_cfar_so.m

%% 程序初始化
clc;clear all;close all;

%% 均匀背景噪声（单目标&多目标）
% shape=[200];
% variance=200;
% noise_db=20;
% noise_p=10.^(noise_db./10);
% show_out=0;
% [ xc ] = env_uniform(variance,  shape, noise_db,show_out);

% 单目标
% SNR1=15;    signal1_p=10.^(SNR1./10).*noise_p;
% xc(1,90)=signal1_p;
% 
% 多目标
% SNR1=20;signal1_p=10.^(SNR1./10).*noise_p;
% xc(1,90)=signal1_p;
% SNR2=15;signal2_p=10.^(SNR2./10).*noise_p;
% xc(1,102)=signal2_p;


%% 杂波边缘背景噪声（单目标&多目标）
shape=[100,200];
variance=200;
noise_db=[20,30];
noise_p=10.^(noise_db./10);
show_out=0;
[ xc ] = env_edge(variance,  shape, noise_db,show_out);

% 多目标
SNR1=15;signal1_p=10.^(SNR1./10).*noise_p(1,end);
SNR2=12;signal2_p=10.^(SNR2./10).*noise_p(1,end);
SNR3=8;signal3_p=10.^(SNR3./10).*noise_p(1,end);
SNR4=5;signal4_p=10.^(SNR4./10).*noise_p(1,end);

loc1=randi([43,44],1,1);
xc(1,loc1)=signal1_p;
loc2=randi([46,48],1,1);
xc(1,loc2)=signal3_p;
loc3=randi([50,53],1,1);
xc(1,loc3)=signal2_p;
loc4=randi([55,58],1,1);
xc(1,loc4)=signal1_p;
% loc5=randi([90,93],1,1);    %接近杂波区，但是依然在
% xc(1,loc5)=signal1_p;
% xc(1,loc5)=signal2_p;     %这种杂波边缘有目标的情况，还需要优化
% xc(1,loc5)=signal3_p;

loc6=randi([102,108],1,1);
xc(1,loc6)=signal3_p;

%% 算法结果&图谱显示
N=36;
pro_N=2;
PAD=10^(-4);
[ index, XT ] = cfar_so( abs(xc), N, pro_N, PAD);

figure;
plot(10.*log(abs(xc))./log(10)),hold on;
plot(index,10.*log(abs(XT))./log(10)),hold on;

 

图8 SO-CFAR效果图

 

 

cfar_os.m　　有序统计恒虚警算法实现

优点：多目标检测性能很好；

缺点：杂波边缘区域虚警概率提高；

function [ index, XT ] = cfar_os( xc, N, k, pro_N, PAD)
%   假设回波服从高斯分布
%   alpha赋值有些问题，一个比较复杂的高次函数

%% 计算alpha
% syms alpha PFA;
% PFA(alpha)=gamma(N-1).*gamma(N-k+alpha-1)./gamma(N-k-1)./gamma(N+alpha-1);
% [alpha,~,~]=solve(PFA(alpha)==PAD,'ReturnConditions', true) ;

alpha=N.*(PAD.^(-1./N)-1);

index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
XT=zeros(1,length(index));

for i=index
    cell_left=xc(1,i-N/2-pro_N/2:i-pro_N/2-1);
    cell_right=xc(1,i+pro_N/2+1:i+N/2+pro_N/2);
    cell_all=cat(2,cell_left,cell_right);
    cell_sort=sort(cell_all);
    
    Z=cell_sort(1,k);
    
    XT(1,i-N/2-pro_N/2)=Z.*alpha;
end

end

test_cfar_os.m

%% 程序初始化
clc;clear all;close all;

%% 均匀背景噪声（单目标&多目标）
% shape=[200];
% variance=200;
% noise_db=20;
% noise_p=10.^(noise_db./10);
% show_out=0;
% [ xc ] = env_uniform(variance,  shape, noise_db,show_out);

% 单目标
% SNR1=15;    signal1_p=10.^(SNR1./10).*noise_p;
% xc(1,90)=signal1_p;
% 
% 多目标
% SNR1=20;signal1_p=10.^(SNR1./10).*noise_p;
% xc(1,90)=signal1_p;
% SNR2=15;signal2_p=10.^(SNR2./10).*noise_p;
% xc(1,102)=signal2_p;


%% 杂波边缘背景噪声（单目标&多目标）
shape=[100,200];
variance=200;
noise_db=[20,30];
noise_p=10.^(noise_db./10);
show_out=0;
[ xc ] = env_edge(variance,  shape, noise_db,show_out);

% 多目标
SNR1=15;signal1_p=10.^(SNR1./10).*noise_p(1,end);
SNR2=12;signal2_p=10.^(SNR2./10).*noise_p(1,end);
SNR3=8;signal3_p=10.^(SNR3./10).*noise_p(1,end);
SNR4=5;signal4_p=10.^(SNR4./10).*noise_p(1,end);

loc1=randi([43,44],1,1);
xc(1,loc1)=signal1_p;
loc2=randi([46,48],1,1);
xc(1,loc2)=signal3_p;
loc3=randi([50,53],1,1);
xc(1,loc3)=signal2_p;
loc4=randi([55,58],1,1);
xc(1,loc4)=signal1_p;
loc5=randi([90,93],1,1);    %接近杂波区，但是依然在
%     xc(1,loc5)=signal1_p;
xc(1,loc5)=signal2_p;     %这种杂波边缘有目标的情况，还需要优化
%     xc(1,loc5)=signal3_p;

loc6=randi([102,108],1,1);
xc(1,loc6)=signal3_p;

%% 算法结果&图谱显示
N=36;
pro_N=2;
k=2.*N./4;
PAD=10^(-4);
[ index, XT ] = cfar_os( abs(xc), N, k, pro_N, PAD);

figure;
plot(10.*log(abs(xc))./log(10)),hold on;
plot(index,10.*log(abs(XT))./log(10)),hold on;

 

图9 OS-CAFR效果图

 

cfar_tc.m　　杂波图恒虚警算法实现

优点：通过提高累积帧数，杂波背景噪声功率趋近于理论值

缺点：对于慢移动目标有自遮掩效果

function [ index, XT, xc_tpn ] = cfar_tc( xc, xc_tp, N, pro_N, PAD)
%   假设回波服从高斯分布
%   
xc_tc=zeros(1,length(xc));
alpha=N.*(PAD.^(-1./N)-1);
xc_tc(1,2:end-1)=xc_tp(1,2:end-1)./3+xc(1,1:end-2)./3+xc(1,3:end)./3;
xc_tpn=xc_tp(1,1:end)./2+xc(1,1:end)./2;
xc_tc(1,1)=xc_tpn(1,1);
xc_tc(1,end)=xc_tpn(1,end);

index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
XT=zeros(1,length(index));

for i=index
    cell_left=xc_tc(1,i-N/2-pro_N/2:i-pro_N/2-1);
    cell_right=xc_tc(1,i+pro_N/2+1:i+N/2+pro_N/2);
    Z=(sum(cell_left)+sum(cell_right))./N;
    
    XT(1,i-N/2-pro_N/2)=Z.*alpha;
end

end

test_cfar_tc.m

%% 程序初始化
clc;clear all;close all;

%% 均匀背景噪声（单目标&多目标）
% shape=[200];
% variance=200;
% noise_db=20;
% noise_p=10.^(noise_db./10);
% show_out=0;
% [ xc ] = env_uniform(variance,  shape, noise_db,show_out);

% 单目标
% SNR1=15;    signal1_p=10.^(SNR1./10).*noise_p;
% xc(1,90)=signal1_p;
% 
% 多目标
% SNR1=20;signal1_p=10.^(SNR1./10).*noise_p;
% xc(1,90)=signal1_p;
% SNR2=15;signal2_p=10.^(SNR2./10).*noise_p;
% xc(1,102)=signal2_p;


%% 杂波边缘背景噪声（单目标&多目标）
shape=[100,200];
variance=200;
noise_db=[20,30];
noise_p=10.^(noise_db./10);
show_out=0;
% [ xc ] = env_edge(variance,  shape, noise_db, show_out);

% 多目标
SNR1=15;signal1_p=10.^(SNR1./10).*noise_p(1,end);
SNR2=12;signal2_p=10.^(SNR2./10).*noise_p(1,end);
SNR3=8;signal3_p=10.^(SNR3./10).*noise_p(1,end);
SNR4=5;signal4_p=10.^(SNR4./10).*noise_p(1,end);

% loc1=randi([43,44],1,1);
% xc(1,loc1)=signal1_p;
% loc2=randi([46,48],1,1);
% xc(1,loc2)=signal3_p;
% loc3=randi([50,53],1,1);
% xc(1,loc3)=signal2_p;
% loc4=randi([55,58],1,1);
% xc(1,loc4)=signal1_p;
% loc5=randi([90,93],1,1);    %接近杂波区，但是依然在
% xc(1,loc5)=signal1_p;
% xc(1,loc5)=signal2_p;     %这种杂波边缘有目标的情况，还需要优化
% xc(1,loc5)=signal3_p;

% loc6=randi([102,108],1,1);
% xc(1,loc6)=signal3_p;

%% 迭代算法实现（杂波图）
N=36;
u=0.001;
PAD=10^(-4);
alpha=N.*(PAD.^(-1./N)-1);
xtp=zeros(1,shape(end));
d0=0;
d1=[0];
for i=1:1:4000
    [ xc ] = env_edge(variance,  shape, noise_db,show_out);
    
    loc1=randi([43,44],1,1);
    xc(1,loc1)=signal1_p;
    loc2=randi([46,48],1,1);
    xc(1,loc2)=signal3_p;
    loc3=randi([50,53],1,1);
    xc(1,loc3)=signal2_p;
    loc4=randi([55,58],1,1);
    xc(1,loc4)=signal1_p;
    
    loc6=randi([102,108],1,1);
    xc(1,loc6)=signal3_p;
    
    xtp=xtp.*(1-u)+xc.*u;
    d0=d0.*(1-u)+u;
    d1=cat(2,d1,d0);
    
    if mod(i,100)==0
        %% 图谱显示
        figure;
        plot(10.*log(abs(xc))./log(10)),hold on;
        plot(10.*log(abs(xtp.*alpha./d0))./log(10)),hold on;
    end

end

 

图10 杂波图CFAR效果图

 

五，总结

 　　理想情况下的CFAR算法：

　　①多目标检测性能很好，几乎互不影响；

　　②杂波边缘区域，突变很快，虚警率稳定；

　　③慢移动目标自遮掩改善；

　　④算法实用，复杂度适中，理论可实现；

 

六，源程序完善[添加环境杂波模拟程序-一维距离元]　　2019.08.02

env_uniform.m文件

function [ xc ] = env_uniform(variance,  shape, power_db, show_out)
%UNIFORM_EN 此处显示有关此函数的摘要
%   此处显示详细说明
if (nargin==3)
    show_out=0;
end

c=10^(power_db/10);     % 这里是幅度——功率
xc=c + random('Normal',0,variance,1,shape); 
if show_out==1
    figure;plot(10.*log(abs(xc))./log(10));
end

end

env_edge.m文件

function [ xc ] = env_edge(variance,  shape, power_db, show_out)
%UNIFORM_EN 此处显示有关此函数的摘要
%   此处显示详细说明
if (nargin==3)
    show_out=0;
end

c=10.^(power_db./10);               % 这里是幅度——功率
xc=random('Normal',0,variance,1,shape(1,end)); 
xc(1,1:end)=xc(1,1:end)+c(1,1);
index=1;
for i=1:length(power_db)
    xc(1,index:shape(1,i))=xc(1,index:shape(1,i)).*c(1,i)./c(1,1);
    index=shape(1,i)+1;
end


if show_out==1
    figure;plot(20.*log(abs(xc))./log(10));
end

end