胡曉寧 汪丙南* 向茂生 王鐘斌
①(中國科學(xué)院空天信息創(chuàng)新研究院 北京 100190)
②(中國科學(xué)院微波成像技術(shù)國家級重點實驗室 北京 100190)
③(中國科學(xué)院大學(xué) 北京 100049)
數(shù)字高程模型(Digital Elevation Model,DEM)在地形測繪、水文、農(nóng)業(yè)等諸多領(lǐng)域均有著廣泛的應(yīng)用[1–3]。利用干涉合成孔徑雷達(dá)(Interferometric Synthetic Aperture Radar,InSAR)技術(shù)反演DEM具有全天時、全天候、高效率、高精度等優(yōu)勢,是目前常用的DEM獲取手段之一[4]。在地形起伏劇烈的區(qū)域,InSAR干涉條紋變得十分密集,特別是星載平臺。在如橫斷山脈等地形突變區(qū)域,甚至?xí)霈F(xiàn)干涉相位模糊現(xiàn)象。密集的干涉條紋不僅增加了相位展開的難度,并且會在解纏后的干涉相位中引入更多的數(shù)據(jù)處理誤差,進(jìn)而影響高程反演的精度。
通過引入外源DEM可以降低干涉條紋的密度,已有一些文獻(xiàn)對該方法進(jìn)行了研究[5–8]。考慮到中國已有覆蓋全國的地形數(shù)據(jù),以及已有的各國研究團(tuán)隊通過各種方法獲取的一些全球數(shù)字高程數(shù)據(jù),如利用InSAR技術(shù)獲取的SRTM(Shuttle Radar Topography Mission)DEM[9],利用光學(xué)成像的立體像對技術(shù)獲取的ASTER GDEM(The Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model)[10],以及通過對SRTM進(jìn)行再處理,并將其與精細(xì)化后的ASTER GDEM高程相融合而生成的NASADEM[11]。上述數(shù)字高程數(shù)據(jù)可通過互聯(lián)網(wǎng)下載獲取,為在InSAR處理中引入外源DEM提供了條件。
后向投影(BackProjection,BP)算法是一種完全時域的無誤差成像算法,具有聚焦效果好、相位保持精度高、運動誤差補(bǔ)償精度高且算法簡單、適用于各種工作模式以及任意運動軌跡的雷達(dá)等優(yōu)點,在合成孔徑雷達(dá)(Synthetic Aperture Radar,SAR)成像中得到了廣泛的應(yīng)用[12]。文獻(xiàn)[13]對基于BP算法的InSAR方法開展了相關(guān)研究,理論分析了基于BP算法和基于頻域成像算法的InSAR方法的異同,但沒有進(jìn)行實驗驗證。文獻(xiàn)[14]研究了毫米波的后向投影InSAR成像方法,并進(jìn)行了實驗驗證,但在實測InSAR數(shù)據(jù)處理實驗中僅獲取了干涉相位,未進(jìn)行后續(xù)的高程反演研究。后向投影成像算法也被用于獲取干涉處理所需的SAR圖像,以簡化時序干涉處理的流程[15]??傮w而言,針對基于BP算法的干涉處理目前研究較少,尤其是針對理論模型和實測InSAR數(shù)據(jù)的實驗還缺乏深入的研究[16–18]。
針對上述問題,為了降低干涉相位的條紋密度,本文將外源DEM與BP算法相結(jié)合,提出了一種基于DEM輔助后向投影模型的InSAR高程反演方法。該方法利用BP算法獲得聚焦的SAR圖像,并在成像過程中引入輔助DEM,然后建立基于BP成像算法的InSAR高程反演模型,實現(xiàn)高精度的DEM反演。外源DEM輔助以及BP算法的引入,不僅可以降低干涉條紋的密度,減少干涉相位的纏繞,多數(shù)情況下甚至可以避免相位纏繞,同時圖像對無需干涉配準(zhǔn),簡化了干涉處理的流程。仿真實驗與機(jī)載InSAR系統(tǒng)獲取的實測數(shù)據(jù)處理結(jié)果證明了算法的有效性。
傳統(tǒng)InSAR處理中多采用頻域算法成像,獲取的干涉相位中包含平地相位和高程相位。平地效應(yīng)會大大增加干涉條紋的密集程度,去平地操作雖然可以降低干涉條紋的密度,但在地形起伏較大的區(qū)域干涉條紋仍然十分密集,甚至出現(xiàn)相位模糊的問題,進(jìn)而影響高程反演。在本文方法中,外源DEM的輔助和BP算法的引入大大降低了干涉條紋的密度,并且多數(shù)情況下避免了相位纏繞。本節(jié)將從成像處理和干涉處理兩個環(huán)節(jié)對該方法進(jìn)行詳細(xì)介紹,并給出具體的流程。
BP成像算法的基本原理是對距離壓縮后的數(shù)據(jù)精確補(bǔ)償投影面上的采樣點到天線相位中心的回波延時相位,并逐點逐脈沖地進(jìn)行相干累積,從而得到聚焦圖像[19]。
主要包括以下步驟:
(1)根據(jù)成像區(qū)域確定成像空間,對成像空間進(jìn)行網(wǎng)格劃分,并將輔助DEM插值到劃分的網(wǎng)格中作為BP成像過程中的投影面。
(2)對原始回波數(shù)據(jù)進(jìn)行距離壓縮。
(3)計算投影面上的采樣點到天線相位中心的距離延時,進(jìn)行多普勒相位補(bǔ)償,對合成孔徑內(nèi)每個天線相位中心處相位補(bǔ)償后的數(shù)據(jù)進(jìn)行相干累加,以實現(xiàn)該采樣點的方位向聚焦。
(4)遍歷成像空間中的每一個采樣點,得到聚焦成像的結(jié)果。
BP算法不存在雷達(dá)和目標(biāo)的斜距近似誤差,成像信號模型十分精確,相位保持精度也很高,并且BP算法根據(jù)雷達(dá)軌跡信息進(jìn)行運動誤差補(bǔ)償,具有很高的補(bǔ)償精度,為SAR圖像在后續(xù)的干涉處理中獲取高精度的DEM提供了基礎(chǔ);另外BP算法的成像過程處于空間直角坐標(biāo)系下,易于引入外源輔助DEM,并且可以避免復(fù)雜的圖像拼接和地理編碼過程。這也是選擇BP算法進(jìn)行成像處理的原因。
由于BP算法的方位向聚焦步驟中補(bǔ)償?shù)亩嗥绽障辔豁棸司W(wǎng)格點到雷達(dá)平臺的距離歷史信息,而這正是傳統(tǒng)InSAR中干涉相位的來源,因此基于BP成像算法的干涉幾何模型有別于傳統(tǒng)InSAR幾何模型。考慮到實際情況下輔助DEM存在地形起伏,本文建立了基于BP算法曲面投影下的InSAR高程反演模型,下面進(jìn)行具體的分析。
成像過程中引入的輔助DEM為粗精度DEM,與目標(biāo)的真實高程之間存在誤差,因此目標(biāo)在投影曲面上的投影位置將偏離目標(biāo)的真實位置,并且由于BP成像過程中網(wǎng)格的離散化采樣,目標(biāo)的投影位置也會偏離成像網(wǎng)格的采樣點位置,因此目標(biāo)真實的散射中心與成像過程中投影的網(wǎng)格點分別到雷達(dá)天線的斜距存在偏差,進(jìn)而導(dǎo)致BP成像方位向聚焦過程中多普勒相位補(bǔ)償存在誤差。該誤差正是基于時域BP算法獲取的SAR圖像進(jìn)行干涉處理時干涉相位的來源,幾何模型如圖1所示。
圖1 基于BP成像算法的干涉模型Fig.1 The interferometric model based on BP imaging algorithm
圖1中z軸為高度向,y軸指向雷達(dá)平臺的飛行方向,x軸由右手坐標(biāo)系確定。A1和A2分別代表主、副天線相位中心的位置,下視角為φ,基線長度B,基線與水平方向之間的夾角為α。由于目標(biāo)P的真實高度與輔助DEM之間存在高度偏差,主副天線將P分別投影到Q1,Q2,而成像網(wǎng)格中對應(yīng)的離散采樣點為C。天線相位中心A1的高程為H。點C,Q1和Q2的高程值分別為投影曲面上的采樣點C與目標(biāo)點P的真實位置之間的高程偏差為Δh,θ為與垂直方向的夾角。
根據(jù)余弦定理可得
在傳統(tǒng)InSAR的干涉相位模型中,P和Q1兩點間的干涉相位差與其高程差有關(guān),根據(jù)高程相位的表達(dá)式[20],該干涉相位差如式(6)所示。
根據(jù)圖1中的幾何關(guān)系,Q1和C點之間的干涉相位差可以表示為引起的平地相位和引起的高程相位之和的形式,如式(7)所示。
分析式(8),后向投影InSAR高程反演模型得到的干涉相位由兩部分組成:Δh的高度差引起的干涉相位差對應(yīng)的斜距變化引起的干涉相位差Δφr。在實際情況中,Δφr遠(yuǎn)小于Δφh,可以忽略不計(以3.2節(jié)的實驗參數(shù)為例,根據(jù)式(8)計算高程誤差為10m時對應(yīng)的Δφr與Δφh分別為–0.0166rad和–1.7rad)。因此最終得到后向投影干涉模型下高程反演公式為
由式(9)可以看出,本文基于后向投影算法的InSAR高程反演模型獲得的干涉相位去除了成像參考面對應(yīng)的地形相位,僅包含輔助DEM的高程誤差對應(yīng)的高程相位,因此無需去平地即可獲得條紋密度大大降低的干涉相位圖,并且與基于頻域成像算法的傳統(tǒng)InSAR高程反演方法中得到的去平地后的干涉相位相比,其條紋更加稀疏。
下面分析本文方法中相位纏繞的發(fā)生條件。圖2中P1,P2為干涉相位發(fā)生2π纏繞的兩個點目標(biāo),h1,he1分別為P1的真實高程值和輔助DEM對應(yīng)點處的高程值,h2,he2分別為P2的真實高程值和輔助DEM對應(yīng)點處的高程值,Δh1和Δh2則分別為P1和P2輔助DEM對應(yīng)點處的高程誤差。根據(jù)式(9)計算基于后向投影模型得到的干涉相位變化2π,即出現(xiàn)相位纏繞現(xiàn)象時對應(yīng)的輔助DEM高程誤差變化值,用Δh2π表示,如式(10)所示。Δh2π對應(yīng)圖2中的Δh2?Δh1。根據(jù)高程模糊度的定義以及圖2中的幾何關(guān)系,推導(dǎo)本文基于輔助DEM的后向投影InSAR高程反演模型的高度模糊數(shù)計算公式,如式(11)所示。與傳統(tǒng)InSAR方法中的高度模糊數(shù)相比,本文方法的高度模糊數(shù)更大,即干涉條紋更加稀疏,且更不易發(fā)生干涉相位纏繞。
圖2 高程模糊度幾何示意圖Fig.2 Geometric description of the height of ambiguity
以3.2節(jié)的X波段機(jī)載SAR系統(tǒng)參數(shù)為例,根據(jù)式(10)計算得到Δh2π大于100m,即輔助DEM的高程誤差變化大于100m時,干涉相位才會纏繞,而對于目前可以獲取的外源DEM而言,其精度遠(yuǎn)高于該值,因此即使在地形起伏較大的區(qū)域也不會出現(xiàn)相位纏繞的現(xiàn)象,從而避免了相位解纏處理。在實際情況中,考慮到基線誤差會引入隨距離向變化的相位誤差,當(dāng)輔助DEM高程誤差與基線誤差導(dǎo)致的干涉相位小于2π時,基于后向投影算法的In-SAR模型得到的干涉相位不會出現(xiàn)相位纏繞,避免相位解纏過程。
基于DEM輔助后向投影模型的InSAR高程反演方法在簡化干涉處理流程的同時,減少了解纏過程中的相位損失,提高了DEM反演的精度。綜上所述,本文所提的高程反演方法具有精確的成像信號模型、運動誤差補(bǔ)償模型以及很高的相位保持能力。
本文給出了一種基于DEM輔助后向投影模型的InSAR高程反演方法,處理流程如圖3所示,并與傳統(tǒng)基于頻域算法的InSAR處理流程進(jìn)行對比。傳統(tǒng)InSAR采用頻域成像算法得到聚焦的SAR圖像后,需要進(jìn)行圖像配準(zhǔn)、去平地、干涉相位濾波、相位解纏等處理,然后根據(jù)高程反演模型得到DEM重建結(jié)果。而基于DEM輔助后向投影模型的InSAR高程反演方法簡化了數(shù)據(jù)處理流程,無需圖像配準(zhǔn)以及相位解纏操作,減少了數(shù)據(jù)處理過程中的相位損失,為獲取高精度DEM提供了保障。并且基于后向投影的InSAR模型位于地理空間坐標(biāo)系下,可以避免后續(xù)復(fù)雜的圖像拼接和三維定位等處理。具體步驟如下:
圖3 InSAR處理流程Fig.3 InSAR processing chain
(1)BP成像。獲取實驗場景的外源DEM數(shù)據(jù),并插值作為成像的參考面,利用時域的后向投影算法獲取聚焦的SAR圖像。
(2)圖像干涉。主副天線獲取的復(fù)圖像直接共軛相乘得到后向投影InSAR模型下的干涉相位。
(3)相位濾波。為了降低相位噪聲對高程反演精度的影響,對干涉相位進(jìn)行濾波處理。
(4)高程反演。根據(jù)式(9),實現(xiàn)干涉相位到高程的轉(zhuǎn)換。由于式(9)計算的是BP算法成像過程中參考面相對真實地形的高程偏差,需要將其加回到成像過程中采用的輔助DEM,即可得到最終反演的DEM結(jié)果。
利用仿真實驗驗證基于DEM輔助后向投影模型的InSAR高程反演方法的有效性及優(yōu)越性。仿真實驗的參數(shù)如表1所示。實驗?zāi)M地形起伏較大的區(qū)域,其高程設(shè)置如圖4(a)所示。根據(jù)InSAR回波信號模型,仿真生成主副天線的回波信號,并加入–30dB的高斯白噪聲。實驗過程中分別采用了本文所提的基于后向投影的InSAR方法與傳統(tǒng)InSAR方法對回波信號進(jìn)行成像和干涉處理。
表1 仿真參數(shù)Tab.1 Simulation parameters
首先利用后向投影InSAR高程反演方法,引入存在高程誤差的輔助DEM作為BP算法成像的參考面并生成聚焦的SAR圖像。干涉相位由主輔圖像對直接共軛相乘即可得到,無需圖像配準(zhǔn),如圖4(b)所示。根據(jù)2.2節(jié)的分析,基于后向投影算法的In-SAR高程反演模型獲得的干涉相位已去除了成像過程中大部分的地形相位,因此干涉條紋密度大大降低。同時由于輔助DEM的高程誤差較小,殘余的地形相位數(shù)值較小不足以引起干涉相位的纏繞,如圖4(b)所示,避免了相位解纏操作。
圖4 仿真實驗結(jié)果圖Fig.4 Simulation results
然后根據(jù)傳統(tǒng)InSAR高程反演方法,采用距離多普勒(RD)算法成像,對得到的復(fù)圖像對進(jìn)行配準(zhǔn)、去平地以及相位濾波處理,獲取的干涉相位如圖4(c)所示。由于仿真區(qū)域地形起伏較大,因此干涉條紋十分密集,并且在高程劇烈變化的區(qū)域出現(xiàn)了條紋模糊的問題,如圖4(c)中紅色矩形框區(qū)域所示,采用常規(guī)的干涉處理將無法得到該區(qū)域正確的相位解纏結(jié)果。
針對圖4(c)中的相位模糊問題,本文在傳統(tǒng)InSAR方法中進(jìn)一步引入DEM輔助相位解纏處理。將根據(jù)SAR圖像生成的原始干涉相位與基于輔助DEM模擬的干涉相位進(jìn)行差分處理,得到殘余干涉相位,如圖4(d)所示。該干涉相位去除了輔助DEM對應(yīng)的地形相位,大大降低了干涉條紋的密度。將殘余干涉相位與模擬的干涉相位求和即可得到原始干涉相位圖的解纏結(jié)果。
由圖4(b)和圖4(d)可以看出,輔助DEM的引入避免了條紋密集區(qū)域的相位模糊問題,并且可以得到正確的相位解纏結(jié)果。但是在傳統(tǒng)InSAR方法中引入輔助DEM需要模擬DEM對應(yīng)的干涉相位,并需要對模擬相位與原始干涉相位進(jìn)行配準(zhǔn),增加了干涉處理流程的復(fù)雜度。而本文所提后向投影In-SAR高程反演方法在解決相位模糊問題的同時簡化了干涉處理流程,因此具有更大的優(yōu)勢。
在生成InSAR回波信號時,仿真的實驗場景中設(shè)置了25個強(qiáng)散射點作為地面控制點,用于統(tǒng)計高程反演的誤差。本文所提后向投影InSAR方法的實驗結(jié)果如圖5(a)所示,25個控制點處的高程反演誤差的均值為–0.0326m,標(biāo)準(zhǔn)差為0.2510m。
圖5 兩種方法的高程反演誤差Fig.5 The elevation inversion errors of two methods
進(jìn)一步對比傳統(tǒng)InSAR高程反演方法,采用RD算法成像并在相位解纏過程中引入DEM輔助。統(tǒng)計控制點處的高程反演誤差,如圖5(b)所示,均值為0.0349m,標(biāo)準(zhǔn)差為0.3484m。由于BP算法的成像模型更加精確,相位保持精度更高,因此后向投影InSAR高程反演精度略優(yōu)于基于RD成像的傳統(tǒng)InSAR方法,即本文所提簡化的基于輔助DEM后向投影InSAR高程反演方法可以實現(xiàn)高精度DEM的反演。
仿真實驗驗證了本文方法在反演高精度DEM的同時,無需圖像配準(zhǔn)操作,大大降低干涉條紋的密度,甚至避免了相位解纏步驟,簡化了干涉處理的流程。在地形起伏劇烈的區(qū)域,基于后向投影成像模型的InSAR高程反演方法通過輔助DEM的引入,可以有效地避免相位模糊問題。
(1)平地區(qū)域?qū)嶒灲Y(jié)果
在本節(jié)中,采用第2節(jié)提出的基于后向投影的InSAR高程反演方法對實測機(jī)載InSAR數(shù)據(jù)進(jìn)行處理。實驗所用的機(jī)載InSAR數(shù)據(jù)是在中國內(nèi)蒙古獲取的X波段雙天線數(shù)據(jù),高程測量標(biāo)稱精度1m,具體參數(shù)如表2所示。
表2 X波段機(jī)載InSAR系統(tǒng)參數(shù)Tab.2 X-band airborne InSAR system parameters
實驗中選取的成像區(qū)域大小為5.0km×1.6km,引入30m分辨率的SRTM DEM作為輔助DEM,利用BP算法獲取聚焦的SAR圖像,并與RD算法得到的SAR圖像對比。選取成像場景中的強(qiáng)散射點目標(biāo),對該目標(biāo)鄰近區(qū)域進(jìn)行16倍FFT插值,圖6(a)和圖6(b)分別是插值后地面強(qiáng)散射點處BP算法和RD算法的主輔圖像對成像結(jié)果,圖6(c)和圖6(d)分別為目標(biāo)沿距離向的能量分布??梢钥闯?,BP算法得到的主輔圖像已配準(zhǔn),而RD算法得到的主輔圖像對中目標(biāo)位置在距離向偏移了59個單元。插值結(jié)果證明了本文所提后向投影InSAR高程反演模型可以實現(xiàn)自配準(zhǔn),避免圖像配準(zhǔn)過程。
圖6 強(qiáng)散射目標(biāo)點成像結(jié)果Fig.6 The focused images of intense scatterer
BP算法得到的主圖像幅度圖如圖7(a)所示。實驗過程中在該成像區(qū)域布放的地面控制點位置信息如圖7(a)中紅色三角形標(biāo)注。由于主輔圖像已實現(xiàn)配準(zhǔn),因此直接共軛相乘即可得到干涉相位,如圖7(b)所示。根據(jù)第2節(jié)的理論分析,基于DEM輔助后向投影模型的InSAR高程反演方法由于在成像過程中引入了外源DEM,可直接得到去地形相位后的干涉相位,該干涉相位由輔助DEM成像參考面與真實地形間高程偏差引起,因此殘余干涉相位數(shù)值較小,如圖7(b)所示,變化范圍小于2π/5,避免了相位解纏處理步驟。相位濾波后直接進(jìn)行高程反演,輔助DEM的高程偏差如圖7(c)所示。將計算的偏差與輔助DEM的高程值相加,最終獲得成像區(qū)域的高程反演結(jié)果,如圖7(d)所示。
圖7 機(jī)載InSAR數(shù)據(jù)處理結(jié)果圖Fig.7 Airborne InSAR data processing results
利用實驗中布放的部分地面控制點(未參與干涉定標(biāo))作為檢查點,對傳統(tǒng)方法與本文方法高程反演誤差進(jìn)行對比驗證,如表3所示。對比結(jié)果表明,傳統(tǒng)InSAR成像和干涉處理結(jié)果能夠達(dá)到1m標(biāo)稱高程精度要求(高程精度0.85m),然而利用本文提出的DEM輔助干涉處理算法可以達(dá)到更高的處理精度(高程精度0.65m)。由于后向投影成像和干涉處理幾何,是一種完全時域的無誤差成像算法,在成像信號模型、運動誤差補(bǔ)償模型上理論最為精確,相位保持精度更高;同時在外源DEM輔助下,避免配準(zhǔn)和相位解纏過程,減少處理過程中的相位損失。因此,本文方法處理精度優(yōu)于傳統(tǒng)InSAR成像和干涉處理算法,可以獲得高精度的DEM反演結(jié)果,且具有簡化的干涉處理流程。
表3 地面檢查點處高程反演誤差Tab.3 The elevation inversion errors of ground detection points
(2)丘陵區(qū)域?qū)嶒灲Y(jié)果
本文另外選取了丘陵區(qū)域的實驗數(shù)據(jù),進(jìn)一步驗證所提的基于DEM輔助后向投影模型的In-SAR高程反演方法。BP算法得到的主輔圖像對直接共軛相乘得到干涉相位如圖8(a)所示。雖然實驗區(qū)域地形起伏較大,但由于輔助DEM的引入,基于后向投影的InSAR模型在成像過程去除了絕大多數(shù)的地形相位,圖8(a)中的干涉相位范圍小于π,無需相位解纏,簡化了后續(xù)的干涉處理。而傳統(tǒng)InSAR方法得到的干涉相位如圖8(b)所示,由于平地效應(yīng)與地形起伏,圖8(b)中存在3個纏繞的干涉條紋,與圖8(a)相比后續(xù)干涉處理步驟更加復(fù)雜。圖8(c)是利用BP算法得到的幅度圖。根據(jù)DEM輔助下基于后向投影的InSAR模型得到的干涉相位及式(9)反演輔助DEM的偏差,并加回到輔助DEM,最終得到該成像區(qū)域的DEM如圖8(d)所示。由丘陵區(qū)域?qū)嶒灁?shù)據(jù)的處理結(jié)果可以看出,本文基于輔助DEM與BP成像算法的InSAR方法可以有效地減少干涉相位的纏繞,避免相位解纏處理,驗證了所提基于后向投影的InSAR高程反演模型的有效性。
圖8 機(jī)載InSAR山區(qū)數(shù)據(jù)處理結(jié)果圖Fig.8 Airborne InSAR data processing results of mountainous area
本文提出了一種基于DEM輔助后向投影模型的InSAR高程反演方法。該方法利用BP算法運動補(bǔ)償精度高和保相性能好的優(yōu)勢,同時在BP成像過程中引入外源DEM輔助,降低干涉條紋的密度,在多數(shù)情況下還實現(xiàn)了干涉處理流程的簡化—無需圖像配準(zhǔn)與相位解纏操作。本文最后設(shè)計了仿真實驗,并對X波段雙天線機(jī)載In-SAR數(shù)據(jù)進(jìn)行了處理,通過所提方法中簡化的干涉處理流程實現(xiàn)了高精度DEM的反演,驗證了該方法的有效性。然而,在基于BP算法的干涉模型中還存在一些問題,例如實際的InSAR系統(tǒng)中天線相位中心存在測量誤差,在高分辨成像時需要考慮對方位向配準(zhǔn)產(chǎn)生的影響,以及本文模型的高程反演誤差理論分析等,這將是今后的研究方向。