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

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

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

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

1 引言

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

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

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

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

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

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

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

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

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

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

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

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

透視收縮現象是在雷達衛星天線中心到達地物頂部和底部的斜距之差小于其所對應的地距之差的情形之下發生的。如圖3所示透視收縮現象進一步發展，當傳感器天線中心到地物頂部的斜距小于或等于其到地物底部的斜距時，則會在雷達影像中出現頂底倒置和疊掩現象[23]。在陡峭的地形區域，當迎坡角大于入射角時，SAR的側視特性會使同一分辨率單元內多個信號具有相同的距離和多普勒頻率，這就產生了疊掩。當地形的背坡角大于(90°–入射角)時，地形的背面沒有回波信號就產生了陰影。總之，影像所固有的透視收縮、頂底倒置、疊掩和陰影等幾何畸變會給干涉測量數據處理和解譯帶來極大的不便，當地形起伏較大時，這些幾何畸變甚至會嚴重影響干涉測量結果的精度。

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

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

預處理之后的DEM，在有效區域加權融合的權值選擇尤為重要。在本文的處理流程中，干涉圖經過了NL-InSAR處理，得到了無偏估計的相干系數，將有效區域利用相干系數進行加權，如式(6)所示，整個處理流程總結為如圖6所示。

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

其中，γi為相干系數，N為待融合的DEM的數量，在本文中為2，當影像對多于2時，該方法依然可以擴展。

4 實驗及分析

4.1 實驗數據集

本實驗采用TerraSAR-X/TanDEM-X雙基觀測模式數據，即兩顆衛星同時經過觀測區域，其中一顆衛星發射電磁波并接收地面的反射回波，另一顆衛星同時接收地面的反射回波，如圖7(a)所示。實驗地區位于北京西部地區如圖7(b)所示，影像內包括密集建筑區和山區，散射特性復雜，包含較多陰影和疊掩區域。其中升軌干涉影像對于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，數據統計信息如表2所示。

4.2 基于NL-InSAR的單基線DEM

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

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

表2 數據集Tab.2 Datasets

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

4.3 升降軌融合DEM精度分析

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

表3 缺失區域統計表Tab.3 Statistics table of invalids

表4 與SRTM DEM對比高程殘差統計表(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；在山區精度不如平地，但是融合之后山地比融合前精度有所提升。由于SRTM DEM本身精度有限，參考文獻[10]采用2個地面控制點對獲取的DEM進行精度評價，本文采用于2016年實地測量獲取的4個地面控制點對DEM進行精度分析，如圖11中紅色圓點所示，地面控制點位于地形較為平緩的區域，地面控制點的測量采用靜態測量精度達到毫米級的第3代GNSS (Global Navigation Satelllite System)接收機。精度情況如表5所示，融合后平地DEM精度為6.09 m，而升軌DEM為6.74 m，降軌DEM為6.67 m。在山區由于地面控制點數量有限，不能進行精確地對比，但是可以預測融合之后山區DEM精度也會有所提高。

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

5 結論

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

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