Gamma分布在海面目標檢測中的應用*

周 圍， 朱 勇， 杜玉晗， 謝俊好

(1. 中國電子科技集團公司第三十八研究所， 安徽合肥 230088；2. 哈爾濱工業大學電子與信息工程學院， 黑龍江哈爾濱 150001)

0 引言

雷達最重要的任務之一是目標檢測。固定門限檢測是一種經典的檢測方法，它利用檢測背景噪聲、雜波水平的先驗信息計算固定的檢測門限，若回波信號的幅度超過門限則判定為目標。該方法在先驗信息完全已知的均勻雜波背景下具有較好的效果，而在時變的雜波環境中檢測性能會急劇下降。對海雷達在海面目標檢測過程中不可避免地要接收到來自海面的后向散射信號，通常我們將這種信號稱為海雜波。由于海雜波功率水平在實際情況中是時變的，通常難以設置一個合適的門限保證良好的檢測概率和較低的虛警率，因此固定門限檢測器幾乎難以直接應用[1]。

為了在對目標正確檢測的基礎上保持恒定的虛警率，恒虛警(CFAR)處理技術在現代雷達檢測系統中獲得廣泛應用[2-6]，其核心思想是利用參考回波樣本對未知檢測背景的信息進行估計并自適應地設定檢測門限，在特定雜波、噪聲背景下實現目標的恒虛警率自動檢測。為了選擇、設計合適的檢測器及工作參數，需要對雜波背景進行統計特性分析。研究海雜波幅度統計特性主要利用海雜波的回波包絡或者功率進行非相干統計建模，對海面目標檢測算法的設計具有重要意義。在低分辨率雷達中，海雜波可以用散斑(Speckle)模型來建模，其幅度概率密度函數服從瑞利(Rayleigh)分布。隨著雷達距離分辨率的提高，海雜波幅度分布函數出現“拖尾”現象，并逐漸偏離瑞利分布。近年來，隨著對海雜波形成機理及物理含義的進一步理解，復合高斯模型在高分辨海雜波統計建模上得到廣泛認可與深入研究發展。學者們提出了多種復合分布模型，如K分布[7]、K+N分布[8]和帕累托分布[9]等，并利用大量的實測數據進行驗證。

需要注意的是，在實際情況中，檢測器類型和參考窗長度需要根據雜波背景和實際干擾目標特性來決定，雖然近年來CFAR檢測得到廣泛研究，但最優、通用的CFAR檢測器是幾乎無法獲得的[10]，因此在設計檢測方法、制定檢測準則時必須仔細考慮實際檢測背景。同時，雖然經過長期的研究分析，學者們對于海雜波的物理機理已有相對成熟的結論，但由于其散射機制的復雜性以及雜波特性對實際雷達系統參數、外界環境的依賴性，在實際應用中很難找到一個與海雜波實際情況完全符合的模型[11]，目前人們對海雜波特性的研究與建模仍然主要依靠實際測量值與觀察結果。

為豐富海雜波建模方式、進一步提升海雜波中目標恒虛警探測能力，本文結合大量實測數據處理結果，詳細分析了利用Gamma分布對海雜波幅度統計特性建模擬合的可行性，同時介紹了兩種基于Gamma分布海雜波背景下的CFAR檢測器，通過理論分析和計算機仿真驗證了它們在均勻雜波背景和非均勻雜波背景下的檢測性能。結果表明，海雜波的幅度分布可以用Gamma分布進行建模，且提出的兩類CFAR檢測器在各類典型雜波背景下具有良好的目標檢測能力和魯棒的虛警控制能力。

1 Gamma分布與擬合結果

1.1 Gamma分布

Gamma分布的概率密度函數如下所示:

(1)

式中，x為雷達回波幅度，α和β分別為形狀參數和尺度參數，Γ(·)為Gamma函數。相應地，Gamma分布的累計概率密度函數(CDF)可以表示為

(2)

其中γ(s,x)為低階不完全Gamma函數，定義為

(3)

Gamma分布隨機變量的n階矩可以表示為

(4)

作為一種雙參數概率密度函數，Gamma分布族的形狀可以通過改變其中一個參數(例如形狀參數α)來進行靈活調整，如圖1所示。這意味著Gamma分布可以用于描述不同海面環境下的雜波統計特性。由圖1可以看出，當形狀參數較小時Gamma分布概率密度函數出現明顯拖尾，反之則概率密度函數趨于“平坦”。

圖1 不同形狀參數下的Gamma分布概率密度函數

1.2 實測海雜波數據統計分析

本節中，共統計分析來自3個實測海雜波數據庫的樣本，它們分別來自3次大型外場試驗。其中，2個數據庫來源于某X波段對海ISAR試驗雷達(以下簡稱“試驗雷達”)，分別采集于中國山東省青島市和威海市，另外一個數據庫來源于加拿大McMaster大學研制的IPIX相參全極化雷達(以下簡稱“IPIX雷達”)。

1.2.1 試驗雷達實測數據統計結果

本節對實測海雜波的幅度進行概率密度函數模型的擬合優度檢驗(Good-of-Fitness Test)，待測試概率密度函數(PDF)包括瑞利分布、Gamma分布、對數正態分布、K分布、K+N分布、韋布爾分布和帕累托分布。在本例中，瑞利分布、Gamma分布、K分布、K+N分布和帕累托分布的分布參數由矩估計方法獲得，對數正態和韋布爾分布的參數可通過最大似然估計方法(MLE)來獲得[12]。

圖2(a)給出了試驗雷達在青島采集的實測數據的直方圖擬合結果。圖中橫坐標為歸一化幅度，縱坐標為互補累計分布函數(1-CDF)。考慮到K+N分布是K分布的廣義形式，因此在圖2(a)中省略了K分布的擬合結果。從圖2(a)中可以直觀看出，實測海雜波的幅度統計分布可以近似由Gamma分布來擬合。

為了獲得更為精確、定量的擬合結果，通常采用χ2檢驗來進行擬合優度檢驗。χ2檢驗統計量可以表示為

(5)

式中，Nbin為分組數，Ai為觀測頻數，Ti為理論頻數。

表1給出了試驗雷達在青島采集的海雜波數據的參數估計與χ2檢驗統計結果。在這批數據來自6個距離單元，每個距離單元包含3 000個海雜波時域回波樣本。從表1的統計結果可以看出，Gamma分布的χ2檢驗統計結果最小，即Gamma分布在這三者中具有最優的擬合性能。

表1 海雜波參數估計與擬合優度檢驗結果(青島數據)

續表

圖2(b)給出了試驗雷達于山東威海采集的海雜波數據統計結果。這批數據一共包含了來自30個距離單元的海雜波時域回波，每個距離單元共包含16 384個樣本。從圖中可以看出，Gamma分布可以很好地對海雜波幅度進行建模。

(a) 青島實測數據擬合結果

考慮到已有研究證明K+N分布和帕累托分布對實際海雜波的擬合效果優于對數正態分布和韋布爾分布[9,13]，圖3給出了利用Gamma分布、 K+N分布和帕累托分布對在威海采集的海雜波數據的擬合優度結果。圖3的結果表明,Gamma分布具有相對更小的χ2值，即更優的擬合效果。

圖3 Gamma，K+N和帕累托分布擬合優度檢驗結果(威海數據)

1.2.2IPIX雷達實測數據統計結果

本節統計分析IPIX雷達于1993年在Dartmouth采集的實測數據，其雷達參數與數據相關信息見表2。

表2 IPIX雷達參數

為了盡量避免不同的外界環境對統計結果帶來的影響，提高分析結果的代表性，本例中主要研究數據文件#310和#311。這兩批數據文件的采集時間和對應的海況天氣與試驗雷達的數據類似。每個數據文件包含14個距離單元，每個距離單元包含217個時域樣本。其中，由于第6至第9個距離單元中含有合作目標的回波信息，因此在雜波統計時將這幾個距離單元的數據剔除。

圖4給出了IPIX雷達實測數據的直方圖擬合結果，其中橫坐標為歸一化幅度，縱坐標為互補累計分布函數。從圖中可以看出，Gamma分布、K+N分布和帕累托分布皆可以描述真實海雜波的幅度統計模型。表3和表4給出了利用上述3種分布對這兩批數據建模的統計分析與擬合優度檢驗結果。結果表明，Gamma分布相比于K+N分布具有相對更小的χ2值；經計算，帕累托分布的χ2值遠大于Gamma分布與K+N分布，因此在表中將其省略。因此，利用Gamma分布對實際海雜波幅度分布進行統計建模是合理的，這也為后續的目標檢測奠定了基礎。

圖4 IPIX雷達實測數據擬合結果

表3 海雜2檢驗結果(IPIX #310數據)

表4 海雜波χ2檢驗結果(IPIX #311數據)

2 Gamma分布海雜波下的目標檢測

2.1 固定門限檢測

當檢測背景樣本服從如式(1)所示的Gamma分布時，則固定門限檢測器的虛警率可以表示為

式中，T為檢測門限。

將e-t進行冪級數展開，式(6)可以改寫為

(7)

幾種典型形狀參數下的虛警率結果如圖5所示。從圖中可以看出，當門限值給定時，其虛警率對形狀參數仍有依賴性。這意味在實際工作中，雜波背景的形狀參數需要作為先驗信息提供給檢測器。事實上，幾乎在所有雷達系統工作之前，都需要事先對監視海域進行觀測，獲得并分析大量的雜波數據來對雷達系統各項信號處理技術方法進行評估。在雷達系統實際工作時可以事先收集足夠的真實海雜波數據，利用Gamma分布對海雜波進行建模分析并將估計出的形狀參數作為先驗信息；這種處理方式也類似于雜波圖法[14]。此外，當監視區域距離范圍跨度較大時，可以將監視區域按照距離或者估計出的形狀參數分割為若干個子區域[15]，在進行目標檢測時可以根據各個子區域的特點給檢測器配置不同的參數。為不失一般性，本文后續分析皆假設形狀參數先驗已知。

圖5 固定門限檢測器在不同形狀參數下的虛警率

固定門限檢測器在不同形狀參數下的檢測概率如圖6所示。由圖6易知，這類檢測器在不同形狀參數下的檢測性能差異較大。同時，由于實際環境下海雜波的時變特性，這類檢測器幾乎僅存在理論上的可用性，無法適應現代雷達系統的應用。

圖6 固定門限檢測器在不同形狀參數下的檢測概率

2.2 基于最大似然估計的CFAR檢測器

2.2.1 檢測器設計

CFAR檢測器通常采用滑窗技術，其基本流程框圖如圖7所示。經過檢波后的雷達回波樣本序貫存儲在延遲線上，根據具體用途將其分為CUT、參考單元和保護單元。CUT為感興趣的待檢測單元；參考單元中的回波樣本用來估計雜波的平均水平；為了避免CUT中目標的能量泄漏到參考窗內而影響雜波水平估計，在CUT的兩側通常放置若干個保護單元。自適應檢測門限可以表示為雜波水平估計值與乘積因子的乘積，當CUT中樣本的強度超過該門限，則判定目標存在。

圖7 CFAR檢測器基本流程圖

考慮到實測雷達海雜波幅度分布可以用Gamma分布進行建模擬合，本節介紹一種基于最大似然估計的Gamma分布海雜波背景下的CFAR檢測器，即MLG-CFAR[5]。

MLG-CFAR檢測器利用MLE方法對雜波平均功率進行估計。假設參考窗中的樣本是獨立同分布的Gamma隨機變量X1，…，XNr，其中Nr為參考窗的長度，則參考單元樣本幅度的似然函數可以寫為

L(α,β|X1,…,XNr)=

(8)

考慮似然函數L恒為正，因此可以在上式的兩端同時取自然對數，從而得到

lnL(α,β|X1,…,XNr)=

-Nrαlnβ-NrlnΓ(α)+

(9)

尺度參數β的最大似然估計可以通過解如下所示的一元方程獲得:

(10)

由式(10)可以得到尺度參數的估計值為

(11)

mp=Γ(α+p)βp/Γ(α)

(12)

則雜波平均功率的MLE(即MLG估計器)可以表示為

(13)

式中，E[·]為求數學期望，α為已知的形狀參數。易證MLG估計器是對Gamma分布雜波平均功率的一種漸進有效的估計器。

2.2.2 在Gamma背景下的檢測性能

(14)

通過隨機變量的參數變換，可獲得z的概率密度函數為

(15)

MLG-CFAR檢測器對雜波背景的平均功率估計值及其概率密度函數已由式(12)和式(15)給出，則該檢測器在Gamma分布下虛警率表達式為

(16)

式中，TMLG為歸一化乘積因子。式(16)中不包含β，且由于形狀參數α為已知的先驗信息，所以MLG-CFAR的虛警率僅與歸一化乘積因子TMLG有關，這表明MLG-CFAR檢測器對尺度參數β是恒虛警的。

圖8(a)和(b)分別給出了上述檢測器在均勻雜波背景下、不同信雜比時的檢測概率曲線和信雜比為15 dB時的接收機工作曲線(ROC)。本例中，設檢測背景為各向同性的Gamma分布雜波，α=2，β=100。上述各類檢測器的參考窗總長度都為16，左、右兩個保護窗長度都為3，給定期望虛警率為Pfa=10-6。為了提高仿真結果的精度，本例中針對每一組參數都進行了105次獨立蒙特卡洛實驗。由圖8中的結果可以看出，MLG-CFAR相比于上述各類CFAR檢測器都具有最優的檢測性能。

(a) 檢測概率

圖9給出了幾類檢測器在非均勻背景下的檢測性能。非均勻雜波背景下主要包含了雜波邊緣和多目標兩種情況。在雜波邊緣主要關注CFAR檢測器的虛警率控制能力，在多目標情況下主要關注CFAR檢測器對主目標的檢測能力。圖9(a)給出了上述幾種CFAR檢測器在雜波邊緣的虛警率變化曲線。其中，檢測場景共包含48個距離單元，其中第25個距離單元為雜波與瑞利噪聲區域的交界位置。本例中，雜波幅度服從形狀參數α=2、尺度參數β=100的Gamma分布，且雜噪比為10 dB。各類CFAR檢測器的參考窗總長度都為16，左、右兩個保護窗長度都為3，期望虛警率設為10-3。由仿真結果可以看出，在雜波邊緣GO-CFAR具有最優的虛警率控制能力；雖然MLG-CFAR的虛警率控制能力并不讓人很滿意，但是其性能與CA-CFAR和OS-CFAR相當。

為了更有針對性地對檢測性能進行分析與對比，本文結合實際雷達回波，主要考慮了3種多目標情形的假設：1)參考窗內有一個干擾目標，其干雜比(ICR)為10 dB；2)參考窗內有一個干擾目標，該目標具有與CUT內主目標相同的功率；3)參考窗內有5個干雜比為10 dB的干擾目標。為了更好地驗證多目標情況下的“遮蓋效應”，避免引入由于目標RCS起伏帶來的額外干擾項，這幾類目標的幅度都設為非起伏模型。

圖9(b)至(d)分別給出了各類CFAR檢測器在上述3種多目標情況下的檢測性能。在圖9 (b)中，參考窗內存在一個干擾目標，其干雜比為10 dB。由圖中檢測概率曲線可以明顯觀察到，由于“遮蓋效應”的存在，CA-CFAR和GO-CFAR的檢測性能有所下降。OS-CFAR在估計雜波平均功率時進行了排序處理，因此該檢測器在多目標的情況具有較大優勢。本例中，OS-CFAR檢測器中的樣本次序指數設為12，即在參考窗總長度為16時最多可容忍4個干擾目標。考慮到本仿真中僅有一個干擾目標，因此OS-CFAR具有較好的魯棒性。隨著主目標信雜比的逐漸增大，CA-CFAR和OS-CFAR的檢測性能逐漸趨近。相比于其他幾種CFAR算法，當主目標信雜比比較高時，MLG-CFAR具有較大的優勢；當主目標信雜比相對較低時，MLG-CFAR的檢測性能明顯優于CA-CFAR和GO-CFAR，僅略遜于OS-CFAR。

圖9(c)給出了干擾目標的功率與主目標的功率相等時的檢測概率曲線。此時，OS-CFAR具有良好的檢測性能，而CA-CFAR和GO-CFAR由于嚴重的“遮蓋效應”，在進行目標檢測時出現飽和狀態，即無法正確檢測主目標；這種情況可以通過適當增大參考窗長度來緩解。雖然此時MLG-CFAR的檢測性能不盡如人意，但依舊具有一定的魯棒性，不會出現檢測飽和的極端情況。

圖9(d)考慮了參考窗內同時出現5個不起伏的干擾目標的情況，每一個干擾目標的干雜比都是10 dB。考慮到此時干擾目標數目為5個，已經超過了OS-CFAR的容忍范圍，因此OS-CFAR的檢測性能急劇下降。由于GO-CFAR的選大特性，依然會受到“遮蓋效應”的影響而難以達到理想的檢測性能，但由于其雜波平均功率是Nr/2個樣本的算術平均值，干擾目標帶來的影響得到一定的削弱，所以檢測性能相對于OS-CFAR有一定改善。CA-CFAR通過對Nr個雜波樣本的功率求算數平均來獲得雜波平均功率的估計值，樣本數目比GO-CFAR多一倍，因此能夠進一步削弱由干擾目標帶來的影響，獲得相對于GO-CFAR較優的檢測性能也是符合預期的。MLG-CFAR檢測器在上述幾種CFAR檢測器中獲得最優的檢測性能。

(a) 雜波邊緣的虛警率(α=2，β=100)

2.3 加權幅度循環迭代CFAR檢測器

2.3.1 檢測器設計

在非均勻雜波背景下影響CFAR檢測能力的主要原因是此時參考樣本不再滿足獨立同分布的假設，導致檢測器對雜波平均功率的估計出現偏差，進而造成虛警率的上升和檢測概率的下降。如果一個檢測器在估計雜波平均功率時能夠最大程度地削弱異常值、尖峰或者潛在干擾目標的影響，那么可以預測該檢測器在非均勻雜波背景下能具有良好的魯棒性與適應性。基于上述思路，本節介紹了一種基于加權幅度循環迭代的CFAR檢測器(WAI-CFAR)[2]。

WAI-CFAR檢測器對雜波平均水平的估計是一個循環迭代過程，它主要包括兩個步驟：排序和加權平均。假設參考窗中的樣本是獨立同分布的Gamma隨機變量數據集:

(17)

則這些樣本經過升序排列后為

{X(1),X(2), …,X(Nr)}

(18)

X(1)=ηwX(1)+(1-ηw)X(Nr)

(19)

式中，ηw為加權因子且ηw∈[0,1]。將由式(19)獲得的處理結果放回剔除幅度最大和最小的樣本后的數據集，可獲得經過一次循環迭代處理之后的新數據集，即

(20)

對式(20)給出的新數據集繼續進行升序排列，可以獲得升序之后新樣本序列，即

{X′(1),X′(2), …,X′(Nr-1)}

(21)

(22)

WAI-CFAR檢測器中的加權系數ηw可以根據實際需要靈活選取。當ηw=0.5時，一次WAI處理輸出的結果為兩個輸入樣本的幅度算術平均；當ηw>0.5時，幅度較小的樣本在最終輸出結果中貢獻較大；反之，則幅度較大的樣本具有更大的貢獻。本文主要分析了WAI-CFAR檢測器在5種典型的加權系數下的性能，如表5所示。

表5 5種典型的加權系數

2.3.2 在Gamma背景下的檢測性能

圖10給出了WAI-CFAR檢測器在不同Gamma分布參數和不同加權系數下的虛警概率曲線。為不失一般性，本例中檢測器的參考窗總長度為16，左、右兩個保護窗長度都為3。為了更好地體現WAI-CFAR在不同尺度參數下的通用性，仿真分析中主要考慮了β=100,150,1 000,1 500和10 000五種情況。為了保證蒙特卡洛仿真的精度，針對每一種參數組合都進行了109次獨立的蒙特卡洛試驗。

圖10(a)給出了當α=0.5時WAI-CFAR檢測器在不同尺度參數時的虛警率曲線，其中5種典型加權系數由表4給出。由圖10(a)可以看出，當加權系數給定時，WAI-CFAR檢測器在不同尺度參數的Gamma分布雜波中的虛警率曲線幾乎是各自重疊在一起，這說明此時的虛警率僅隨乘積因子的變化而變化，與尺度參數β的取值無關。考慮到蒙特卡洛仿真的精度有限，經過局部放大雖然可以觀察到不同尺度參數下的虛警率曲線并不完全相同，但其微小差異可以忽略不計。

圖10(b)和(c)分別給出了當α為2和10時的虛警率曲線。仿真結果證明，盡管WAI-CFAR檢測器在Gamma分布雜波中的虛警率對形狀參數α有依賴性，但該檢測器在形狀參數給定時對尺度參數具有恒虛警特性。

(a) 形狀參數為0.5

圖11(a)和(b)分別給出了上述幾種檢測器在期望虛警率為Pfa=10-6時的檢測概率以及當目標信雜比為16 dB時的接收機工作曲線。本例通過蒙特卡洛方法獲得上述定量結果，同時為了提高仿真結果的精度，對于每種參數下的檢測概率估計都進行了105次獨立試驗。由結果中可以看出，CA-CFAR檢測器和加權系數為參數2和參數3的WAI-CFAR檢測器的檢測性能相當，且優于OS-CFAR和GO-CFAR檢測器。

圖11 CA,GO,OS和WAI-CRAR在均勻Gamma背景下的性能(α=2，β=100)

圖12(a)給出了上述幾種CFAR檢測器在雜波邊緣的虛警率變化曲線。為保證仿真結果對比的客觀性，本例中檢測場景仿真參數與圖9一致。由仿真結果可以看出，CA-CFAR的虛警率在雜波邊緣急劇上升，OS-CFAR的虛警率控制能力也不盡如人意，而GO-CFAR具有最優的虛警率控制能力。此時，WAI-CFAR檢測器的虛警率變化情況與加權系數的選擇有關。在加權系數為參數3，4和5時，WAI-CFAR的虛警率控制能力較差；當加權系數為參數2時，WAI-CFAR的虛警率在雜波邊緣有一定上升且高于GO-CFAR，但是明顯低于CA-CFAR和OS-CFAR檢測器；WAI-CFAR在加權系數為參數1時的虛警率控制能力整體上優于其他幾種方法。

在圖12(b)中，參考窗內存在一個干雜比為10 dB干擾目標。由結果中可以看出，由于“遮蓋效應”GO-CFAR的檢測性能急劇下降。由于WAI-CFAR檢測器可以利用加權系數進行靈活配置，當加權系數為參數2到參數5時性能較優，當目標信雜比較高時甚至優于OS-CFAR檢測器。

在圖12(c)中，參考單元中包含一個與主目標的功率相等的干擾目標。此時，OS-CFAR在估計雜波平均水平時通過樣本幅度排序將干擾目標的信息剔除，在干雜比較大的情況下也具有良好的檢測性能。由于嚴重的“遮蓋效應”，CA-CFAR檢測器的檢測性能急劇下降。GO-CFAR和加權系數為參數1的WAI-CFAR檢測出現飽和狀態，無法正確檢測主目標。WAI-CFAR在加權系數為參數4或參數5時具有較為魯棒的性能，其定量結果與OS-CFAR檢測器相當。

圖12(d)考慮了參考窗內同時出現5個干雜比為10 dB的干擾目標的情況。此時，由于干擾目標數目超過了OS-CFAR檢測器的抗干擾能力的極限，OS-CFAR的檢測性能急劇下降，只有當主目標的信雜比充分大(例如31 dB)時，才能獲得大于90%的檢測概率；嚴重的“遮蓋效應”同樣導致CA-CFAR和GO-CFAR的檢測性能難以令人滿意。當加權系數為參數4或參數5時，參考窗中幅度較小的樣本對雜波平均水平的估計值貢獻較大，此時即使干擾目標的干雜比較大，但經過加權平均后對最終的輸出統計量不會產生嚴重影響，因此此時WAI-CFAR檢測器具有最優的檢測性能。

圖12 CA,GO,OS和WAI-CFAR在非均勻Gamma背景下的檢測性能

3 結束語

本文首先利用在中國青島、威海采集的X波段海雜波實測數據以及加拿大McMaster大學提供的X波段雷達海雜波實測數據進行統計分析，證明了海雜波幅度分布可以Gamma分布來進行建模，并給出了直觀的直方圖擬合結果和定量的擬合優度檢驗結果。隨后介紹了兩種適用于Gamma分布雜波中的CFAR檢測器，即MLG-CFAR和WAI-CFAR。理論分析和仿真結果表明,這兩類檢測器在Gamma分布形狀參數已知的前提下對尺度參數具有恒虛警特性，且在均勻雜波背景和非均勻雜波背景下相比于幾類經典CFAR算法具有更優的檢測性能。