廣域監視模式下新的雜波加噪聲譜密度矩陣估計方法

閆 賀 鄭明潔 李 飛 艾加秋②

①(中國科學院電子學研究所 北京 100190)

②(中國科學院研究生院 北京 100039)

1 引言

對寬廣的區域存在的動目標進行檢測、定位、跟蹤在民用和軍用領域均有迫切的需求，在民用領域，可以實現重點區域的交通流量控制；在軍用領域，可以提供戰區的部隊調動情況，為戰場指揮官提供重要的決策信息。

廣域監視模式是一種通過天線的電掃描或機械掃描來實現對大范圍的、感興趣的區域中存在的運動目標進行全天時、全天候監視的工作模式。這種工作模式相對于條帶 SAR/GMTI模式有自己的特點：可以多次從多個角度照射到運動目標，提高運動目標的檢測概率，使得對目標運動軌跡的描述成為可能[1]。美國的 AN/APY_3(JSTARS),AN/APY-8(Lynx),TeSAR,HiSAR,ARL-M，德國的PAMIR，英國的ASTOR，以色列的EL/M-2055，北約的 SOSTAR-X等雷達系統都具有廣域監視的工作模式，這些雷達系統對地面動目標進行檢測時面臨的最大問題都是如何抑制較強的地雜波。現在已經知道美國的JSTARS系統采用的是雜波抑制干涉方法(CSI)[2]而德國的 PAMIR系統采用的是Scan-MTI算法[3-5]。CSI算法是一種通道對消的算法，而 Scan-MTI算法是一種空時處理算法，兩種算法在處理實際數據中均取得了較好的結果。

本文對廣域監視模式進行了研究，安排如下：在第2節，證明了廣域監視模式下，當滿足一定的條件時，每個多普勒單元的雜波譜密度矩陣的秩均為1；第3節在第2節的基礎上，推導了一種新的雜波加噪聲譜密度矩陣估計方法，并分析提出：通過第i-1和第i+1個多普勒單元的信息來估計第i個多普勒單元的雜波加噪聲譜密度矩陣結合特征矢量子空間投影的方法可以大大減弱目標污染樣本的影響；第4節仿真了兩種波束指向下場景的回波，并采用上述方法和常規譜密度矩陣求逆方法對場景回波進行處理，顯示了本文所提方法在處理目標污染樣本時的優越性；在第5節中得出結論。

2 廣域監視模式下雜波特性分析

Ender[6,7]指出在SAR-GMTI模式下，采用1DT算法時[8]，如果方位向PRF滿足Nyquist采樣定律，則每個多普勒單元對應的雜波譜密度矩陣的秩為1(文獻[9-11]進一步討論了該算法在SAR-GMTI模式中的應用)，但該結論的前提是方位向時間帶寬積必須比較大。在廣域監視模式下，波束在每個掃描波位上停留的時間非常短，方位向時間帶寬積不可能滿足該條件，因此不能就此認為廣域監視模式下每個多普勒單元的雜波譜密度矩陣的秩也為 1。下面對廣域監視模式下每個多普勒單元的雜波譜密度矩陣的秩進行研究，并且假設方位向PRF滿足Nyquist采樣定律。

與SAR成像模式下的回波類似，廣域監視模式下的每個通道的回波仍然可以看成系統的沖激函數和隨時間變化的地面后向散射系數之間的卷積，即

ci(t)表示第i個通道(相對于參考通道)的回波。通過空間幾何關系(如圖1所示)可知第i個通道經過距離壓縮后的沖激函數可以表示為

圖1 廣域掃描時的幾何關系圖

式(2)中的Ta表示方位向發射脈沖重復時間，N表示在一個掃描波位上發射的脈沖數。由空間幾何關系可知，地面靜止點和載機的實時距離可以近似表示為

式中V表示載機飛行的速度，Ri0表示地面靜止點和載機的初始距離，Ri(t,V)表示地面靜止點和載機的實時距離。從式(2)可以看到，系統的沖激函數由 3部分組成，其中第3部分是由天線的雙程加權和通道的位置差異(相對參考通道)而導致回波在幅度和相位上的變化，由天線的知識可知

式(4)中，Di(θ(t))表示天線的雙程加權，ri表示第i個通道相對參考通道的距離差。通常廣域監視系統在設計的時候，一個掃描波位的駐留時間是非常短的，ai(θ(t))的變化很小，可以認為是不變的，這里改用ai(θ)表示。

對式(2)進行傅里葉變換可得

根據方位向多普勒頻率和方位角度的一一對應關系，ai(θ)可以寫成

其中fθ表示θ對應的多普勒頻率。將式(6)代入式(5)，并將結果寫成向量的形式(這是對所有通道的聯合表達式)，可知

在進行方位向FFT的處理中，根據每個波位上的方位向采樣點數N對系統PRF進行分割，這樣相鄰兩個頻率采樣點的間隔為PRF/N，根據方位向多普勒頻率和方位角度的一一對應關系可知這實際上是將場景在方位向上進行分割，可以把方位向分割的每個塊產生的回波等效成位于該塊中心的點產生的回波，這種等效處理在方位向采樣點越多的情況下越精確[12]。利用這種等效處理，可得方位向FFT變換后第k個多普勒單元對應的表達式應該為

其中c(fk)表示與多普勒頻率fk有關的復常數。需要說明的是為了簡化公式推導，上述考慮的是不加權的方位向FFT變換，實際采用1DT算法處理時使用的是方位向加權FFT變換，這樣處理可以大大抑制來自其他多普勒單元的雜波，減輕空域濾波負擔，但同時會導致該多普勒單元雜波譜展寬。

根據Ender在文獻[13]中的結論(假設雜波在空間上的分布是平穩的)，第k個多普勒單元對應的雜波譜密度矩陣為

其中CP(fk)表示雜波的功率譜，γ(fk)表示與多普勒頻率fk有關的復常數。通過式(9)可知第k個多普勒單元對應的雜波譜密度矩陣秩為1，如果讓k遍歷所有的多普勒單元，可得到結論：在廣域掃描時，如果方位向PRF滿足Nyquist采樣定律且每個掃描波位上的脈沖數足夠多，則任意多普勒單元對應的雜波譜密度矩陣的秩均為 1。根據此特點，我們可以采用特征矢量子空間投影方法來計算濾波器的權系數[4]。

如果把噪聲分量(設噪聲功率為σ2，且噪聲是白噪聲)也考慮進去的話，并認為雜波和噪聲之間是不相關的，可以得到第k個多普勒單元的雜波加噪聲譜密度矩陣。

其中σ2表示噪聲經過距離壓縮和方位 FFT處理后的功率。

3 新的雜波加噪聲譜密度矩陣估計方法

根據第2節中的結論可以得到一種新的譜密度矩陣估計方法：通過某個多普勒單元的雜波加噪聲譜密度矩陣來估計另一個多普勒單元的雜波加噪聲譜密度矩陣，下面給出推導過程。

假設通過最大似然估計方法從回波數據(RD域)中得到第k個多普勒單元的雜波加噪聲譜密度矩陣為

其中Nr為待處理的距離單元數，Z(nr,fk)表示第r個距離單元第k個多普勒單元的回波矢量。

如果想通過第k個多普勒單元的雜波加噪聲譜密度矩陣來估計第j個多普勒單元的雜波加噪聲譜密度矩陣，根據第2節中的推導可知R(fk)和R(fj)可寫成如下形式：

則需要在第k個多普勒單元的雜波加噪聲譜密度矩陣上補償相應的幅度和相位因子。由文獻[13]可知，如果存在這樣的轉換矩陣

其中“./”表示向量或矩陣的點除運算。也就是說它滿足

那么我們可以得到(因為γ(fk)≈γ(fj))

則通過第k個多普勒單元估計的第j個多普勒單元的雜波加噪聲譜密度矩陣為

也就是說對通過最大似然方法估計的第k個多普勒單元的雜波加噪聲譜密度矩陣乘上相應的轉換矩陣，再對噪聲功率分量進行適當的調整，就可以近似得到第j個多普勒單元的雜波加噪聲譜密度矩陣。整個信號處理流程如圖2所示。

圖2 新的雜波加噪聲譜密度矩陣估計方法流程圖

另外，在上述處理過程中，需要知道天線方向圖的頻率域表達式，它是從回波數據中估計得到的[1]。

(fk)表示的是第i個通道第k個多普勒單元對應的天線方向圖的頻域表達式的值。

對于 SAR/GMTI系統廣域監視模式，由于距離分辨率比一般的預警雷達高的多，在進行空時處理時，會有充足的距離單元樣本來估計雜波加噪聲譜密度矩陣。但如果對每個距離單元都單獨進行雜波加噪聲譜密度矩陣的估計會使運算量急劇增加，所以通常對廣域掃描數據作如下的簡化處理：選取一段距離單元，并用該段內所有的距離單元作為樣本來估計一個雜波加噪聲譜密度矩陣，然后將它用于該段內全部距離單元的雜波抑制[3-6]。文獻[12]指出如果含有目標導向矢量的信號混入到訓練樣本中相對于功率的非均勻性對自適應濾波器性能的影響更大。據此本文提出新的雜波加噪聲譜密度矩陣估計算法：通過第i-1個和第i+1多普勒單元的信息來估計第i個多普勒單元的雜波加噪聲譜密度矩陣，結合式(11)-式(17)，可知其處理過程如下：

其中轉換矩陣為

需要說明的是，我們并沒有采用奇異檢測算法剔除被目標污染的樣本[14]，而是直接用所有待測的距離單元作為樣本，通過大量的樣本的平均作用來減弱目標污染樣本的影響，接著通過第i-1個和第i+ 1 個多普勒單元的信息來估計第i個多普勒單元的雜波加噪聲譜密度矩陣，這樣處理的好處是：如果用第i個多普勒單元的數據來估計自身的雜波加噪聲譜密度矩陣，則存在于第i個多普勒單元的運動目標會導致自適應濾波輸出的雜波剩余增加而目標信號減弱，引起虛警率提高而目標檢測性能下降，但通過第i-1個和第i+1個多普勒單元的信息來估計第i個多普勒單元的雜波加噪聲譜密度矩陣可以大大減弱存在第i個多普勒單元的運動目標的影響，同時這相當于把用來估計雜波加噪聲譜密度矩陣的樣本數增加了一倍。

在此基礎上，根據第2節的結論，在方位向PRF滿足Nyqiust采樣定律的時候，可以通過特征矢量子空間投影的方法來設計第i個多普勒單元的雜波抑制濾波器矩陣[4]。

其中v(fi)表示R′(fi)最大特征值對應的特征向量。從理論上講，該濾波器可以將雜波抑制到噪聲水平。

則用于該段所有距離單元的第i個多普勒單元的濾波器可以表示為

其中b(θs)表示天線指向對應的波束形成器。

4 算法仿真及分析

為了驗證上述雜波加噪聲譜密度矩陣估計方法的有效性，本文在計算機上進行了仿真實驗，系統的仿真參數和PAMIR系統基本類似，見表1所示。這里選取了波束指向為0°(正側視)和15°兩種情況下場景的回波(假設地面后向散射系數服從高斯分布)，在每個距離單元上按照一個多普勒單元對應一個方向角，在每個方位角上布置1個點目標，所以在每個距離單元上都布置了128個點，并截取了第1000-1299距離單元共300個距離單元的回波數據。

表1 仿真系統參數

在兩種情況下的場景回波中都加入了幅值相同的6個動目標信息，且都位于波束指向的中心，經過距離壓縮和方位向 FFT后的SNR≈20 dB,CNR≈30 dB，他們的速度和所處距離單元見表2。

表2 回波中加入的動目標信息

為了表述方便，這里把通過第i-1個和第i+1個多普勒單元的信息估計第i個多普勒單元的雜波加噪聲譜密度矩陣并聯合特征矢量子空間投影的方法，稱為算法1，把直接通過第i個多普勒單元的所有距離單元來估計自身的雜波加噪聲譜密度矩陣并對該矩陣直接求逆的方法，稱為算法 2，將波束指向0°時的場景，稱為場景1，將波束指向15°時的場景，稱為場景2，在處理仿真回波時采用的是1DT處理算法。

把5個通道接收到的回波數據進行距離壓縮和方位向FFT變換到RD域，圖3(a)，圖3(b)分別表示了通道3接收到的場景1，場景2的RD域回波功率圖，從圖中可以看出6個動目標均被雜波所淹沒。

圖4(a)、圖4(b)分別表示場景1和場景2的回波經過算法1進行雜波抑制后的功率圖，從兩幅圖中可以看出雜波得到了有效的抑制，6個動目標信息均凸顯出來，從而表明了算法1的有效性。

為了顯示算法 1在處理目標污染樣本時的優勢，這里將算法1和算法2的處理結果進行了對比，如圖5，圖6所示。分別選取波束指向0°時第1個目標的距離向切面圖和波束指向15°時第5個目標的距離向切面圖，并考察了樣本數目依次選取 300,100,50時的情況。從圖5(a)，圖5(b)中可以看出，采用算法2進行雜波抑制后，目標的相對功率(相對噪聲功率)隨著樣本數目的減少而降低，在圖5(a)中，目標1的相對功率分別為：18.5 dB(300),14.5 dB(100),12.5 dB(50)，在圖5(b)中，目標5的相對功率分別為：18.8 dB(300),16.2 dB(100),12.9 dB(50)。這是因為當樣本數目比較多時，含有運動目標信息的樣本的影響會因較多樣本數目的平均作用而被大大削弱，但是這種影響會隨著樣本數目的減少而逐漸顯現出來。而從圖6(a)，圖6(b)中可以看出，采用算法1進行雜波抑制后，目標的相對功率幾乎不隨樣本數目的減少而變化，在圖6(a)中，目標1的相對功率均為20 dB，在圖6(b)中，目標5的相對功率均為19 dB，可見算法1有效地避免了目標污染樣本的影響。

圖3 通道3接收到兩個場景的回波功率圖

圖4 采用算法1處理后回波功率圖

圖5 算法2處理時不同樣本情況下目標切面圖

5 結論

本文首先證明了廣域監視模式下進行雜波抑制時，如果方位向采樣點數足夠多，在方位向PRF滿足Nyqiust采樣定律時，每個多普勒單元的雜波譜密度矩陣的秩為1，然后推導了通過第k個多普勒單元的雜波加噪聲譜密度矩陣來估計第j個多普勒單元的雜波加噪聲譜密度矩陣的計算公式，在此基礎上提出通過第i-1個和第i+1個多普勒單元的信息來估計第i個多普勒單元的雜波加噪聲譜密度矩陣的方法，聯合特征矢量子空間投影可以大大減弱目標污染樣本對空時處理性能的影響，最后結合PAMIR系統的參數，驗證了上述方法的有效性。

圖6 算法1處理時不同樣本情況下目標切面圖

[1]Cerutti-Maori D,Klare J,Brenner A R,et al..Wide-area traffic monitoring with the SAR/GMIT system PAMIR[J].IEEE Transactions on Geoscience and Remote Sensing,2008,46(10)：3019-3030.

[2]Entzminger J N,Senior J R,Fowler C A,et al..JointSTARS and GMTI：past,present and future[J].IEEE Transactions on Aerospace and Electronic Systems,1999,35(2)：748-761.

[3]Cerutti-Maori D,Gierull C H,and Ender J H G.Experimental verification of SAR-GMTI improvement through antenna switching[J].IEEE Transactions on Geoscience and Remote Sensing,2010,48(4)：2066-2075.

[4]Wilden H and Brenner A R.The SAR/GMTI airborne radar PAMIR：technology and performance[C].IEEE MTT-S,Anaheim,California,May 23-28,2010：534-537.

[5]Ender J H G,Berens P,Brenner A R,et al..Multi channel SAR/MTI system development at FGAN：from AER to PAMIR[C].European Conference on Synthetic Aperture Radar,Toronto,2002：1697-1701.

[6]Ender J H G.Space-time processing for multichannel synthetic aperture radar[J].Electronics＆Communication Engineering Journal,1999,11(1)：29-38.

[7]Ender J H G.Detection and estimation of moving target signals by multi-channel SAR[C].European Conference on Synthetic Aperture Radar,K?nigswinter,Germany,1996：411-417.

[8]Ward J.Space-time Adaptive Processing for Airborne Radar[M].NASA STI/Recon Technical Report N,1994：109-119.

[9]魏俊,孫進平,袁運能,等.一種基于頻域 STAP處理的多通道 SAR-GMTI算法[J].系統工程與電子技術,2008,30(1)：1-5.Wei Jun,Sun Jin-ping,Yuan Yun-neng,et al..New multi-channel SAR-GMTI frequency STAP algorithm[J].Systems Engineering and Electronics,2008,30(1)：1-5.

[10]齊維孔,禹衛東.多載頻星載MIMO-SAR系統動目標檢測技術研究[J].電子與信息學報,2010,32(10)：2365-2370.Qi Wei-kong and Yu Wei-dong.Research on ground moving target indication for spaceborne MIMO-SAR of multiple carrier-frequencies[J].Journal of Electronics＆Information Technology,2010,32(10)：2365-2370.

[11]張曉光,劉祥峰,付琨.基于 MIMO-SAR體制的空頻域自適應動目標檢測技術研究[J].電子與信息學報,2010,32(8)：1997-2001.Zhang Xiao-guang,Liu Xiang-feng,and Fu Kun.Study on ground moving target indication of multiple input multiple output space-borne SAR[J].Journal of Electronics＆Information Technology,2010,32(8)：1997-2001.

[12]吳迪,朱岱寅,朱兆達.基于改進 FRACTA 算法的多通道SAR 動目標檢測技術[J].電子與信息學報,2010,32(9)：2201-2207.Wu Di,Zhu Dai-yin,and Zhu Zhao-da.Moving target detection for multi-channel SAR based on improved FRACTA algorithm[J].Journal of Electronics＆Information Technology,2010,32(9)：2201-2207.

[13]Ender J H G.Subspace transformation techniques applied to multi-channel SAR/MTI[C].IEEE International Geoscience and Remote Sensing Symposium,Hamburg,Germany,July 1999：38-40.

[14]吳建新.機載相控陣雷達STAP及目標參數估計方法研究[D].[博士論文],西安電子科技大學,2009.