基于壓縮感知的下視三維SAR成像新方法

李學仕 孫光才 徐 剛 邢孟道 梁 毅

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

1 引言

傳統的合成孔徑雷達(Synthetic Aperture Radar,SAR)只能得到場景的 2維像，現實中，各種各樣的應用需求都需要雷達具有下視3維SAR成像能力，如特殊地形的勘探、測繪、城市規劃以及環境破壞的影響分析等[1,2]。傳統的單基雷達無法工作于下視，這是由于飛行航跡左右兩側的點多普勒特性一致，跨航向方向沒辦法得到分辨。為了獲取這一維的分辨，通常需在跨航向維擺放一個線性陣列或等效陣列。然而，跨航向維的分辨率與線陣長度密切相關，這就要求陣列天線長度不能太短。由于陣元間距必須小于1/2波長，才能保證天線方向圖不會出現柵瓣，因此在孔徑長度不能太短的情況下，陣元數就會非常多。然而在實際應用中，由于線性陣列天線設計的成本和實現的復雜度等問題，使得我們無法使用這樣大型的陣列。

近年來，壓縮感知(Compressed Sensing,CS)的概念被提出，它使得信息理解和獲取等方面發生了革命性的變化，各國掀起了 CS理論在信號處理應用等領域的研究熱潮[3,4]。由 CS理論可知，在滿足有限等距性質(Restricted Isometry Condition,RIC)準則的情況下，通過構造有效的稀疏表征基利用求解最小l1-范數正則化問題，目標就可以從有限的觀測樣本中以極大的概率得到完美恢復，并且具有良好的抗噪性能[5-7]?；谶@些特性，CS技術已經推廣到與SAR相關的諸多領域，文獻[8]針對在逆合成孔徑雷達(Inverse Synthetic Aperture Radar,ISAR)系統中，方位采樣缺損的情況下，通過構造方位稀疏的 ISAR成像模型，建立稀疏表征域及基空間表達式，通過優化方法來求稀疏解，該方法對實測數據的處理獲得了良好的結果。文獻[9]針對條帶體制下的SAR成像，提出基于場景方位向小波稀疏表示的壓縮感知成像方法，實測數據成像結果表明所提算法具有較好的有效性。文獻[10]把CS技術運用到3維層析SAR中，通過實測數據處理驗證了CS方法在3維層析SAR中應用的合理性。從這些研究成果中可以看出，CS在ISAR成像、SAR成像、3維SAR層析等方面具有獨特的優勢。本文正是基于 CS的思想，在采用跨航向維稀疏陣列配置下，研究下視3維SAR成像新方法。

本文內容安排如下：第1節為引言；第2節建立了下視3維SAR成像幾何模型，進行了雙通道等效單通道處理的公式推導，并討論了多發多收線陣天線的工作形式以及下視3維SAR成像場景稀疏性分析；第3節介紹了基于CS理論的下視3維SAR成像方法；第4節給出了整個成像處理流程；第 5節對 CS理論強調要滿足的兩個條件進行了詳細的分析；第6節通過仿真數據處理，驗證了本文提出算法的有效性；第7節對全文進行了總結。

2 下視3維SAR成像幾何模型

下視3維SAR成像幾何模型如圖1所示，載機平臺高度為h，沿航向速度為v，設載機位于O點時的時刻為沿航向慢時間tm的起點，線陣天線陣元沿跨航向稀疏分布，由上一節分析知，在實際應用中由于線陣天線設計的成本及復雜度等問題，要求陣元在跨航向維采用稀疏布陣的方式。

2.1 雙通道等效單通道處理

這里下視3維SAR成像采用多發多收模式，根據等效相位中心原理，當收發分置的兩個天線之間基線長度相對于作用距離很短時，可以等效為兩者連線的中心作自發自收的單通道工作模式，但需要補償一個常數相位[11]，這樣處理可以避免由于發射天線、接收天線之間存在一定的水平基線而導致回波信號斜距歷程表現為雙根號形式。

圖1 下視3維SAR成像幾何模型

如圖1所示，經過tm時間后，任意一發射天線的位置坐標可以表示為(xT,vtm,h)，任意一接收天線的位置坐標可以表示為(xR,vtm,h)，則場景中點目標P的收發斜距分別為

收發天線的等效相位中心到P點的距離為

則等效之前與等效之后收發雙程的距離差為

從式(5)可以看出，補償的相位依賴于場景中不同的點目標位置，但通常取載機正下方的點對應的垂直距離來代替R0，由此引入的誤差通常是可以忽略的，如線陣天線長度為8 m，發射信號波長為0.03 m，載機平臺高度為10 km，跨航向波束寬度為20o，則由此導致的相位誤差小于π/100，這個誤差在成像處理中是可以忽略的，因而可以對成像場景進行統一的相位補償。

2.2 多發多收線陣天線的信號形式

采用等效相位中心原理，雖然可以簡化系統回波的收發模型，但是并不能增加天線有效陣元的數量。為了基于有限的陣元，獲得更多的等效陣元數，這里采用多發多收模式，在少量陣元數情況下可以獲得更多的有效陣元。多發多收線陣天線有兩種工作方式，如圖2所示，所有陣元都同時具有發射和接收功能。第1種如圖2(a)所示，發射天線不同時刻順序發射信號，全部接收天線同時接收回波信號，當所有發射天線依次執行完一次發射時稱為一個周期。第2種如圖2(b)所示，雙發射天線同時發射正交或近似正交的兩路信號，比如正負調頻信號，在接收端通過校正處理分離正交的兩路信號，從而區分出不同發射陣元發射的信號，同樣，當所有發射天線都執行完一次發射時為一個周期。為了獲得較多的等效相位中心，這兩種工作方式，都需要把一個周期內的等效相位中心校正到某個參考時刻處。從信號處理的角度可以知道，不同發射時刻的信號等效認為是存在一個時間的延時，將兩維的相位中心投影到某參考時刻可以認為是在航向維做時間的平移，可以通過在航向維對應的頻域乘線性相位操作來完成。經過處理后，就得到一個排列相對緊密的線陣，如圖2中的參考位置。此時可以認為相對緊密的陣列在對地面目標進行觀測。

圖2 多發多收線陣天線信號形式示意圖

2.3 下視3維SAR成像稀疏性分析

下視3維SAR成像場景沿航向的一個切面如圖3所示，在跨航向維，除極少數連續的成片目標外，大多數目標與跨航向維分辨率相比，都可視為點目標，而且，成片的目標如綠色植被等通常認為是連續的背景，可以將其視為噪聲處理。下視3維SAR成像實際是對場景目標的輪廓進行成像，這樣對每個高度單元而言，如圖中的虛線所示，場景都滿足稀疏性，上述因素決定下視3維SAR場景是可以被稀疏重建的。

3 下視3維SAR成像算法

經過等效相位中心處理，就可以認為收發雙通道做單通道自發自收。下面將主要分析稀疏等效相位中心陣列配置下的下視3維成像算法。

3.1多通道航向維處理

經過第2節的處理后，可以認為是在參考位置處一個排列相對緊密的線陣在做自發自收操作，這樣tm時刻第n個等效相位中心的位置坐標可以表示為(xn,vtm,h)，其中xn為第n個等效相位中心在跨航向的位置，則等效相位中心到P點的距離可以表示為

圖3 高度-跨航向2維場景平面

對每一個陣元而言，其錄取的信號可以認為是正側視情況下的回波數據。針對每個陣元的多次回波信號，可以采用傳統的成像處理算法對其進行成像處理。由于傳統成像算法比較成熟，加之本文的重點是跨航向維的稀疏處理，這里對每個陣元的航向維成像不做更多的敘述。這里假設每個陣元的信號都經過RD算法處理。經過成像處理后，各陣元信號可以表示為

其中B為發射信號帶寬，Ba為航向維帶寬。經過航向維處理之后，下面將分析高度和跨航向維處理。

3.2 高度和跨航向2維處理

將式(8)變換到高度維頻域得

其中Ar表示發射信號窗函數ar(·)的頻域形式，fr為高度維頻率，fc為發射信號載頻，將式(9)變換到跨航向維頻域可以得到兩維傅里葉表達式。由于跨航向維采用稀疏布陣的方式，所以對式(9)進行跨航向維的傅里葉變換要利用傅里葉變換的定義，設傅里葉變換矩陣為ψ，其表達式形式為

將式(9)變換到2維頻域后，進行高度和跨航向維的徙動校正、二次脈沖壓縮后，變換到2維時域進行跨航向維dechirp處理，然后進行跨航向維傅里葉變換即完成了成像處理過程，其中徙動校正函數和二次脈沖壓縮函數分別為

其中φ為點目標到陣元的斜視角，Rs為選定的參考距離。下面把式(9)之后的操作寫成矩陣表示形式為

其中x表示成像結果，ψ表示跨航向維傅里葉變換矩陣，ω1～ω6分別表示徙動校正，二次脈沖壓縮，高度維逆傅里葉變換，跨航向維逆傅里葉變換，跨航向維dechirp處理，跨航向維傅里葉變換操作對應的矩陣形式。實際處理時，是把x與S表示為長列向量的形式，這樣上述的操作都可以寫成矩陣左乘S的形式，如式(14)：

其中Ψ表示ψ對應的操作，Ψ的矩陣形式是以ψ為對角元素的一個大矩陣，Ω1～Ω6表示ω1～ω6對應的操作，同樣都是以矩陣ω的元素構成的大矩陣，并用A表示上述所有的操作，則式(14)可以寫成下面的表達式：

這里，S可以看作是觀測信號，x看作為要恢復的原始場景，這樣把式(15)寫成下面的形式：

其中Φ表示A的反操作，由CS理論知，x可以通過l0-范數最小化進行重建。

這是一個NP難題，一般地，采用l1-范數來逼近l0-范數，從而獲得稀疏解[12,13]。此時問題變為

在高度不確定的系統中，要獲得有效解，CS理論強調要滿足兩個條件[6,7]：第一，場景要滿足稀疏性，由第2節分析知，這一點是滿足的。第二，觀測矩陣要滿足RIC準則。觀測矩陣的RIC特性與等效相位中心的擺放有關，因此我們可以通過設置收發天線的位置來得到較好的結果，這將在第4節中單獨進行討論。于是式(18)的求解可以轉化為優化問題求解，構造代價函數為

式(19)中，第1項為觀測誤差，第2項為點目標特性增強，系數α主要平衡兩者之間的關系，其值的大小可以通過統計參數估計得到，具體計算方法可以參見文獻[6]。式(19)的求解需要數值迭代算法，由文獻[14]知，基于常規Hess矩陣的牛頓或者擬牛頓算法，其求解結果一般較差，因而本文采用基于更新Hess矩陣的擬牛頓算法對上式進行求解[14]。

為避免l1-范數在零處不可導問題，對式(19)進行近似得

其中ξ為一個很小的正數。對式(20)進行求梯度得

其中η為迭代步長因子，式(22)迭代終止的條件為δ為一個很小的正數。

4 下視3維SAR成像算法流程

經過第3節的處理，就可以對下視3維目標進行SAR成像了。本文給出方法的處理流程如圖4所示。對于多通道收發系統而言，首先是等效相位中心處理，經過處理后多通道收發系統可以近似為是各陣元獨立收發。然后對每個陣元進行航向維處理，航向維的處理可以看作是傳統的正側視情況，因而可以采用RD算法進行實現，需要注意的是在航向維的頻域要進行相位補償，將一個周期內的等效相位中心校正到同一條直線上，從而得到一個排列緊密的線陣。對每個陣元成像處理后，接下來就是跨航向的處理，進行高度和跨航向維成像，這兩維的處理相當于是考慮彎曲校正的dechirp處理，由于跨航向維采用稀疏布陣方式，因而在這兩維結合稀疏信號處理方法進行成像，最終通過正則化問題求解得到下視3維SAR圖像。

圖4 下視3維SAR成像處理流程

5 陣列配置與RIC準則

根據稀疏信號處理理論可以知道，下視3維場景能否被正確重建，取決于兩個因素[6,7]：第一，場景的稀疏度，第二，觀測矩陣Φ是否滿足RIC準則。

在前面分析知，下視3維SAR場景是滿足稀疏特性的。下面主要分析觀測矩陣Φ的RIC特性。RIC準則最早由Candes和Tao提出，當觀測矩陣Φ滿足RIC條件時，則觀測矩陣的列可以近似看作為一個正交的系統。為詳細分析陣元分布對觀測矩陣RIC條件的影響，先從簡單的情況分析，就是不考慮彎曲校正的情況，這時Φ就是一個逆傅里葉變換(IDFT)矩陣，維數為N×M。當考慮徙動校正的情況時，要把觀測數據和原始場景數據分別寫成列向量的形式，這時Φ就為一個大矩陣，其維數為O(N×L)×O(M×L)，其中L為高度向采樣點數，兩者對于陣元分布對觀測矩陣的 RIC影響是相似的，而且當觀測矩陣維數較小時，更便于觀察仔細。

陣元分布如圖5所示，進行等效后虛擬陣元位置如圖6所示，不考慮彎曲情況下，此時得到觀測矩陣Φ的列自相關矩陣如圖7所示，觀測矩陣的每一列表示點目標的回波，觀測矩陣的列自相關矩陣表示兩點之間的相關性，從理論上講兩點相鄰越近則相關性越大，如兩個點在一個分辨單元之內，兩點相鄰越遠則相關性越小，如兩個點在一個分辨單元之外。從仿真中也可以看出矩陣Φ的列自相關矩陣很接近一單位陣，這說明觀測矩陣Φ滿足RIC特性[15]。

圖5 等效前陣元擺放位置

圖6 等效后虛擬陣元位置

圖7 觀測矩陣Φ的列自相關矩陣

由于不考慮徙動校正的情況，沿跨航向維會帶來幾何形變，同一高度單元的點聚焦之后分布在不同的高度單元。因而要進行徙動校正，這樣同一高度單元的點成像后仍在同一條直線上。這時對觀測矩陣進行RIC測試結果如圖8所示，從局部放大圖可以看出，矩陣Φ的列近似滿足正交性，說明矩陣Φ仍滿足RIC特性。圖8之所以沒有像圖7一樣，這是把左乘、右乘、點乘矩陣都寫成左乘矩陣形式所決定的。

6 仿真數據處理結果與分析

為驗證本文算法的有效性，下面進行點目標仿真驗證，具體仿真參數如表1所示。

本仿真中采用的陣元分布如圖5所示，場景中點目標分布位置如圖9所示。共9個點，分布在3個高度上，地平面上分布4個點，關于中心對稱成十字擺放，跨航向維兩點相距200 m，沿航向兩點相距400 m。高度為30 m的地方分布一個點，在航向和跨航向2維平面的中心位置處。高度為60 m同樣分布4個點，其在航向和跨航向2維平面的位置與地平面上的點分布相同。

圖8 觀測矩陣Φ的列自相關矩陣

表1 仿真參數設置

采用本文的處理方法，即圖4的處理流程，成像結果如圖10所示，成像后點目標位置與原始場景中點目標位置基本相同，從成像效果可以看出本文算法是有效的。為進一步驗證本文所提算法，對沿航向位置為零處進行切片處理，并對直接采用dechirp算法進行成像，與本文采用CS方法在相同條件下進行比較，其中圖11為直接采用dechirp算法成像結果，圖12為采用本文算法成像結果。圖13(a)為圖11中A1點的等高線圖，圖13(b)為圖12中A2點的等高線圖，從圖13(a)可以看出直接采用dechirp算法成像時，跨航向維旁瓣非常嚴重，圖中旁瓣的形狀是由于在時寬-帶寬積較小的情況下，對跨航向維旁瓣進行徙動校正時發生錯誤，導致旁瓣與主瓣落在不同的高度單元內[16]，而采用本文提出的算法，如圖13(b)所示，旁瓣得到有效的抑制，等高線圖很規則，主瓣旁瓣能夠明顯分開，說明成像效果良好，以上進一步驗證了本文提出算法的有效性。

7 結束語

本文研究了基于CS思想的下視3維SAR成像算法，解決了垂直航向維實孔徑天線成像對陣元間距的限制問題。文中首先給出了雙通道等效單通道處理需補償的相位項，同時介紹了多發多收模式下線陣天線的工作形式，在此基礎上本文對下視3維SAR成像處理算法進行了研究，并給出了整個成像處理過程的流程圖。文中還對 CS處理需滿足的兩個條件進行了詳細的分析，驗證了CS方法在下視3維SAR成像處理中應用的合理性。最后，通過仿真數據處理驗證了本文所提算法的有效性及正確性。

圖9 點目標分布位置

圖10 點目標成像結果

圖11 直接采用dechirp算法成像結果

圖12 本文成像算法

圖13 點目標等高線圖

[1]Klare J.A new airborne radar for 3D imagingsimulation study of ARTINO[C].European Conference on Synthetic Aperture Radar,Germany,May 16-18,2006：1-4.

[2]Fery O and Meier E.3-D Time-domain SAR imaging of a forest using airborne multibaseline data at L- and P-bands[J].IEEE Transactions on Geoscience and Remote Sensing,2011,49(10)：3660-3664.

[3]Potter L C,Ertin E,Parker J T,et al..Sparsity and compressed sensing in radar imaging[J].Proceedings of the IEEE,2010,98(6)：1006-1020.

[4]Zhang Y X,Sun J P,Tian J H,et al..Compressed sensing SAR imaging with real data[C].International Congress on Image and Signal Processing,Yantai,China,Oct.16-18,2010：2026-2029.

[5]Ji S H,Xue Y,and Carin L.Bayesian compressive sensing[J].IEEE Transactions on Signal Processing,2008,56(6)：2346-2356.

[6]Candes E J and Tao T.Near-optimal signal recovery from random projections：universal encoding strategies?[J].IEEE Transactions on Information Theory,2006,52(12)：5406-5425.

[7]Candes E J and Tao T.Decoding by linear programming[J].IEEE Transactions on Information Theory,2005,51(12)：4203-4215.

[8]Li Jun,Xing Meng-dao,and Wu Shun-jun.Application of compressed sensing in sparse aperture imaging of radar[C].Asian Pacific Conference on Synthetic Aperture Radar Proceedings,Xi,an,China,Oct.26-30,2009：651-655.

[9]王偉偉,廖桂生,吳孫勇.基于小波稀疏表示的壓縮感知SAR成像算法研究[J].電子與信息學報,2011,33(6)：1440-1446.Wang Wei-wei,Liao Gui-sheng,and Wu Sun-yong.A compressive sensing imaging approach based on wavelet sparse representation[J].JournalofElectronics＆Information Technology,2011,33(6)：1440-1446.

[10]Zhu X X and Bamler R.Tomographic SAR inversion byL1-norm regularization-the compressive sensing approach[J].IEEE Transactions on Geoscience and Remote Sensing,2010,48(10)：3839-3846.

[11]梁毅.調頻連續波SAR信號處理[D].[博士論文],西安電子科技大學,2009.Liang Yi.Signal processing of FMCW synthetic aperture radar data[D].[Ph.D.dissertation],Xidian University,2009.

[12]Chen S,Donoho D L,and Saunders M.Atomic decomposition by basis pursuit[J].SIAM Journal on Scientific Computing,1998,20(1)：33-61.

[13]Donoho D L.Compressed sensing[J].IEEE Transactions on Information Theory,2006,52(4)：1289-1306.

[14]Cetin M and Kard W C.Feature-enhanced synthetic aperture radar image formation based on nonquadratic regularization[J].IEEE Transactions on Image Processing,2001,10(4)：623-631.

[15]Alonso M T,López-Dekker P,and Mallorquí J J.A novel strategy for radar imaging based on compressive sensing[J].IEEE Transactions on Geoscience and Remote Sensing,2010,48(12)：4285-4295.

[16]Moreira A,Mittermayer J,and Scheiber R.Extend chirp scaling algorithm for air- and spaceborne SAR data processing in stripmap and ScanSAR imaging modes[J].IEEE Transactions on Geoscience and Remote Sensing,1996,34(5)：1123-1136.