孫希龍 余安喜 董 臻 梁甸農(nóng)

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

1 引言

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

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

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

2 差分SAR層析基本原理

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

圖1 SAR成像示意圖

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

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

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

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

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

圖2 差分SAR層析數(shù)據(jù)獲取模式

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

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

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

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

其中

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

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

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

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

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

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

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

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

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

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

4 實驗結(jié)果

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

表1 ENVISAT-ASAR系統(tǒng)參數(shù)

圖3 球場的SAR圖像

圖4 觀測數(shù)據(jù)分布情況

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

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

5 結(jié)論

圖5 高程分辨能力仿真圖

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

圖7 球場的3維模型

圖8 實測數(shù)據(jù)處理結(jié)果

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

[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層析技術(shù)的研究與應(yīng)用[J]. 遙感技術(shù)與應(yīng)用, 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]. 電子與信息學(xué)報，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.