秦小芳 張華春 張 衡② 王 宇 劉華有②

①(中國(guó)科學(xué)院電子學(xué)研究所 北京 100190)

②(中國(guó)科學(xué)院大學(xué) 北京 100049)

1 引言

數(shù)字高程模型(Digital Elevation Model, DEM)是描述地球表面形狀的3維數(shù)字模型，由一系列包含有地理平面坐標(biāo)和高程的數(shù)據(jù)集組成[1]，在科學(xué)研究、經(jīng)濟(jì)建設(shè)和軍事領(lǐng)域都有著重要的應(yīng)用價(jià)值。某些特定的應(yīng)用場(chǎng)合如地震形變提取、地形監(jiān)測(cè)，高分辨率高精度DEM顯得尤為重要，但是其提取通常非常復(fù)雜。干涉合成孔徑雷達(dá)(Interferometric Synthetic Aperture Radar, InSAR)由于其全天時(shí)、全天候的工作特性，且作為一種主動(dòng)式傳感器，成為精確獲取DEM的成熟技術(shù)之一。目前覆蓋全球的DEM數(shù)據(jù)非常有限，應(yīng)用最廣泛的是SRTM[2]，該數(shù)據(jù)集覆蓋了全球陸地面積的80%以上，在美國(guó)區(qū)域公開(kāi)的分辨率為30 m，其他區(qū)域的分辨率為90 m。在2010年，德國(guó)宇航局(DLR)發(fā)射的TanDEM-X衛(wèi)星與2007年發(fā)射的TerraSAR-X衛(wèi)星組成了真正意義上的全球第1個(gè)星載雙基SAR系統(tǒng)[3,4]，其首要目的是獲取全球高精度數(shù)字高程模型。其標(biāo)準(zhǔn)產(chǎn)品是符合HRTI-3標(biāo)準(zhǔn)的DEM，分辨率為12 m，相對(duì)高程精度為2 m，絕對(duì)高程精度為10 m[5–7]。在特定的區(qū)域，DLR也制作了符合HRTI-4標(biāo)準(zhǔn)分辨率為6 m的DEM[8]。

在DEM的獲取過(guò)程中非常關(guān)鍵的一步是相位解纏，然而由于相位噪聲的存在，高精度相位解纏變得非常困難，有時(shí)甚至出現(xiàn)嚴(yán)重的錯(cuò)誤，這嚴(yán)重影響了DEM的準(zhǔn)確度和精度[9]，文獻(xiàn)[10]分析了多視數(shù)與相位梯度和噪聲梯度的關(guān)系，采用多視迭代的方法降低相位噪聲，得到較高分辨率和精度的DEM，文獻(xiàn)[11]也采用多次迭代的方法，并將升降軌數(shù)據(jù)利用基線信息進(jìn)行融合，在山區(qū)和稀疏城區(qū)得到較高分辨率和精度的DEM。但是采用多次迭代的方式，重復(fù)步驟較多，且處理誤差在迭代的過(guò)程中會(huì)產(chǎn)生傳遞，不利于大范圍快速的DEM生成。非局部干涉(NonLocal Interferometric SAR,NL-InSAR)相位濾波是近幾年發(fā)展比較迅速的一種可以精確估計(jì)干涉條紋同時(shí)降低噪聲的濾波方法，文獻(xiàn)[12]將NL-InSAR應(yīng)用在干涉圖估計(jì)中，取得了較好的降噪效果，同時(shí)細(xì)節(jié)信息得到保持。文獻(xiàn)[8]采用非局部濾波得到分辨率為6 m的高精度DEM。本文采用NL-InSAR方法，抑制相位噪聲，只需進(jìn)行一次多視處理，避免了多次迭代時(shí)的誤差傳遞問(wèn)題，得到分辨率為4 m的DEM；同時(shí)，利用外部DEM輔助進(jìn)行干涉處理，可以提高干涉處理的可靠性、可解性和精度[13,14]，從降低干涉圖噪聲和干涉條紋頻率兩個(gè)方面降低相位解纏的難度，最終得到較高分辨率和較高精度的DEM。另外NLInSAR可以利用較快的并行計(jì)算方案，在較快時(shí)間內(nèi)完成高質(zhì)量濾波處理。

在山區(qū)和地表起伏較大的復(fù)雜地形區(qū)域，雷達(dá)的側(cè)視成像造成疊掩和陰影[15]等數(shù)據(jù)缺失區(qū)域，影響了生成的DEM質(zhì)量。升降軌的視線是對(duì)稱的，因此疊掩和陰影區(qū)域的缺失數(shù)據(jù)可以通過(guò)融合這兩種互補(bǔ)的數(shù)據(jù)進(jìn)行恢復(fù)。已有研究大多采用加權(quán)平均法進(jìn)行DEM的融合，Carrasco等人基于相干性融合了升降軌ERS DEM數(shù)據(jù)，首先設(shè)定相干性閾值，對(duì)于待確定的像素點(diǎn)，如果兩幅圖像在該點(diǎn)的相干性都高于閾值，就采用加權(quán)平均的方法計(jì)算，否則選擇相干性高的DEM估計(jì)作為該點(diǎn)的高程[16]。文獻(xiàn)[17]基于相干系數(shù)融合多側(cè)視方向DEM，并在融合過(guò)程中借助陰影和疊掩區(qū)識(shí)別方法、陰影去除方法以及相干系數(shù)圖，提高了DEM精度。但同時(shí)利用多頻和多基線的干涉結(jié)果時(shí)，單一的選擇權(quán)重并不合理，文獻(xiàn)[18]基于擴(kuò)展的卡爾曼濾波器融合多頻和多基線DEM制作了高質(zhì)量的DEM。針對(duì)TerraSAR-X/TanDEM-X雙基升降軌數(shù)據(jù)融合，近幾年也有一些研究工作，文獻(xiàn)[19]采用理論高度誤差圖計(jì)算權(quán)重融合了雙基SAR DEM，然而隨著數(shù)據(jù)源的增多，會(huì)增加升降軌融合的計(jì)算量，而文獻(xiàn)[11]利用基線信息計(jì)算權(quán)值融合升降軌DEM，在稀疏城區(qū)得到的結(jié)果較好，但在建筑物過(guò)于密集的區(qū)域精度較前者偏低。本文在NL-InSAR的基礎(chǔ)上，提出一種基于相干系數(shù)的升降軌DEM融合方法，先將升降軌DEM中疊掩陰影等數(shù)據(jù)缺失區(qū)域互相恢復(fù)，然后利用NL-InSAR無(wú)偏估計(jì)的相干系數(shù)進(jìn)行加權(quán)融合，得到精度較高的DEM。

2 基于NL-InSAR的快速高分辨率高精度DEM重建技術(shù)

兩幅配準(zhǔn)的單視復(fù)圖像(Single Look Complex,SLC)每一個(gè)分辨單元可以用3個(gè)參數(shù)來(lái)描述：散射系數(shù)、真實(shí)的干涉相位、相干系數(shù)。假設(shè)z和是兩幅配準(zhǔn)圖像上對(duì)應(yīng)的像素點(diǎn)，根據(jù)Goodman模型[20]，z與z′服從如下分布：

其中，R與R′是散射系數(shù)，β是真實(shí)相位，D是相干系數(shù)。

非局部相位濾波[12]對(duì)3個(gè)參數(shù)R,β,D聯(lián)合估計(jì)，其結(jié)果一般會(huì)優(yōu)于僅利用相位信息的濾波方法，SAR成像過(guò)程表明相位突變的地方通常伴隨著幅度或者相干性的變化，因此不同的參數(shù)相互支撐彼此的估計(jì)。其核心是基于圖像分塊間的相似性質(zhì)去除噪聲，采用較大的搜索窗搜索相似的像素點(diǎn)，相似的像素不再局限于某個(gè)局部區(qū)域，當(dāng)前像素的估計(jì)值由圖像中與它具有相似鄰域結(jié)構(gòu)的像素加權(quán)得到。實(shí)際處理過(guò)程中通常分為3個(gè)步驟[21]，如圖1(a)所示，首先根據(jù)相似性準(zhǔn)則搜索出相似塊，塊之間相似性通常定義為在每個(gè)對(duì)應(yīng)像素的相似度的總和為所選取的像素的相似性；選出相似塊之后，這些塊在第2步中結(jié)合起來(lái)估計(jì)雷達(dá)特征，如相干系數(shù)等；最后將第2步中的估計(jì)值投影到雷達(dá)圖像空間。

假設(shè)兩幅圖像的散射特性相同，即R=R′，這種假設(shè)在相干性較高的區(qū)域是成立的，同時(shí)令為去平地后的相位，則有：

NL-InSAR使用權(quán)重最大似然估計(jì)(Weighted Maximum Likelihood Estimation, WMLE)來(lái)估計(jì)權(quán)重，在分塊s中，將參數(shù)估計(jì)記作：

在文獻(xiàn)[22]中Seymour和Cumming區(qū)別于傳統(tǒng)的最大似然估計(jì)假設(shè)，令R=R′，將公式擴(kuò)展到WMLE中，令Θs=(Rs，βs，Ds)以及Ot=(At，At′，φt)。

由于噪聲的存在，基于局部小窗口的相干性估計(jì)不僅會(huì)受到噪聲的影響，而且會(huì)產(chǎn)生一定的偏差。而NL-InSAR濾波從非局部的像素塊估計(jì)干涉圖參數(shù)，對(duì)于相干性的估計(jì)更加準(zhǔn)確。圖2對(duì)比了NL-InSAR與傳統(tǒng)boxcar濾波器的相干性估計(jì)，圖中可以看出，NL-InSAR估計(jì)的相干系數(shù)更接近真實(shí)值，在相干性較低的區(qū)域，傳統(tǒng)濾波器過(guò)估計(jì)現(xiàn)象尤其明顯。在文獻(xiàn)[8]中提到基于NL-InSAR濾波估計(jì)得到的相干系數(shù)相當(dāng)于13×13的boxcar濾波器估計(jì)的相干系數(shù)。

3 基于相干系數(shù)的升降軌DEM融合技術(shù)

透視收縮現(xiàn)象是在雷達(dá)衛(wèi)星天線中心到達(dá)地物頂部和底部的斜距之差小于其所對(duì)應(yīng)的地距之差的情形之下發(fā)生的。如圖3所示透視收縮現(xiàn)象進(jìn)一步發(fā)展，當(dāng)傳感器天線中心到地物頂部的斜距小于或等于其到地物底部的斜距時(shí)，則會(huì)在雷達(dá)影像中出現(xiàn)頂?shù)椎怪煤童B掩現(xiàn)象[23]。在陡峭的地形區(qū)域，當(dāng)迎坡角大于入射角時(shí)，SAR的側(cè)視特性會(huì)使同一分辨率單元內(nèi)多個(gè)信號(hào)具有相同的距離和多普勒頻率，這就產(chǎn)生了疊掩。當(dāng)?shù)匦蔚谋称陆谴笥?90°–入射角)時(shí)，地形的背面沒(méi)有回波信號(hào)就產(chǎn)生了陰影?？傊?，影像所固有的透視收縮、頂?shù)椎怪?、疊掩和陰影等幾何畸變會(huì)給干涉測(cè)量數(shù)據(jù)處理和解譯帶來(lái)極大的不便，當(dāng)?shù)匦纹鸱^大時(shí)，這些幾何畸變甚至?xí)?yán)重影響干涉測(cè)量結(jié)果的精度。

單基線DEM的獲取過(guò)程采用非局部濾波和外部DEM差分處理保證了相位的可靠性，但是對(duì)于側(cè)視成像本身帶來(lái)的缺失數(shù)據(jù)無(wú)法恢復(fù)。星載SAR傳感器的平臺(tái)通常是太陽(yáng)同步軌道衛(wèi)星，它可以從不同的軌道方向，對(duì)同一地區(qū)進(jìn)行觀測(cè)，獲取升軌和降軌干涉像對(duì)，如圖4所示，構(gòu)成升降軌InSAR。若觀測(cè)時(shí)側(cè)視方向保持不變，則雷達(dá)波束的照射方向相反，升軌和降軌干涉圖將為同一地物提供相反入射方向的干涉相位，可以有效填補(bǔ)疊掩和陰影等地區(qū)的缺失數(shù)據(jù)。

在進(jìn)行融合之前，先對(duì)升降軌DEM進(jìn)行預(yù)處理，如圖5所示。先提取公共的區(qū)域，將疊掩陰影區(qū)域進(jìn)行掩膜處理，如果疊掩和陰影區(qū)域在另外一個(gè)DEM中處于非疊掩和陰影區(qū)域并且在相干性閾值之上，則直接用另一個(gè)DEM的高度值取代；如果兩者都為疊掩和陰影區(qū)域或者相干性都處于閾值之下則仍設(shè)為無(wú)效值；如果不為疊掩陰影等區(qū)域并且相干性在閾值之上，則保持不變，如表1所示，升降軌DEM處理時(shí)按照此邏輯處理表進(jìn)行預(yù)處理。

預(yù)處理之后的DEM，在有效區(qū)域加權(quán)融合的權(quán)值選擇尤為重要。在本文的處理流程中，干涉圖經(jīng)過(guò)了NL-InSAR處理，得到了無(wú)偏估計(jì)的相干系數(shù)，將有效區(qū)域利用相干系數(shù)進(jìn)行加權(quán)，如式(6)所示，整個(gè)處理流程總結(jié)為如圖6所示。

表1 升降軌DEM融合處理表Tab.1 Logic table for ascending and descending pass TanDEM-X raw DEMs fusion

其中，γi為相干系數(shù)，N為待融合的DEM的數(shù)量，在本文中為2，當(dāng)影像對(duì)多于2時(shí)，該方法依然可以擴(kuò)展。

4 實(shí)驗(yàn)及分析

4.1 實(shí)驗(yàn)數(shù)據(jù)集

本實(shí)驗(yàn)采用TerraSAR-X/TanDEM-X雙基觀測(cè)模式數(shù)據(jù)，即兩顆衛(wèi)星同時(shí)經(jīng)過(guò)觀測(cè)區(qū)域，其中一顆衛(wèi)星發(fā)射電磁波并接收地面的反射回波，另一顆衛(wèi)星同時(shí)接收地面的反射回波，如圖7(a)所示。實(shí)驗(yàn)地區(qū)位于北京西部地區(qū)如圖7(b)所示，影像內(nèi)包括密集建筑區(qū)和山區(qū)，散射特性復(fù)雜，包含較多陰影和疊掩區(qū)域。其中升軌干涉影像對(duì)于2014年8月19日獲取，垂直基線104.27 m，模糊高度為73.31 m，方位向和距離向的分辨率為1.36 m和1.98 m，入射角為43.34°，因此地距分辨率為2 m；降軌干涉影像對(duì)于2014年4月17日獲取，垂直基線89.22 m，模糊高度為93.21 m，方位向和距離向的分辨率為1.36 m和2.05 m，入射角為45.86°，地距分辨率為2 m，數(shù)據(jù)統(tǒng)計(jì)信息如表2所示。

4.2 基于NL-InSAR的單基線DEM

首先采用NL-InSAR濾波，并與目前在InSAR中應(yīng)用較多的Goldstein濾波進(jìn)行比較。得到的濾波效果圖如圖8所示，圖8(a)是NL濾波的結(jié)果，圖8(b)是Goldstein濾波的結(jié)果，從中可以明顯看出，NLInSAR濾波效果較好，在有效抑制噪聲的同時(shí)對(duì)于細(xì)節(jié)特征的保持更為豐富；而基于窗口的Goldstein濾波含噪更多，在人造建筑區(qū)域如道路和建筑的細(xì)節(jié)損失較為明顯。

同時(shí)，NL-InSAR濾波對(duì)于相干性的估計(jì)更加準(zhǔn)確，選取相干系數(shù)圖部分結(jié)果如圖9所示，圖9(a)采用NL-InSAR估計(jì)得到的相干系數(shù)圖，圖9(b)是Goldstein濾波估計(jì)的相干系數(shù)圖。圖像所在區(qū)域中存在相干性比較高的建筑物、失相干區(qū)域的水域、相干系數(shù)較低的草地，在非局部濾波的相干系數(shù)圖可以明顯看出這些差別，而Goldstein濾波之后的相干系數(shù)整體偏高，同時(shí)從相干系數(shù)統(tǒng)計(jì)直方圖可以看出，Goldstein濾波器得到的相干系數(shù)較集中分布在高相干性區(qū)域，NL-InSAR得到的相干系數(shù)更加符合所在區(qū)域的相干性分布特征，Goldstein濾波器相干系數(shù)存在過(guò)估計(jì)的現(xiàn)象，非局部濾波產(chǎn)生無(wú)偏估計(jì)的相干系數(shù)的同時(shí)保持了圖像的邊緣細(xì)節(jié)。

表2 數(shù)據(jù)集Tab.2 Datasets

為了降低相位解纏的難度，減少相位解纏誤差，利用外部DEM進(jìn)行差分處理，削弱區(qū)域地形引起的密集條紋。根據(jù)SAR系統(tǒng)參數(shù)、軌道參數(shù)等，利用外部DEM仿真SAR圖像，將仿真SAR圖像與實(shí)際SAR圖像配準(zhǔn)，建立外部DEM與實(shí)際SAR圖像的對(duì)應(yīng)關(guān)系[13]。再利用外部DEM仿真干涉圖，將原始干涉圖與仿真干涉圖進(jìn)行差分干涉處理得到纏繞的差分干涉圖，解纏差分干涉圖，并加上原來(lái)的仿真干涉圖，得到原始干涉圖的解纏相位。最后根據(jù)InSAR幾何關(guān)系模型反演高度，再進(jìn)行地理編碼得到升降軌DEM。采用文獻(xiàn)[10]中多視迭代的方法，本文也進(jìn)行了多視迭代的實(shí)驗(yàn)，結(jié)果如圖10(b)所示，可以看到采用NL-InSAR相位濾波得到的結(jié)果細(xì)節(jié)更加清晰，且只進(jìn)行一次多視處理，處理過(guò)程快速簡(jiǎn)單。

4.3 升降軌融合DEM精度分析

在地形復(fù)雜區(qū)域，數(shù)據(jù)受疊掩和陰影影響的無(wú)效區(qū)域較多，因此利用升軌和降軌分別重建的DEM互相融合從而恢復(fù)無(wú)效區(qū)域。升降軌數(shù)據(jù)處理分別得到各自的DEM，將其地理編碼到WGS84坐標(biāo)系下，根據(jù)經(jīng)緯度信息提取出公共部分，如圖11所示。為了進(jìn)一步提升升降軌DEM配準(zhǔn)精度，同樣將升降軌幅度進(jìn)行地理編碼，將地理編碼后的升降軌幅度進(jìn)行配準(zhǔn)，使用最小二乘得到擬合的多項(xiàng)式系數(shù)，使用該多項(xiàng)式系數(shù)實(shí)現(xiàn)對(duì)升降軌DEM的配準(zhǔn)。升軌的無(wú)效點(diǎn)數(shù)占總共的4.93%，降軌數(shù)據(jù)占4.52%，升軌數(shù)據(jù)無(wú)效區(qū)域相對(duì)于降軌數(shù)據(jù)較高的原因是由于其入射角低于降軌數(shù)據(jù)，因此實(shí)驗(yàn)區(qū)域有更多的像素點(diǎn)由于坡度角超過(guò)入射角而受到疊掩陰影等的影響。將升降軌的無(wú)效區(qū)域互相恢復(fù)之后，得到如圖12、圖13所示融合后的部分DEM結(jié)果，圖中可以明顯看到無(wú)效區(qū)域的減少。統(tǒng)計(jì)結(jié)果如表3中所示，融合之后無(wú)效區(qū)域比例減少到1.34%。

表3 缺失區(qū)域統(tǒng)計(jì)表Tab.3 Statistics table of invalids

表4 與SRTM DEM對(duì)比高程殘差統(tǒng)計(jì)表(SRTM)Tab.4 Comparison of height difference with respect to SRTM DEM

為了對(duì)所獲取的DEM進(jìn)行精度分析，采用SRTM和地面控制點(diǎn)進(jìn)行對(duì)比。與SRTM進(jìn)行對(duì)比如表4所示，在平地融合之后均方根誤差(Root Mean Square Error, RMSE)為6.61 m，而升軌DEM為7.37 m，降軌DEM為6.81 m；在山區(qū)精度不如平地，但是融合之后山地比融合前精度有所提升。由于SRTM DEM本身精度有限，參考文獻(xiàn)[10]采用2個(gè)地面控制點(diǎn)對(duì)獲取的DEM進(jìn)行精度評(píng)價(jià)，本文采用于2016年實(shí)地測(cè)量獲取的4個(gè)地面控制點(diǎn)對(duì)DEM進(jìn)行精度分析，如圖11中紅色圓點(diǎn)所示，地面控制點(diǎn)位于地形較為平緩的區(qū)域，地面控制點(diǎn)的測(cè)量采用靜態(tài)測(cè)量精度達(dá)到毫米級(jí)的第3代GNSS (Global Navigation Satelllite System)接收機(jī)。精度情況如表5所示，融合后平地DEM精度為6.09 m，而升軌DEM為6.74 m，降軌DEM為6.67 m。在山區(qū)由于地面控制點(diǎn)數(shù)量有限，不能進(jìn)行精確地對(duì)比，但是可以預(yù)測(cè)融合之后山區(qū)DEM精度也會(huì)有所提高。

從統(tǒng)計(jì)結(jié)果來(lái)看，融合后的DEM精度有一定的提升。但仍然可能存在一定的誤差。分析其原因，數(shù)據(jù)的一部分區(qū)域北京城區(qū)高樓林立，建筑物分布非常密集，散射特性復(fù)雜，NL-InSAR在建筑物密集城區(qū)的適用性需要進(jìn)一步的研究；另外升降軌干涉影像對(duì)獲取的時(shí)間間隔較長(zhǎng)，且位于不同的季節(jié)，山區(qū)植被隨季節(jié)變化導(dǎo)致散射特性會(huì)產(chǎn)生差異，因此也有可能產(chǎn)生誤差。

5 結(jié)論

本文基于TerraSAR-X/TanDEM-X雙基SAR干涉數(shù)據(jù)，采用非局部濾波濾除相位噪聲，并利用外部DEM進(jìn)行差分處理，有效減少了相位解纏誤差，相比于多視迭代處理，這種方法處理更為便捷，并且有效防止了處理過(guò)程中的誤差傳遞?；谏弟墧?shù)據(jù)，將疊掩和陰影等無(wú)效區(qū)域進(jìn)行恢復(fù)，數(shù)據(jù)無(wú)效區(qū)域的比例相比于融合前的升軌4.93%、降軌4.52%，融合后降低到1.34%。利用相干系數(shù)進(jìn)行融合，最終得到分辨率為4 m、精度為6.09 m的高分辨率高精度DEM。

表5 高程精度Tab.5 Accuracy of different DEMs