改進的基于特征子空間的SAR圖像射頻干擾抑制算法

周春暉 李 飛 李 寧 鄭慧芳 王翔宇②

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

②(中國科學院大學 北京 100049)

1 引言

合成孔徑雷達(Synthetic Aperture Radar,SAR)是一種微波遙感設備，在遙感領域具有舉足輕重的地位。由于其具備優良的全天時、全天候、寬測繪帶、遠程成像等能力，已被廣泛應用于軍事和民用等諸多領域，如地質勘探、地形測繪、災害評估、海洋應用、軍事偵察、農林檢測[1,2]。頻率越低，信號覆蓋和滲透性能越好，越有利于環境資源的檢測，以及軍事監視。然而，低頻段內同時存在移動通信信號、電視信號、廣播信號等其它信號[3]，它們會對SAR信號產生干擾，即所謂的射頻干擾(Radio Frequency Interference, RFI)。在SAR圖像中，RFI通常沿距離向以密集的雨滴或明亮的條紋的形式出現，會對有用信號產生許多不良影響，如降低信干噪比(Signal to Interference plus Noise Ratio, SINR)，淹沒潛在的感興趣目標，嚴重降低SAR成像質量，在高分辨率情況下該影響尤其嚴重。因此，射頻干擾檢測和抑制對SAR圖像應用具有重要的意義。

到目前為止，國內外學者已經提出了眾多方法來抑制RFI，這些方法通常可以分為兩類：參數化方法和非參數化方法。參數化方法，通常是基于將RFI看作是一系列正弦信號疊加的假設，通過特定的數學方法[4,5]來估計正弦干擾信號的振幅、頻率和相位，然后從混有干擾的回波信號中減去估計出的正弦信號，進而實現干擾抑制。該類方法包括最小均方法(Least Mean Square, LMS)[6,7]，最大似然估計法[8]等。在最優的情況下，利用參數化方法，可以準確估計出RFI信號模型參數，實現其與假設模型參數的完全匹配。然而，在實際情況下，射頻干擾的數學模型相當復雜，或多或少會出現模型失配，產生模型參數估計誤差，進而導致RFI估計不準確。此外，參數化方法需要額外的參數估計過程，大大降低了信號處理效率。其中，LMS自適應濾波法是一種典型的參數化方法，通過使干擾信號的方向上能量最小化來減少RFI。但是，該方法不能滿足SAR實時成像的要求，同時可能會導致大量有用信號的損失。非參數方法，可以在沒有任何先驗知識或參數模型的情況下實現RFI抑制。多數非參數化方法都可以利用原始數據自適應地實現RFI抑制，其中最簡單最典型的方法就是陷波法[9]。陷波法首先檢測包含干擾的信號頻譜，辨別出RFI信號的尖峰脈沖，然后通過頻譜置零的方式來實現RFI抑制。陷波法可以簡單快速地實現RFI抑制，但是也會給SAR成像結果帶來不利影響，例如信號頻譜的不連續，距離向分辨率的降低，信噪比(Signalto-Noise Ratio, SNR)的損失，等等。子空間濾波法將RFI投影到干擾子空間[10]，然后從原始信號中去除投影部分，進而實現RFI抑制。然而，該方法對整塊數據進行處理，會損失部分有用信號。一些非參數方法可以基于成像后的復數據來實現，子帶譜對消(Subband Spectral Cancelation, SSC)法是其中典型的一種[11]。SSC法不需要進行干擾檢測和估計，通過SAR圖像不同距離子帶頻譜對消來實現RFI抑制。但是如果頻譜不具有對稱性，就會導致有用信號的大量損失。

本文提出了一種應用于SAR圖像的改進的特征子空間射頻干擾抑制算法。首先，進行初步的RFI檢測，判斷數據中是否存在干擾。若存在干擾，則進行后續的RFI檢測，在頻域和時域分別判斷干擾所在的區域。最后，在前兩步的基礎上，利用傳統的特征子空間法來去除原始數據中的干擾信號。仿真和實測實驗的處理結果表明，本文所提出的方法處理效率更高，效果更好。

本文將圍繞改進的特征子空間法展開論述。第2節介紹了射頻干擾的特點。第3節詳細介紹了我們提出的方法。第4節給出了實驗結果和相應的分析。最后，第5節給出本文的結論。

2 RFI的特性

在SAR系統中，包含噪聲和RFI的目標信號可用如下的數學模型表示

其中S表示有用信號，I代表RFI信號[12,13]，N代表噪聲信號。

由于有用信號和噪聲信號均可視為白噪聲信號，因此可以直接將上述公式表示為

圖1為原始數據(機載L波段SAR系統)的距離向頻譜，其中圖1(a)為完整頻譜圖，圖1(b)和圖1(c)分別為圖1(a)中的紅框(b)和紅框(c)內的部分的放大。沿距離向的亮條紋即為RFI，這些RFI各不相同，有的在帶內，有的在帶外，有的具有周期性，有的僅出現在某些回波內。然而，它們有一個共同特征––窄帶性。

對實測數據進行分析，可以看出回波信號的時域波形受到窄帶射頻干擾的嚴重影響。一方面，由于SAR天線與射頻干擾源的相對位置的變化和脈沖重復時間的差異，很難確定窄帶射頻干擾的絕對功率。另一方面，對多個目標回波的疊加導致回波的波形在時域上是無序的，因此很難僅靠觀察時域波形來判斷是否存在干擾[14]。

事實上，由于射頻信號的窄帶特性，其能量高度集中在頻率域，所以我們可以利用該特點進行頻域檢測，這樣可使干擾的識別和判斷更容易。

圖2是由圖1中的數據沿方位向疊加而得到的。通過對實測數據的統計分析可知，SAR數據中窄帶RFI的帶寬分布在0.1 MHz到10 MHz的范圍內，這要比SAR系統的帶寬小得多。與射頻干擾的帶寬相比，可以認為SAR系統的工作帶寬比前者大得多，所以通常在后續處理中將射頻干擾視為窄帶干擾信號。

在該L波段機載SAR實測數據中，主要干擾所在頻段為1335–1345 MHz，產生原因可能為航空無線電導航、無線電定位、衛星無線電導航(地對空)。周期干擾所在頻段為1377–1390 MHz，產生原因應該為無線電定位。

3 算法介紹

為了實現干擾抑制，首先需要檢測和識別RFI信號。傳統的特征子空間法只是在原始數據的頻域進行粗略檢測，以確定RFI在該數據中是否存在，之后對存在RFI的數據進行后續處理。而事實上，干擾在頻域具有窄帶特性，在時域也具有區域局部性，對整塊數據進行特征值分解處理會帶來有用信號的損失。同時，特征值分解處理非常耗時，對整塊數據進行處理效率很低。本文提出了改進的特征子空間法，主要分為3部分：初步RFI檢測，分別在距離頻域和方位時域進行特定的RFI檢測，RFI抑制。其中，初步RFI檢測是對整塊數據中是否存在干擾進行檢測，在存在干擾的情況下進行后續處理；距離頻域RFI檢測利用了支持向量回歸(Support Vector Regression, SVR)算法，得到干擾所在的距離頻域范圍；方位時域RFI檢測則是通過向量RFI_detect得到干擾所在的方位時域范圍；最后，只對存在干擾的數據區域進行特征值分解處理，將得到的處理結果與不含干擾的數據部分合并，進行后續成像。該方法的主要步驟流程圖如圖3所示。該方法的細節部分將在后面繼續討論。應該注意的是，在本文中，SAR圖像的橫向代表方位向，縱向代表距離向，即一行代表一個距離門，一列代表一條距離線。

3.1 初步RFI檢測

為了避免無意義的RFI抑制處理，需要進行初步RFI檢測，判斷數據中是否存在RFI。首先，對原始數據X沿距離向進行快速傅里葉變換(Fast Fourier Transform, FFT)，得到原始距離頻域方位時域數據XF。然后，計算XF每一行(沿方位向)幅度值的平均值，得到一個列向量XF_Mean，如下：

其中，k代表距離門序號，K代表距離門總數，XF由原始數據X經距離向FFT得到。可以從XF_Mean中獲取距離頻域信息。通過對比XF_Mean中各個元素的幅度值，可以很容易地實現初步RFI檢測。

3.2 分別在距離頻域和方位時域進行特定RFI檢測

在進行初步RFI檢測之后，并且確定該原始距離頻域方位時域數據XF存在干擾的情況下，需要分別在距離頻域和方位時域進行特定RFI檢測。

3.2.1 距離頻域RFI檢測首先需要進行特定的距離頻域RFI檢測，即需要檢測干擾所在的距離向頻譜范圍。可以利用多種方法來進行頻域RFI檢測，例如SVR等方法[15,16]。下面簡單介紹一下SVR方法。

假設給定訓練數據｛(x1，y1)，···，(x?，y?)｝?χ×R ，其中χ表示輸入空間，訓練數據即為圖2所用到的數據。在ε-SV回歸中，我們的目標是尋找一個擬合函數f(x)，使得對于任意輸入xi，都能實現函數輸出f(xi)與真實輸出yi相比最多有的偏差，同時使擬合函數的結構風險最小。

首先對線性函數f的條件進行描述，該函數具有以下形式：

其中，〈·，·〉代表在中進行點積。為使式(5)中的函數結構風險最小，我們要尋找一個盡可能小的權重。為此，可以使范數最小化，即‖w‖2=〈w，w〉。

我們可以將該問題轉化為凸優化問題：

其中常數C＞0，它決定了函數f的結構風險與最大允許偏差ε之間的權衡。這對應于處理一個-不敏感損失函數|ξ|ε，如下：

我們可以利用核函數將SVR應用于非線性的情況下。常用的核函數包括多項式核、高斯核(Radial Basis Function, RBF)、拉普拉斯核、Sigmoid核等[17]。

在這里，我們選擇高斯核進行處理。輸入數據為圖2的結果(向量數據)，利用采用高斯核的SVR算法處理，得到RBF模型(也是向量數據)。圖2的結果與得到的RBF模型相減，可以得到擬合誤差，公式如下：

其中，eF代表擬合誤差，RS代表圖2的結果(向量數據)，RBF代表得到的RBF模型(向量數據)。

為擬合誤差eF設定閾值，高于閾值的判斷為干擾部分，即可判斷出RFI所在的距離向頻譜范圍。

在確定出干擾所在的頻譜范圍之后，可以將數據XF直接在距離頻域方位時域空間上劃分為兩部分，XF1和XF2。其中，XF1包含干擾，XF2不含干擾。

3.2.2 方位時域RFI檢測此外，還需要特定的時域RFI檢測(針對包含干擾的數據XF1)。可以依次檢測XF1中每條回波數據中是否存在干擾，以便進行后續RFI抑制處理。具體來說，對數據XF1沿距離向求幅度平均值，得到行向量dL；對數據XF沿距離向求幅度平均值，得到行向量dT。可以定義參數向量Ie如下：

其中Ie代表判斷每條回波數據是否存在干擾的參數向量，dL代表對數據XF1沿距離向求幅度平均值而得到的行向量，dT代表對數據XF沿距離向求幅度平均值而得到的行向量，Th代表設定的的閾值。

之后，可以得到向量dRFI

其中，dRFI用于表示每條回波數據中是否存在RFI,Ie代表判斷每條回波數據是否存在RFI, sign代表用于判斷參數符號的符號函數(相關參數大于0時，結果為1；相關參數等于0時，結果為0；相關參數小于0時，結果為–1)。

然后，綜合頻域和時域檢測的結果，可以確定出干擾在數據XF1中存在的范圍，得到距離頻域方位時域數據XF1t和XF2t，其中XF1t包含干擾，XF2t不含干擾。分別對XF1t與XF2t沿距離向進行IFFT，得到2維時域數據X1t與X2t(由于IFFT為線性變換，X1t與X2t之和為原始2維時域數據X)，其中X1t包含干擾，X2t不含干擾。

3.3 RFI抑制

在經過上述處理之后，我們成功完成了RFI檢測，然后利用傳統的特征子空間法[18]對包含干擾的2維時域數據X1t進行后續的RFI抑制處理。

整個處理流程包括數據分塊，相關矩陣估計，特征值分解，干擾子空間構建，數據投影和數據重排。

假設xm代表第m個距離向采樣點(m=1，2，···，M)，代表某條回波數據向量，其中(·)T表示轉置。

首先，將包含干擾的2維時域數據X1t分成K個L維子向量，其中K=M-L+1。第k(k=1，2，···，K)個子向量定義為：

然后，可以構建一個L×K的數據矩陣，即包含干擾的2維時域數據X1t，如下：

我們可以利用前面的公式來得到干擾的估計協方差相關矩陣，如下：

如果下面的不等式成立：

其中r(r＜L)代表主特征值的個數。可以利用主特征值對應的特征向量，來構建特征子空間，如下：

從原始數據之中去除掉求得的干擾數據，即可得到不含干擾的數據，有如下形式：

之后，我們可以構建一個新的不含干擾的數據矩陣：

3.4 算法總結與說明

(1) 距離頻域檢測部分：

(a) 距離頻域RFI檢測(即SVR)的輸入數據為圖2的結果(向量數據)，利用采用高斯核的SVR算法處理，得到RBF模型(也是向量數據)。

(c) 在確定出干擾所在的頻譜范圍之后，可以將數據XF直接在距離頻域方位時域空間上劃分為兩部分，XF1和XF2。其中，XF1包含干擾，XF2不含干擾。

(2) 方位時域檢測部分：

(a)對XF1進行方位時域檢測，進一步縮小包含干擾的數據范圍，得到包含干擾的部分XF1t，以及不包含干擾的部分XF2t(XF1t與XF2t之和為XF)。

(b) 分別對XF1t與XF2t沿距離向進行IFFT，得到2維時域數據X1t與X2t(由于IFFT為線性變換，X1t與X2t之和為原始2維時域數據X)。

(3) RFI抑制部分：

(a) 對包含干擾的2維時域數據X1t進行第3.3節的RFI抑制(即利用傳統特征子空間法處理)，得到2維時域數據X1ts。

(b) 最后，將X1ts與X2t重新組成數據矩陣X′ ，進行后續正常的成像。

(4) 需要特別說明的是：

(1) FFT和IFFT均為線性變換，先求和后變換與先變換后求和對處理結果沒有影響。

(2)XF,XF1,XF2,XF1t與XF2t均為距離頻域方位時域數據，X,X1t,X2t,X1ts與X′ 均為2維時域數據。

4 實驗結果與分析

在該節中，分別利用仿真實驗和機載SAR實測數據進行驗證，實驗結果證明了本文所提算法在SAR圖像干擾抑制中的有效性。

4.1 仿真實驗

為了初步證明提出的改進方法的有效性，首先進行了仿真實驗。仿真實驗的主要系統參數如表1所示。在該仿真實驗中，添加了干信比為40 dB的單頻干擾，然后分別利用傳統的特征子空間法和改進的特征子空間法進行處理。傳統方法和改進方法的處理結果分別如圖4和圖5所示，它們相應的處理結果的關鍵指標如表2所示。

表1 仿真實驗的主要系統參數Tab. 1 Main system parameters for experiment

通過對比圖4和圖5的點目標仿真結果，可以看出與傳統的特征子空間法相比，改進的特征子空間法在距離向上旁瓣降低，顯著提升成像質量。由表2可知，與傳統的特征子空間法相比，改進的特征子空間法在距離向上PSLR降低了約1.5 dB，ISLR降低了約1.7 dB。經改進的特征子空間法處理后得到的PSLR和ISLR都與理想點目標的指標基本一致，進而證明了改進算法的有效性。

4.2 實測實驗

利用實測數據對改進的特征子空間法進一步進行驗證。我們用到的L波段機載SAR實驗數據由中國科學院電子學研究所提供[19]。成像場景是海岸線區域。機載實驗系統的主要系統參數如表3所示。

完整的原始SAR圖像如圖6所示。我們可以從中看出整幅圖像明顯受到了RFI的影響，嚴重影響了后續的應用。

表2 距離向/方位向關鍵指標的仿真實驗結果(dB)Tab. 2 Simulation results of key indicators in range /azimuth direction (dB)

表3 機載實驗系統的主要系統參數Tab. 3 Main system parameters for experiment

為了進一步證明提出的改進方法的有效性，利用實測數據進行了對比實驗。首先選取圖6中由紅色方框標記的區域“A” ，該區域明顯受到了RFI的影響；然后我們分別利用傳統特征子空間法和改進的特征子空間法來對該區域(原始數據)進行處理，處理結果如圖7所示。圖7(a)為未進行RFI抑制的原始圖像。我們可以從中看出很多亮條紋，這些亮條紋是由RFI生成的，沿距離向分布，嚴重影響成像質量。圖7(b)為利用傳統特征子空間法進行RFI抑制后的圖像。可以看出RFI成功得到抑制，然而也同時丟失了很多細節信息，特別是海岸線區域。圖7(c)為利用改進的特征子空間法進行RFI抑制后的圖像。可以從中看出RFI成功得到抑制，同時與傳統方法相比，圖像中明顯保留了更多的細節信息。換句話說，圖7(c)的成像質量明顯優于圖7(a)和圖7(b)。利用相同的方法對圖6中的區域“B” 進行處理，可得圖8，對比結果可以得到與圖7相同的結論。

在運算效率方面，與傳統特征子空間法相比，改進的特征子空間法可以顯著提升運算效率。其原因是改進方法增加了預處理過程，選出包含RFI的距離頻域方位時域數據，進而將后續處理范圍限定在相對較小的數據范圍內。

從理論上看，特征子空間處理中，處理的是一個2維矩陣。需要將每一個列向量分成若干段(有重疊)，然后分別進行處理，最后還要還原回矩陣，運算量十分大。改進的特征子空間法只在檢測出的存在干擾的部分進行特征值分解處理，可以節省大量時間。而檢測算法的運算效率都很高。其中，初步RFI檢測和時域RFI檢測都只有幾行代碼。而頻域RFI檢測(SVR檢測)處理的對象是一個向量，運算量很小，與后續干擾抑制處理相比，也可以忽略不計。

從實驗結果來看，在我們的實驗中，CPU型號：Intel(R) Xeon(R) CPU E5405 @ 2.00 GHz 2.00 GHz。SVR檢測時間：12.871571 s。傳統特征子空間法時間：40366.133685s。改進的特征子空間法時間：10442.629862 s。處理數據的尺寸為6144×32768。

理論和實驗結果均可說明，新方法比傳統方法運算效率更高。

5 結束語

本文針對SAR圖像RFI抑制，提出了一種改進的特征子空間法。在特征值分解處理之前添加了預處理模塊，即在距離頻域和方位時域分別進行RFI檢測。與傳統特征子空間法相比，改進的特征子空間法可以在更短的時間內得到更好的成像效果，即同時提高了成像效果和運算效率。改進的特征子空間法可以廣泛應用于工作在低波段的SAR系統中，以便處理RFI問題。

[1]Goriachkin O V. Azimuth resolution of spaceborne P, VHF-band SAR[J].IEEE Geoscience and Remote Sensing Letters, 2004, 1(4): 251–254. DOI: 10.1109/LGRS.2004.833777.

[2]鄧云凱, 趙鳳軍, 王宇. 星載SAR技術的發展趨勢及應用淺析[J]. 雷達學報, 2012, 1(1): 1–10. DOI: 10.3724/SP.J.1300.2012.20015.Deng Yun-kai, Zhao Feng-jun, and Wang Yu. Brief analysis on the development and application of spaceborne SAR[J].Journal of Radars, 2012, 1(1): 1–10. DOI: 10.3724/SP.J.1300.2012.20015.

[3]Meyer F J, Nicoll J B, and Doulgeris A P. Correction and characterization of radio frequency interference signatures in L-band synthetic aperture radar data[J].IEEE Transactions on Geoscience and Remote Sensing, 2013, 51(10): 4961–24972. DOI: 10.1109/TGRS.2013.2252469.

[4]Nguyen L and Soumekh M. Suppression of radio frequency inteference (RFI) for synchronous impulse reconstruction ultra-wideband radar[C]. Proceedings of SPIE 5808,Algorithms for Synthetic Aperture Radar Imagery XII,Orlando, Florida, USA, 2005: 178–184.

[5]Ojowu Jr O O and Li Jian. RFI suppression for synchronous impulse reconstruction UWB radar using RELAX[J].International Journal of Remote Sensing Applications, 2013, 3(1): 33–46.

[6]Lord R T and Inggs M R. Efficient RFI suppression in SAR using LMS adaptive filter integrated with range/Doppler algorithm[J].Electronics Letters, 1999, 35(8): 629–630. DOI:10.1049/el:19990437.

[7]Luo X, Ulander L M H, Askne J,et al.. RFI suppression in ultra-wideband SAR systems using LMS filters in frequency domain[J].Electronics Letters, 2001, 37(4): 241–243. DOI:10.1049/el:20010153.

[8]Ulug B. An Algorithm for Sinusoidal Interference Reduction Using Iterative Maximum Likelihood Estimation Techniques[M]. 1992.

[9]Buckreuss S. Filtering interferences from p-band SAR data[C]. Proceedings of EUSAR’98 Conference, Friedrichshafen, Germany, 1998.

[10]Ruan Hang, Ye Wei, Yin Can-bin,et al.. Wide band noise interference suppression for SAR with dechirping and eigensubspace filtering[C]. Proceedings of 2010 International Conference on Intelligent Control and Information Processing, Dalian, China, 2010: 39–42.

[11]Feng Jin, Zheng Hui-fang, Deng Yun-kai,et al.. Application of subband spectral cancellation for SAR narrow-band interference suppression[J].IEEE Geoscience and Remote Sensing Letters, 2012, 9(2): 190–193. DOI: 10.1109/LGRS.2011.2163150.

[12]Zhou Feng, Tao Ming-liang, Bai Xue-ru,et al.. Narrowband interference suppression for SAR based on independent component analysis[J].IEEE Transactions onGeoscience and Remote Sensing, 2013, 51(10): 4952–4960.DOI: 10.1109/TGRS.2013.2244605.

[13]Zhou Feng and Tao Ming-liang. Research on methods for narrow-band interference suppression in synthetic aperture radar data[J].IEEE Journal of Selected Topics in Applied Earth Observations and Remote Sensing, 2015, 8(7):3476–3485. DOI: 10.1109/JSTARS.2015.2431916.

[14]鄭慧芳, 楊淋, 馮錦. SAR窄帶干擾抑制的子帶子空間濾波技術研究[J]. 電子與信息學報, 2013, 35(12): 2836–2842. DOI:10.3724/SP.J.1146.2013.00201.Zheng Hui-fang, Yang Lin, and Feng Jin. Research on the subband subspace filtering for narrow band interference suppression in SAR[J].Journal of Electronics&Information Technology, 2013, 35(12): 2836–2842. DOI:10.3724/SP.J.1146.2013.00201.

[15]Smola A J and Sch?lkopf B. A tutorial on support vector regression[J].Statistics and Computing, 2004, 14(3):199–222. DOI: 10.1023/B:STCO.0000035301.49549.88.

[16]Huang Guang-bin, Zhou Hong-ming, Ding Xiao-jian,et al..Extreme learning machine for regression and multiclass classification[J].IEEE Transactions on Systems,Man,and Cybernetics,Part B(Cybernetics), 2012, 42(2): 513–529.DOI: 10.1109/TSMCB.2011.2168604.

[17]周志華. 機器學習[M]. 北京: 清華大學出版社, 2016: 128–139.Zhou Zhi-hua. Machine Learning[M]. Beijing: Tsinghua University Press, 2016: 128–139.

[18]Zhou Feng, Wu Ren-biao, Xing Meng-dao,et al..Eigensubspace-based filtering with application in narrowband interference suppression for SAR[J].IEEE Geoscience and Remote Sensing Letters, 2007, 4(1): 75–79. DOI:10.1109/LGRS.2006.887033.

[19]Li Ning, Wang R, Deng Yun-kai,et al.. Autofocus correction of residual RCM for VHR SAR sensors with light-small aircraft[J].IEEE Transactions on Geoscience and Remote Sensing, 2017, 55(1): 441–452. DOI:10.1109/TGRS.2016.2608423.