天繪二號雙星InSAR成像與DSM生成技術(shù)

向建冰，呂孝雷，付希凱，薛飛揚，云 燁，葉 宇，何 可

1. 中國科學(xué)院空天信息創(chuàng)新研究院，北京 100194； 2. 中國科學(xué)院大學(xué)電子電氣與通信工程學(xué)院，北京 100049； 3. 中國科學(xué)院空天信息創(chuàng)新研究院空間信息處理與應(yīng)用系統(tǒng)技術(shù)重點實驗室，北京 100190； 4. 中國資源衛(wèi)星應(yīng)用中心，北京 100094

合成孔徑雷達(synthetic aperture radar,SAR)是利用微波獲取地面信息的主動式傳感器，它不受天氣條件限制，能夠穿透雨雪霧，全天候、全天時獲取數(shù)據(jù)，成為了高分辨率對地觀測和全球資源管理不可或缺的手段之一，已廣泛用于測繪、監(jiān)測和水文學(xué)等領(lǐng)域[1]。星載合成孔徑雷達干涉技術(shù)(InSAR)是對星載SAR圖像進行干涉處理的一種技術(shù)，它可以生成數(shù)字地表模型圖(DSM)、監(jiān)測地形形變和實現(xiàn)三維重建，是全球大范圍地形測繪最有效的觀測技術(shù)之一。目前世界上能夠用于InSAR處理的衛(wèi)星主要有加拿大的RADARSAT系列，意大利的COSMO-SkyMed以及德國的TerraSAR-X/TanDEM-X系統(tǒng)[2]。文獻[3—4]分別用ERS-1/2數(shù)據(jù)，TanDEM-X系統(tǒng)的數(shù)據(jù)利用InSAR技術(shù)實現(xiàn)了地形提取。而目前國產(chǎn)的SAR影像主要來源于高分三號，但高分三號重訪周期長且軌道誤差較大，時間去相干和大氣去相干嚴(yán)重[5-6]，大部分的高分?jǐn)?shù)據(jù)不滿足干涉測量條件，國內(nèi)對于InSAR數(shù)據(jù)需求迫切。天繪二號是國內(nèi)首個編隊飛行的X波段合成孔徑雷達衛(wèi)星，相比于傳統(tǒng)InSAR系統(tǒng)，編隊飛行能夠克服時間去相干和大氣去相干等相干源，獲得高精度的干涉對，對于提升我國遙感觀測能力具有重大意義[7]。

衛(wèi)星編隊是指衛(wèi)星在非常接近的軌道面上協(xié)同飛行，通過調(diào)整衛(wèi)星飛行軌道的偏心率、近地點角和升交點赤經(jīng)等，可以使雙星具有3種靈活工作模式[8]，分別為跟飛模式、雙站模式和雙站交替模式，其中雙站模式是其中一顆衛(wèi)星作為脈沖發(fā)射衛(wèi)星，兩顆衛(wèi)星同時接收回波信號，即一發(fā)雙收模式。與傳統(tǒng)重軌干涉相比，雙站模式具有0 s基線的獨特優(yōu)勢，能夠消去時間去相干源和大氣去相干源[9]，極大提高了有效干涉數(shù)據(jù)的質(zhì)量和數(shù)量。

在雙星SAR中，目標(biāo)的距離歷程與收發(fā)衛(wèi)星均有關(guān)，其距離歷程表現(xiàn)為雙根式，使得無法準(zhǔn)確解析駐定相位點，難用表達式表示回波信號的頻域形式，而經(jīng)典成像算法(如ω-k算法、CS算法、RD算法)的導(dǎo)出均基于回波信號的二維頻域表達式。近年來，很多學(xué)者對于雙站模式下的成像算法做了相關(guān)研究，文獻[10]提出了適用于任何構(gòu)型雙基SAR的時域成像方法，這種方法可以產(chǎn)生最佳的圖像質(zhì)量但同時也帶來了巨大的計算復(fù)雜性。與時域成像方法相比，頻域成像算法能極大地提高處理效率，其中求解點目標(biāo)的二維頻譜是算法推導(dǎo)中重要的一步。文獻[11—13]將平行構(gòu)型的雙基SAR斜距歷程等效成單基SAR的形式，進而可以運用傳統(tǒng)的單基SAR成像方法，例如Range-Doppler(RD)算法處理雙基數(shù)據(jù)。文獻[14]提出Loffeld 's bistatic formula (LBF)方法對收發(fā)相位歷程在各自駐相點處二階泰勒展開，得到任意構(gòu)形下雙基SAR點目標(biāo)的二維近似頻譜。文獻[15—16]基于LBF方法并利用二維逆變標(biāo)傅里葉變換(inverse scaled fast Fourier transform,ISFT)和chirp scaling (CS)算法聚焦雙基數(shù)據(jù)。此外，文獻[17]采用級數(shù)反演法(method of series reversion,MSR)來求解任意構(gòu)形雙基地SAR點目標(biāo)的二維頻譜。其中單站等效法可近似解析雙星SAR二維頻譜，同時后續(xù)處理可沿用傳統(tǒng)單基成像算法，具有高效便利的優(yōu)勢，適合實際工程中雙星SAR條帶模式成像。

本文針對天繪二號干涉成像處理及三維重建進行了研究和試驗，首先介紹了基于雙曲等效的雙星干涉成像擴展線性調(diào)頻變標(biāo)(extended chirp scaling,ECS)算法，并引入預(yù)濾波處理混合基線帶來的圖像間相干性降低和干涉相位誤差等問題。然后闡述了干涉處理方法，設(shè)計了用InSAR技術(shù)重建DSM的技術(shù)流程。最后用天繪二號原始數(shù)據(jù)進行了干涉成像試驗和干涉處理，并對成像結(jié)果和干涉處理結(jié)果進行了分析，驗證了成像算法的聚焦性和保相性，以及干涉處理和三維重建的能力。

1 雙星InSAR干涉成像

天繪二號SAR工作模式為條帶模式，這種模式可以覆蓋較大的幅寬，有利于完成大范圍測繪任務(wù)[9]。雙站模式下的天繪二號的軌道構(gòu)型與TanDEM-X/TerraSAR-X系統(tǒng)類似[7]，近似為一個螺旋式軌道，在一次成像時間內(nèi)，兩顆衛(wèi)星可看作平移不變構(gòu)型，即兩顆衛(wèi)星速度一致且基線保持不變。在雙星成像算法中，雙曲等效算法將單站等效算法擴展到了斜視情況并適用于平移不變構(gòu)型[13]，可以便利地求解平移不變構(gòu)型下雙基SAR點目標(biāo)的二維近似頻譜。因此筆者將雙曲等效法和ECS成像算法結(jié)合來處理雙星SAR數(shù)據(jù)。主輔星相對于目標(biāo)的斜距歷程可以表示為

(1)

式中，RT(t)和RR(t)分別表示主星和輔星相對目標(biāo)的斜距歷程；tc為場景中心時刻；t為方位向時間；Rm、Vm、φm分別表示雙曲等效后的斜距、速度和斜視角；Rx,c、vx和φx,c(x=T,R)分別表示主輔星相對于目標(biāo)點的斜距、速度和斜視角。

等效參數(shù)的計算如式(2)—式(4)所示的雙曲等效方法完全等效了雙星斜距歷程的常數(shù)、線性和二次項，同時也補償了部分三次及以上項[13]

Rm=(RT,c+RR,c)/2

(2)

(3)

(4)

將上述等效參數(shù)代入單站ECS成像算法中即可以實現(xiàn)輔星ECS聚焦成像，同時由于主星是自發(fā)自收，其成像方法可采用單星成像方法。

編隊飛行的主輔星之間不可避免地存在沿航向基線和垂直航向基線，稱之為混合基線?；旌匣€將引起主輔星回波數(shù)據(jù)間方位向多普勒譜偏移和距離波束譜偏移，導(dǎo)致主輔SAR圖像間相干性降低和干涉相位誤差等問題。

針對垂直航向基線引起的距離波束譜偏移導(dǎo)致的主輔圖像間相干性降低和干涉相位誤差問題，本文采用距離預(yù)濾波技術(shù)濾除不重疊的距離波束譜，即根據(jù)反射率譜移動量，截取公共的距離波數(shù)譜，再對具有公共距離波數(shù)譜的信號進行距離壓縮處理。垂直航向基線導(dǎo)致距離波數(shù)譜偏移的解釋如下：SAR天線接收的回波數(shù)據(jù)可以看作是地面反射率與入射信號的卷積，從頻域上看是地面反射率譜與信號譜的乘積，其表達式如下

(5)

(6)

式中，R(ω)為地面反射率的傅里葉變換；ω0為雷達中心頻率；θ1和θ2分別為兩顆衛(wèi)星的下視角；θ=(θ1+θ2)/2；β為地面坡度；W(ω)為SAR系統(tǒng)發(fā)射和接收的帶通濾波器，其帶寬為B。不同下視角之間的地面反射率譜的偏移量為

(7)

式中，Δθ為視角差；f0為雷達載頻。不重疊譜段(非相干譜段)相當(dāng)于重疊譜段(相干譜段)的噪聲，從而降低了相干性。因此必須切除不重疊譜段同時保留重合譜段。

針對沿航向基線引起的多普勒譜偏移導(dǎo)致的主輔圖像間相干性降低和干涉相位誤差問題，本文采用方位預(yù)濾波技術(shù)濾除不重疊的方位多普勒譜。具體實現(xiàn)方法如下：根據(jù)主輔衛(wèi)星回波數(shù)據(jù)的多普勒譜中心頻率和多普勒帶寬確定出重疊的多普勒譜，并濾除非重疊部分的多普勒譜。

圖1 方位預(yù)濾波多普勒譜段Fig.1 Diagram of azimuth pre-filtering Doppler spectrum

2 DSM生成方法

InSAR利用干涉相位信息，對同一區(qū)域的兩幅SAR單視復(fù)圖像進行干涉處理，生成DSM數(shù)據(jù)。InSAR處理要求兩幅SAR圖像具有較好的相干性，數(shù)字地表模型圖(DSM)為干涉處理結(jié)果之一[18-19]，它反映了圖像區(qū)域的高度信息。雙星模式下InSAR測量幾何如圖3所示，兩衛(wèi)星同時照射同一區(qū)域，去除了時間去相干性和大氣去相干等去相干源，產(chǎn)生的圖像對能直接用于干涉處理。

圖2 干涉成像算法流程Fig.2 The flowchart interferometry imaging algorithm

圖3 InSAR測量幾何Fig.3 The geometry of InSAR

圖3中P點為A星、B星照射的同一目標(biāo)點,P點的高程為h，S1、S2分別表示A星和B星在照射P點時的位置,R1、R2分別表示A星和B星到P點的斜距，B表示A星和B星之間的基線矢量，ΔB表示平行基線，α表示B與水平方向的夾角，θ為A星照射點的下視角，根據(jù)圖2所示的幾何關(guān)系可以把P點的高程表示為

h=H-R1·

(8)

由于干涉測量需要同一地區(qū)的兩幅相干性很高的SAR圖像,同時干涉相位不能直接從復(fù)圖像中得到，實現(xiàn)DSM生成需要通過復(fù)圖像配準(zhǔn)、干涉圖生成、干涉相位濾波、相位解纏、基線估計及DSM重建等步驟。因此本文設(shè)計了如圖4所示的干涉處理流程來完成DSM生成。

圖4 DSM生成流程Fig.4 The flowchart of DSM generation

復(fù)數(shù)圖像配準(zhǔn)是干涉處理第1步，包括兩幅單視復(fù)圖像之間的偏移值多項式計算和復(fù)影像重采樣[20-22]，兩幅影像配準(zhǔn)精度要優(yōu)于0.1像元，然后對配準(zhǔn)好的復(fù)圖像進行干涉處理，得到干涉條紋圖。配準(zhǔn)后干涉所形成的相位條紋密度一般較大，不利于后續(xù)相位濾波和解纏處理，通過基線估計計算條紋頻率來去除平地相位得到去平地后的干涉條紋圖。為了抑制干涉相位圖中的噪聲，提高相位解纏的效率和精度，本文采用自適應(yīng)濾波對相位進行降噪處理[23-24]。在得到解纏相位圖后，通過外部控制點精化基線重建DSM，其中基線是干涉處理中的一個重要參數(shù)，基線的長度和指向直接影響著圖像的相干性、測高靈敏度等，它也是影響DSM重建的精度的最重要誤差源之一[25-27]，在干涉處理中先通過軌道法和快速傅里葉變換法估計基線來去除平地相位[28-30]，后續(xù)通過引入外部控制點精化基線。

其中基于軌道法的基線估計是從基線的空間幾何關(guān)系出發(fā)，利用已知的衛(wèi)星軌道狀態(tài)矢量來獲得基線分量的一類方法，得到的基線只能是真實基線的一個近似解?；诳焖俑道锶~變換法的基線估計主要利用干涉圖中蘊含的信息如干涉條紋、干涉相位差等來估算基線分量。在軌道信息不完整或者精度較差的情況下，采用基于快速傅里葉變換的基線估計方式可以很好克服這個影響。相位與基線的關(guān)系滿足

(9)

代入衛(wèi)星高度H=R1cosθ得

(10)

則有

(11)

令

則有

(12)

式中，rmin和rmax對應(yīng)邊緣像元的斜距；kmin和kmax分別是干涉條紋相對應(yīng)的函數(shù)譜，可從如下解纏相位的傅里葉變換中得

基線精估計是基于解纏相位、對應(yīng)地面高度、基線之間的關(guān)系建立起來的方程，通過解這個方程可以求解出基線信息。對于地面點Ti的干涉相位觀測值為

式中，φi為Ti的解纏相位；φc為常數(shù)相位；rmi和rsi分別為地面點Ti至主影像和輔影像的斜距。斜距關(guān)系為

式中，Bx、By、rmix和rmiy分別為基線距rmi的x、y分量，則有

(13)

式中，θi為地面點Ti的視角；H為雷達高度；hi為地面點高度；rg為地球半徑。由于參數(shù)Bx、By、H和φc的變化率很小且呈線性。當(dāng)已知至少7個地面控制點時，可以解出未知參數(shù)Bx0、Bx1、By0、By1、H0、H1和φc。地面外部控制點的精度對最終的結(jié)果有一定的影響，并且控制點的數(shù)量也會有一定的影響，通過合適地選擇控制點，能夠精確地算出基線信息，而且該基線結(jié)果還能消除由于初始基線不準(zhǔn)確引起的DSM相對偏移和傾斜。

在生成DSM時，由于雙星編隊下輔星成像幾何與主星自發(fā)自收的情況不同，需要結(jié)合干涉相位方程、斜距方程和多普勒方程求解三維定位，此處可采用雙基等效相位中心法或聯(lián)合主圖像信息的等效相位中心法[31]。

3 試驗結(jié)果與分析

3.1 成像試驗

本文所選試驗區(qū)為某試驗場山地區(qū)域，A星和B星工作模式為雙站模式，其中B星為主星，A星為輔星，兩星均用同一波位觀測。對原始回波解壓縮后用基于雙曲等效的雙星干涉ECS成像算法處理，圖5為該地區(qū)的成像結(jié)果。

圖5 跟飛模式下雙星成像幅度Fig.5 Amplitude figure of two satellites imaging in follow-flight mode

在A星成像幅度圖的紅色矩形框圈出的區(qū)域中找到一個角反射器，如圖6(a)所示，對它插值128倍，結(jié)果如圖6(b)所示，可見輔星聚焦質(zhì)量很好。

圖6 角反射器成像結(jié)果Fig.6 Imaging result of corner reflector

3.2 干涉處理與DSM生成

選取3.1節(jié)中的雙星模式下的復(fù)數(shù)SAR圖像對進行干涉處理，干涉處理技術(shù)流程如圖4所示。首先需要對復(fù)數(shù)SAR圖像對進行裁剪和配準(zhǔn)，配準(zhǔn)后距離向配準(zhǔn)精度為0.051像元，方位向為0.02像元，滿足干涉處理要求。相干性是影響地形提取的重要因素，圖7是截取區(qū)域配準(zhǔn)后的相干性圖，圖中紅色區(qū)域為相干系很好的區(qū)域，綠色和藍色區(qū)域表示相干性較差區(qū)域，可以看出SAR圖像對整體相干性非常好，雙星模式下的雙星SAR圖像具有很好的相干性。

圖8顯示了干涉處理結(jié)果，其中圖8(a)是干涉條紋圖，對其進行去平地處理，得到如圖8(b)所示的地形高程變化引起的相位變化。數(shù)據(jù)所對應(yīng)的基線為556.628 1 m。然后，用最小費用流法(MCF)對濾波后的去平地后相位進行相纏，得到如圖8(c)所示的解纏相位圖。在求解三維定位后，采用精化后的基線生成數(shù)字高程圖，并轉(zhuǎn)換到地圖投影坐標(biāo)系下得到圖9中地理編碼后的RAWDSM。

圖7 相干性示意Fig.7 Coherence illustration

圖8 干涉條紋與解纏相位Fig.8 Interference fringes and unwrapped phase

圖9(a)為利用試驗場山地區(qū)域SAR圖像數(shù)據(jù)進行干涉處理和DSM重建后得到的該地區(qū)的數(shù)字高程圖，圖9(b)為同一地區(qū)的SRTM DSM，可以看到，本文處理生成的DEM與SRTM的DSM相比，整體地形輪廓與變化趨勢及細(xì)節(jié)幾乎一致，驗證了成像算法的保相性和數(shù)據(jù)的干涉測量以及三維重建能力。

圖9 數(shù)字地表模型圖(DSM)Fig.9 Digital surface model (DSM)

4 總 結(jié)

本文對雙星模式下的InSAR干涉成像和三維重建進行了研究和試驗。首先介紹了基于雙曲等效的雙星干涉ECS成像算法，然后闡述了DSM重建原理與技術(shù)流程。本文選取了天繪二號雙星模式下的某試驗場附近區(qū)域數(shù)據(jù)進行了成像試驗和DSM重建試驗，成像試驗結(jié)果驗證了干涉成像算法的良好聚焦性，配準(zhǔn)后的數(shù)據(jù)的相干性圖反映了雙星模式下的天繪二號數(shù)據(jù)的高相干性，最后干涉處理與DSM重建結(jié)果驗證了干涉成像算法的保相性和天繪二號的良好InSAR能力。