深入分析:恒虚警率检测算法之Switch-CFAR

  转发和使用请注明来源,以下为本人精心整理,还请尊重本人劳动成果与产权!由于本人现有知识和能力有限,如存在错误之处请指正!下面为正文内容:

  1.S-CFAR检测算法(Switching,开关CFAR)

   S-CFAR检测算法是通过参考单元内数据确定筛选门限,将小于筛选门限的参考单元划分至S0,将大于筛选门限的参考单元划分至S1,当S0的单元数目n0大于抗干扰数目时,可用S0来估计杂波背景功率水平。

  实现的步骤具体可以分为三步:

  (1)根据CA的思想求解出门限水平T=αZ,其中α为衰减系数(通过人为选择)。

  (2)根据(1)中求解得到的门限水平T分别与参考窗内每个单元采样值xi进行比较,如果xi>T则将xi放入S1序列,否则将xi放入S0序列中。

  (3)如果S0序列中单元数目n0超过整数门限NT,那么采用S0中的样本估计阈值,反之则利用全体参考单元的样本估计阈值。参考单元中是否存在目标的判决规则如下:

x>β0n0Σxi=S0xi,n0>NT 

x>β12NΣi=12Nxi,n0<=NT

  其中,S-CFAR检测算法的检测性能与αβ0β1NT有关。S-CFAR的设计流程中可以简化设计使β0=β1=β,其余参数的设计如下:(1)设计干扰目标数目(最大可容纳的干扰目标数量)为Ni,整数门限NT满足关系NT=2NNi1.(2)在选定NT以后,对每个区间的α值,根据Pfa式可以计算出来门限系数β

  改进的SOS-CFAR算法中,更新了判决规则为如下,其中k为选择序号。

x>β0x(k)xi,n0>NT  

x>β12NΣi=12Nxi,n0<=NT 

  比较S-CFAR与SOS-CFAR两种算法的区别主要在于更新了n0>NT时的判决规则(门限水平估计准则),这两种算法的框图可以如下图表示(图片绘制较粗糙,但可以反映原理)。

 

 

  为什么会出现上述两种判决规则,如何理解n0与整数门限NT间的关系呢?我的理解是,如果背景杂波占据了Nt/N比例(这个比例比较大)的参考窗,那么可以用这部分的参考单元即S0的采样值作为功率水平的估计(估计方式可以采用CA,OS等);反之如果背景杂波占据的比例非常小,其强度可以忽略不计,则可以用参考窗内全体采样值的平均进行估计(准确地说此时参考窗内被强杂波占据)。  

  理解判决规则后自然需要推导出S-CFAR算法的虚警概率PFA,经典版本S-CFAR的虚警概率具体形式如下:

P(N,NT,α,β0,β1)=11+σn0=0Ntm=M0min(M,n0)n=0n0mk=0m(Mm)(2NMn0m)(n0mn)(mk)(1)k+nV(k,n)(1+σI)M+11+σn0=Nt+12Nm=M0min(M,n0)n=0n0mk=0m(Mm)(2NMn0m)(n0mn)(mk)(1)k+nW(k,n)(1+σI)m

  其中,W(k,n)的具体形式如下:

W(k,n)=ξk,nn0[1ξk,NσI(1+σI)(ξk,n+ρk,n+(1+σ)1)]m(ρk,n+(1+σ)1)(ξk,n+ρk,n+(1+σ)1)n0ξk,n=n0β01(k+n)α0ρk,n=α[Mm+k1+σI+2NMn0+m+n]

   其中,V(k,n)的具体形式如下:

V(k,n)=ξk,n2N[1ξk,NσI(1+σI)(ξk,n+ρk,n+(1+σ)1)]M(ρk,n+(1+σ)1)(ξk,n+ρk,n+(1+σ)1)2Nξk,n=2Nβ11(2Nn0+k+n+α)0ρk,n=α[Mm+k1+σI+2NMn0+m+n]

  上述虚警概率和相关公式主要参考了“Constant false-alarm rate algorithm based on test cell information”等论文(可在谷歌学术上搜索得到),具体公式推导由于篇幅问题不再展开。从虚警概率公式可见,该公式与先前的恒虚警率检测算法的虚警概率公式(多为关于门限系数α的闭合函数形式)不同,由多个变量(包含衰减系数α,整数门限Nt,门限系数β,功率系数σ)共同作用,因此采用二分法求解门限系数显得困难,故多篇论文中采用了交叉验证的方式来求解合适的门限系数β

  设置功率系数σ1σ为0,参考窗长度2N=64N=32,设置整数门限Nt=32,设置衰减系数α为[0.1,1]并且步进为0.1,设置门限系数ββ1为[10:50]并且步进为0.1,代入上述的虚警概率公式中求解并取对数形式,可以得到虚警概率关于不同门限系数和衰减系数的曲线如下图所示.可以看到当衰减系数从1降低至0.1时,S-CFAR的虚警概率曲线逐渐偏离CA-CFAR的虚警概率曲线,这意味着要实现相同的检测性能时S-CFAR需要花销比CA-CFAR更大的门限系数代价,同时在均匀背景下为实现相同的性能S-CFAR需要有更多的CFAR损失;而衰减系数更大时其逼近CA,这不利于多目标背景下的目标检测。因此在衰减系数的选择上需要进行折中选择,这里可以选择门限系数α=0.3.

 

 

 

  选择完合适的衰减系数后需要进一步确定整数门限Nt,整数门限的选择影响着最大可容纳干扰数目(Ni=nNNt1)。从图中可见,随着整数门限的下降,虚警概率Pfa关于门限系数β的曲线会逐渐偏离Nt=32时的曲线。整数门限越大,意味着参考窗内可容纳干扰数目越小,对多目标背景下目标检测越不利。因此,选择合适的整数门限,对多目标背景检测是有意义的。这里可以选择整数门限Nt=17

 

   综上分析,这里可以选择的参数为:虚警概率pfa=105,衰减系数α=0.3,整数门限Nt=17,门限系数β=15.03.

  以上程序的代码如下所示。

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
%% S-CFAR的虚警概率计算
clc;clear;close all;
 
%% 参数设置
sigma = 0;
sigmaI = 0;
alpha_array = 0.1:0.1:1.0; %衰减系数
beta_array = 10:0.1:50;
M = 0;
N = 32;
Nt = 17;
Pfa = zeros(length(alpha_array),length(beta_array));
 
%% 计算虚警概率
for index = 1:length(alpha_array)
    alpha = alpha_array(index);
    for index1 = 1:length(beta_array)
        beta0 = beta_array(index1);
        Qn0 = 0;
        for n0 = 0:Nt
            M0 = max(0,n0 - 2 * N + M);
            for m = M0:min(n0,M)
                for n = 0:(n0-m)
                    for k = 0:m
                        rou = alpha * ((M - m + k) / (1 + sigmaI) + 2 * N - M - n0 + m + n);
                        phy =  2 * N / beta0  - (2 * N - n0 + k + n) * alpha;
                        if(phy < 0)
                            break;
                        else
                           V = phy^(2*N) * (1-(phy*sigmaI/(1+sigmaI)/(phy+rou+(1+sigma)^(-1)))^(-M)) / (rou+(1+sigma)^(-1)) / (phy+rou+(1+sigma)^(-1))^(2*N);
                           Qn0 = Qn0 + 1 / (1+sigma) * nchoosek(M,m) * nchoosek(2*N-M,n0-m) * nchoosek(n0-m,n) * nchoosek(m,k) * (-1)^(k+n) / (1+sigmaI)^M * V;
                        end
                    end
                end
            end
        end
 
        Pn0 = 0;
        for n0 = Nt+1:2*N
            M0 = max(0,n0 - 2 * N + M);
            for m = M0:min(M,n0)
                for n = 0:n0-m
                    for k = 0:m
                        phy = n0 / beta0 - (k + n) * alpha;
                        rou = alpha * ((M - m + k) / (1 + sigmaI) + 2* N - M - n0 + m + n);
                        if phy < 0
                            break;
                        else
                           W = phy^n0 * (1 - (phy*sigmaI/(1+sigmaI))/(phy+rou+(1+sigma)^(-1)))^(-m)/(rou +(1+sigma)^(-1))/(phy + rou + (1+sigma)^(-1))^n0;
                           Pn0 = Pn0 + 1 / (1+sigma) * nchoosek(M,m) * nchoosek(2*N-M,n0-m) * nchoosek(n0-m,n) * nchoosek(m,k) * (-1)^(k+n) * W / (1+sigmaI)^m;
                        end
                    end
                end
            end
        end
        Pfa(index,index1) = Pn0 + Qn0;
    end
end
Pfa_ca = (1+beta_array/2/N).^(-2*N);
%% 绘制曲线
plot(beta_array,log10(Pfa_ca),'Color','#FFD700','LineWidth',2,'Marker','<','MarkerIndices',2:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(1,:)),'Color','#DC143C','LineWidth',2,'Marker','+','MarkerIndices',4:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(2,:)),'Color','#FF1493','LineWidth',2,'Marker','*','MarkerIndices',6:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(3,:)),'Color','#8B008B','LineWidth',2,'Marker','v','MarkerIndices',8:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(4,:)),'Color','#483D8B','LineWidth',2,'Marker','^','MarkerIndices',10:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(5,:)),'Color','#0000FF','LineWidth',2,'Marker','o','MarkerIndices',12:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(6,:)),'Color','#6495ED','LineWidth',2,'Marker','d','MarkerIndices',14:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(7,:)),'Color','#00BFFF','LineWidth',2,'Marker','p','MarkerIndices',16:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(8,:)),'Color','#00CED1','LineWidth',2,'Marker','s','MarkerIndices',18:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(9,:)),'Color','#00FA9A','LineWidth',2,'Marker','x','MarkerIndices',20:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(10,:)),'Color','#DAA520','LineWidth',2,'Marker','>','MarkerIndices',22:4:length(beta_array));hold on;
h = legend('CA','\alpha=0.1','\alpha=0.2','\alpha=0.3','\alpha=0.4','\alpha=0.5','\alpha=0.6','\alpha=0.7','\alpha=0.8','\alpha=0.9','\alpha=1.0','NumColumns',2);
set(h,'edgecolor','none');
xlabel('\beta');
ylabel('log_{10}Pfa');
grid minor;

 

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
%% S-CFAR的虚警概率计算
clc;clear;close all;
 
%% 参数设置
sigma = 0;
sigmaI = 0;
alpha = 0.3; %衰减系数
beta_array = 5:0.1:50;
M = 0;
N = 32;
Nt_array = 12:5:32;
Pfa = zeros(length(Nt_array),length(beta_array));
 
%% 计算虚警概率
for index = 1:length(Nt_array)
    Nt = Nt_array(index);
    for index1 = 1:length(beta_array)
        beta0 = beta_array(index1);
        Qn0 = 0;
        for n0 = 0:Nt
            M0 = max(0,n0 - 2 * N + M);
            for m = M0:min(n0,M)
                for n = 0:(n0-m)
                    for k = 0:m
                        rou = alpha * ((M - m + k) / (1 + sigmaI) + 2 * N - M - n0 + m + n);
                        phy =  2 * N / beta0  - (2 * N - n0 + k + n) * alpha;
                        if(phy < 0)
                            break;
                        else
                           V = phy^(2*N) * (1-(phy*sigmaI/(1+sigmaI)/(phy+rou+(1+sigma)^(-1)))^(-M)) / (rou+(1+sigma)^(-1)) / (phy+rou+(1+sigma)^(-1))^(2*N);
                           Qn0 = Qn0 + 1 / (1+sigma) * nchoosek(M,m) * nchoosek(2*N-M,n0-m) * nchoosek(n0-m,n) * nchoosek(m,k) * (-1)^(k+n) / (1+sigmaI)^M * V;
                        end
                    end
                end
            end
        end
 
        Pn0 = 0;
        for n0 = Nt+1:2*N
            M0 = max(0,n0 - 2 * N + M);
            for m = M0:min(M,n0)
                for n = 0:n0-m
                    for k = 0:m
                        phy = n0 / beta0 - (k + n) * alpha;
                        rou = alpha * ((M - m + k) / (1 + sigmaI) + 2* N - M - n0 + m + n);
                        if phy < 0
                            break;
                        else
                           W = phy^n0 * (1 - (phy*sigmaI/(1+sigmaI))/(phy+rou+(1+sigma)^(-1)))^(-m)/(rou +(1+sigma)^(-1))/(phy + rou + (1+sigma)^(-1))^n0;
                           Pn0 = Pn0 + 1 / (1+sigma) * nchoosek(M,m) * nchoosek(2*N-M,n0-m) * nchoosek(n0-m,n) * nchoosek(m,k) * (-1)^(k+n) * W / (1+sigmaI)^m;
                        end
                    end
                end
            end
        end
        Pfa(index,index1) = Pn0 + Qn0;
    end
end
 
%% 绘制曲线
plot(beta_array,log10(Pfa(1,:)),'Color','#DC143C','LineWidth',2,'Marker','+','MarkerIndices',2:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(2,:)),'Color','#FF1493','LineWidth',2,'Marker','*','MarkerIndices',4:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(3,:)),'Color','#8B008B','LineWidth',2,'Marker','v','MarkerIndices',6:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(4,:)),'Color','#483D8B','LineWidth',2,'Marker','^','MarkerIndices',8:4:length(beta_array));hold on;
plot(beta_array,log10(Pfa(5,:)),'Color','#0000FF','LineWidth',2,'Marker','o','MarkerIndices',10:4:length(beta_array));hold on;
h = legend('Nt=12','Nt=17','Nt=22','Nt=27','Nt=32','NumColumns',1);
set(h,'edgecolor','none');
xlabel('\beta');
ylabel('log_{10}Pfa');
grid minor;

  

 

  将S-CFAR在邻近单目标背景下进行仿真,随机产生在200个距离单元内产生信杂比分别为15dB、20dB的邻近目标,经典均值类ML-CFAR和有序统计类OS-CFAR以及S-CFAR的检测结果如下图所示,OS-CFAR和S-CFAR成功检测除了邻近目标,具有较好的多目标检测性能。

 

   具体代码如下所示。

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
clc;clear;close all;
%% 均匀背景
% % lambda = 1; %指数分布参数
% % dB = 20; %背景杂波功率
% % SCR = 15; %信杂比
% % num = [1,1024]; %距离单元设置
% % [x,target] = uniform_background(lambda,dB,SCR,num); %产生均匀背景
 
%% 多目标背景
lambda = 1; %指数分布参数
dB = 20; %背景杂波功率
SCR1 = 15; %信杂比
SCR2 = 20; %信杂比
num = [1,200]; %距离单元设置
localRule = [20,3,10]; %距离规则
[x,target] = multi_background(lambda,dB,SCR1,SCR2,num,localRule); %产生邻近单目标背景
 
%% CFAR参数
Pfa = 1e-6; %虚警概率
NSlide = 32;%滑动单元大小
Pro_cell = 4; %保护单元大小
 
%% 绘图参数
len = num(2);
plotNum = 2;
markersize = 5;
LineWidth = 1.5;
figure(1);
plot(target(:,1),form_Power_P2DB(target(:,2)),'s','Color',[1 0 1],'MarkerSize',markersize+2,'LineWidth',LineWidth),hold on;
plot(form_Power_P2DB(x),'Color','[0.5 0.5 1]','LineWidth',LineWidth),hold on;
 
%% 恒虚警率检测算法测试
Nt = 17;beta = 15.03;
rate = 0.75; %选择序号
alpha(1,1) = ca_threhold(Pfa,NSlide); %CA-CFAR标称因子计算
alpha(2,1) = go_threhold(Pfa,NSlide); %GO-CFAR标称因子计算
alpha(3,1) = so_threhold(Pfa,NSlide); %SO-CFAR标称因子计算
alpha(4,1) = os_threhold(Pfa,NSlide,rate); %OS-CFAR标称因子计算
 
[result_S] = func_cfar_s(x,beta,NSlide,Nt,Pro_cell); %求解S-CFAR门限
[result_CA] = func_cfar_ca(x,alpha(1,1),NSlide,Pro_cell); %求解CA-CFAR门限
[result_GO] = func_cfar_go(x,alpha(2,1),NSlide,Pro_cell); %求解GO-CFAR门限
[result_SO] = func_cfar_so(x,alpha(3,1),NSlide,Pro_cell); %求解SO-CFAR门限
[result_OS] = func_cfar_os(x,alpha(4,1),NSlide,Pro_cell,rate); %求解OS-CFAR门限
%% 绘制检测曲线
[T_opt] = T_opt_solve(Pfa,lambda,dB,num); %求解理想门限
plot(form_Power_P2DB(T_opt),'Color','#3CB371','Marker','+','MarkerIndices',2:plotNum:len,'LineWidth',LineWidth);
hold on;
plot(form_Power_P2DB(result_CA{1,2}),'Color','#DC143C','Marker','*','MarkerIndices',4:plotNum:len,'LineWidth',LineWidth);
hold on;
plot(form_Power_P2DB(result_GO{1,2}),'Color','#ADFF2F','Marker','v','MarkerIndices',6:plotNum:len,'markersize',markersize,'LineWidth',LineWidth);
hold on;
plot(form_Power_P2DB(result_SO{1,2}),'Color','#8A2BE2','Marker','^','MarkerIndices',8:plotNum:len,'markersize',markersize,'LineWidth',LineWidth);
hold on;
plot(form_Power_P2DB(result_OS{1,2}),'Color','#0000CD','Marker','o','MarkerIndices',10:plotNum:len,'markersize',markersize,'LineWidth',LineWidth);
hold on;
plot(form_Power_P2DB(result_S{1,2}),'Color','#7CFC00','Marker','p','MarkerIndices',12:plotNum:len,'markersize',markersize,'LineWidth',LineWidth);
hold on
grid minor;
xlabel('\fontname{宋体}距离单元');
ylabel('\fontname{宋体}功率值\fontname{Times New Roman}(dB)');
ylim([0 50]);
h = legend('TARGET','CLUTTER','OPT','CA-CFAR','GO-CFAR','SO-CFAR','OS-CFAR','S-CFAR','Location','SouthEast','NumColumns',1);
set(h,'edgecolor','none')

 

 

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
%% 本程序主要实现S-CFAR恒虚警率检测算法的函数形式
function result = func_cfar_s(x,beta,NSlide,Nt,Pro_cell)
    %x:原始杂波数据
    %alpha:标称因子
    %Nslide:滑动窗大小
    %Pro_cell:保护单元大小
    %Nt:整数门限
 
    persistent left;
    persistent right;
    persistent Half_Slide;
    persistent Half_Pro_cell;
    persistent len;
     
    if isempty(left)
        left = 1 + Half_Pro_cell + Half_Slide; %设置左边界
        right = length(x) - Half_Pro_cell - Half_Slide; %设置右边界
        Half_Slide = NSlide / 2; %半滑动窗
        Half_Pro_cell = Pro_cell / 2; %一侧保护单元长度
        len = length(x); %杂波单元
    end
 
     
    T = zeros(1,len); %CMLD检测阈值
    target = java.util.LinkedList; %利用Java链表来实现目标的存储
 
    for i = 1:left - 1 %左边界
        cell_right = x(1,i + Half_Pro_cell + 1 : i + Half_Slide * 2 + Half_Pro_cell);
        S0 = cell_right(cell_right <= mean(cell_right));
        S1 = cell_right(cell_right > mean(cell_right));
        if length(S0) > Nt
            T(1,i) = mean(S0) * beta;
        else
            T(1,i) = mean(cell_right) * beta;
        end
        if T(1,i) < x(i)
            target.add(i); %加入目标
        end
    end
     
    for i = left:right %中间区域
        cell_left = x(1,i - Half_Pro_cell - Half_Slide : i - Half_Pro_cell - 1);
        cell_right = x(1,i + Half_Pro_cell + 1 : i + Half_Pro_cell + Half_Slide);
        cell = [cell_left,cell_right];
        S0 = cell(cell <= mean(cell));
        S1 = cell(cell > mean(cell));
        if length(S0) > Nt
            T(1,i) = mean(S0) * beta;
        else
            T(1,i) = mean(cell) * beta;
        end
        if T(1,i) < x(i)
            target.add(i); %加入目标
        end
    end
     
    for i = right + 1 : len
        cell_left = x(1,i - Half_Pro_cell - Half_Slide * 2 : i - Half_Pro_cell - 1);
        S0 = cell_left(cell_left <= mean(cell_left));
        S1 = cell_left(cell_left > mean(cell_left));
        if length(S0) > Nt
            T(1,i) = mean(S0) * beta;
        else
            T(1,i) = mean(cell_left) * beta;
        end
        if T(1,i) < x(i)
            target.add(i);
        end
    end
     
    result = {'CFAR_S',T,target}; %采用字典类型
end
     
    

  

 

 

 

 

 

 

 

 

posted @   信海  阅读(1487)  评论(0编辑  收藏  举报
编辑推荐:
· Linux系列:如何用 C#调用 C方法造成内存泄露
· AI与.NET技术实操系列(二):开始使用ML.NET
· 记一次.NET内存居高不下排查解决与启示
· 探究高空视频全景AR技术的实现原理
· 理解Rust引用及其生命周期标识(上)
阅读排行:
· DeepSeek 开源周回顾「GitHub 热点速览」
· 物流快递公司核心技术能力-地址解析分单基础技术分享
· .NET 10首个预览版发布:重大改进与新特性概览!
· AI与.NET技术实操系列(二):开始使用ML.NET
· 单线程的Redis速度为什么快?
点击右上角即可分享
微信分享提示