一種差分SAR層析高分辨成像方法

孫希龍 余安喜 董 臻 梁甸農

(國防科技大學電子科學與工程學院 長沙 410073)

1 引言

差分 SAR層析成像技術(Differential SAR Tomography)由Lombardini[1]于 2005年提出。該技術利用從不同空間、時間位置對同一目標的多次觀測，在保持傳統SAR方位向合成孔徑的同時，在斜距垂向(normal-to-slant-range direction, nsr direction)-斜距形變速率向 2維平面內形成合成孔徑，使SAR具備了方位向-距離向-高程向-形變速率向的分辨能力，所以差分SAR層析成像也被稱為4維SAR成像[2]。

目前，差分SAR層析的觀測數據通常由SAR系統的重軌任務獲得。由于當前在軌運行的SAR系統多為單雷達(如 Envisat)或雙雷達系統(如Tandem-X)，并且這些系統在設計之初并未考慮差分SAR層析的應用需求，所以對同一目標的重軌觀測次數較少，并且軌道間的空間和時間間隔非均勻。此時，差分SAR層析的觀測數據為空間-時間2維全孔徑數據的稀疏采樣，基于傅里葉變換的成像方法無法得到令人滿意的聚焦效果，旁瓣高且分辨率低。基于空間譜估計的差分SAR層析成像方法[1]雖然具有較好的高度向-形變速率向分辨能力，但是卻需要使用獨立同分布的多視數據估計觀測的協方差矩陣，從而降低了方位向-距離向分辨率。城區建筑的4維成像是差分SAR層析技術的重要應用領域之一[3,4]。此種應用中，強散射元的尺寸小、密度高、散射特性不滿足各態歷經性[5,6]，需要獨立同分布多視數據的SAR層析成像方法并不適用。基于奇異值分解(Singular Value Decomposition, SVD)的成像方法[7]和基于Backus-Gilbert的成像方法[8]雖然能夠保持方位向-距離向分辨率，但是它們的高度向-形變速率向分辨率依然較低。利用稀疏的觀測數據，如何在保持方位向-距離向分辨率的同時實現高度向-形變速率向的高分辨成像是差分 SAR層析研究中的熱點和難點問題。

對于城區建筑目標，同一 SAR方位向-距離向分辨單元內通常僅存在少數個強散射元[5]，差分SAR層析的觀測數據可以模型化為復正弦信號和的形式。鑒于此，本文將 RELAX 算法[9]引入到差分SAR層析成像中，基于確定性點目標模型，提出了一種適用于城區建筑目標的差分 SAR層析成像方法。該方法對多軌單視數據進行處理，在保持方位向-距離向分辨率的同時實現了高度向-形變速率向的高分辨成像。在仿真數據和Envisat-ASAR實測數據處理實驗中，本文方法的性能得到了驗證。

2 差分SAR層析基本原理

為了實現對目標的4維成像，差分SAR層析除了需保持傳統SAR在方位向的合成孔徑外，還需要在斜距垂向(nsr direction)和斜距形變速率向上形成合成孔徑，這主要通過SAR系統的重軌飛行任務來實現。假設具有P個雷達的SAR系統從不同的空間和時間位置對同一目標進行了L次觀測，圖1給出了一次航過中的SAR成像示意圖，為了方便，假設編隊衛星沿nsr向分布。經過方位向-距離向2維壓縮后可得P×L幅SAR單視復圖像。首先選取其中一幅圖像作為主圖像，在通過圖像配準和去斜操作之后，圖像序列中某一方位向-距離向分辨單元的復序列可表示為[1,2]

圖1 SAR成像示意圖

其中s為nsr向高度，[ -smax,smax]為目標的nsr向跨度，v為斜距向形變速率， [-vmax,vmax]為形變速率的跨度，a(s,v)為目標的雷達散射特性函數，λ為雷達波長，r為主圖像中心斜距，bp,l為垂直基線，tp,l為第(p,l)幅單視復圖像相對于主圖像的時間間隔。令ξp,l= 2bp,l/λr,ηp,l=- 2tp,l/λ，則式(1)可記為如下形式：

由式(2)知差分SAR層析的觀測數據為同一分辨單元內目標的雷達散射特性函數在高度和形變速率向的2維聯合譜。

圖2給出了3種差分SAR層析數據獲取模式。圖2(a)為理想的全孔徑數據獲取模式，觀測在垂直基線-時間平面內均勻分布并且滿足奈奎斯特采樣定理。在這種理想的全孔徑模式下，對觀測數據進行2維傅里葉變換就可得到較好的聚焦效果，分辨率由垂直基線總跨度和時間總跨度決定。但是當前在軌運行的SAR系統多為雙雷達(如Tandem-X)或單雷達系統(如Envisat)，這些系統在設計之初并未考慮差分SAR層析的應用需求，對同一目標的重軌觀測次數少，并且軌道間的空間和時間間隔非均勻，所以理想的全孔徑數據在實際中是無法獲得的。圖2(b)和圖2(c)為兩種實際中的數據獲取模式，圖2(b)對應于雙雷達SAR系統，圖2(c)對應于單雷達SAR系統。不難看出，實際中獲取的數據為理想全孔徑數據的稀疏采樣。此時對觀測數據進行2維傅里葉變換將無法得到理想的聚焦效果，必須尋求其他的成像方法。

3 基于RELAX的差分SAR層析成像

利用差分SAR層析對城區建筑目標進行4維成像是當前的研究熱點。城區場景電磁散射主要包含如下幾種主要成分(如圖1所示)：(1)傾斜面造成的單次反射，如傾斜屋角或屋頂等；(2)垂直墻體和地面構成的2面角造成的2次反射；(3)復雜結構造成的3次反射；(4)粗糙表面造成的漫散射。前3種散射的強度較強，散射元尺寸小，通常可以用點目標模型表示[3]。由粗糙表面引起的漫散射的散射強度較弱，通常被作為雜波。同一方位向-距離向分辨單元內包含的由前3種散射造成的強散射元個數由建筑的結構復雜度決定。文獻[5,10]的研究成果表明，對于高分辨率和中等分辨率的SAR系統而言，同一方位向-距離向分辨單元內的強散射元個數通常不大于4。

圖2 差分SAR層析數據獲取模式

假設在某一方位向-距離向分辨單元內共有K個強散射元，則式(2)可記為如下形式：

其中ep,l為對時間去相關、大氣相位和雜波等進行補償之后的殘余噪聲[1]。由式(3)可知，觀測數據g(p,l)為多個復正弦信號的疊加，差分SAR層析的高度向-形變速率向成像可等效為復正弦信號幅度參數a(sk,vk)和頻譜參數sk,vk的估計問題。

RELAX 算法是在CLEAN 算法基礎上發展起來的一種基于非線性最小方差(NLS)準則的高分辨譜估計算法。RELAX算法對噪聲不做任何限定性假設，具有較強的適應性和魯棒性。鑒于此，本文使用RELAX算法對各強散射元的未知參量a(sk,vk),sk,vk進行估計。令

將式(3)寫成矩陣形式為

其中

hk為H的列，⊙代表Hadamard積。通過最小化如下的NLS代價函數對各強散射元的參數進行估計：

求解式(12)是一個非線性多維最優化問題，需要進行全局極值的多維搜索，計算量相當大，并可能無法求出最優解。RELAX算法基于信號分離機理，用迭代的方法將極為復雜的多維搜索問題簡化進行求解。在給出基于RELAX的差分SAR層析成像迭代算法之前，先做下述準備，令

假設除第k個強散射元外的其他強散射元已被估計出，則根據式(11)和式(12)，第k個強散射元的參量估計為

下面給出基于RELAX的差分SAR層析的具體迭代步驟：

步驟 1 假設K=1。通過式(14)和式(15)求得

步驟 2 假設K=2。將求得的代入式(13)計算得到g2，通過式(14)和式(15)求得；將求得的代入式(13)計算得到g1，進而通過式(14)和式(15)重新求得。重復步驟2直至滿足內部收斂條件。具體的內部收斂條件將在后面給出。

步驟 3 假設K=3。將求得的代入式(13)計算得到g3，通過式(14)和式(15)求得；將求得的代入式(13)計算得到g1，進而通過式(14)和式(15)重新求得；將求得的代入式(13)計算得到g2，進而通過式(14)和式(15)重新求得。重復步驟3直至滿足內部收斂條件。

步驟 4 重復上述相同的步驟直至K等于實際的強散射元個數或者預先設定的強散射元個數。

上述每一步驟中，通常以鄰近兩次迭代中由式(10)所示的NLS代價函數的變化值是否小于某個門限εNLS作為內部收斂的判據。εNLS越小，則表示兩次迭代中得到的各參量估計值變化越小，迭代算法越收斂，估計精度越高。但是隨著εNLS的減小，上述迭代算法的計算負擔也會相應增加，所以需要在估計精度和計算負擔之間取得折中。通常εNLS取為10-3即可得到較好的效果，在確保了估計精度的基礎上又不會使迭代算法的計算負擔過重。上述迭代算法中的另一個問題是強散射元個數K的估計問題，具體估計方法可參見文獻[11, 12]。在本文后續的實驗中，使用最小化GAIC(Generalized Akaike Information Criterion)代價函數的強散射元個數估計方法。

上面給出的基于RELAX的差分SAR層析成像方法是一種參數化的譜估計方法，汲取了CLEAN算法的思想，采用信號分離機理，基于符合城區目標實際散射的點散射目標模型，及時剔去已估計出的散射元的信號成份，盡量消除各散射元之間的相互影響，利用反復迭代逼近給定的數學模型，從而實現各參量的精確估計。與傅里葉變換和SVD等非參數化方法相比，該方法避免了點擴展函數(Point Spread Function, PSF)間的相互影響，特別是由于稀疏采樣導致的高旁瓣問題，從而能夠得到較好的分辨能力。但是，這并不意味著該方法具有無限的分辨能力，當兩個散射元的距離過小時，信號的保持性和剔除性降低，散射元間的相互影響增大，估計性能下降。在后續的仿真實驗中將對本文方法和SVD方法的分辨能力進行比較。

4 實驗結果

為了對本文方法的性能進行驗證，下面將給出仿真數據和實測數據的處理結果。實測數據處理中使用了Envisat自2003年10月至2007年11月間獲得的20軌ASAR數據，目標選定為位于美國加利福尼亞州的安那罕天使球場(Angel Stadium of Anaheim，北緯 33°48'0.11''，西經 117°52'58.85'')。表1為Envisat-ASAR的系統參數。圖3為球場的SAR圖像。圖4給出了以2005年6月15日獲得的圖像為主圖像時觀測的分布情況。因為Envisat為單雷達 SAR系統，所以在該種情況下P= 2 0,L=1。

表1 ENVISAT-ASAR系統參數

圖3 球場的SAR圖像

圖4 觀測數據分布情況

(1)仿真數據處理結果 利用實測數據的系統參數和軌道參數進行仿真實驗。實驗中設定兩個強散射元，散射元的散射強度，加性噪聲e的每個分量都服從零均值復高斯分布 C N(0,σ2)且相互獨立，散射元的信噪比 S NR1= S NR2=SNR3= 1 0 dB，其中第i個強散射元的信噪比定義為。首先，設定兩個散射元的形變速率相同，nsr向間隔Δs在0 m至40 m之間變化。圖5給出了不同Δs情況下本文方法和SVD方法得到的散射元nsr向位置。SVD方法中的散射元位置是通過取聚焦結果的峰值點位置得到的。圖5中的虛線給出了兩個散射元的真實位置；十字形點給出了兩個散射元無法分辨時的結果；圓形點給出了可分辨情況下的兩個散射元位置的估計結果。觀察圖5可知：當Δs大于15 m時，本文方法就可得到兩個散射元nsr位置較精確的估計；當Δs大于24 m時，SVD方法才可將兩個散射元較精確的分辨開來。接下來，設定兩個散射元的 nsr向位置相同，形變速率差Δv在0 cm/Yr至 5 cm/Yr之間變化。圖 6給出了不同Δv情況下本文方法和SVD方法得到的散射元形變速率值。當Δv大于 0.6 cm/Yr時，本文方法就可得到兩個散射元形變速率較精確的估計結果；當Δv大于1 cm/Yr時，SVD方法才可將兩個散射元較精確的分辨開來。上面的實驗證明本文方法具有明顯優于SVD方法的分辨能力。在本文給出的實際采樣和典型噪聲水平情況下，與SVD方法相比，本文方法的分辨能力可提高約40%。

(2)實測數據處理結果 下面使用本文方法對安那罕天使球場的實測數據進行處理，處理中內部收斂條件εNLS為 1 0-3，使用 GAIC(Generalized Akaike Information Criterion)方法對強散射元個數進行估計。圖3中的星號標出了進行差分SAR層析處理的4個方位向-距離向分辨單元。為了對處理結果進行評估，圖7給出了球場的3維模型。該球場頂棚的最大高度為42 m，距離向遠端露天看臺的高度為24 m。圖8為處理結果，結果中nsr向高度和斜距向形變速率已經根據雷達入射角轉換到了垂直于水平面的高度和形變速率上。圖8(a)為第1個分辨單元的處理結果。該分辨單元位于球場之外，僅存在一個由地面復雜結構引起的強散射元。設定該強散射元的高度和形變速率為零，其他各分辨單元的處理結果都依據此進行校正。圖8(b)為第2個分辨單元的處理結果。(0 m, 0cm/Yr)處的強散射元由地面復雜結構引起，(28 m, 0 cm/Yr)處的強散射元由球場前墻體的復雜結構引起。兩個強散射元的形變速率都為零，這說明球場的前墻體相對于地面未發生相對形變。圖8(c)為第3個分辨單元的處理結果。(0 m, 0 cm/Yr)處的強散射元同樣由地面復雜結構引起。(43 m, -1.8cm/Yr)處的強散射元由球場的金屬頂棚引起，其高度與體育場頂棚的實際高度吻合的較好，其形變應該由頂棚的老化等因素造成。圖8(d)為第4個分辨單元的處理結果，該分辨單元內存在兩個強散射元，分別位于(4 m, 0.9cm/Yr)和(18 m, -0.2 cm/Yr)處。由于該分辨單元跨過體育場遠端看臺區，所以兩個強散射元應該由座椅的復雜結構引起，強散射元的形變應該由座椅的老化和損壞等因素造成。

5 結論

圖5 高程分辨能力仿真圖

圖6 形變速率分辨能力仿真圖

圖7 球場的3維模型

圖8 實測數據處理結果

在城區建筑目標的差分SAR層析成像中，同一方位向-距離向分辨單元內通常僅包含少數個可由點目標模型表示的強散射元。鑒于此，本文提出了基于RELAX算法的差分SAR層析成像方法。該方法對多軌單視數據進行處理，可在保持方位向-距離向分辨率的同時實現高度向-形變速率向的高分辨成像。在利用實測數據系統參數和軌道參數進行的仿真實驗中，對本文方法和 SVD方法的高度向-形變速率向分辨能力進行了分析。在利用 Envisat-ASAR實測數據對安那罕天使球場進行的差分SAR層析處理中，本文方法得到了較好的結果，獲得了美國安那罕天使球場的高度和形變信息。

[1]Lombardini F. Differential tomography: a new framework for SAR interferometry[J].IEEE Transactions on Geoscience and Remote Sensing, 2005, 43(1): 37-44.

[2]Gianfranco F, Diego R, and Serafino F. Four-dimensional SAR imaging for height estimation and monitoring of single and double scatterers[J].IEEE Transactions on Geoscience and Remote Sensing, 2009, 47(1): 224-237.

[3]Reale D, Gianfranco F, Pauciullo A,et al.. Advanced techniques and new high resolution SAR sensors for monitoring urban areas[C]. Processing of IEEE International Geoscience and Remote Sensing Symposium, Honolulu Hawaii, USA, 2010: 1800-1803.

[4]張紅, 江凱, 王超, 等.SAR層析技術的研究與應用[J]. 遙感技術與應用, 2010, 25(3): 282-287.Zhang H, J iang K, Wang C,et al.. The current status of SAR tomography[J].Remote Sensing Technology and Application,2010, 25(3): 282-287.

[5]Zhu X X and Bamler R. Very high resolution spaceborne SAR tomography in urban environment[J].IEEE Transactions on Geoscience and Remote Sensing, 2010,48(12): 4296-4308.

[6]Reale D, Fornaro G, and Pauciull A. Tomographic imaging and monitoring of buildings with very high resolution SAR data[J].IEEE Geoscience and Remote Sensing Letters, 2011,8(4): 661-665.

[7]Serafino F, Soldovieri F, and Lombardini F. Singular value decomposition applied to 4D SAR imaging[C]. Processing of IEEE International Geoscience and Remote Sensing Symposium, Seoul, Korea, 2005: 2701-2704.

[8]任笑真, 楊汝良. 一種基于逆問題的差分干涉SAR層析成像方法[J]. 電子與信息學報，2010, 32(3): 582-586.Ren X Z and Yang R L. An inverse problem based approach for differential SAR tomography imaging[J].Journal of Electronics&Information Technology, 2010, 32(3):582-586.

[9]Li J and Stoica P. Efficient mixed-spectrum estimation with application to target feature extraction[J].IEEE Transactions on Signal Processing, 1996, 44(2): 281-295.

[10]Budillon A, Evangelista A, and Schirinzi G. Threedimensional SAR focusing from multipass signals using compressive sampling[J].IEEE Transactions on Geoscience and Remote Sensing, 2011, 49(1): 488-499.

[11]Zhu X X, Adam N, and Brcic R. Space-borne high resolution SAR tomography: experiments in urban environment using TS-X data[C]. Processing of Joint Urban Remote Sensing Event, Shanghai, China, 2009: 1-8.

[12]Reale D, Pauciullo A, Fornaro G,et al.. A scatterers detection scheme in SAR tomography for reconstruction and monitoring of individual buildings[C]. Processing of Joint Urban Remote Sensing Event, Munich, Germany, 2011:249-252.