雷達無線電系列（二）經典CFAR算法圖文解析與實現（matlab）

時間 2019-11-13

標籤雷達無線電系列經典 cfar 算法圖文解析實現 matlab 欄目 MATLAB 简体版

原文原文鏈接

一，CFRA基礎知識介紹算法

簡介

　　恆虛警檢測技術是指雷達系統在保持虛警機率恆定條件下對接收機輸出的信號與噪聲做判別以肯定目標信號是否存在的技術。

前提

　　因爲接收機輸出端中確定存有噪聲(包括大氣噪聲、人爲噪聲、內部噪聲和雜波等)，而信號通常是疊加在噪聲上的。這就須要在接收機輸出的噪聲或信號加噪聲條件下，採用檢測技術判別是否有目標信號。

偏差機率

　　任何形式的判決必然存在着兩種偏差機率：發現機率和虛警機率。當接收機輸出端存在目標回波信號，而判決時判爲有目標的機率爲Pd，判爲無目標的機率爲1-Pad。當接收機輸出端只有噪聲時，而判爲有目標的機率爲Pfa。因爲噪聲是隨機變量，其特徵可用機率密度函數表示，所以信號加噪聲也是一隨機變量 [1]

具體過程

　　恆虛警檢測器首先對輸入的噪聲進行處理後肯定一個門限，將此門限與輸入端信號相比，如輸入端信號超過了此門限，則判爲有目標，不然，判爲無目標。

方法

　　爲了使系統的虛警機率保持恆定，此門限是隨輸入噪聲變化而進行快速的自適應調整的，噪聲處理方法是隨噪聲的不一樣分佈而異，於是恆虛警檢測技術包括恆虛警處理技術和目標檢測技術兩大部分。恆虛警處理技術包括快門限、慢門限處理技術、目標檢測包括似然比檢測、二進制檢測、序貫檢測和非參量檢測等 [1]。

二，經典CFAR算法解析與性能圖譜展現網絡

　　①均值類CFRAdom

　　核心思想是經過對參考窗內採樣數據取平均來估計背景功率。ide

　　CA-CFAR、GO-CFAR、SO-CFAR算法這三個是最經典的均值類CFAR算法，後續有其變形log-CFAR算法可是原理是相同的。函數

圖1 均值類CFAR檢測器原理圖性能

　　②統計有序CFAR優化

　　核心思想：經過對參考窗內的數據由小到大排序選取其中第K個數值假設其爲雜波背景噪聲。spa

圖2 OS-CFAR檢測器原理圖3d

　　③自適應CFARcode

　　針對不一樣的雜波選用不一樣的決策方法來進行。

圖3 HCE-CFAR檢測器原理圖

圖4 VI-CFAR檢測器原理圖

　　④雜波圖CFAR

　　經典的空域雜波圖CFAR算法，假設每一個距離元和方位角上的平方率值聽從獨立同分布，經過對時間維上的數據進行遺忘迭代來擬合雜波背景噪聲。

圖5 雜波圖CFAR檢測器原理圖

三，CFRA發展示狀

　　均值類-統計有序類-自適應-雜波圖-變換域-壓縮感知-神經網絡

四，經典CFAR算法實現源碼（matlab）

cfar_ac.m　　單元平均恆虛警算法實現

優勢：損失率最少的一種算法；

缺點：多目標遮掩，雜波邊緣性能也欠佳；

 1 function [ index, XT ] = cfar_ac( xc, N, pro_N, PAD)
 2 %   假設回波服從高斯分佈
 3 %   
 4 alpha=N.*(PAD.^(-1./N)-1);
 5 index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
 6 XT=zeros(1,length(index));
 7 
 8 for i=index
 9     cell_left=xc(1,i-N/2-pro_N/2:i-pro_N/2-1);
10     cell_right=xc(1,i+pro_N/2+1:i+N/2+pro_N/2);
11     Z=(sum(cell_left)+sum(cell_right))./N;
12     
13     XT(1,i-N/2-pro_N/2)=Z.*alpha;
14 end
15 
16 end

View Code

test_cfar_ca.m

 1 %% 程序初始化
 2 clc;clear all;close all;
 3 
 4 %% 均勻背景噪聲（單目標&多目標）
 5 % shape=[200];
 6 % variance=200;
 7 % noise_db=20;
 8 % noise_p=10.^(noise_db./10);
 9 % show_out=0;
10 % [ xc ] = env_uniform(variance,  shape, noise_db,show_out);
11 
12 % 單目標
13 % SNR1=15;    signal1_p=10.^(SNR1./10).*noise_p;
14 % xc(1,90)=signal1_p;
15 % 
16 % 多目標
17 % SNR1=20;signal1_p=10.^(SNR1./10).*noise_p;
18 % xc(1,90)=signal1_p;
19 % SNR2=15;signal2_p=10.^(SNR2./10).*noise_p;
20 % xc(1,102)=signal2_p;
21 
22 
23 %% 雜波邊緣背景噪聲（單目標&多目標）
24 shape=[100,200];
25 variance=200;
26 noise_db=[20,30];
27 noise_p=10.^(noise_db./10);
28 show_out=0;
29 [ xc ] = env_edge(variance,  shape, noise_db,show_out);
30 
31 % 多目標
32 SNR1=15;signal1_p=10.^(SNR1./10).*noise_p(1,end);
33 SNR2=12;signal2_p=10.^(SNR2./10).*noise_p(1,end);
34 SNR3=8;signal3_p=10.^(SNR3./10).*noise_p(1,end);
35 SNR4=5;signal4_p=10.^(SNR4./10).*noise_p(1,end);
36 
37 loc1=randi([43,44],1,1);
38 xc(1,loc1)=signal1_p;
39 loc2=randi([46,48],1,1);
40 xc(1,loc2)=signal3_p;
41 loc3=randi([50,53],1,1);
42 xc(1,loc3)=signal2_p;
43 loc4=randi([55,58],1,1);
44 xc(1,loc4)=signal1_p;
45 loc5=randi([90,93],1,1);    %接近雜波區，可是依然在
46 %     xc(1,loc5)=signal1_p;
47 xc(1,loc5)=signal2_p;     %這種雜波邊緣有目標的狀況，還須要優化
48 %     xc(1,loc5)=signal3_p;
49 
50 loc6=randi([102,108],1,1);
51 xc(1,loc6)=signal3_p;
52 
53 %% 算法結果&圖譜顯示
54 N=36;
55 pro_N=2;
56 PAD=10^(-4);
57 [ index, XT ] = cfar_ac( xc, N, pro_N, PAD);
58 
59 figure;
60 plot(10.*log(abs(xc))./log(10)),hold on;
61 plot(index,10.*log(abs(XT))./log(10)),hold on;

View Code

圖6 CA-CFAR效果圖

cfar_go.m　　最大選擇恆虛警算法實現

優勢：雜波邊緣區域虛警機率下降

缺點：多目標遮掩

 1 function [ index, XT ] = cfar_go( xc, N, pro_N, PAD)
 2 %   假設回波服從高斯分佈
 3 %   alpha賦值有些問題，一個比較複雜的高次函數
 4 alpha=N.*(PAD.^(-1./N)-1);
 5 
 6 index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
 7 XT=zeros(1,length(index));
 8 
 9 for i=index
10     cell_left=xc(1,i-N/2-pro_N/2:i-pro_N/2-1);
11     cell_right=xc(1,i+pro_N/2+1:i+N/2+pro_N/2);
12     Z=max([mean(cell_left),mean(cell_right)]);
13     
14     XT(1,i-N/2-pro_N/2)=Z.*alpha;
15 end
16 
17 end

View Code

test_cfar_go.m

 1 %% 程序初始化
 2 clc;clear all;close all;
 3 
 4 %% 均勻背景噪聲（單目標&多目標）
 5 % shape=[200];
 6 % variance=200;
 7 % noise_db=20;
 8 % noise_p=10.^(noise_db./10);
 9 % show_out=0;
10 % [ xc ] = env_uniform(variance,  shape, noise_db,show_out);
11 
12 % 單目標
13 % SNR1=15;    signal1_p=10.^(SNR1./10).*noise_p;
14 % xc(1,90)=signal1_p;
15 % 
16 % 多目標
17 % SNR1=20;signal1_p=10.^(SNR1./10).*noise_p;
18 % xc(1,90)=signal1_p;
19 % SNR2=15;signal2_p=10.^(SNR2./10).*noise_p;
20 % xc(1,102)=signal2_p;
21 
22 
23 %% 雜波邊緣背景噪聲（單目標&多目標）
24 shape=[100,200];
25 variance=200;
26 noise_db=[20,30];
27 noise_p=10.^(noise_db./10);
28 show_out=0;
29 [ xc ] = env_edge(variance,  shape, noise_db,show_out);
30 
31 % 多目標
32 SNR1=15;signal1_p=10.^(SNR1./10).*noise_p(1,end);
33 SNR2=12;signal2_p=10.^(SNR2./10).*noise_p(1,end);
34 SNR3=8;signal3_p=10.^(SNR3./10).*noise_p(1,end);
35 SNR4=5;signal4_p=10.^(SNR4./10).*noise_p(1,end);
36 
37 loc1=randi([43,44],1,1);
38 xc(1,loc1)=signal1_p;
39 loc2=randi([46,48],1,1);
40 xc(1,loc2)=signal3_p;
41 loc3=randi([50,53],1,1);
42 xc(1,loc3)=signal2_p;
43 loc4=randi([55,58],1,1);
44 xc(1,loc4)=signal1_p;
45 % loc5=randi([90,93],1,1);    %接近雜波區，可是依然在
46 % xc(1,loc5)=signal1_p;
47 % xc(1,loc5)=signal2_p;     %這種雜波邊緣有目標的狀況，還須要優化
48 % xc(1,loc5)=signal3_p;
49 
50 loc6=randi([102,108],1,1);
51 xc(1,loc6)=signal3_p;
52 
53 %% 算法結果&圖譜顯示
54 N=36;
55 pro_N=2;
56 PAD=10^(-4);
57 [ index, XT ] = cfar_go( abs(xc), N, pro_N, PAD);
58 
59 figure;
60 plot(10.*log(abs(xc))./log(10)),hold on;
61 plot(index,10.*log(abs(XT))./log(10)),hold on;

View Code

圖7 GO-CFAR效果圖

cfar_so.m　　最小選擇恆虛警算法實現

優勢：多目標效果有改進；

缺點：雜波邊緣區域虛警機率提高；

 1 function [ index, XT ] = cfar_so( xc, N, pro_N, PAD)
 2 %   假設回波服從高斯分佈
 3 %   alpha賦值有些問題，一個比較複雜的高次函數
 4 alpha=N.*(PAD.^(-1./N)-1);
 5 
 6 index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
 7 XT=zeros(1,length(index));
 8 
 9 for i=index
10     cell_left=xc(1,i-N/2-pro_N/2:i-pro_N/2-1);
11     cell_right=xc(1,i+pro_N/2+1:i+N/2+pro_N/2);
12     Z=min([mean(cell_left),mean(cell_right)]);
13     
14     XT(1,i-N/2-pro_N/2)=Z.*alpha;
15 end
16 
17 end

View Code

test_cfar_so.m

 1 %% 程序初始化
 2 clc;clear all;close all;
 3 
 4 %% 均勻背景噪聲（單目標&多目標）
 5 % shape=[200];
 6 % variance=200;
 7 % noise_db=20;
 8 % noise_p=10.^(noise_db./10);
 9 % show_out=0;
10 % [ xc ] = env_uniform(variance,  shape, noise_db,show_out);
11 
12 % 單目標
13 % SNR1=15;    signal1_p=10.^(SNR1./10).*noise_p;
14 % xc(1,90)=signal1_p;
15 % 
16 % 多目標
17 % SNR1=20;signal1_p=10.^(SNR1./10).*noise_p;
18 % xc(1,90)=signal1_p;
19 % SNR2=15;signal2_p=10.^(SNR2./10).*noise_p;
20 % xc(1,102)=signal2_p;
21 
22 
23 %% 雜波邊緣背景噪聲（單目標&多目標）
24 shape=[100,200];
25 variance=200;
26 noise_db=[20,30];
27 noise_p=10.^(noise_db./10);
28 show_out=0;
29 [ xc ] = env_edge(variance,  shape, noise_db,show_out);
30 
31 % 多目標
32 SNR1=15;signal1_p=10.^(SNR1./10).*noise_p(1,end);
33 SNR2=12;signal2_p=10.^(SNR2./10).*noise_p(1,end);
34 SNR3=8;signal3_p=10.^(SNR3./10).*noise_p(1,end);
35 SNR4=5;signal4_p=10.^(SNR4./10).*noise_p(1,end);
36 
37 loc1=randi([43,44],1,1);
38 xc(1,loc1)=signal1_p;
39 loc2=randi([46,48],1,1);
40 xc(1,loc2)=signal3_p;
41 loc3=randi([50,53],1,1);
42 xc(1,loc3)=signal2_p;
43 loc4=randi([55,58],1,1);
44 xc(1,loc4)=signal1_p;
45 % loc5=randi([90,93],1,1);    %接近雜波區，可是依然在
46 % xc(1,loc5)=signal1_p;
47 % xc(1,loc5)=signal2_p;     %這種雜波邊緣有目標的狀況，還須要優化
48 % xc(1,loc5)=signal3_p;
49 
50 loc6=randi([102,108],1,1);
51 xc(1,loc6)=signal3_p;
52 
53 %% 算法結果&圖譜顯示
54 N=36;
55 pro_N=2;
56 PAD=10^(-4);
57 [ index, XT ] = cfar_so( abs(xc), N, pro_N, PAD);
58 
59 figure;
60 plot(10.*log(abs(xc))./log(10)),hold on;
61 plot(index,10.*log(abs(XT))./log(10)),hold on;

View Code

圖8 SO-CFAR效果圖

cfar_os.m　　有序統計恆虛警算法實現

優勢：多目標檢測性能很好；

缺點：雜波邊緣區域虛警機率提升；

 1 function [ index, XT ] = cfar_os( xc, N, k, pro_N, PAD)
 2 %   假設回波服從高斯分佈
 3 %   alpha賦值有些問題，一個比較複雜的高次函數
 4 
 5 %% 計算alpha
 6 % syms alpha PFA;
 7 % PFA(alpha)=gamma(N-1).*gamma(N-k+alpha-1)./gamma(N-k-1)./gamma(N+alpha-1);
 8 % [alpha,~,~]=solve(PFA(alpha)==PAD,'ReturnConditions', true) ;
 9 
10 alpha=N.*(PAD.^(-1./N)-1);
11 
12 index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
13 XT=zeros(1,length(index));
14 
15 for i=index
16     cell_left=xc(1,i-N/2-pro_N/2:i-pro_N/2-1);
17     cell_right=xc(1,i+pro_N/2+1:i+N/2+pro_N/2);
18     cell_all=cat(2,cell_left,cell_right);
19     cell_sort=sort(cell_all);
20     
21     Z=cell_sort(1,k);
22     
23     XT(1,i-N/2-pro_N/2)=Z.*alpha;
24 end
25 
26 end

View Code

test_cfar_os.m

 1 %% 程序初始化
 2 clc;clear all;close all;
 3 
 4 %% 均勻背景噪聲（單目標&多目標）
 5 % shape=[200];
 6 % variance=200;
 7 % noise_db=20;
 8 % noise_p=10.^(noise_db./10);
 9 % show_out=0;
10 % [ xc ] = env_uniform(variance,  shape, noise_db,show_out);
11 
12 % 單目標
13 % SNR1=15;    signal1_p=10.^(SNR1./10).*noise_p;
14 % xc(1,90)=signal1_p;
15 % 
16 % 多目標
17 % SNR1=20;signal1_p=10.^(SNR1./10).*noise_p;
18 % xc(1,90)=signal1_p;
19 % SNR2=15;signal2_p=10.^(SNR2./10).*noise_p;
20 % xc(1,102)=signal2_p;
21 
22 
23 %% 雜波邊緣背景噪聲（單目標&多目標）
24 shape=[100,200];
25 variance=200;
26 noise_db=[20,30];
27 noise_p=10.^(noise_db./10);
28 show_out=0;
29 [ xc ] = env_edge(variance,  shape, noise_db,show_out);
30 
31 % 多目標
32 SNR1=15;signal1_p=10.^(SNR1./10).*noise_p(1,end);
33 SNR2=12;signal2_p=10.^(SNR2./10).*noise_p(1,end);
34 SNR3=8;signal3_p=10.^(SNR3./10).*noise_p(1,end);
35 SNR4=5;signal4_p=10.^(SNR4./10).*noise_p(1,end);
36 
37 loc1=randi([43,44],1,1);
38 xc(1,loc1)=signal1_p;
39 loc2=randi([46,48],1,1);
40 xc(1,loc2)=signal3_p;
41 loc3=randi([50,53],1,1);
42 xc(1,loc3)=signal2_p;
43 loc4=randi([55,58],1,1);
44 xc(1,loc4)=signal1_p;
45 loc5=randi([90,93],1,1);    %接近雜波區，可是依然在
46 %     xc(1,loc5)=signal1_p;
47 xc(1,loc5)=signal2_p;     %這種雜波邊緣有目標的狀況，還須要優化
48 %     xc(1,loc5)=signal3_p;
49 
50 loc6=randi([102,108],1,1);
51 xc(1,loc6)=signal3_p;
52 
53 %% 算法結果&圖譜顯示
54 N=36;
55 pro_N=2;
56 k=2.*N./4;
57 PAD=10^(-4);
58 [ index, XT ] = cfar_os( abs(xc), N, k, pro_N, PAD);
59 
60 figure;
61 plot(10.*log(abs(xc))./log(10)),hold on;
62 plot(index,10.*log(abs(XT))./log(10)),hold on;

View Code

圖9 OS-CAFR效果圖

cfar_tc.m　　雜波圖恆虛警算法實現

優勢：經過提升累積幀數，雜波背景噪聲功率趨近於理論值

缺點：對於慢移動目標有自遮掩效果

 1 function [ index, XT, xc_tpn ] = cfar_tc( xc, xc_tp, N, pro_N, PAD)
 2 %   假設回波服從高斯分佈
 3 %   
 4 xc_tc=zeros(1,length(xc));
 5 alpha=N.*(PAD.^(-1./N)-1);
 6 xc_tc(1,2:end-1)=xc_tp(1,2:end-1)./3+xc(1,1:end-2)./3+xc(1,3:end)./3;
 7 xc_tpn=xc_tp(1,1:end)./2+xc(1,1:end)./2;
 8 xc_tc(1,1)=xc_tpn(1,1);
 9 xc_tc(1,end)=xc_tpn(1,end);
10 
11 index=1+N/2+pro_N/2:length(xc)-N/2-pro_N/2;
12 XT=zeros(1,length(index));
13 
14 for i=index
15     cell_left=xc_tc(1,i-N/2-pro_N/2:i-pro_N/2-1);
16     cell_right=xc_tc(1,i+pro_N/2+1:i+N/2+pro_N/2);
17     Z=(sum(cell_left)+sum(cell_right))./N;
18     
19     XT(1,i-N/2-pro_N/2)=Z.*alpha;
20 end
21 
22 end

View Code

test_cfar_tc.m

 1 %% 程序初始化
 2 clc;clear all;close all;
 3 
 4 %% 均勻背景噪聲（單目標&多目標）
 5 % shape=[200];
 6 % variance=200;
 7 % noise_db=20;
 8 % noise_p=10.^(noise_db./10);
 9 % show_out=0;
10 % [ xc ] = env_uniform(variance,  shape, noise_db,show_out);
11 
12 % 單目標
13 % SNR1=15;    signal1_p=10.^(SNR1./10).*noise_p;
14 % xc(1,90)=signal1_p;
15 % 
16 % 多目標
17 % SNR1=20;signal1_p=10.^(SNR1./10).*noise_p;
18 % xc(1,90)=signal1_p;
19 % SNR2=15;signal2_p=10.^(SNR2./10).*noise_p;
20 % xc(1,102)=signal2_p;
21 
22 
23 %% 雜波邊緣背景噪聲（單目標&多目標）
24 shape=[100,200];
25 variance=200;
26 noise_db=[20,30];
27 noise_p=10.^(noise_db./10);
28 show_out=0;
29 % [ xc ] = env_edge(variance,  shape, noise_db, show_out);
30 
31 % 多目標
32 SNR1=15;signal1_p=10.^(SNR1./10).*noise_p(1,end);
33 SNR2=12;signal2_p=10.^(SNR2./10).*noise_p(1,end);
34 SNR3=8;signal3_p=10.^(SNR3./10).*noise_p(1,end);
35 SNR4=5;signal4_p=10.^(SNR4./10).*noise_p(1,end);
36 
37 % loc1=randi([43,44],1,1);
38 % xc(1,loc1)=signal1_p;
39 % loc2=randi([46,48],1,1);
40 % xc(1,loc2)=signal3_p;
41 % loc3=randi([50,53],1,1);
42 % xc(1,loc3)=signal2_p;
43 % loc4=randi([55,58],1,1);
44 % xc(1,loc4)=signal1_p;
45 % loc5=randi([90,93],1,1);    %接近雜波區，可是依然在
46 % xc(1,loc5)=signal1_p;
47 % xc(1,loc5)=signal2_p;     %這種雜波邊緣有目標的狀況，還須要優化
48 % xc(1,loc5)=signal3_p;
49 
50 % loc6=randi([102,108],1,1);
51 % xc(1,loc6)=signal3_p;
52 
53 %% 迭代算法實現（雜波圖）
54 N=36;
55 u=0.001;
56 PAD=10^(-4);
57 alpha=N.*(PAD.^(-1./N)-1);
58 xtp=zeros(1,shape(end));
59 d0=0;
60 d1=[0];
61 for i=1:1:4000
62     [ xc ] = env_edge(variance,  shape, noise_db,show_out);
63     
64     loc1=randi([43,44],1,1);
65     xc(1,loc1)=signal1_p;
66     loc2=randi([46,48],1,1);
67     xc(1,loc2)=signal3_p;
68     loc3=randi([50,53],1,1);
69     xc(1,loc3)=signal2_p;
70     loc4=randi([55,58],1,1);
71     xc(1,loc4)=signal1_p;
72     
73     loc6=randi([102,108],1,1);
74     xc(1,loc6)=signal3_p;
75     
76     xtp=xtp.*(1-u)+xc.*u;
77     d0=d0.*(1-u)+u;
78     d1=cat(2,d1,d0);
79     
80     if mod(i,100)==0
81         %% 圖譜顯示
82         figure;
83         plot(10.*log(abs(xc))./log(10)),hold on;
84         plot(10.*log(abs(xtp.*alpha./d0))./log(10)),hold on;
85     end
86 
87 end

View Code

圖10 雜波圖CFAR效果圖

五，總結

　　理想狀況下的CFAR算法：

　　①多目標檢測性能很好，幾乎互不影響；

　　②雜波邊緣區域，突變很快，虛警率穩定；

　　③慢移動目標自遮掩改善；

　　④算法實用，複雜度適中，理論可實現；

六，源程序完善[添加環境雜波模擬程序-一維距離元]　　2019.08.02

env_uniform.m文件

 1 function [ xc ] = env_uniform(variance,  shape, power_db, show_out)
 2 %UNIFORM_EN 此處顯示有關此函數的摘要
 3 %   此處顯示詳細說明
 4 if (nargin==3)
 5     show_out=0;
 6 end
 7 
 8 c=10^(power_db/10);     % 這裏是幅度——功率
 9 xc=c + random('Normal',0,variance,1,shape); 
10 if show_out==1
11     figure;plot(10.*log(abs(xc))./log(10));
12 end
13 
14 end

View Code


env_edge.m文件

 1 function [ xc ] = env_edge(variance,  shape, power_db, show_out)
 2 %UNIFORM_EN 此處顯示有關此函數的摘要
 3 %   此處顯示詳細說明
 4 if (nargin==3)
 5     show_out=0;
 6 end
 7 
 8 c=10.^(power_db./10);               % 這裏是幅度——功率
 9 xc=random('Normal',0,variance,1,shape(1,end)); 
10 xc(1,1:end)=xc(1,1:end)+c(1,1);
11 index=1;
12 for i=1:length(power_db)
13     xc(1,index:shape(1,i))=xc(1,index:shape(1,i)).*c(1,i)./c(1,1);
14     index=shape(1,i)+1;
15 end
16 
17 
18 if show_out==1
19     figure;plot(20.*log(abs(xc))./log(10));
20 end
21 
22 end