一種高分辨率多基線InSAR三維層析成像方法

柳祥樂，何向成，楊汝良

(1.湖北汽車工業(yè)學(xué)院電子信息系， 湖北十堰442002)(2.中國航天二院第二十五研究所， 北京100854;3.中國科學(xué)院電子學(xué)研究所， 北京100190)

0 引言

合成孔徑雷達(dá)(SAR)通過采用距離向脈沖壓縮、方位向合成孔徑技術(shù)來實(shí)現(xiàn)高分辨率的地形成像。SAR所成的像是二維像，是三維分布的散射體在二維平面上的一個(gè)投影。干涉合成孔徑雷達(dá)(InSAR)有兩副天線，一副天線發(fā)射雷達(dá)波，兩副天線同時(shí)接收回波，兩副天線接收的回波具有一定的相干性，經(jīng)干涉處理后便可得到能反應(yīng)地形高度的三維成像［1］。然而，InSAR所成的三維像只是地表的高度像，是一個(gè)面圖，并不是反映散射體空間分布的立體三維像。由于只有兩副天線，InSAR對高度向上分布的散射體分辨力很差，難以實(shí)現(xiàn)真正的立體三維成像。

有兩種方案可以獲得高度向的分辨能力，一種是利用多通道的極化信息，這就是極化干涉測量(PoIn-SAR)，另一種稱為多基線 InSAR(Multi-baseline In-SAR)，又稱為“多基線層析成像SAR”(Multi-baseline Tomography SAR)。

多基線InSAR是傳統(tǒng)InSAR的擴(kuò)展，傳統(tǒng)的In-SAR系統(tǒng)只有一條基線，而多基線InSAR在垂直于視線的方向依次增加多副SAR天線，或用重復(fù)航過的方式對同一地區(qū)成像，相當(dāng)于在垂直于視線的方向上合成了一個(gè)較大的孔徑，因而可以獲得高度向上的分辨力。當(dāng)用波長較長的電磁波對諸如冰、土壤、森林植被等半透明的媒質(zhì)進(jìn)行照射成像時(shí)，利用其高度向的分辨力，能夠得到任何一個(gè)斷層的影像，這一點(diǎn)類似于CT和核磁共振(NMR)的斷層攝影技術(shù)，因此，多基線InSAR此時(shí)也被稱為層析成像SAR(Tomography SAR)。多基線InSAR的示意圖可用圖1來所示，圖中共有K條航線或K副天線，它們沿著垂直于視線方向的圓弧排成一列。

圖1 多基線InSAR示意圖

早期的多基線InSAR研究多在實(shí)驗(yàn)室進(jìn)行，一般是三天線系統(tǒng)，雖然基線的數(shù)量很少，但在植被研究、干涉和層析成像方面也獲得了一定的應(yīng)用［2-3］。德國宇航局(DLR)首次成功地使用其機(jī)載E-SAR系統(tǒng)進(jìn)行了層析成像飛行實(shí)驗(yàn)［4］，對Oberpfaffenhofen地區(qū)進(jìn)行了真正意義上的層析成像，他們采用的是載機(jī)重復(fù)航過的方式，得到14個(gè)平行的航線(13條基線)，實(shí)現(xiàn)了高度維的聚焦成像。在星載情況下，用小衛(wèi)星編隊(duì)飛行的“分布式SAR”也能夠形成層析成像所需要的構(gòu)形。雖然提出的時(shí)間不長，但多基線層析成像在地質(zhì)學(xué)、林學(xué)、考古學(xué)及冰川研究中都展示出了巨大的潛力。

多基線InSAR三維成像處理的一般步驟是:先對觀測區(qū)進(jìn)行常規(guī)的二維SAR成像，得到同一區(qū)域的多幅高分辨率圖像，然后對多幅圖像進(jìn)行精確配準(zhǔn)，最后對配準(zhǔn)后的多幅圖像進(jìn)行高度維成像聚焦。二維SAR成像技術(shù)目前已經(jīng)比較成熟，多基線InSAR三維層析成像最主要的難點(diǎn)是高度維的成像問題，目前主要高度維成像算法是采用DFT聚焦技術(shù)，由于實(shí)際多基線InSAR系統(tǒng)中基線的數(shù)量(對應(yīng)于飛行航線數(shù))比較少，在滿足空間采樣定理的條件下，使高度維的合成孔徑長度較短，受限于瑞利界，高度維的分辨率較低，而靠增多基線條數(shù)來提高分辨率對多基線InSAR系統(tǒng)來說又非常不現(xiàn)實(shí)。為提高高度維的成像分辨率，本文將陣列信號處理技術(shù)引入到多基線InSAR高度維成像中來。仿真實(shí)驗(yàn)表明，這種成像方法可以突破瑞利界的限制，在保持基線條數(shù)較少的情況下大大提高高度維的分辨率，實(shí)現(xiàn)超分辨成像。

1 多基線InSAR三維層析成像之高度維成像的陣列信號模型

如上所述，基于DFT算法的多基線InSAR三維成像的基本過程是:首先，對各個(gè)天線所得的數(shù)據(jù)進(jìn)行常規(guī)的SAR成像處理，得到多幅二維SAR復(fù)圖像，對這些復(fù)圖像使用一定的方法進(jìn)行高精度的配準(zhǔn)，使各個(gè)同名點(diǎn)對齊;然后，對配準(zhǔn)后的復(fù)圖像沿高度維逐點(diǎn)進(jìn)行傅里葉變換，將數(shù)據(jù)變換到空頻域即可得到高度維上的影像［4］。三維成像處理的步驟如圖2所示。

多基線InSAR的方位向和距離向的成像與常規(guī)的二維SAR成像沒有區(qū)別，其特殊之處在于高度維聚焦成像。如前所述，通常受限于實(shí)際條件，多基線In-SAR的基線條數(shù)非常少，而基線長度受到空間奈奎斯特采樣定理約束，不能太長，于是在高度維總的合成孔徑長度受限，高度維分辨率較低。從實(shí)際情況看，基線條數(shù)通常少于20，采樣數(shù)據(jù)序列的長度很短，而與之相比，距離向和方位聚焦成像時(shí)，滿足采樣定理約束的數(shù)據(jù)點(diǎn)常常能達(dá)到數(shù)千個(gè)。在有限的數(shù)據(jù)集下獲得更高的分辨率，可以使用陣列信號處理技術(shù)，為此首先需要建立一個(gè)基于陣列信號處理的多基線InSAR高度維成像模型。為集中研究高度維的聚焦問題，本文假定SAR系統(tǒng)的二維點(diǎn)擴(kuò)展函數(shù)是理想的δ函數(shù)(實(shí)際上是sinc函數(shù))，在方位向、距離向已經(jīng)實(shí)現(xiàn)了理想的聚焦。

圖2 多基線InSAR三維成像處理原理框圖

圖1中，r表示斜距向，即視線方向，n表示垂直于視線的方向。通常各個(gè)SAR天線都沿著垂直于視線方向的一個(gè)圓弧上排列，但由于各條天線之間的間距都遠(yuǎn)小于它們到成像目標(biāo)區(qū)的距離，因此可以近似地認(rèn)為K副天線都排列在一條直線上，形成一個(gè)天線陣，如圖3所示。

圖3 簡化的多基線InSAR成像幾何(距離/高度維截面)

當(dāng)天線尺寸和基線長度遠(yuǎn)遠(yuǎn)小于天線相位中心到成像區(qū)域的距離時(shí)，多基線InSAR的陣列信號模型可以表達(dá)為

這是一個(gè)多視模型，變量n表示第n視，共有N視。M是距離/方位平面上某點(diǎn)對應(yīng)的高度線上的散射點(diǎn)的個(gè)數(shù)，本文中假定散射點(diǎn)的個(gè)數(shù)是已知的。對確定的n，y(n)和v(n)都是K維的復(fù)矢量，y(n)表示接收數(shù)據(jù)矢量，v(n)表示加性噪聲。sm(n)是高度向上第m個(gè)散射點(diǎn)的復(fù)反射系數(shù)，可以認(rèn)為它是一個(gè)窄帶信源。am是第m個(gè)源的方向矢量，在均勻陣列的情況下

式中:上標(biāo)T表示轉(zhuǎn)置;φm是對應(yīng)于全長基線(即在陣列中最外側(cè)的兩個(gè)天線間的基線長度)的干涉相位，根據(jù)InSAR的模型，它可以表示為

式中:θm是參考天線(如第1個(gè))對第 m個(gè)源的入射角;α是天線傾角;B是全基線長度;λ是雷達(dá)波長。

am也可以表示成另外一種形式，設(shè)以第1個(gè)天線為參考天線，對第m個(gè)目標(biāo)第2個(gè)天線與第1個(gè)天線間的干涉相位為

式中:d是各個(gè)天線間的間距，即基本基線的長度。在均勻陣列的情況下，各個(gè)基線對應(yīng)的干涉相位為

φmk=(k-1)ωm，(k=2，3，…，K)

于是方向矢量

am=〔1，ejωm，ej2ωm，…，ej(K-1)ωm〕T

于是模型(1)可以表示成簡潔的矩陣形式

式中:y(n)=［y1(n)，y2(n)，…，yK(n)］T(n=1，2，…，N)是一個(gè) K×N的接收數(shù)據(jù)矩陣;s(n)=［s1(n)，s2(n)，…，sM(n)］T(n=1，2，…，N)是 K×N 的信源矩陣;v(n)=［v1(n)，v2(n)，…，vk(n)］T(n=1，2，…，N)是 K×N的噪聲矩陣;A(ω)是方向矩陣

當(dāng)陣列均勻時(shí)它是一個(gè)Vandermonde矩陣。

上面提到的變量n表示的是第n視，共有N個(gè)視數(shù)據(jù)，這里的視數(shù)是指“N個(gè)單視”圖像數(shù)據(jù)，相當(dāng)于陣列信號處理中的快拍數(shù)。多基線InSAR超分辨率高度維成像模型要求在常規(guī)的SAR二維成像過程中存儲(chǔ)N幅單視復(fù)圖像。這可以在距離向壓縮完成后，方位向處理之前加上預(yù)置濾波器，然后再進(jìn)行各個(gè)單視的成像處理，同時(shí)存儲(chǔ)多個(gè)單視復(fù)圖像，一種典型的處理方法如圖4所示［5］。

圖4 獲得多個(gè)單視復(fù)圖像處理

2 高分辨率多基線InSAR高度維成像算法

建立了陣列信號處理模型后，就可以使用相關(guān)的信號處理方法進(jìn)行高度維成像了。二維SAR成像時(shí)，方位、距離坐標(biāo)相同而高度不同的各個(gè)點(diǎn)是層疊在一起的，在高度維上沒有分辨能力。而在多基線InSAR里面，多個(gè)SAR天線在空間中形成了一個(gè)陣列，高度不同的各個(gè)散射點(diǎn)回波到達(dá)陣列的角度互不相同，通過這些不同的角度可以將層疊在一起的各個(gè)散射點(diǎn)分辨出來，從而形成高度維的分辨能力，這就是陣列信號處理實(shí)現(xiàn)多基線InSAR高度維成像的基本原理。多基線InSAR高度維成像分辨率取決于對散射點(diǎn)回波波達(dá)角(DOA)差異的分辨率能力，采用現(xiàn)代信號處理技術(shù)，對DOA的估計(jì)能夠達(dá)到極高的精度，因而多基線InSAR可以得到極高的高度維分辨率。

DOA估計(jì)有很多種方法，本文中著重討論多重信號分類(MUSIC)算法，將它應(yīng)用于多基線InSAR高度維成像中，并將這種算法與業(yè)界已有的傳統(tǒng)成像算法進(jìn)行比較。

2.1 MUSIC算法高度維成像原理

現(xiàn)代信號處理方法的一個(gè)核心思想是充分利用信號信息的結(jié)構(gòu)特征，而自相關(guān)矩陣或自協(xié)方差矩陣往往含有豐富的信息特征。對于上節(jié)所建的多基線In-SAR高度維成像模型，觀測數(shù)據(jù)的自相關(guān)矩陣為Ry=E{y(n)yH(n)}，上標(biāo)H表示共軛轉(zhuǎn)置。當(dāng)航線數(shù)為K時(shí)，它是一個(gè)K×K維的矩陣?？蓪ζ溥M(jìn)行特征值分解，有

式中:λi(i=1，2，…，K)為特征值，并且 λ1＞λ2＞…＞λK;U是特征向量矩陣。在加性高斯白噪聲的條件下，有M個(gè)特征值大于某個(gè)閾值，稱為主特征值，或信號特征值;剩下的K-M個(gè)小于該閾值，稱為次特征值，或噪聲特征值。在很多情況下，主特征值和次特征值差異明顯，主特征值要明顯大于次特征值，并且在大信噪比時(shí)，主特征值的個(gè)數(shù)一般與信源(散射點(diǎn))的個(gè)數(shù)相等。將特征向量矩陣U按相應(yīng)的特征值大小分塊，有

S是矩陣U的前M個(gè)主特征值對應(yīng)的特征向量組成的矩陣;G是后K-M個(gè)次特征值對應(yīng)的特征向量組成的矩陣?？梢宰C明［6］

將式(3)代入到式(8)得

即，在M個(gè)散射點(diǎn)來波方向滿足aH(ω)G=0，而在其他方向aH(ω)G≠0。

將式(9)改寫為標(biāo)量形式，定義函數(shù)

作為相應(yīng)點(diǎn)目標(biāo)的高度維像。

由式(9)可知，在有目標(biāo)存在的角度方向上，P(ω)出現(xiàn)峰值，而在沒有目標(biāo)的位置上，P(ω)迅速衰減，形成了高度維的點(diǎn)目標(biāo)像。讓a(ω)中的ω以一定的角度間隔進(jìn)行搜索，就可以得到全部的目標(biāo)高度像。

值得注意的是，對成像而言除了分辨率以外，還要關(guān)心旁瓣問題。根據(jù)式(10)MUSIC法的P(ω)在目標(biāo)存在的方向上出現(xiàn)峰值，在偏離目標(biāo)方向時(shí)P(ω)會(huì)迅速衰減，但最后并不衰減到零，而是取一個(gè)有限值。對成像而言，此有限值是一種平均化的背景亮度，它會(huì)使圖像的對比度下降，成像時(shí)可以減去這個(gè)值。

2.2 高度維成像處理步驟

將高度維高分辨率成像處理步驟總結(jié)如下:

步驟1 常規(guī)的二維SAR成像，可以采用RD，CS，ECS等算法，要求校正距離徙動(dòng)。成像后所有的圖像都要精確配準(zhǔn)，配準(zhǔn)可以采用相干系數(shù)法、平均波動(dòng)函數(shù)法、最大干信比法、最小殘余點(diǎn)數(shù)法等［7］。值得注意的是成像處理器必須存儲(chǔ)每一個(gè)單視的復(fù)圖像，而不是存儲(chǔ)多視相加后的多視復(fù)圖像，這一點(diǎn)容易從頻域多視處理器上實(shí)現(xiàn)。

步驟2 精配準(zhǔn)后的K幅N個(gè)單視SAR復(fù)圖像上的每一點(diǎn)，都形成一個(gè)K×N維的觀測數(shù)據(jù)矩陣y(n)。

步驟3 求y(n)自相關(guān)矩陣Ry=E{y(n)yH(n)}。

步驟4 根據(jù)式(10)求一個(gè)方位、高度坐標(biāo)位置的高度維像，可以通過基本干涉相位角變化搜索來實(shí)現(xiàn)。成像過程中要減最小值并做歸一化處理。

3 數(shù)字仿真實(shí)驗(yàn)

為了驗(yàn)證陣列信號處理進(jìn)行高度維成像的性能，本文進(jìn)行了數(shù)字仿真實(shí)驗(yàn)。實(shí)驗(yàn)中假定天線等間距排列，形成均勻線陣(ULA)，若不均勻，需要將數(shù)據(jù)插值到均勻柵格再處理。數(shù)據(jù)的模擬是根據(jù)式(2)的模型進(jìn)行的，這里設(shè)噪聲為零均值的、獨(dú)立同分布(IID)的加性高斯白噪聲，而且各個(gè)視之間也是相互獨(dú)立的，N視的噪聲向量可以通過隨機(jī)數(shù)發(fā)生器運(yùn)行N次來求得。目標(biāo)信源s(n)可以是各種信源，如高斯信源。在不同的視，它們的值是不同的，為了減小計(jì)算量，本文假定信源與視無關(guān)，它只是一個(gè)確定性的矢量，即假設(shè)s(n)=［τ1，τ2，…，τM］T，τi(i=1，2，…，M)是各個(gè)散射點(diǎn)的散射強(qiáng)度。這一假設(shè)并不影響算法的性質(zhì)，卻可以減小仿真的復(fù)雜度。通過上述方法不難得到模擬的觀測數(shù)據(jù)矩陣y(n)，它是一個(gè)K×N的矩陣。

MUSIC算法都需要用到觀測數(shù)據(jù)的自相關(guān)矩陣Ry=E{y(n)yH(n)}，實(shí)際數(shù)據(jù)是有限的數(shù)據(jù)，因此只能對Ry進(jìn)行估計(jì)，本文用式(11)估計(jì)自相關(guān)矩陣

圖5所示的是成像區(qū)一個(gè)透視切面，即一個(gè)方位/高度截面，水平方向是方位向，垂直方向是高度向。沿方位向每隔一個(gè)方位分辨單元排列一個(gè)散射點(diǎn)，共10個(gè)散射點(diǎn)，各點(diǎn)在高度向呈階梯排布，高度間隔為5 m，10個(gè)點(diǎn)的高度值為0 m～45 m，如圖5b)所示。實(shí)驗(yàn)所用的數(shù)據(jù)為:天線K=8根，N=10個(gè)單視，白噪聲的方差為=1，信源功率(i=1，2，…，10)，信噪比/=10 dB。

圖5 三維層析成像區(qū)中的一個(gè)透視切面

圖6是使用傳統(tǒng)的傅里葉變換法和使用MUSIC算法的成像結(jié)果。在圖中，散射點(diǎn)高度換算成了干涉相位角的大小。

圖6 十個(gè)散射點(diǎn)的成像結(jié)果

為細(xì)致觀察成像結(jié)果，取出方位向第6列像，如圖7所示。

圖7 傳統(tǒng)高度維成像算法與MUSIC成像算法結(jié)果比較(K=8，N=10，信噪比 SNR=10 dB)

從圖6和圖7可以看出，MUSIC算法分辨率遠(yuǎn)高于傳統(tǒng)的DFT方法，而且MUSIC法的峰值旁瓣很低，可低于-30 dB，而傳統(tǒng)的DFT的方法第一旁瓣約為-13.2 dB(未加窗)。

再考查分辨率情況，以圖7b)中的3 dB寬度作為成像算法的干涉相位角分辨率，使用插值法可以求出兩種方法的基本干涉相位角的分辨率分別為:傳統(tǒng)DFT法為 φb≈0.715 4 rad=41°，MUSIC 法為 φm≈0.058 9 rad=3.4°，相對于傳統(tǒng)DFT法，提高約12倍，分辨率的改善十分顯著。

可算得傳統(tǒng)DFT法分辨率約為4.55 m，而MUSIC法約為0.38 m，為超分辨，遠(yuǎn)優(yōu)于DFT法。

最后再考查積分旁瓣比情況，積分旁瓣比可以通過對圖7b)一維沖擊響應(yīng)功率(幅度平方)進(jìn)行積分得到，定義為［8］

式中:Pmain為主瓣功率;Ptotal為總功率，分子為旁瓣的總能量。

將主要成像性能參數(shù)總結(jié)如表1所示，表中Pmain計(jì)算時(shí)把兩個(gè)相鄰零點(diǎn)之間的部分作為主瓣(這時(shí)算出來的ISLR比用3 dB寬度作為主瓣時(shí)的ISLR典型值要低)。

表1 主要成像性能參數(shù)對比(基線數(shù) K=8，視 N=10，SNR=10 dB)

4 結(jié)束語

多基線InSAR三維層析成像的關(guān)鍵是高度維成像。由于實(shí)際條件的限制，多基線InSAR系統(tǒng)的基線數(shù)目通常很少，數(shù)據(jù)序列的長度很短，采用常規(guī)的對觀測數(shù)據(jù)直接進(jìn)行傅里葉變換的方法難以達(dá)到較高的高度維分辨率。本文建立了一個(gè)基于陣列信號處理的成像模型，在此模型的基礎(chǔ)上，使用陣列信號中的DOA估計(jì)算法進(jìn)行多基線InSAR高度維成像，高度維分辨率大為提高，同時(shí)成像旁瓣也很低。

除本文提到的算法外，陣列信號處理中還有許多優(yōu)良的DOA估計(jì)算法，本文沒有一一提及，使用這類方法，可以在基線較少的情況下達(dá)到很高的高度維分辨率，或者在一定的分辨率要求下只使用很少的基線數(shù)。減少一根基線意味著載機(jī)少飛行一次，或少一顆小衛(wèi)星，使得系統(tǒng)簡化，成本降低，因此本文提出的成像方法對經(jīng)濟(jì)地解決多基線InSAR系統(tǒng)平臺的實(shí)際困難具有一定的指導(dǎo)意義。

［1］ 王 超，張 紅，劉 智.星載合成孔徑雷達(dá)干涉測量［M］.北京:科學(xué)出版社，2002.Wang Chao，Zhang Hong，Liu Zhi.Spaceborne synthetic aperture radar interferometry［M］.Beijing:Science Press，2002.

［2］ Treuhaft R，Moghaddam M，van Zyl J，et al.Estimating vegetation and surface topographic parameters from multibase radar interferometry［C］//Geoscience and Remote Sensing Symposium.Lincoln，NE:IEEE Press，1996:978-980.

［3］ Pasquali P，Prati C，Rocca F，et al.A 3-D SAR experiment with EMSL data［C］//Geoscience and Remote Sensing Symposium.Firenze:IEEE Press，1995:784-786.

［4］ Reigberand A，Moreira A.First demonstration of airborne SAR tomography using multibaseline L-band data［J］.IEEE Transactions on Geoscience and Remote Sensing，2000，38(5):2142-2152.

［5］ 張澄波.綜合孔徑雷達(dá)——原理、系統(tǒng)分析與應(yīng)用［M］.北京:科學(xué)出版社，1989.Zhang Chengbo.Synthetic aperture radar:principle，systems analysis and application［M］.Beijing:Science Press，1989.

［6］ 張賢達(dá).現(xiàn)代信號處理［M］.2版.北京:清華大學(xué)出版社，2002.Zhang Xianda.Modern signal processing［M］.2nd ed.Beijing:Tsinghua University Press，2002.

［7］ 柳祥樂，左艷軍，楊汝良.基于最少殘余點(diǎn)數(shù)的InSAR復(fù)圖像配準(zhǔn)［J］. 現(xiàn)代雷達(dá)，2007，29(3):37-41.Liu Xiangle，Zuo Yanjun，Yang Ruliang.InSAR SLC based on minimal residues［J］.Modern Radar，2007，29(3):37-41.

［8］ 美 Ian G.Cumming，F(xiàn)rank H.Wong.合成孔徑雷達(dá)成像——算法與實(shí)現(xiàn)［M］.洪 文，胡東輝，等譯.北京:電子工業(yè)出版社，2007.Ian G Cumming，F(xiàn)rank H.Wong.Digital processing of synthetic aperture radar data:algorithms and implementation［M］.Hong Wen，Hu Donghui，et al，tran.Beijing:Publishing House of Electronics Industry，2007.