利用稀疏非負矩陣分解的大轉角SAR成像方法

許 然，李亞超，邢孟道

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

合成孔徑雷達（Synthetic Aperture Radar，SAR）通過發射寬帶脈沖信號并進行方位相干積累，實現二維（距離向和方位向）高分辨成像，在地表測繪和目標識別等領域得到了廣泛應用.為了使SAR圖像包含豐富的特征信息和保證目標的檢測識別概率，最直接的途徑就是提高成像分辨率.近年來寬帶射頻半導體技術的發展使得SAR能夠取得亞米級，甚至厘米級的距離分辨率.而與此同時，為了在方位向獲得與距離向同等級別的分辨能力，則需要比常規SAR更寬角度的方位積累，即進行大轉角SAR成像（Wide Angle SAR，WASAR）[1].

在進行大轉角SAR成像的情況下，成像處理將變得困難.首先，回波數據的二維頻域支撐區將是占據很大角度范圍的同心圓環扇形區域，而不再用矩形網格近似；其次，大的相干積累角度意味著散射點在距離向和方位向的徙動現象嚴重，需要加以補償；此外，由于積累角度很大，必須考慮目標的各向散射異性，即各個散射點的回波能量分布在長方位孔徑上，但并不一致.這些問題使得常規小轉角下的理想假設不再成立，因此不能直接使用如極坐標格式算法（PFA）等基于傅里葉變換的成像算法.針對大轉角SAR成像問題，常用方法是將整個成像間隔分為若干個有重疊區段的子孔徑進行處理，這樣每塊子數據便滿足小轉角的近似假設.文獻[2]認為子塊數據滿足線性等均勻采樣條件，散射點未發生走動，可直接利用快速傅氏變換算法（FFT）進行成像，然后再進行旋轉融合獲得大轉角綜合圖像，但在高分辨率大目標情況下，不徙動條件難以滿足；文獻[1，3]用逆問題對成像處理進行描述，構建帶有正則化稀疏約束的最優化問題，對每塊子孔徑數據分別迭代求解成像再取最大值融合，能夠一定程度地增強二維分辨能力.但該方法會產生虛假點，需要構建復雜的成像投影算子，并且對計算機存儲和計算性能要求很高；文獻[4]針對WASAR中二維信號缺失的問題，提出了一種自適應幅相估計的方法實現超分辨成像，但其譜估計的精度在很大程度上依賴于較高的信噪比.

非負矩陣分解（Non－negative Matrix Factorization，NMF）是在所有元素非負性條件下的一種矩陣分解方法，可以獲得更直觀和更符合人腦認知的加性基表示，其作為一種穩健的特征降維分析工具，被廣泛應用于信號處理和模式識別等領域[5－7].在遙感探測方面，NMF也得到了一定的關注[8－9]，但并沒有在雷達成像領域的應用.筆者提出了一種基于稀疏NMF的大轉角SAR成像算法，首先將全孔徑劃分為互有重疊的子孔徑成像段以滿足短孔徑假設條件，然后采用PFA算法獲得目標在各個不同方位角的子孔徑圖像，再通過構造帶有稀疏增強的正則項約束的NMF目標函數，將子圖像映射到NMF子空間，完成全孔徑圖像綜合與稀疏特性增強.該方法不需要構造成像逆算子，采用改進的乘性迭代算法保證收斂性與非負性，能夠在有效提高圖像信噪（雜）比的同時，獲得特征信息更清晰的目標圖像.

1 大轉角雷達成像模型

圖1 聚束SAR模式示意圖

提高剛體目標的成像時間，可以獲得更豐富的散射信息，有利于目標的檢測和識別.例如，對地面目標可采用高分辨聚束SAR成像模式（如圖1所示），其本質是通過增加相干積累角以提高方位分辨率，并獲得多角度雷達散射截面（RCS）信息.在常規短孔徑近似下，目標區域被認為是由若干各向同性的散射點組成的，不同點的散射系數僅是其空域位置的函數；而在寬角度成像模式下，這種近似將不再合理.因為通常在實際情況下，大多數場景點目標是各向異性的，即散射系數還與方位觀測角存在耦合關系.令雷達發射脈沖為p（t），對于場景中某一點p，其在方位角θ處的回波可表示為

其中，Rp（θ）為目標點與雷達相位中心在方位θ處的斜距，σ（xp，yp；θ）是目標點的散射系數，c是光速，（xp，yp）是點p在斜距平面內的投影坐標.實際上，目標點與雷達視線的擦地角也隨θ變化，對散射系數產生額外的微小調制，但在小場景遠場假設下，這種變化可以忽略.如果雷達發射線性調頻信號為

其中，rect（·）表示矩形窗函數；Tp為脈沖寬度；γ＝B Tp，為調頻率，B為信號帶寬.則雷達接收到的場景總回波在經過去載頻、匹配濾波和運動補償等預處理后，可以表示為

其中，V 代 表雷達照射區域；Kx＝4πfsinθc，為距離向波數；Ky＝4πfcosθc，為方位向波數；θ∈[－Θ2，Θ2]，表示方位觀測角變化范圍；f∈ [fc－B 2，fc＋B 2]，為信號頻率.常規短孔徑成像模式下，信號的二維波數譜體現為角度支撐不大的扇區（通常不超過10°），給定Θ便可以直接利用PFA算法完成極坐標到直角坐標的插值去除距離和方位的耦合，再采用IFFT進行成像.然而在大轉角情況下，這種二維耦合現象嚴重，波數譜占據了很大的角度區間，其支撐區接近于同心圓環（如圖2所示），有效插值區域存在上限，進而出現頻譜利用率較低的情況，使得直接使用PFA有一定的角度限制[10]；反向傳播（Back－Projection，BP）算法雖然沒有使用約束條件，但全孔徑相干積分會使能量平均到整個成像區間，使得弱散射點被抑制，抗噪性能不強，而且全孔徑BP算法運算量較大.為了解決這些問題，通常是犧牲方位分辨率，將整個全孔徑劃分為若干子孔徑以保證PFA性能，同時降低各向異性效應，然后分別成像，獲得不同觀測角下的子圖像，再最終完成旋轉融合.子孔徑數據可以表示為

圖2 WASAR二維波數域

其中，N表示子數據總塊數，Ω表示子孔徑方位積累角，Δθ表示子孔徑間的角度差.若令Δθ＜Ω，則相鄰子孔徑間有Ω－Δθ大小的角度重疊.由于PFA算法會損失一定的頻譜，重疊可以防止子孔徑間的積累角度出現較大斷裂，提高譜利用率.其次，這種劃分能平滑子圖像間的散射系數變化，同時保護較短方位響應范圍的弱散射點，防止其相位歷程被切斷而造成能量弱化和分辨能力降低.選擇子孔徑的長度時，應在滿足PFA使用條件和保證方位分辨率的情況下，盡量縮短孔徑長度以減小各向異性對成像質量的影響.

接下來對各個子孔徑數據分別使用PFA算法進行成像，即將sn（f，θ）進行極坐標到直角坐標系的插值.通常是分別沿著距離維和方位維完成一維插值重采樣，然后再對插值后矩形網格數據進行二維IFFT變換得到子圖像.為了抑制副瓣可以在該步驟中加入窗函數.子圖像代表了目標在不同方位觀測角處的SAR圖像，要獲得全孔徑特征信息，則應對子圖像旋轉到統一坐標系配準后再進行全孔徑圖像綜合.最簡單也最容易理解的方法是，對所有配準子圖像直接求和.這種方式可以直觀地反映目標的特征信息，但弱散射點的能量不能夠得到很好的體現，尤其在低信噪比情況下目標和噪聲都進行了非相干疊加，使得圖像質量不高；文獻[1]采用了基于一般似然比檢驗（GLRT）的方法，將綜合圖像中的像素灰度定義為所有子圖像中該像素位置處的灰度最大值，這種方法獲得的目標圖像特征信息更清晰，但抗噪性仍然不高，且破壞了目標原始的散射能量分布，這對于目標散射信息提取和目標識別十分不利.因此，筆者提出了一種基于稀疏非負矩陣分解的全孔徑圖像綜合方法，在圖像域通過加入稀疏約束的迭代優化，融合各個子孔徑下的目標散射信息，抑制噪聲（雜波）的同時，獲取散射特征信息增強的全孔徑SAR圖像，提高圖像的可解譯能力.

2 基于稀疏NMF的全孔徑圖像綜合

非負矩陣分解是以所有矩陣元素均非負作為約束條件的矩陣分解方法，其問題可以描述為：已知一個非負矩陣V∈Rm×n，求非負矩陣W ∈Rm×r和非負矩陣H ∈Rr×n，使得

其中，r滿足（m＋n）r＜mn.NMF可看做是將原始數據V映射到低維特征子空間中，在減少數據量的同時，獲得基于數據自身的特征基向量的線性表示，而這種表示的基向量本身和權值系數也滿足非負條件.NMF算法得到的基向量組W具有一定的正交性和稀疏性，即體現了原始數據的局部特征.那么由非負的線性組合來近似表達V，則體現了整體由局部組成的直觀認知，具有明確的物理意義.W中基向量的個數r直接決定了特征空間的維數.根據NMF算法的理論，當r＝1時，原始數據V中各個樣本均可由子空間的惟一基向量加權近似表示，那么該基向量應盡可能包含原始樣本集的完整特征.如果令各個配準后的子孔徑圖像作為NMF的特征樣本輸入V，迭代得到惟一基W作為輸出的綜合全孔徑圖像，則該圖像就包含了目標所有方位觀測角下的散射特性.而在NMF算法中加入稀疏約束，可以增強分解結果本身具有的稀疏性，提高散射點的分辨能力，并且抑制旁瓣、噪聲和雜波.

目前，已經有很多關于稀疏NMF的研究：文獻[11]實現了基于NMF的稀疏編碼，但其針對的是分解系數矩陣而不是基矩陣；文獻[12]采用了一種定量控制基矩陣或系數矩陣稀疏度的NMF算法，采用l1范數約束，通過梯度下降法進行迭代求解，但每步迭代都需要復雜的稀疏投影操作，且無法保證非負性.相比較于l2范數和l1范數，加以lk范數正則項（k＜1）約束的解具有更稀疏的結構[13].考慮點散射模型下，大轉角剛體目標的雷達圖像由強散射點組成這一稀疏特性，可采用歐氏距離作為NMF算法的代價函數，并向其中加入對基向量（也就是全孔徑圖像）的lk范數懲罰（0＜k＜1），可得到如下的最優化問題：

其中，V ＝ [v ，v，v，… ， v]，v＝vec（I），I是由sn（f，θ）經過PFA算法得到的第n幅子圖像的灰度矩

123Nnnn陣，vec（·）表示矩陣按列堆疊成列向量操作；W 是 NMF輸出的全孔徑綜合圖像幅度向量，其長度M＝Nr×Ncr，Nr和Ncr分別為距離向和方位向采樣點數.式（6）中的第1項 V －WH用來控制NMF分解結果與原始數據的擬合誤差，αW為稀疏約束正則項，α是懲罰因子，·k代表矩陣元lk范數，定義為大小決定了W的稀疏程度，k越小，則解中非零元素個數越少，也就越稀疏.α和k作為正則化參數，一同起著平衡分解誤差與稀疏程度的作用，應根據實際應用場合進行適當選取.當信噪比較低或雜波較強時，需合理采用較大的α和較小的k以增強圖像的稀疏度，抑制噪聲和雜波的影響.

最常用的NMF算法是基于交替迭代的梯度下降算法，即按照下式進行迭代：

利用矩陣跡函數的性質和求導法則，得到F（W，H）分別對Wia和Haj進行偏導數，即

Lee和Seung令式（7）中的μia和δaj隨著W和H的即時迭代結果進行更新，得到一種乘性迭代規則[14]，該算法易于實現，同時保證了非負性和收斂性能.在此基礎上，根據具體的代價函數形式，得到修改后的乘性迭代規則.令

其中，ε是一個很小的正數，例如可以取ε＝10－7，以避免分母為零的情況出現.將式（9）和式（10）代入式（7），可得迭代規則為

與其他經典NMF算法一樣，在每次W迭代后需對其基向量進行歸一化操作，這里選擇l2范數歸一化：

NMF對初始值并不敏感，通常采用非負隨機矩陣作為初始即可快速收斂.將迭代終止條件設置為

其中，δ＞0是個很小的門限常數.按照式（11）進行迭代，直至滿足式（12）.此時的Wopt經過向量轉換為矩陣的逆操作vec－1（Wopt），即得到稀疏特征增強后的全孔徑綜合圖像.文中方法的全孔徑圖像綜合過程如圖3所示.值得指出的是，子孔徑分割處理和PFA算法不可避免地會降低方位分辨率.但實際上這種綜合圖像的分辨率并不比全孔徑圖像差很多，因為通常目標的方位散射響應范圍有限，而且存在遮擋等情況，這樣成像的分辨率存在上限，不會隨著積累角度無限制增加.而通過對圖像特征的稀疏增強，可以抑制副瓣與雜波噪聲，這在一定程度上改善了散射點的分辨能力.此外，文中方法應用范圍不僅局限于WASAR，還可應用于其他成像場合，例如多基地雷達協同成像[15]和大轉角ISAR成像等，前者通過方位向數據拼接形成長角度孔徑來提高圖像分辨率，后者對非機動或合作目標采用長時間跟蹤成像模式，可以通過運動補償等效為WASAR.這些應用本質上都屬于寬角度積累問題，均可用文中方法進行成像處理.

圖3 全孔徑圖像綜合示意圖

3 實驗結果與分析

為了驗證文中方法的有效性，采用Backhoe公開數據[16]進行大轉角圖像重構實驗.該數據方位合成角為110°，信號頻率范圍為7～13GHz，人為地對原始數據加以－5dB的加性高斯白噪聲.通過計算該數據的二維頻域支撐區范圍，可知全孔徑PFA算法的有效插值區域過小，無法直接適用.圖4（a）給出了全孔徑BP算法的成像結果.可以看出，噪聲與目標能量相當，極大地影響了其特征辨識.將全孔徑分為29個Δθ＝2.9°，每個轉角Θ＝28°的子孔徑，然后采用PFA算法加漢明窗獲得29幅子圖像，對子圖像旋轉配準后，分別使用直接求和疊加方法和文獻[1]中的GLRT方法，得到的綜合全孔徑圖像分別如圖4（b）和圖4（c）所示.與圖4（a）相比較，這兩種子孔徑PFA方法成像結果分辨率稍差，但重疊子孔徑處理避免了散射能量的全局平均，尤其對弱散射點保護較好，目標信息更為豐富；但圖4（b）中的目標特征在噪聲背景下仍不夠清晰，圖4（c）中目標散射特性雖然較好，但噪聲水平依然很高.

圖4（d）給出了采用文中方法迭代得到的綜合全孔徑圖像（α＝0.6，k＝0.8）.可以看出，圖4（d）中目標結構輪廓清晰，分辨能力得到了提高，噪聲和旁瓣被很好地抑制，目標的稀疏特征得到增強.比較這3種子孔徑綜合圖像所對應的灰度向量幅度（圖5（a）～（c）所示），GLRT方法得到的圖像雖然特征清晰，但噪聲也被放大，散射能量分布已被破壞；文中方法在保留散射能量分布情況的同時，令目標圖像的特征信息和信噪比量[12]，得到稀疏度sp＝0.903 8，并將該數值設置為稀疏參數，使用文獻[12]中的稀疏約束NMF方法（NMFSC），得到綜合圖像和相應的灰度向量分別如圖4（e）和圖5（d）所示.與文中方法結果相比，相同稀疏度下，該方法得到的全孔徑綜合圖像噪聲水平較高，目標散射信息的表述能力稍差.為了定量比較圖像質量，分別對不同方法下的WASAR成像結果進行了圖像熵和目標－背景雜波比的衡量比較，其結果如表1所示.可以看出，文中方法結果的熵值最小，TCR最高，表明其目標表征能力最強.另外，筆者還對幾種方法在PC平臺上的運行時間進行了統計對比，文中方法運算時間雖比圖像旋轉求和及GLRT方法的耗時，但仍處于一個數量級，對比最終圖像質量，這種運算時間上的增加是可以接受的.而與全孔徑BP算法和NMFSC方法相比較，文中方法又具有圖像質量和運算時間上的雙重優勢，體現出了良好的性能.

表1 不同方法下的全孔徑圖像質量對比

圖4 不同方法的全孔徑圖像

圖5 不同方法全孔徑圖像灰度矢量

4 結束語

提出了一種基于稀疏非負矩陣分解的WASAR成像方法.該方法首先將全孔徑分為若干重疊的子孔徑，分別進行PFA成像，得到若干不同視角下的子圖像，然后構造帶有稀疏增強正則項的NMF優化問題，利用改進的乘性迭代規則，對旋轉配準后的子圖像在NMF子空間進行融合，獲得最終的綜合全孔徑圖像.Backhoe數據的處理結果證明了文中方法有很強的目標特征增強能力和去噪能力.文中方法實現簡單，但正則化參數的自動選取仍需要進一步研究.

[1] Moses R L，Potter L C，Cetin M.Wide Angle SAR Imaging[C]//Proceedings of SPIE：5427，Algorithms for Synthetic Aperture Radar ImageryⅪ.Bellingham：SPIE，2004：164－175.

[2] Moore T G.A FFT Based Algorithm for the Formation of Wide－angle ISAR Images Using EIP [C]//IEEE International Radar Conference.Piscataway：IEEE，1995：392－395.

[3] Stojanovic I，Cetin M，Karl W C.Joint Space Aspect Reconstruction of Wide－angle SAR Exploiting Sparsity[C]//Proceedings of SPIE：6470，Algorithm for Synthetic Aperture Radar ImageryⅩⅤ.Bellingham：SPIE，2008：697005.

[4] 劉啟，洪文，譚維賢，等.寬角合成孔徑雷達二維缺失數據自適應幅相估計成像方法 [J].電子與信息學報，2012，34（3）：616－621.Liu Qi，Hong Wen，Tan Weixian，et al.Adaptive Tuning Missing－data Amplitude and Phase Estimation Method in Wide Angle SAR [J].Journal of Electronics ＆Information Technology，2012，34（3）：616－621.

[5] 殷海青，劉紅衛.一種基于L1稀疏正則化與非負矩陣分解的盲源信號分離新算法 [J].西安電子科技大學學報，2010，37（5）：835－841.Yin Haiqing，Liu Hongwei.New Blind Source Separation Algorithm Based on L1Sparse Regularization and Nonnegative Matrix Factorization[J].Journal of Xidian University，2010，37（5）：835－841.

[6] Lee D D，Seung H S.Learning the Parts of Objects by Non－negative Matrix Factorization[J].Nature，1999，401（21）：788－791.

[7] Li Y，Ngom A.A New Kernel Non－negative Matrix Factorization and Its Application in Microarray Data Analysis[C]//IEEE Symposium on Computational Intelligence in Bioinformatics and Computational Biology.Piscataway：IEEE，2012：371－378.

[8] Guo R，Zhang L，Xing M，et al.Polarimetric SAR Image Fusion Using Nonnegative Matrix Factorisation and Improved－RGB Model[J].Electronics Letters，2010，46（20）：1399－1401.

[9] Yokoya N，Yairi T，Iwasaki A.Coupled Nonnegative Matrix Factorization Unmixing for Hyperspectral and Multispectral Data Fusion[J].IEEE Transactions on Geoscience and Remote Sensing，2012，50（2）：528－537.

[10] Zhou F，Bai X，Xing M，et al.Analysis for Wide－angle Radar Imaging[J].IET Radar，Sonar and Navigation，2011，5（4）：449－457.

[11] Eggert J，Korner E.Sparse Coding and NMF [C]//IEEE International Conference on Neural Networks.Piscataway：IEEE，2004：2529－2533.

[12] Hoyer P O.Non－negative Matrix Factorization with Sparseness Constraints[J].Journal of Machine Learning Research，2004，4（9）：1457－1469.

[13] Cetin M.Feature－enhanced Synthetic Aperture Radar Imaging[D].Boston：Boston University，2001.

[14] Lee D D，Seung H S.Algorithms for Non－negative Matrix Factorization [C]//Neural Information Processing Systems Conference.Cambridge：Massachusetts Institute of Technology Press，2000：556－562.

[15] Bucciarelli M，Pastina D.Multi－grazing ISAR for Side－view Imaging with Improved Cross－range Resolution [C]//IEEE Radar Conference.Piscataway：IEEE，2011：939－944.

[16] Air Force Research Laboratory.Backhoe Data Sample and Visual Challenge Problem[EB/OL].[2012－08－18].Https：//www.sdms.afrl.af.mil/.