TerraSAR-X/TanDEM-X升降軌雙基干涉模式獲取DEM方法研究

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

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

預(yù)處理之后的DEM，在有效區(qū)域加權(quán)融合的權(quán)值選擇尤為重要。在本文的處理流程中，干涉圖經(jīng)過了NL-InSAR處理，得到了無偏估計的相干系數(shù)，將有效區(qū)域利用相干系數(shù)進行加權(quán)，如式(6)所示，整個處理流程總結(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)影像對多于2時，該方法依然可以擴展。

4 實驗及分析

4.1 實驗數(shù)據(jù)集

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

4.2 基于NL-InSAR的單基線DEM

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

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

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

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

4.3 升降軌融合DEM精度分析

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

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

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

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

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

5 結(jié)論

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

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