TerraSAR-X/TanDEM-X獲取高精度數字高程模型技術研究

杜亞男， 馮光財， 李志偉， 朱建軍， 彭星

中南大學地球科學與信息物理工程學院雷達遙感研究室， 長沙 410083

?

TerraSAR-X/TanDEM-X獲取高精度數字高程模型技術研究

杜亞男， 馮光財*， 李志偉， 朱建軍， 彭星

中南大學地球科學與信息物理工程學院雷達遙感研究室， 長沙 410083

以雙星系統(tǒng)(TerraSAR-X/TanDEM-X)下的bistatic數據模式為例研究了差分干涉獲取高精度DEM產品的融合算法和技術流程.針對不同觀測幾何條件下(升降軌，不同入射角)觀測數據的畸變和缺失，提出一種迭代的顧及垂直基線、陰影和疊影的數據融合新方法重建高分辨率高精度的數字高程模型，并對TanDEM-X融合的DEM在不同地物屬性特征下(山區(qū)及高樓林立的城區(qū))的精度進行定量分析.本文采用了兩對覆蓋珠海、澳門區(qū)域的升降軌TerraSAR-X/TanDEM-X干涉對進行融合處理，并通過收集的高精度Lidar數據進行精度比較分析.此外，本文還定量分析了TanDEM-X的DEM對常規(guī)DInSAR技術的改進，并與SRTM、ASTER的結果進行對比.結果表明：提出的升降軌融合方法較單一軌道平臺能夠較好地改正或減弱由于幾何畸變引起的高程信息缺失或錯誤，通過與Lidar數據的對比發(fā)現TanDEM-X的融合DEM在山區(qū)的精度較高，其殘差的標準差為3.5 m，較單一軌道(升、降軌)分別降低8%和22%，能夠通過迭代的方法獲取高分辨率(可達2～5 m)、高精度的地形信息；而針對城區(qū)密集建筑物的復雜地形來說，融合的DEM的精度稍低，殘差的標準差為11.8 m，但較單一軌道(升、降軌)來說有較大改進，其殘差標準差分別降低了28%和22%；而在分布較為稀疏的居民區(qū)，融合的DEM能夠較Lidar數據獲取更好的建筑物高度及外形信息，此時的殘差標準差可達5 m.同時，TanDEM-X的融合DEM作為外部DEM能夠較SRTM和ASTER來說更好地去除地形信息，尤其在山區(qū)及高程建筑密集分布的城區(qū)，從而利于后續(xù)的相位解纏和形變信息的精確獲取和解譯，為更高精度的時序InSAR形變監(jiān)測提供有利條件.

數字高程模型； TerraSAR-X/TanDEM-X； 幾何畸變； 升降軌融合

1 引言

數字高程模型(DEM)是用一組有序數值陣列形式表示地面高程的一種實體地面模型，它是地學領域用來表示地形特征的重要表達形式.它作為一種基礎的地理信息數據已經在我們日常生活以及國防領域有著非常廣泛的應用，如地形制導的軍事打擊、差分干涉雷達形變監(jiān)測、地質勘測、各種線路(高鐵、城際、高速公路、電力線路等)的規(guī)劃、城區(qū)地籍變化監(jiān)測等.然而目前能覆蓋全球的DEM數據還非常有限，一種是SRTM，2000年由美國太空總署(NASA)和國防部國家測繪局(NIMA)聯(lián)合發(fā)射的“奮進”號航天飛機測量得到的，覆蓋了全球陸地面積80%以上，在美國區(qū)域公開的DEM分辨率為30 m，全球其他區(qū)域公開的DEM分辨率為90 m (Farr & Kobrick, 2000; Feng et al., 2015a)；另外2011年10月美國空間局(NASA)和日本經濟產業(yè)省(METI)發(fā)布了ASTER GDEM V2.0，它覆蓋了全球陸地面積的99%以上范圍，其空間分辨率可達30 m，雖然其在空間分辨率和高程精度方面都較V1.0有較大提升，但是精度仍然非常有限(Tachikawa et al., 2011; Feng et al., 2015b).當然，除了以上兩類全球覆蓋DEM外，德國宇航局(DLR)在2011年也公開了其X波段的SRTM，其分辨率較前兩類來說有著進一步的提高(約25 m)，但局限于網狀不規(guī)則覆蓋(Walker et al., 2007).值得注意的是SRTM和ASTER都是基于十多年前的影像數據為基礎生成的DEM，因此很多區(qū)域現實性都較差，無法反映最新的地形特征.

雖然通過航空攝影測量、機載Lidar、地面GPS和全站儀等方法也可以獲取局部的高分辨率數字高程模型，但是這些方法不僅耗時耗力而且容易受到天氣、環(huán)境等眾多因素的干擾，因此本文重點關注利用星載InSAR技術獲取大范圍DEM信息.為了獲取更高分辨率、更高精度和更新的DEM，國內外InSAR學者進行了眾多有益的嘗試，并提出了許多算法和技術流程.首先，廖明生和林輝(2003)、肖金群等(2012)都利用歐空局ERS1/2串行飛行克服時間去相關的影響(時間間隔1天)，基于傳統(tǒng)SRTM的輔助，成功獲取了研究區(qū)域內更高分辨率和更高精度的DEM；另外Ferretti等(2012)也嘗試利用時序InSAR的思想來克服時間去相關和大氣延遲等誤差的影響，從而重建高精度的DEM；近年來隨著SAR系統(tǒng)硬件的提升，SAR影像的空間分辨率越來越高(1 m左右)和重返周期也越來越小(3天)，因此通過融合不同傳感器、升降軌和不同入射角的DEM測量數據也成為了一種生成高分辨率和高精度DEM的技術手段(Jiang et al., 2014; Ashutosh et al., 2013; Feng et al., 2015c).雖然以上各種方法均能在一定條件上獲取較好的DEM，但是其監(jiān)測方式都為重復軌道干涉測量，仍然容易受到時間、空間上的去相關以及大氣延遲和軌道誤差的影響.此外，受到各SAR系統(tǒng)參數，性能以及覆蓋范圍和影像數量等多方面制約，以上各種方法并不能保證在任何區(qū)域均適用.鑒于此，德國宇航局(DLR)2010年發(fā)射TanDEM-X衛(wèi)星與之前2007年發(fā)射的TerraSAR-X衛(wèi)星進行編隊飛行，組成了一個雙星分布式SAR系統(tǒng).目的是通過這種雙星編隊的組合方式獲取覆蓋全球的高精度數字高程模型(Krieger et al., 2007)，從而為我們提供了一個嶄新的、更加廣闊的、可靠的DEM制作研究平臺.

TerraSAR-X/TanDEM-X(TSX/TDX)相比于傳統(tǒng)InSAR系統(tǒng)最主要特點就是能夠獲取“0”時間基線和高精度軌道的干涉對，能夠很好地克服各種失相關、大氣延遲和軌道誤差等引起的相位噪聲(Krieger et al., 2007).經過近幾年的發(fā)展，圍繞TanDEM-X和TerraSAR-X獲取高精度的DEM也開展了一些研究工作，主要分為兩方面：1)TSX/TDX數據處理與DEM融合.現在常用的DEM融合算法主要針對升降軌、不同傳感器InSAR和光學DEM等數據融合，而融合的關鍵便是定權，主要是通過與相干性、多視視數相關的概率密度函數積分得到的干涉相位標準差以及高程模糊度來進行定權(Jiang et al., 2014; Rossi et al., 2013a; Deo et al., 2014)，它是傳統(tǒng)星載重復軌道InSAR系統(tǒng)生成DEM技術的延續(xù)，然而該算法隨著數據源的增多會急劇增加融合的計算量，不利于大范圍的DEM制作.更重要的是由于TSX/TDX系統(tǒng)的特殊性，其時間去相關、大氣噪聲等因素影響幾乎可以忽略不計，因此在該情況下再使用傳統(tǒng)概率密度函數積分的方法不但不能顯著提高融合的精度，反而增加了DEM融合計算的復雜度和計算量，嚴重制約DEM融合的批處理能力.另外，現有的升降軌融合方法并沒有對失相干處(如水域)及兩個軌道中均處于幾何畸變處的空值進行處理，這也使得融合后的DEM不完整(Jiang et al., 2014; Rossi et al., 2012; Deo et al., 2014)；2) TSX/TDX 生成的DEM精度分析，現有的研究區(qū)域多針對地形陡峭的山區(qū)或者建筑物分布較為稀疏的居民區(qū)進行精度分析(Rossi et al., 2013; Deo et al., 2013; Brautigam et al., 2012).而對于地形特征復雜的區(qū)域(如沿海區(qū)域的水域、幾何畸變較為嚴重的高樓林立的城區(qū)等)的數據融合以及精度評定等問題仍然缺少定量實驗和研究.同時值得注意的是，已有的數據融合方法多基于TSX/TDX Stripmap模式下的數據，其空間分辨率為HRTI3標準下的12 m(Rossi et al., 2012; Deo et al., 2013; Brautigam et al., 2012).而這對于城區(qū)的地形信息解讀還遠遠不夠，因此如何在僅用Stripmap格式下的數據獲取更高分辨率(原始分辨率，高達2 m)和高精度的DEM還需要進行系統(tǒng)的研究和嘗試.

本文主要針對目前TSX/TDX干涉生成中存在的DEM融合，精度評定問題和缺陷進行系統(tǒng)研究，提出一種新的適用于大范圍DEM制作的升降軌融合算法和流程對雷達側視成像中的幾何畸變處的信息進行相應的恢復，以及通過不同地物特征分析DEM的精度水平，并定量分析了該DEM對DInSAR相位干涉的改進效果.具體地，我們選擇我國珠海和澳門為試驗區(qū)域，利用TSX/TDX得到的DEM與外部Lidar獲得的高精度參考DEM進行對比，并對其精度和潛在制約因素進行分析，為TanDEM-X全球覆蓋DEM的制作和應用提供參考.同時，將迭代獲取的DEM應用于高分影像的DInSAR技術，利于后續(xù)的相位解纏和形變信息的解譯，為更高精度的時序InSAR形變監(jiān)測提供有利條件.

2 TSX/TDX提取DEM原理及誤差分析

2.1 TSX/TDX雷達衛(wèi)星系統(tǒng)

TanDEM-X是德國宇航局于2010年發(fā)射的，是TerraSAR-X的姊妹衛(wèi)星.TSX/TDX采用HELIX衛(wèi)星運行方式并進行編隊飛行形成一個單軌雙天線系統(tǒng)對地球進行觀測從而得到3D地表數字高程模型(Krieger et al., 2007).該系統(tǒng)下兩顆衛(wèi)星相隔很近(<400 m)，克服了重復軌道干涉時的時間去相關，同時TSX/TDX的基線直接由星載GPS測量得到，且可以實時靈活地調整和量測(Krieger et al., 2007)，其精度較常規(guī)SAR衛(wèi)星(厘米級)的基線來說更為精確，達到毫米級.此系統(tǒng)下可以進行三種模式的觀測： 1) 靜態(tài)模式(monostatic mode)，TSX和 TDX互相獨立； 2) 收發(fā)分置模式(bistatic模式)，即一發(fā)雙收； 3) 交互收發(fā)分置模式(alternating monostatic模式)，通過每次脈沖后改變發(fā)射器來最小化時間去相關.此外，除了交叉軌道下的地形監(jiān)測，TSX/TDX的收發(fā)配置模式下還能進行順軌干涉測量，可以用來識別地表運動的物體，比如交通、艦船監(jiān)測應用，海洋學中海浪波譜分析等.該雙星系統(tǒng)獲取的DEM精度較以往任何機載、星載獲取的DEM都要高很多，其垂直方向的相對精度在局部坡度角<20%時候達到2 m，>20%時達到4 m，絕對精度也優(yōu)于10 m，為全球DEM測量開創(chuàng)了新紀元(Krieger et al., 2007; Tachikawa et al., 2011).

2.2 基本原理

本文采用的SAR數據為TDX一發(fā)雙收模式(bistatic)下獲取的，如圖1a，該模式下信號由其中一個衛(wèi)星發(fā)射，兩個衛(wèi)星接受，通過對兩個回波信號相位差進行分析來獲取地形信息.同時，bistatic模式下的SAR影像對的時間基線幾乎為0，能夠最大限度地克服大氣、時間去相關和軌道對干涉相位的影響，利于干涉相位的解譯從而獲得精確的地形信息.圖1b為TSX/TDX雷達衛(wèi)星系統(tǒng)的干涉原理圖.

其中，TSX為TerraSAR-X衛(wèi)星，TDX為TanDEM-X衛(wèi)星，B為兩者之間的基線長度，BX為基線B在水平方向上的分量，BZ為基線B在豎直方向上的分量，P為地表的任意一點，θ為主影像的入射角，α為垂直基線與水平方向的夾角，r1，r2分別為地面點P分別到兩SAR傳感器的距離，h為P點的高程，H為衛(wèi)星到地面的距離.由其幾何關系可得

圖1 TSX/TDX雷達衛(wèi)星干涉測量原理Fig.1 Theory map of interferometry of TSX/TDX

(1)

h=H-r1·cosθ,

(2)

(3)

(4)

其中，λ為波長，由于TSX/TDX的bistatic模式為單軌雙天線系統(tǒng)，K=1.由式(2,4)可知，影響測高精度的主要因素有：衛(wèi)星的斜距、基線長、平臺高度等.

2.3 誤差分析

由于X波段波長較短(3.1 cm)，對于地形起伏較大的山區(qū)及城區(qū)來說，地形相位將在干涉圖上產生較密的條紋，這將會影響到后面的相位解纏甚至造成解纏失敗.因此，引入外部DEM來去掉大部分的地形相位從而舒緩差分干涉圖中的條紋，再通過差分相位進行解纏來恢復殘余的地形信息.本文中用到的TSX/TDX下的bistatic模式的差分干涉相位表達式如下：

φ=φtop_err+φorb+φaps+φnoise,

(5)

主要包含以下幾部分：殘余地形相位φtop_err、軌道殘差相位φorb、大氣相位φaps及噪聲相位φnoise.由于該雙星系統(tǒng)的時間基線幾乎為0，因此可認為主從影像成像時刻的大氣分布情況一致，因此φaps可忽略不計；TSX/TDX系統(tǒng)有著較常規(guī)SAR更精確的軌道信息，同時，我們還可以通過多項式擬合的方式來去除殘余軌道相位(Fengetal., 2012;Zhouetal., 2013;Fengetal., 2015a)，因此軌道殘差相位φorb也可以忽略不計，從而得到僅含有地形信息的干涉相位.由2.2節(jié)可知影響測高精度的主要因素為平臺高度H, 水平基線BX，豎直基線BZ，入射角θ，斜距r，由式(4)可知：

(6)

其中：

(7)

由以上精度分析可知，TanDEM-XInSAR干涉測高的精度與基線成反比，基線越長其測高精度越高；而斜距的精度受到硬件設備的影響，如衛(wèi)星采樣時鐘抖動及衛(wèi)星鐘的不確定性，該項誤差與其產生的高程誤差為同一數量級，對高程精度影響較小，可以忽略不計.

3 升降軌數據迭代融合算法提出和DEM技術流程

傳統(tǒng)DInSAR技術雖然能夠通過干涉相位獲取地面的數字高程模型，但是單一軌道的干涉對容易受到高程模糊度、幾何畸變、陰影和疊影等因素的影響而無法獲取某些區(qū)域的地形信息.本文采用升降軌結合的方式，通過不同飛行方式所具有的不同幾何關系和參數，如入射角、垂直基線、高程模糊度等，來對因單一軌道產生的陰影、疊影等區(qū)域的地形信息進行恢復和補充.目前大部分升降軌融合的方法多考慮與相干性、多視視數相關的干涉相位的概率密度函數，繼而得到最后的權重因子(Jiangetal., 2014;Deoetal., 2014).但是TSX/TDX系統(tǒng)的時間基線幾乎為0，因此其相干性除了在幾何畸變處外，不管是在山區(qū)、還是在平坦地區(qū)均很好，多視視數越大，其相干性與干涉相位的標準差之間的相干性越小(Toughetal., 1995).因此，本實驗中沒有將相干性作為定權因子，而是通過提出一種升降軌數據迭代融合的簡潔算法來實現大范圍DEM數據制作處理.同時，由2.3節(jié)可知，不同垂直基線對應著不同的高程敏感性(高程模糊度)，因此，本文除了考慮疊影、陰影在升降軌中的位置信息外，還引入了垂直基線用來定權，從而得到更加精確的DEM模型.

升降軌數據融合的技術流程如圖2所示，首先通過對單獨的升降軌的CoSSC數據對(CoSSC為DLR通過TanDEM-XProcessor(ITP)數據處理軟件進行了一定的輻射校正、精確配準后的產品(Nestoretal., 2010))進行常規(guī)的干涉處理，由于CoSSC數據已經進行了精確配準，因此只需要選擇一定的多視視數進行干涉得到干涉圖，接著引入外部90m分辨率的SRTM數據來削弱區(qū)域地形引起的密集條紋.此處，由于SRTM所得到的高程為正高(orthometicheight)，因此需要將高程異常(geoidalheight)去除使其轉換至WGS84坐標系下，一般而言在局部區(qū)域這個高程異常影響較小可以忽略.然后，對得到的差分干涉圖進行濾波和相位解纏，得到殘余的地形信息，值得注意的是，由于第一次引入的外部SRTM分辨率較低，從而導致多視視數不能過小(過小的視數會導致外部SRTM無法配準至SAR坐標系下，同時在高樓林立處因外部DEM分辨率不夠而導致相位解纏誤差)，而最后得到的DEM的分辨率也并不理想.為了獲得更為豐富的DEM細節(jié)信息，本文通過迭代的思想將第一次處理得到的最終DEM作為第二次的外部DEM，同時減小多視視數，再重復差分干涉處理從而得到更高分辨率的DEM，該操作能夠很好地緩解因單視下高樓區(qū)的相位解纏誤差.此外，本文對TanDEM-X獲取的DEM中的水系也進行了一定的處理，通過借鑒PS-InSAR選點的思想，將振幅和相干性雙向閾值對水系的輪廓進行一定的判定，記錄其位置信息，并對其統(tǒng)一賦值(通過計算其周邊像素來獲得平均值)達到光滑水面的作用，從而得到更為精確可靠的數字高程模型.另外不同幾何條件下獲取的SAR影像在地理編碼的時候有著一定的偏差，因此，還需要將升降軌得到的高程圖統(tǒng)一配準到相同的坐標系統(tǒng)下.本文通過在反編碼以后的強度圖中選擇一定數目(20個)的GCP控制點，再通過最小二乘的思想對其進行一次多項式擬合，得到影像對之間的偏移量關系式.

圖2 升降軌數據融合流程圖Fig.2 Flowchart of the fusion of ascending and descending data

對配準后的升降軌數據獲取的DEM進行一定的融合處理，本文首先分別對升降軌外部DEM進行疊影、陰影區(qū)域探測：當坡度角為正的時候，如果坡度角大于等于入射角，此時在雷達影像上將產生疊影，同時也需要考慮坡度角小于入射角時被疊影影響的區(qū)域；當坡度角為負的時，如果坡度角的絕對值大于等于入射角時，雷達影像上表現為陰影，同時，當坡度角絕對值小于入射角時，某些區(qū)域也受到陰影的影響.接著分別記錄探測得到的疊影、陰影區(qū)的位置信息為Mi，此處的位置矩陣中僅含有0、1兩類元素，以疊影的情況為例，如果升軌中該像素位于疊影區(qū)域(記為0)而降軌中該對應像素不受陰影或者疊影的影響，此時其位置信息為1.然后，計算升降軌干涉對中每個像素上的垂直基線，記為Bi，取其絕對值.當升降軌影像多余2對時(N對TSX/TDX影像)，該方法仍然可以進行擴展，只要依次探測出各個干涉對中的陰影、疊影位置信息，并記錄，同時仍然按照垂直基線的絕對值來進行定權，得到公式(8)，其適應性和擴展性強：

(8)

其中i代表一個平臺(升降軌、入射角等不同的TSX/TDX像對)的像對，而N為總共的平臺個數，Hi為每個平臺的DEM值，Hfusion為融合后的DEM.最后，本文在i=2的情況下，將垂直基線、疊影及陰影組成的權重因子對TXS/TDX的DEM進行融合來獲取高精度數字高程模型.

4 實驗及分析

4.1 實驗數據處理

本文選取珠海、澳門等地作為實驗區(qū)域，選取了兩個不同時間點的升降軌TSX/TDX影像對作為本文的實驗數據，詳細數據信息如表1所示.干涉對的覆蓋范圍如圖3所示，每景影像的覆蓋面積約為33 km×56 km.該實驗區(qū)域最高海拔約為450 m，高程大于30 m的平均地形坡度為16°，地勢總體而言較為平坦，但是地貌和地物特征較為復雜，除了植被茂盛的山區(qū)外，還有較多的農田、水系及高樓林立的城市，因而有利于分析在不同地形和地物特征下TanDEM-X的數字高程模型的精度水平.為了驗證本文算法獲得的DEM的精度和可靠性，還獲得了覆蓋澳門及部分珠海區(qū)域的高精度Lidar數據得到的數字高程模型，其空間分辨率為5 m, 垂直精度高達1 m.

圖3 研究區(qū)域及影像覆蓋黑色實線框分別代表了兩個影像時間點的TSX/TDX干涉對，其中一對升軌數據，兩對降軌數據；藍色實線框代表了本次選取的升降軌重復區(qū)域.Fig.3 The research area and coverage of SAR imagesThe black solid rectangles represent three pairs of TSX/TDX, including an ascending pair and two descending pairs; the blue solid rectangles represent the overlap area of ascending and descending which will be used in the following parts.

表1 本文選取的TerraSAR-X/TanDEM-X數據參數

現以其中一對CoSSC數據為例，介紹本實驗中的計算參數和流程：首先對已經配準好的數據進行常規(guī)的差分干涉處理.為了抑制相位噪聲，對干涉圖進行方位向距離向4×4的多視處理(多視后分辨率約10 m)，再引入外部SRTM(分辨率為90 m)進行差分處理，得到去除基本地形的差分干涉圖.接著，根據生成的差分干涉圖進行相應的判斷，如果信噪比較高的話，可以直接進行相位解纏工作，如果信噪比偏低，建議進行適當濾波處理；然后通過人工屏蔽的方法避免或者減弱低相干區(qū)域(海域及孤立的島嶼)對相位解纏的影響，使其不參與解纏計算，對濾波后的差分干涉圖中相干性高于閾值(本文選擇0.45)的點采用最小費用流(MCF)方法進行解纏，此處，為了后續(xù)將相對相位轉為絕對相位，本文相位解纏的時候選擇的參考對象是整景影像.此外，本文還通過對幅度、相干性設置一定的閾值來獲取影像中水系的輪廓.最后，將絕對相位轉化為高程信息，并反編碼到地球坐標系下，再與重采樣后的外部SRTM的高程信息相加得到最后的高程信息，此時的DEM的空間分辨率約為10 m.由于本次所收集的Lidar數據的空間分辨率為5 m，因此還需要降低干涉多視視數來獲得更高分辨率的DEM.在第一次獲取的DEM的基礎上，進行方位向距離向2×2的多視，從而得到最終升降軌地距向方位向分辨率3.5 m×4.4 m的高分DEM.當升降軌數據都單獨按照流程圖2處理完后，對升降軌進行配準.由于本次獲取的DEM與收集的Lidar數據分辨率不一致，為了進行更加合理的對比分析，本文采用一定的插值方法(此處采用了雙線性插值)將升降軌TanDEM-X的DEM重采樣至Lidar的空間分辨率，并采用第三節(jié)中提出的融合方法進行融合得到最終的高精度的融合DEM產品.

4.2 實驗結果及DEM精度分析

圖4為我們采用第3節(jié)所述的方法將升降軌重復覆蓋區(qū)域(如圖3中藍色方框所示)進行相應的融合得到的DEM，其中圖4a為升軌TSX/TDX得到的DEM，圖4b為降軌得到的DEM，圖4c為升降軌融合得到的DEM，圖4(d,e,f)為黑色框內的放大圖，從圖中可以看出單一軌道獲得的DEM還存在一定的空值(疊影、陰影區(qū)及低相干區(qū))，但不同幾何條件下的空值位置基本不同，因此通過融合后的DEM能夠有效地恢復出幾何畸變處的信息，如圖5f所示.為了能夠對TanDEM-X的DEM進行一定的精度分析，本文選取了含有高精度Lidar數據覆蓋的澳門及周邊區(qū)，Lidar數據的具體情況如圖5.從圖中，我們可以清晰地看到起伏的山區(qū)、密集的澳門城區(qū)及新區(qū)，圖中的兩個矩形框分別代表了不同的地物，其中藍色代表了山區(qū)、位于珠海境內，如圖5a；紅色代表澳門的老城區(qū),如圖5b，老城區(qū)中分布著密集的居民樓、辦公樓，有的樓高達230多米.

圖6a為融合后的DEM與Lidar的對比圖，從左圖中可以看出圖中有三處光滑區(qū)，均為水域，由于其相干性較低，通過振幅、相干性閾值探測出其位置信息，并賦于周圍像素平均值，同時也不參與后面的結果驗證.圖中的細節(jié)信息保留的較好，如山區(qū)的棱角較分明，陰影、疊影等地方的細節(jié)信息得到較好地恢復，因此，TanDEM-X獲得的DEM在城區(qū)的信息較為豐富，如右上角處的民房及城鎮(zhèn).此外，由于Lidar數據的獲取時間為2010年，在TSX/TDX獲取時間之前，因此其地物之間會有一定的差異，且其在除澳門以外的城區(qū)均無觀測數據.再者，由于升降軌的TanDEM-X數據間相隔較長，約兩年，這也會使得地物之間存在較大的變化，如植被的生長、房屋的改建等，因此融合后的DEM也會存在一些因時間跨度較大導致的差異.為了定量的進行對比，本文選取了兩個剖面，如圖中的黑色實線m，n，其剖面圖對應著圖6(e, f)，圖中升軌、降軌、Lidar及融合的DEM分別用紅、綠、藍及黑色實線表示，從圖6e中的小圖中可以明顯看到，在陰影或者疊影的地方，融合后的DEM明顯優(yōu)于升降軌單獨的DEM值，更接近Lidar數據.實線n對應的剖面圖如圖6f所示，在海拔較高的地方TanDEM-X的DEM與Lidar符合的較好，但是尾部，由于Lidar影像中平地沒有進行監(jiān)測，因而TanDEM-X的DEM較Lidar值偏高，為房屋區(qū).隨后，本文也對山區(qū)的DEM殘差進行了統(tǒng)計分析，如圖6i所示，其中紅色、綠色及藍色分別為為升軌、降軌及融合后的DEM與Lidar數據差值的統(tǒng)計，服從正態(tài)分布.由圖中可以看到融合后的DEM的均值及均方差均有明顯改善：均值由原來的1.75 m(升軌)、1.26 m(降軌)降到1.12 m，均方差由最初的3.88(升軌)、4.56(降軌)減少到3.57，具體見表2.

圖6(c, d)為澳門老城區(qū)的融合DEM及Lidar對比圖，圖中的綠色平滑區(qū)域為水域，通過以上的方法進行識別，并賦予統(tǒng)一數值.從圖中可以看出TanDEM-X的DEM結果較粗糙，尤其是在高樓林立的城區(qū)，由于雷達的側視成像而產生很多陰影、疊影等問題在高樓的地方表現的較為明顯，因此TanDEM-X獲取的DEM在高樓分布較為密集的地方得到的結果較差.然而，在較矮(高程低于30 m)的居民區(qū)TanDEM-X獲得的DEM的高程精度較高，能夠較好地恢復出房屋的高程.此外，本文還對澳門的新城區(qū)進行了一定的實驗，實驗結果表明對于樓房分布較為稀疏的地區(qū)，TanDEM的DEM能夠較Lidar數據更好地恢復出建筑物的外形及高度，去掉三倍中誤差以外的統(tǒng)計值后其標準差基本上在5 m之內.同樣，選擇兩條剖面q，p，從兩條剖面圖(圖6g,h)我們可以看出高程低于30 m的城區(qū)監(jiān)測精度與Lidar數據較符合，但是對于高樓區(qū)，大多會產生低估的現象，該現象主要是由于城區(qū)高樓之間的陰影、疊影以及相位中心的影響導致相位解纏的時候出現問題，尤其是高樓與平地相接的地方，其相位差可能超過了2π，因而導致解纏錯誤.通過統(tǒng)計其與Lidar的殘差圖時(圖6j)，無論是單一平臺還是融合后的DEM，其標準差均較山區(qū)來說大，但是，我們還是可以看到融合后的DEM無論是在殘差均值還是標準差方面都有很大的改進，具體統(tǒng)計數值見表2.

表2 高程殘差統(tǒng)計對比Table 2 Comparison of height difference

4.3 高精度DEM對于差分干涉的改進

除了對TanDEM獲取的DEM進行外部DEM定量比較外，本文還對其進行差分干涉效果的分析.選用一對覆蓋澳門、珠海等區(qū)域的TerraSAR影像對，主從影像的獲取時間分別為2011年12月23日和2012年01月03日，其時間基線為11天，垂直基線為135 m，主影像的入射角為41°.我們選取與升降軌TanDEM所獲取DEM重合的區(qū)域進行實驗，將融合的DEM作為TerraSAR干涉對的外部DEM.為了能夠比較分析TanDEM獲取的融合DEM對差分干涉的改進，本文還選取了其他兩類外部DEM(SRTM和ASTER)進行對比實驗.為了使得SRTM能夠較好地配準至SAR影像坐標下，本文先選用了方位向距離向為5×5的多視進行常規(guī)的DInSAR處理，差分干涉圖的方位向距離向分辨率分別為9.47 m×10.39 m，其結果如圖7(a—c)所示.從圖中的黃色實線圈內可以看到，對于山區(qū)來說，TanDEM作為外部DEM時，其在差分干涉圖中殘存的條紋最少，明顯優(yōu)于另外兩種外部DEM.為了能夠獲取TanDEM對于城區(qū)的差分干涉的改進，本文選擇了覆蓋澳門新城區(qū)的一小部分影像在5×5視數的基礎上進行方位向距離向2×2的多視迭代，此時的空間分辨率高達3.79 m×4.12 m.同樣采用三種不同的外部DEM，其結果如圖7(d—f)，從圖中可以看出，對于分布較為密集的高樓區(qū)，由于TanDEM對于高于30m的建筑物高程恢復能力較弱，因此還殘留較多的條紋，如連接老城區(qū)的高架橋的周邊，但SRTM和ASTER殘留著較TanDEM更密的條紋；對于分布較稀疏的高層建筑物，如“威尼斯人度假村”及周邊的幾個高層建筑，TanDEM能夠對其進行改善，從而得到精確去掉地形信息的差分相位，更有利于形變信息的提取(如圖7f所示).而SRTM及ASTER作為外部DEM的差分干涉圖上仍然殘留較多地形誤差，其可能原因是由于分辨率過低及配準過程存在的偏移誤差.同樣，對于2×2視下的山區(qū)，如圖7(d—f)中黃圈內的山區(qū)，TanDEM較其他兩者的有著明顯的改進，尤其對于陰影、疊影等細節(jié)地方的地形也去的較為干凈.值得注意的是，本次實驗采用的干涉對的垂直基線僅為100多米，如果對于垂直基線較長的干涉對其改進效果將更加顯著.

圖6 TSX/TDX獲取的DEM與Lidar的對比圖(a—b)選取的山區(qū)的融合DEM及對應的Lidar影像；(c—d)選取的城區(qū)的融合DEM及對應的Lidar影像；(e—f)山區(qū)內m,n斷面圖；(g—h)城區(qū)內q，p斷面圖；(i—j)山區(qū)及城區(qū)的殘差統(tǒng)計圖.Fig.6 Comparison between DEM and Lidar data(a—b) Fusion DEM and Lidar data in mountainous area; (c—d) Fusion DEM and Lidar data in urban areas; (e—f) Comparison of two selected profiles named m and n in mountainous areas; (g—h) Comparison of two selected profiles named p and q in urban areas;( i—j) The height difference map of mountainous and urban areas.

為了能夠定量的分析外部DEM在差分干涉圖中的殘余量，我們分別提取了圖7(a—c)中的山區(qū)及城區(qū)范圍，并統(tǒng)計這些范圍內相干性大于0.8的點的相位，其結果如圖8所示，其中橫軸代表這些點的纏繞相位值，單位為rad，縱坐標表示這些參與統(tǒng)計的點的概率密度.由于本次采用的干涉對的時間基線僅為11天，因此可以假設這之間不存在形變信息，因而差分干涉圖中所剩下的僅為殘余地形、大氣及噪聲.因為我們所圈取的范圍對于三類情況都一樣，其大氣和噪聲的影響都是相同的.因此，我們可以通過其纏繞相位的標準差來判斷殘余地形相位的多少.圖8a為山區(qū)范圍內的纏繞相位統(tǒng)計結果，其中紅色實線為SRTM作為外部 DEM改正后的差分干涉相位的分布，綠色實線和藍色實線分別代表外部DEM為ASTER及TanDEM的相位統(tǒng)計圖，其纏繞相位的標準差依次為：1.4 rad, 1.6 rad及1.4 rad.可以看出外部DEM對差分干涉圖中山區(qū)的改善程度由高至低依次為：TanDEM>SRTM>ASTER.同樣，我們也分別圈取了城區(qū)，并做同樣的統(tǒng)計，其纏繞相位統(tǒng)計結果如圖8b所示.當ASTER作為外部DEM時，其纏繞相位并不穩(wěn)定，如圖8b中的綠色實線的左側有著較大的突變，而SRTM和TanDEM作為外部DEM時其相位統(tǒng)計分布較為正常，其纏繞相位標準差依次為：1.7 rad，1.3 rad及1.1 rad.其結果與山區(qū)范圍內的一樣，TanDEM作為外部DEM對差分干涉相位的改善最大，其次為SRTM，ASTER次之.

圖7 TanDEM-X的高分DEM對差分干涉圖的改進(a—c)分別以SRTM、ASTER及TSX/TDX融合DEM為外部DEM得到的多視視數為4×4的差分干涉圖；(d—f)多視視數為2×2的差分干涉圖.Fig.7 Improvement of differential interferograms with TanDEM-X(a—c) Differential interferograms of 4×4 multilook number in azimuth and range direction respectively with SRTM, ASTER and TanDEM; (d—f) Differential interferograms with 2×2 multilook number in azimuth and range directions respectively.

圖8 各外部DEM作用下的差分干涉圖的纏繞相位統(tǒng)計圖紅色實線代表SRTM為外部DEM，綠色和藍色實線分別代表外部DEM為ASTER和TanDEM的相位統(tǒng)計結果.(a)山區(qū)中相干性高于0.8的點的相位統(tǒng)計結果；(b)城區(qū)中相干性高于0.8的點的相位統(tǒng)計結果.Fig.8 Statistical graph of wrapped phase in differential interferogrames with different external DEMThe red solid line represent the statistical results of SRTM and the green and blue solid lines represent the results of ASTER and TanDEM respectively. (a) The statistical results of mountainous areas with coherence over 0.8; (b) The statistical results of urban areas with coherence over 0.8.

綜上所述，通過本文提出的融合升降軌TanDEM數據算法獲取的高分DEM能夠較好地恢復出因幾何畸變(如陰影、疊影等)而損失的高程信息，同時較SRTM、ASTER能獲得更為豐富的細節(jié)信息.在高樓分布較為稀疏的新城區(qū)TanDEM的DEM精度很高，但是對于高度過高的建筑物(遠遠超過高程模糊度)仍會存在一定的低估現象.而在高樓林立的老城區(qū)，由于高建筑物分布過于密集，幾何畸變較為復雜，易產生較為密集的條紋，這也極大地阻礙了后續(xù)的相位解纏，從而導致城區(qū)的高程反演低估甚至失敗.此外，通過TanDEM獲取的DEM對DInSAR的改進較SRTM、ASTER等外部DEM更為明顯，尤其對于山區(qū)，能夠較好地去除地形信息，方便后續(xù)的形變解譯；同時，在分布較為稀疏的新城區(qū)或者開發(fā)區(qū)，TanDEM作為外部DEM去除地形的效果很好，能夠利于后續(xù)的DInSAR差分相位質量的提高，并為時序InSAR提供更為有利的條件.總的來說，TSX/TDX獲得的DEM不管是空間分辨率還是精度都大大優(yōu)于目前我們常用的SRTM和ASTER，但是城區(qū)的細節(jié)信息監(jiān)測仍有待提高.因此高分的TanDEM數據及小視數甚至單視的干涉處理將成為后續(xù)的研究重點，通過更高分辨率的數據獲取城區(qū)精度更高的DEM，同時也為后續(xù)的時序InSAR形變監(jiān)測提供高分辨率外部DEM輔助.

5 結論

本文通過獲取的兩對覆蓋澳門、珠海區(qū)域的TSX/TDX干涉對進行了一系列的干涉處理，成功獲取了每個干涉對的高精度DEM，同時，通過本文提出的一種結合陰影、疊影、垂直基線的加權融合迭代方法在恢復因側視成像而導致的幾何畸變處的信息的同時得到較HRTI3更高空間分辨率的DEM.實驗結果表明幾何畸變恢復效果較為突出，通過引入的高精度的Lidar數據與TSX/TDX獲得的DEM針對不同地物的精度進行一系列的對比分析，我們發(fā)現：針對山區(qū)來說，TSX/TDX獲取的DEM精度很高，且通過升降軌融合后的DEM精度較單獨軌道的精度要高；針對建筑物分布稀疏的城鎮(zhèn)，其精度也較高；但對于高樓林立的城區(qū)來說，由于建筑物高差大且分布過于密集，從而產生較為嚴重的幾何畸變和解纏誤差，導致部分高建筑物的高度產生低估甚至誤估，其精度較前兩者偏低.與此同時，本文還將TanDEM-X獲取的DEM作為TerraSAR的外部DEM來定量分析，通過統(tǒng)計差分干涉圖中纏繞相位的標準差來判斷其對常規(guī)DInSAR的改善效果，結果明顯優(yōu)于SRTM、ASTER.因此本文的研究將為大范圍高分辨率、高精度DEM的制作起到參考作用，同時也為今后InSAR的形變監(jiān)測起到促進作用.

致謝 本研究所用TerraSAR-X/TanDEM-X數據來源于德國宇航局(DLR)項目(XTI_GEOL6191和LAN0820).感謝澳門特別行政區(qū)地圖繪制地籍局提供的Lidar數據，文中部分圖件使用GMT4.5軟件繪制而成.

Ashutosh B, Rajat S C, Kamal J. 2013. Assimilation of DEMs generated from optical stereo and InSAR pair through data fusion.ScienceResearch, 1(3):39-44.

Bhar dwaj A, Chatterjee R S, Jain K. 2013. Assimilation of DEMs generated from optical stereo and InSAR pair through data fusion.ScienceResearch, 1(3): 39-44.

Brautigam B, Rizzoli P, Martone M, et al. 2012. InSAR and DEM quality monitoring of TanDEM-X. ∥ 2012 IEEE International Geoscience and Remote Sensing Symposium (IGARSS). Munich: IEEE, 5570-5573.

Deo R, Manickam S, Rao Y S, et al. 2013. Evaluation of interferometric SAR DEMs generated using TanDEM-X data. ∥ 2013 IEEE International Geoscience and Remote Sensing Symposium (IGARSS). Melbourne, VIC: IEEE, 2079-2082.

Deo R, Rossi C, Eineder M, et al. 2014. Fusion of ascending and descending pass raw TanDEM-X DEM. ∥ 2014 IEEE International Geoscience and Remote Sensing Symposium (IGARSS). Quebec City, QC: IEEE, 21-24.Farr T G, Kobrick M. 2000. Shuttle radar topography mission produces a wealth of data.Eos,TransactionsAmericanGeophysicalUnion, 81(48): 583-585.Feng G C, Ding X L, Li Z W, et al. 2012. Calibration of an InSAR-derived coseimic deformation map associated with the 2011MW9.0 Tohoku-Oki Earthquake.IEEEGeoscienceandRemoteSensingLetters, 9(2): 302-306.

Feng G C, Li Z W, Shan X J, et al. 2015a. Source parameters of the 2014MW6.1 South Napa earthquake estimated from the Sentinel 1A, COSMO-SkyMed and GPS data.Tectonophysics, 665: 139-146, doi: 10.1016/j.tecto.2015.05.018.

Feng G C, Li Z W, Shan X J. 2015b. Geodetic Model of the April 25, 2015MW7.8 Gorkha Nepal Earthquake andMW7.3 Aftershock Estimated from InSAR and GPS Data.GeophysicalJournalInternational, 203(2)：896-900, doi: 10.1093/gji/ggv335.

Feng G C, Xu B, Shan X J, et al. 2015. Coseismic deformation and source parameters of the 24 September 2013 Awaran, PakistanMW7.7 Earthquake derived from optical Landsat 8 satellite images.ChineseJournalofGeophysics(in Chinese), 58(5): 1634-1644, doi: 10.6038/cjg20150515. Ferretti A, Fumagalli A, Novali F, et al. 2012. DEM reconstruction with SqueeSAR. ∥ 2012 Tyrrhenian Workshop on Advances in Radar and Remote Sensing (TyWRRS). Naples: IEEE, 198-201.

Jiang H J, Zhang L, Wang Y, et al. 2014. Fusion of high-resolution DEMs derived from COSMO-SkyMed and TerraSAR-X InSAR datasets.JournalofGeodesy, 88(6): 587-599.

Krieger G, Moreira A, Fiedler H, et al. 2007. TanDEM-X: A Satellite Formation for High-Resolution SAR Interferometry.IEEETransactionsonGeoscienceandRemoteSensing, 45(11): 3317-33441.

Liao M S, Lin H. 2003. Synthetic Aperture Radar Interferometry: Principle and Signal Processing (in Chinese). Beijing: Surveying and Mapping Press.Nestor Y M, Rossi C, Lachaise M, et al. 2010. Interferometric processing algorithms of TanDEM-X data. IEEE International Geoscience and Remote Sensing Symposium (IGARSS), 3518-3521.

Rossi C, Frita T, Eineder M, et al. 2012. Towards an urban DEM generation with satellite SAR interferometry. ∥ ISPRS International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, XXXIX-B7, 73-78.

Rossi C, Eineder M, Fritz T, et al. 2013. Quality assessment of TanDEM-X raw DEMs oriented to a fusion with CartoSAT-1 DEMs. ∥ EARSeL(Congress), 1-9.

Tachikawa T, Kaku M, Iwasaki A, et al. 2011. ASTER Global Digital Elevation Model Version 2-Summary of Validation Results August 31, 2011. Tough R J A, Blacknell D, Quegan S. 1995. A statistical description of polarimetric and interferometric synthetic aperture radar data. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 449(1937): 567-589.Walker W S, Kellndorfer J M, Pierce L E. 2007. Quality assessment of SRTM C- and X-band interferometric data: Implications for the retrieval of vegetation canopy height.RemoteSensingofEnvironment, 106(4): 428-448.

Xiao J Q, Li Z W, Wang C C, et al. 2012. DEM extraction over mountainous area and its accuracy analysis with phase compensation based differential interferometry.GeomaticsandInformationScienceofWuhanUniversity(in Chinese), 37(3): 334-338.

Zhou H, Feng G C, Li Z W, et al. 2013. The fault slip distribution

of the MyanmarMW6.8 earthquake inferred from InSAR measurements.ChineseJournalofGeophysics(in Chinese),56(9): 3011-3021, doi: 10.6038/cjg20130914.

附中文參考文獻

馮光財, 許兵, 單新建等. 2015. 基于Landsat8光學影像的巴基斯坦AwaranMW7.7地震形變監(jiān)測及參數反演研究. 地球物理學報, 58(5): 1634-1644, doi: 10.6038/cjg20150515.

廖明生, 林輝. 2003. 雷達干涉測量—原理與信號處理基礎. 北京: 測繪出版社.

肖金群, 李志偉, 汪長城等. 2012. 相位補償SAR差分干涉提取山區(qū)DEM算法及精度分析. 武漢大學學報(信息科學版), 37(3): 334-338.

周輝, 馮光財, 李志偉等. 2013. 利用InSAR資料反演緬甸MW6.8地震斷層滑動分布. 地球物理學報, 56(9): 3011-3021, doi: 10.6038/cjg20130914.

(本文編輯 汪海英)

Generation of high precision DEM from TerraSAR-X/TanDEM-X

DU Ya-Nan, FENG Guang-Cai*, LI Zhi-Wei, ZHU Jian-Jun, PENG Xing

SchoolofGeoscienceandInfo-PhysicsEngineering,CentralSouthUniversity,Changsha410083,China

Digital Elevation Model (DEM) derived from TerraSAR-X/TanDEM-X, can overcome the influence of atmospheric delay and surface deformation, and has a much higher resolution and precision compared with traditional InSAR. However, the DEM generated from TSX/TDX is vulnerable to geometric distortion such as layover, shadow, foreshortening etc. Moreover, the precision of TSX/TDX DEM, especially for different surface characteristic, such as mountainous areas, urban buildings and flat areas is still far from clear in application. Therefore, it has a practical significance on the guidance for DEM production derived from TSX/TDX.In this paper, we proposed a procedure which combines both ascending and descending measurements to obtain a high precision and resolution DEM product from the bistatic mode of TSX/TDX. We also use an iterative fusion method considering the perpendicular baseline and the location of layover and shadow as weighing. Because the ascending and descending SAR images have different geometrical conditions, such as incident angle and perpendicular baseline, it can regain the DEM affected by geometric distortion. We also adopted a quality assessment of the fusion DEM (5 meters resolution) for two terrain configurations to demonstrate the precision pattern. The fusion DEM in different kinds of terrain configurations: moderate topography and urban areas with dense high-rising buildings, is compared with DEM collected by Lidar technique in the study area. In addition, we qualitatively and quantitatively analyzed the improvement of fusion DEM in the processing of DInSAR technology compared with the existing external DEM.Two pairs of TSX/TDX data, acquired on 21 November 2011 in the ascending mode and 24 October 2013 in the descending mode, respectively, are used to map a fusion DEM in Zhuhai city and Macau. The results indicate that the fusion DEM with the proposed method can recover more information in the distorted areas compared with the DEM derived from a single platform. In order to validate the precision of the fusion DEM, a Lidar DEM with high resolution and precision of the research area are collected and compared. We find that the standard deviation (RMS) in mountainous areas is 3.57 m with a decrease of 8% and 22% for ascending and descending data, respectively. The standard deviation value in urban areas is 11.85 m which also has a drop of 28% and 22% compared with two single platform. Moreover, the precision of the fusion DEM is high in residential areas with sparse buildings with a standard deviation 5 m. In addition, the TerraSAR-X pairs whose temporal baseline is 11days and perpendicular baseline is 135 m are used to conduct comparison of different kinds of external DEM in DInSAR processing. Moreover, statistical analysis of differential phase (pixel whose coherence is higher than 0.8) in mountainous and urban areas are performed separately for a quantitative analysis. The results show that the fusion DEM used as an external DEM can remove topography phase better compared with SRTM and ASTER data especially in the mountainous area. The standard deviations of differential phases is 1.4 rad, 1.6 rad and 1.4 rad for SRTM, ASTER and TanDEM, respectively. In urban areas with dense high buildings, the standard deviations are 1.7 rad 1.3 rad and 1.1 rad respectively, which means higher resolution and precision of external DEM can improve the quality of DInSAR with high resolution SAR images.This work adopted an iterative fusion method to obtain high resolution and precision DEM. Locations of geometric distortion and perpendicular baseline are regarded as weighting factors to improve the efficiency of fusion procedure by comparing with traditional DEM. Besides, we also qualitatively and quantitatively analyzed the improvement of TSX/TDX DEM in traditional DInSAR, which can help us interpret deformation information and obtain a higher precision of time- series of InSAR.

Digital Elevation； Model； TerraSAR-X/TanDEM-X; Geometric distortion; Fusion of ascending and descending

杜亞男， 馮光財， 李志偉等. 2015. TerraSAR-X/TanDEM-X獲取高精度數字高程模型技術研究.地球物理學報，58(9):3089-3102,

10.6038/cjg20150907.

Du Y N, Feng G C, Li Z W, et al. 2015. Generation of high precision DEM from TerraSAR-X/TanDEM-X.ChineseJ.Geophys. (in Chinese),58(9):3089-3102,doi:10.6038/cjg20150907.

10.6038/cjg20150907

P223

2014-08-14，2015-09-08收修定稿

國家自然科學基金(41104003, 41222027，41574005)的資助.

杜亞男，女，1988年生，博士研究生，主要從事時序InSAR監(jiān)測地表形變.

*通訊作者 馮光財，男，講師.E-mail: fredgps@csu.edu.cn