海雜波背景下雷達目標貝葉斯檢測算法

許述文，王喆祥，水鵬朗

(西安電子科技大學 雷達信號處理國家重點實驗室，陜西 西安 710071)

隨著雷達分辨率的提高，中心極限定理的假設不再滿足，海雜波呈現出明顯的非高斯特性，此時再使用高斯背景下的檢測器會造成明顯的虛警或漏檢[1]。高分辨海雜波通常采用符合高斯模型進行建模，復合高斯模型采用兩個獨立的過程描述非高斯雜波：零均值復高斯過程的散斑分量和非負隨機過程的紋理分量。在一個相干處理間隔內，紋理分量可以看作一個隨機常數，此時復合高斯模型退化為球不變隨機向量模型。通過將紋理分量建模為服從不同分布的隨機變量可以得到不同的雜波模型[2]。當紋理分量分別服從伽馬分布、逆伽馬分布、逆高斯(Inverse Gaussian，IG)分布時，雜波的幅度分布服從K分布[3]、廣義帕累托分布[4]和CG-IG分布[5]。這3種幅度模型經常被用來描述不同分辨率下的海雜波。但是，由于海雜波特性受到如極化方式、擦地角、雷達工作頻率等多種條件的影響，以上3種雜波模型的適用范圍不夠廣，同時實測數據也顯示，這3種雜波模型并不能很好地描述不同條件下的海雜波數據[6]。針對這一問題，文獻[6]提出了使用廣義逆高斯(Generalized Inverse Gaussian，GIG)分布來建模紋理分量，通過一定的參數變化可以得到廣義逆高斯分布的3種極限分布：伽馬分布、逆伽馬分布和逆高斯分布。因此，當紋理分量建模為服從廣義逆高斯分布的隨機變量時，可以得到一種更為通用的雜波模型(Compound-Gaussian model with Generalized Inverse Gaussian texture，CG-GIG)。常用的海雜波幅度模型K分布、廣義帕累托分布和CG-IG分布是CG-GIG分布的3種極限分布，所以CG-GIG分布的泛化能力更強，能夠描述低、中、高分辨率的海雜波，這也是K分布、廣義帕累托分布和CG-IG分布所不具備的。實測數據也顯示，使用廣義逆高斯分布來建模紋理分量可以更好地描述實測海雜波數據[6]。

根據上文的描述，CG-GIG模型可以對海雜波進行更為準確的建模。接下來的問題就是如何設計匹配CG-GIG模型的自適應檢測器。在高斯模型下，Kelly根據廣義似然比準則設計了廣義似然比(Generalized Likelihood RaTio，GLRT)檢測器[7]。但是GLRT檢測器的計算量較大。在文獻[8]中，根據兩步法，首先假設協方差矩陣已知設計檢測器，然后使用參考單元數據估計協方差矩陣，并將估計值代入到檢測器中，得到了自適應匹配濾波(Adaptive Matched Filtering，AMF)檢測器。在復合高斯模型下，CONTE提出了一種歸一化匹配濾波(Normalized Matched Filtering，NMF)檢測器，在雜波拖尾很重或者積累脈沖數很大時，NMF的檢測性能趨近于最優[9]。NMF檢測器并不依賴于紋理的先驗分布，所以在不同的紋理下，NMF檢測器會出現一定程度的性能損失。復合高斯雜波背景下的最優檢測器的結構為匹配濾波器的輸出與一個依賴數據項的門限相比[10]。CG-GIG雜波背景下的最優檢測器為GLRT-GIG檢測器[11]，根據RMB準則，當積累脈沖數為N時，只有當參考單元數目大于2N時，檢測器才能有較好的檢測性能[12]。但是在實際的環境中，由于異常數據以及海尖峰等因素的影響，可選取的參考單元數目往往較少[13]。針對這一問題，學者們普遍將先驗知識融入到檢測器的設計中去，減少檢測器的檢測性能對于參考單元數目的依賴。例如當相干脈沖串具有對稱間隔時，協方差矩陣具有雙重對稱性，即關于主對角線對稱的Hermitian結構和次對角線對稱的斜對稱結構[14]。傳統的協方差矩陣估計方法往往只用到了協方差矩陣的Hermitian結構。文獻[15]研究了均勻高斯雜波背景下的斜對稱檢測器，并且推導了檢測概率與虛警概率的數學表達式。文獻[16]研究了在均勻和部分均勻高斯雜波背景下基于協方差矩陣斜對稱結構的子空間目標檢測問題。文獻[17]研究了在K分布雜波背景下的斜對稱檢測器。除了斜對稱結構外，協方差矩陣的Toeplitz結構以及圓結構也經常被用到檢測器設計中[18]。除了利用協方差矩陣的先驗結構外，還可以將散斑協方差視為一個隨機變量，從而利用貝葉斯框架設計檢測器。文獻[19]指出，復高斯隨機向量協方差矩陣逆的共軛先驗分布為逆復Wishart分布。為了在數學上易于處理，通常將散斑協方差矩陣建模為一個服從逆復Wishart分布的隨機矩陣，它包含自由度和均值矩陣兩個參數，自由度用于衡量隨機矩陣與均值矩陣之間的距離，自由度越大，距離越小。文獻[20]設計了非均勻K分布雜波背景下利用散斑協方差矩陣先驗分布的自適應檢測器。文獻[21]中，作者設計了在多輸入多輸出雷達背景下利用散斑協方差矩陣先驗分布的自適應檢測器。文獻[22]中，在逆高斯紋理下，通過迭代估計自由度與紋理分量，設計了一種不依賴輔助數據的貝葉斯檢測器。以上利用協方差矩陣先驗知識的自適應檢測器在參考單元數目較少的情況下，均具有令人滿意的檢測性能。除此之外，自回歸(Auto Regressive，AR)模型[23-24]等參數模型也經常被當作先驗信息用于檢測器的設計。這些先驗信息的使用能夠大大地降低自適應檢測器對于參考單元數目的依賴，從而提升檢測器的檢測性能。除了使用先驗信息外，近年來一些基于數據降維的檢測方法也被廣泛研究，其好處在于不依賴任何的先驗信息，通過對數據進行轉換，大大減少了所需要估計的未知參數，并且在參考單元數目小于脈沖數時也具有不錯的檢測性能[25-27]。

根據貝葉斯框架，筆者在CG-GIG雜波背景下設計了兩種自適應檢測器。將散斑協方差矩陣建模為服從逆復Wishart分布的隨機矩陣，首先根據廣義似然比準則，通過待檢測單元數據來估計自由度與雜波的紋理分量，設計了一個不依賴于輔助數據的貝葉斯檢測器；然后根據最大后驗準則，通過最大似然估計得到紋理分量的估計值，使用輔助數據與先驗知識聯合估計散斑協方差矩陣，得到了一個依賴于輔助數據與先驗知識的貝葉斯檢測器。實驗結果表明，所提出的兩種檢測器在參考單元數目較少時，性能均優于已存在的檢測器，在不同的參考單元數目下，所提出的依賴于輔助數據與先驗知識的貝葉斯檢測器始終具有最優的檢測性能。

1 問題描述

1.1 檢測問題描述

假設雷達在慢時間維積累了N個脈沖，則待檢測單元的信號可以用一個N維的復向量z=[z(1)，z(2)，…，z(N)]T表示，檢測問題可以描述為一個二元假設檢驗問題[22]：

(1)

其中，H0表示目標不存在的假設；H1表示目標存在的假設；α為目標的復幅度；p=[1，exp(j2πfd)，…，exp(j2π(N-1)fd]T，表示目標的導向矢量，fd為目標的歸一化多普勒頻率；zk表示不包含目標信號的參考單元數據；L為參考單元的數目；c和ck為雜波向量。

1.2 信號模型介紹

在一個相干處理間隔(Coherent Processing Interval，CPI)內，紋理分量可以視為一個隨機常數，此時復合高斯模型退化為球不變隨機向量模型。雜波向量c和ck，k=1，2，…，L，可以使用一個正的隨機變量τ(τk)和一個N維零均值復高斯向量u(uk)的乘積表示，數學公式描述如下：

(2)

其中，u表示散斑分量，是一個均值為0的復高斯向量，其協方差R=E(uuH)；τ表示紋理分量。這里使用廣義逆高斯分布對紋理分量進行建模，能夠包含常用的伽馬分布、逆伽馬分布以及逆高斯分布，能更好地擬合實測雜波數據。廣義逆高斯分布的概率密度函數為

(3)

其中，Kp(·)表示p階第二類修正貝塞爾函數。令λ=(ab)1/2，μ=(b/a)1/2，這里μ被稱作尺度參數。

在文中，將散斑協方差矩陣建模為一個服從復逆Wishart分布的一個隨機矩陣，其概率密度函數為

(4)

其中，Σ是一個正定的均值矩陣；Γ(·)是伽瑪函數；v表示自由度，自由度越大，R離均值矩陣Σ的距離越近。通常情況下v是未知的，所以在這里v被建模為一個均勻離散隨機變量[28]。

根據對信號模型和雜波模型的分析，可以得到在H0假設和H1假設下，z的條件概率密度函數為

(5)

1.3 散斑協方差矩陣估計方法

在實際情況中，散斑的協方差矩陣通常是未知的，需要從參考單元的數據中進行估計。常用的散斑協方差矩陣估計方法為歸一化采樣協方差矩陣估計器(Normalized Sample Covariance Matrix Estimator，NSCME)，可表示為[29]

(6)

NSCME對每個距離單元進行了功率歸一化。當ANMF檢測器使用NSCME后，具有對紋理的CFAR特性，但是無法保證對于協方差矩陣結構的CFAR特性。約束漸進最大似然估計器(Constrained Approximate Maximum LikeIihood Estimator，CAMLE)是通過平衡計算復雜度和估計性能提出的迭代估計器，可以保證自適應相干檢測器對雜波紋理和散斑協方差矩陣的CFAR特性，可以表示為[30]：

(7)

2 貝葉斯檢測器設計

2.1 不依賴輔助數據的貝葉斯檢測器

本節設計的貝葉斯檢測器對散斑協方差矩陣R進行積分，所以并不依賴于輔助數據，使用未知參數α的最大似然估計值、τ的最大后驗估計值以及v的最小均方誤差估計值來代替其真實值。檢測器設計準則如下：

(8)

在H1假設下：

(9)

在H0假設下：

(10)

目標復幅度α的最大似然估計,可以通過最大化式(9)得到，也相當于最小化下式，即

L(α)=(z-αp)HΣ-1(z-αp)。

(11)

令式(11)對α求偏導，并令導數為0，可以得到α的最大似然估計為

(12)

前面介紹到，v被建模為一定區間上的均勻離散隨機變量。其MMSE估計通過求其條件后驗期望得到，根據貝葉斯公式，v的條件后驗分布為

(13)

(14)

其中，

(15)

綜上，檢測判決式(8)可以化簡為

(16)

下面計算在H1假設下，τ的MAPE。觀察式(16)的分子，將含有τ的項提出,可以得到：

(17)

對式(17)取對數，然后對τ求偏導，令導數為0，化簡可以得到：

(18)

其中，

(19)

根據笛卡爾符號準則，方程式(18)有惟一正根，并且一元三次方程有解析解。同理，在H0假設下有方程：

(20)

其中，

(21)

(22)

其中，ξ′是修正檢測門限。

2.2 依賴于輔助數據的貝葉斯檢測器

在上一小節中，因為對散斑協方差矩陣進行了積分，所以設計的檢測器雖然能夠使用協方差矩陣的先驗分布信息，但是沒有辦法使用參考單元的信息。為了能夠同時使用先驗信息與參考單元的信息，這里設計了一種同時利用先驗分布信息與參考單元信息的貝葉斯檢測器。

根據最大后驗廣義似然比準則，檢測判決式可表示為

(23)

目標復幅度α的最大似然估計為

(24)

使用兩步法，首先假設R已知，然后再將R用其估計值代替，就得到了自適應檢測器。

在H1假設下，計算τ的MAPE。在式(23)的分子中，提出與τ有關的項可得到：

(25)

(26)

對式(25)取對數，然后對τ求偏導，并令導數為0，解得τ的MAPE為

(27)

同理在H0假設下，τ的MAPE為

(28)

其中，q0=zHR-1z。

將式(24)、(27)、(28)代入式(23)，就可以得到非完全自適應的檢測器：

(29)

下面融合參考單元與先驗分布信息來估計散斑協方差矩陣，從而得到完全自適應的檢測器。與2.1節類似，在第k個參考單元，自由度v的條件后驗密度函數為

(30)

在第k個參考單元上，自由度v的MMSE為

(31)

其中，

(32)

(33)

(34)

將式(31)、(34)代入式(33)，然后取對數，并對散斑協方差矩陣R求偏導，可以得到

(35)

根據矩陣求導法則，可以得到

(36)

(37)

將式(36)、(37)代入式(35)中，并令結果為0，經過化簡得到

(38)

式(38)可以寫為

(39)

其中，

(40)

(41)

3 性能評估

3.1 BGLRT-GIG檢測器的性能評估

首先使用仿真數據來評估BGLRT-GIG檢測器的性能，將文中所提出的BGLRT-GIG檢測器與GLRT-GIG檢測器[11]和GLRT-R檢測器[20]進行對比，均值矩陣Σ建模為具有一階遲滯相關系數ρ的指數相關型協方差矩陣，Σi，j=ρ|i-j|，1≤i，j≤N。對于海雜波ρ∈[0.90，0.99]，這里ρ=0.95，雜波參數設置為λ=5，p=2，u=Kp(λ)/Kp+1(λ)，目標的歸一化多普勒頻率假設已知，且fd=0.2，自由度參數v的上下界分別設置為vl=N+1 和vu=3N，虛警概率設置為Pfa=10-3，檢測門限通過n次獨立的蒙特卡羅實驗得到。其中n=100/Pfa，信雜比(SCR)定義為

(42)

其中，P=μKp+1(λ)/Kp(λ)，為雜波的平均功率。

圖 1展示了不同參考單元數目下，GLRT-GIG、BGLRT-GIG和GLRT-R檢測器的性能曲線。從圖中可以看出，BGLRT-GIG的檢測性能優于GLRT-R的檢測性能，這是因為GLRT-R檢測器并沒有考慮紋理的先驗分布。當參考單元數目較小時，GLRT-GIG檢測器的性能損失較大，BGLRT-GIG檢測器的性能優于GLRT-GIG檢測器的性能，這是因為參考單元數目隨著參考單元數目的增大，GLRT-GIG檢測器的性能會優于BGLRT-GIG檢測器的性能。且隨著參考單元數目的增加，散斑協方差矩陣的估計逐漸準確，此時只依賴協方差矩陣的先驗信息是不夠的，這也是文中設計依賴輔助數據的貝葉斯檢測器的原因。

圖1 不同參考單元數量L下的檢測性能曲線圖

圖2展示了信雜比為-16 dB時，歸一化多普勒頻率對檢測性能的影響。當目標落在主雜波區時，3種檢測器都出現了一定程度的性能損失，隨著目標距離雜波區越來越遠，檢測器的檢測性能逐漸變好。從圖2可以看出，當參考單元數目為14時，在不同的歸一化多普勒頻率下，文中提出的BGLRT-GIG檢測器性能均優于GLRT-GIG檢測器和GLRT-R檢測器。

圖2 信雜比為-16 dB和L=14時歸一化多普勒頻率和檢測性能曲線圖

圖3展示了自由度對于BGLRT-GIG檢測性能的影響。從圖中可以看出，隨著自由度的增大，檢測性能也逐漸提升，這是因為自由度衡量的真實協方差矩陣與先驗矩陣之間的距離，自由度越大，距離越近，先驗矩陣也就更準確。

圖3 不同自由度下檢測性能曲線圖

圖4展示了虛警概率對檢測器檢測性能的影響。隨著虛警概率的不斷提高，兩個檢測器的檢測性能逐漸變好，在不同的虛警概率下，筆者提出的BGLRT-GIG檢測器均擁有最優的檢測性能。

圖4 不同虛警概率下L=14時的性能曲線圖

接下來，使用南非實測海雜波數據評估兩個檢測器在不同參考單元數目條件下的檢測性能。實驗中采用TFC15-001數據集，該組數據集8 h的平均風速為8.14 m/s，8 h平均浪高為5.31 m，雷達的載頻為 9 GHz，脈沖重復頻率為5 kHz，距離分辨率為15 m。圖5展示了不同雜波模型對于實測數據的擬合情況，參數估計方法使用文獻[6]所提出的方法。從圖中可以看出，采用廣義逆高斯紋理建模的復合高斯模型能夠很好地描述這組實測數據。

圖5 TFC15_001數據幅度分布擬合圖

在實測數據中，將第24個距離單元視為待檢測單元，其余數據視為歷史數據，使用NSCM估計器，每16個距離單元估計一次散斑協方差矩陣，最后取平均，當作先驗矩陣。圖 6展示了在實測數據中BGLRT-GIG檢測器與GLRT-GIG檢測器的檢測概率曲線。從圖中可以看出，當參考單元數目較少時，BGLRT-GIG檢測器的性能較好，當參考單元數目比較充足時，GLRT-GIG檢測器的性能更好。

3.2 BMAP-GIG檢測器的性能評估

針對參考單元數目增多時，GLRT-GIG檢測器優于文中提出的BGLRT-GIG檢測器的問題，筆者又提出了一種依賴參考單元數據和協方差矩陣先驗信息的BMAP-GIG檢測器。本小節主要對BMAP-GIG檢測器的性能進行評估，并將其與GLRT-GIG、BGLRT-GIG檢測器進行對比。

首先使用仿真數據對BMAP-GIG檢測器的性能進行評估，實驗參數設置與上一小節相同。圖7展示了在不同參考單元數目下，BGLRT-GIG、GLRT-GIG和BMAP-GIG檢測器的性能曲線。

圖7 不同參考單元數目的仿真數據的性能曲線圖

從圖7中可以看出，隨著參考單元數目的增加，GLRT-GIG檢測器的性能不斷變好，最終優于BGLRT-GIG檢測器。在不同的參考單元數目下，文中提出的BMAP-GIG檢測器始終具有最優的檢測性能。這是因為BGLRT-GIG檢測器只使用了協方差矩陣的先驗信息，GLRT-GIG檢測器只使用了參考單元數據的信息，而BMAP-GIG檢測器同時考慮了先驗信息與參考單元數據的信息，所以在這3種檢測器中，BMAP-GIG檢測器擁有最優的檢測性能。

下面使用南非實測海雜波數據對這3種檢測器進行性能評估，使用的數據集為TFC15_001。在圖 5中已經展示過廣義逆高斯紋理建模的復合高斯模型對這組數據的擬合結果。 門限使用純雜波數據進行計算。圖 8展示了在實測數據中，不同參考單元數目下，3種檢測器的性能曲線。從圖中可以看出，當參考單元數目較少時，GLRT-GIG檢測器的性能損失較大，性能最差。隨著參考單元數目的增加，GLRT-GIG檢測器的性能逐漸提升，最后會優于BGLRT-GIG檢測器。在不同的參考單元數目下，BMAP-GIG檢測器始終具有最優的檢測性能。

圖8 不同參考單元數目的實測海雜波加仿真目標的性能曲線圖

由于BMAP-GIG檢測器的判決式較為復雜，可使用蒙特卡羅實驗的方法研究BMAP-GIG檢測器的虛警性能。圖 9分別展示了檢測器門限對于歸一化多普勒頻率、雜波功率以及雜波一階相關系數的CFAR特性。從圖中可以看出，BMAP-GIG檢測器的檢測門限是近似獨立于這3個參數的。由于BGLRT-GIG檢測器中，需要進行迭代運算，所以相對于GLRT-GIG檢測器計算復雜度較高；在BMAP-GIG檢測器中，需要估計每個參考單元的自由度，并且還需要迭代估計協方差矩陣。所以在GLRT-GIG、BGLRT-GIG、BMAP-GIG檢測器中，BMAP-GIG檢測器計算復雜度最高，BGLRT-GIG檢測器次之，GLRT-GIG檢測器計算復雜度最低。

圖9 BMAP-GIG檢測器的虛警性能

4 結束語

筆者針對參考單元數目較少時，傳統自適應檢測器性能損失較大的問題，在CG-GIG雜波背景下，根據貝葉斯準則，通過將散斑協方差矩陣建模為服從復逆Wishart分布的隨機矩陣，提出了兩種貝葉斯檢測器。其中BGLRT-GIG檢測器不依賴于參考單元信息，BMAP-GIG檢測器同時使用了先驗信息與參考單元的信息。仿真數據實驗和實測數據實驗的結果表明，當參考單元數目較少時，BGLRT-GIG具有較好的檢測性能，并且優于GLRT-GIG檢測器。由于BMAP-GIG檢測器同時使用了先驗信息與參考單元的信息，所以在不同的參考單元數目下，BMAP-GIG檢測器的性能均優于BGLRT-GIG和GLRT-GIG檢測器。