TerraSAR-X/TanDEM-X获取高精度数字高程模型技术研究

杜亚男， 冯光财， 李志伟， 朱建军， 彭星

中南大学地球科学与信息物理工程学院雷达遥感研究室， 长沙 410083



TerraSAR-X/TanDEM-X获取高精度数字高程模型技术研究

杜亚男， 冯光财*， 李志伟， 朱建军， 彭星

中南大学地球科学与信息物理工程学院雷达遥感研究室， 长沙 410083

以双星系统(TerraSAR-X/TanDEM-X)下的bistatic数据模式为例研究了差分干涉获取高精度DEM产品的融合算法和技术流程.针对不同观测几何条件下(升降轨，不同入射角)观测数据的畸变和缺失，提出一种迭代的顾及垂直基线、阴影和叠影的数据融合新方法重建高分辨率高精度的数字高程模型，并对TanDEM-X融合的DEM在不同地物属性特征下(山区及高楼林立的城区)的精度进行定量分析.本文采用了两对覆盖珠海、澳门区域的升降轨TerraSAR-X/TanDEM-X干涉对进行融合处理，并通过收集的高精度Lidar数据进行精度比较分析.此外，本文还定量分析了TanDEM-X的DEM对常规DInSAR技术的改进，并与SRTM、ASTER的结果进行对比.结果表明：提出的升降轨融合方法较单一轨道平台能够较好地改正或减弱由于几何畸变引起的高程信息缺失或错误，通过与Lidar数据的对比发现TanDEM-X的融合DEM在山区的精度较高，其残差的标准差为3.5 m，较单一轨道(升、降轨)分别降低8%和22%，能够通过迭代的方法获取高分辨率(可达2～5 m)、高精度的地形信息；而针对城区密集建筑物的复杂地形来说，融合的DEM的精度稍低，残差的标准差为11.8 m，但较单一轨道(升、降轨)来说有较大改进，其残差标准差分别降低了28%和22%；而在分布较为稀疏的居民区，融合的DEM能够较Lidar数据获取更好的建筑物高度及外形信息，此时的残差标准差可达5 m.同时，TanDEM-X的融合DEM作为外部DEM能够较SRTM和ASTER来说更好地去除地形信息，尤其在山区及高程建筑密集分布的城区，从而利于后续的相位解缠和形变信息的精确获取和解译，为更高精度的时序InSAR形变监测提供有利条件.

数字高程模型； TerraSAR-X/TanDEM-X； 几何畸变； 升降轨融合

1 引言

数字高程模型(DEM)是用一组有序数值阵列形式表示地面高程的一种实体地面模型，它是地学领域用来表示地形特征的重要表达形式.它作为一种基础的地理信息数据已经在我们日常生活以及国防领域有着非常广泛的应用，如地形制导的军事打击、差分干涉雷达形变监测、地质勘测、各种线路(高铁、城际、高速公路、电力线路等)的规划、城区地籍变化监测等.然而目前能覆盖全球的DEM数据还非常有限，一种是SRTM，2000年由美国太空总署(NASA)和国防部国家测绘局(NIMA)联合发射的“奋进”号航天飞机测量得到的，覆盖了全球陆地面积80%以上，在美国区域公开的DEM分辨率为30 m，全球其他区域公开的DEM分辨率为90 m (Farr & Kobrick, 2000; Feng et al., 2015a)；另外2011年10月美国空间局(NASA)和日本经济产业省(METI)发布了ASTER GDEM V2.0，它覆盖了全球陆地面积的99%以上范围，其空间分辨率可达30 m，虽然其在空间分辨率和高程精度方面都较V1.0有较大提升，但是精度仍然非常有限(Tachikawa et al., 2011; Feng et al., 2015b).当然，除了以上两类全球覆盖DEM外，德国宇航局(DLR)在2011年也公开了其X波段的SRTM，其分辨率较前两类来说有着进一步的提高(约25 m)，但局限于网状不规则覆盖(Walker et al., 2007).值得注意的是SRTM和ASTER都是基于十多年前的影像数据为基础生成的DEM，因此很多区域现实性都较差，无法反映最新的地形特征.

虽然通过航空摄影测量、机载Lidar、地面GPS和全站仪等方法也可以获取局部的高分辨率数字高程模型，但是这些方法不仅耗时耗力而且容易受到天气、环境等众多因素的干扰，因此本文重点关注利用星载InSAR技术获取大范围DEM信息.为了获取更高分辨率、更高精度和更新的DEM，国内外InSAR学者进行了众多有益的尝试，并提出了许多算法和技术流程.首先，廖明生和林辉(2003)、肖金群等(2012)都利用欧空局ERS1/2串行飞行克服时间去相关的影响(时间间隔1天)，基于传统SRTM的辅助，成功获取了研究区域内更高分辨率和更高精度的DEM；另外Ferretti等(2012)也尝试利用时序InSAR的思想来克服时间去相关和大气延迟等误差的影响，从而重建高精度的DEM；近年来随着SAR系统硬件的提升，SAR影像的空间分辨率越来越高(1 m左右)和重返周期也越来越小(3天)，因此通过融合不同传感器、升降轨和不同入射角的DEM测量数据也成为了一种生成高分辨率和高精度DEM的技术手段(Jiang et al., 2014; Ashutosh et al., 2013; Feng et al., 2015c).虽然以上各种方法均能在一定条件上获取较好的DEM，但是其监测方式都为重复轨道干涉测量，仍然容易受到时间、空间上的去相关以及大气延迟和轨道误差的影响.此外，受到各SAR系统参数，性能以及覆盖范围和影像数量等多方面制约，以上各种方法并不能保证在任何区域均适用.鉴于此，德国宇航局(DLR)2010年发射TanDEM-X卫星与之前2007年发射的TerraSAR-X卫星进行编队飞行，组成了一个双星分布式SAR系统.目的是通过这种双星编队的组合方式获取覆盖全球的高精度数字高程模型(Krieger et al., 2007)，从而为我们提供了一个崭新的、更加广阔的、可靠的DEM制作研究平台.

TerraSAR-X/TanDEM-X(TSX/TDX)相比于传统InSAR系统最主要特点就是能够获取“0”时间基线和高精度轨道的干涉对，能够很好地克服各种失相关、大气延迟和轨道误差等引起的相位噪声(Krieger et al., 2007).经过近几年的发展，围绕TanDEM-X和TerraSAR-X获取高精度的DEM也开展了一些研究工作，主要分为两方面：1)TSX/TDX数据处理与DEM融合.现在常用的DEM融合算法主要针对升降轨、不同传感器InSAR和光学DEM等数据融合，而融合的关键便是定权，主要是通过与相干性、多视视数相关的概率密度函数积分得到的干涉相位标准差以及高程模糊度来进行定权(Jiang et al., 2014; Rossi et al., 2013a; Deo et al., 2014)，它是传统星载重复轨道InSAR系统生成DEM技术的延续，然而该算法随着数据源的增多会急剧增加融合的计算量，不利于大范围的DEM制作.更重要的是由于TSX/TDX系统的特殊性，其时间去相关、大气噪声等因素影响几乎可以忽略不计，因此在该情况下再使用传统概率密度函数积分的方法不但不能显著提高融合的精度，反而增加了DEM融合计算的复杂度和计算量，严重制约DEM融合的批处理能力.另外，现有的升降轨融合方法并没有对失相干处(如水域)及两个轨道中均处于几何畸变处的空值进行处理，这也使得融合后的DEM不完整(Jiang et al., 2014; Rossi et al., 2012; Deo et al., 2014)；2) TSX/TDX 生成的DEM精度分析，现有的研究区域多针对地形陡峭的山区或者建筑物分布较为稀疏的居民区进行精度分析(Rossi et al., 2013; Deo et al., 2013; Brautigam et al., 2012).而对于地形特征复杂的区域(如沿海区域的水域、几何畸变较为严重的高楼林立的城区等)的数据融合以及精度评定等问题仍然缺少定量实验和研究.同时值得注意的是，已有的数据融合方法多基于TSX/TDX Stripmap模式下的数据，其空间分辨率为HRTI3标准下的12 m(Rossi et al., 2012; Deo et al., 2013; Brautigam et al., 2012).而这对于城区的地形信息解读还远远不够，因此如何在仅用Stripmap格式下的数据获取更高分辨率(原始分辨率，高达2 m)和高精度的DEM还需要进行系统的研究和尝试.

本文主要针对目前TSX/TDX干涉生成中存在的DEM融合，精度评定问题和缺陷进行系统研究，提出一种新的适用于大范围DEM制作的升降轨融合算法和流程对雷达侧视成像中的几何畸变处的信息进行相应的恢复，以及通过不同地物特征分析DEM的精度水平，并定量分析了该DEM对DInSAR相位干涉的改进效果.具体地，我们选择我国珠海和澳门为试验区域，利用TSX/TDX得到的DEM与外部Lidar获得的高精度参考DEM进行对比，并对其精度和潜在制约因素进行分析，为TanDEM-X全球覆盖DEM的制作和应用提供参考.同时，将迭代获取的DEM应用于高分影像的DInSAR技术，利于后续的相位解缠和形变信息的解译，为更高精度的时序InSAR形变监测提供有利条件.

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

2.1 TSX/TDX雷达卫星系统

TanDEM-X是德国宇航局于2010年发射的，是TerraSAR-X的姊妹卫星.TSX/TDX采用HELIX卫星运行方式并进行编队飞行形成一个单轨双天线系统对地球进行观测从而得到3D地表数字高程模型(Krieger et al., 2007).该系统下两颗卫星相隔很近(<400 m)，克服了重复轨道干涉时的时间去相关，同时TSX/TDX的基线直接由星载GPS测量得到，且可以实时灵活地调整和量测(Krieger et al., 2007)，其精度较常规SAR卫星(厘米级)的基线来说更为精确，达到毫米级.此系统下可以进行三种模式的观测： 1) 静态模式(monostatic mode)，TSX和 TDX互相独立； 2) 收发分置模式(bistatic模式)，即一发双收； 3) 交互收发分置模式(alternating monostatic模式)，通过每次脉冲后改变发射器来最小化时间去相关.此外，除了交叉轨道下的地形监测，TSX/TDX的收发配置模式下还能进行顺轨干涉测量，可以用来识别地表运动的物体，比如交通、舰船监测应用，海洋学中海浪波谱分析等.该双星系统获取的DEM精度较以往任何机载、星载获取的DEM都要高很多，其垂直方向的相对精度在局部坡度角<20%时候达到2 m，>20%时达到4 m，绝对精度也优于10 m，为全球DEM测量开创了新纪元(Krieger et al., 2007; Tachikawa et al., 2011).

2.2 基本原理

本文采用的SAR数据为TDX一发双收模式(bistatic)下获取的，如图1a，该模式下信号由其中一个卫星发射，两个卫星接受，通过对两个回波信号相位差进行分析来获取地形信息.同时，bistatic模式下的SAR影像对的时间基线几乎为0，能够最大限度地克服大气、时间去相关和轨道对干涉相位的影响，利于干涉相位的解译从而获得精确的地形信息.图1b为TSX/TDX雷达卫星系统的干涉原理图.

其中，TSX为TerraSAR-X卫星，TDX为TanDEM-X卫星，B为两者之间的基线长度，BX为基线B在水平方向上的分量，BZ为基线B在竖直方向上的分量，P为地表的任意一点，θ为主影像的入射角，α为垂直基线与水平方向的夹角，r1，r2分别为地面点P分别到两SAR传感器的距离，h为P点的高程，H为卫星到地面的距离.由其几何关系可得

图1 TSX/TDX雷达卫星干涉测量原理Fig.1 Theory map of interferometry of TSX/TDX

(1)

h=H-r1·cosθ,

(2)

(3)

(4)

其中，λ为波长，由于TSX/TDX的bistatic模式为单轨双天线系统，K=1.由式(2,4)可知，影响测高精度的主要因素有：卫星的斜距、基线长、平台高度等.

2.3 误差分析

由于X波段波长较短(3.1 cm)，对于地形起伏较大的山区及城区来说，地形相位将在干涉图上产生较密的条纹，这将会影响到后面的相位解缠甚至造成解缠失败.因此，引入外部DEM来去掉大部分的地形相位从而舒缓差分干涉图中的条纹，再通过差分相位进行解缠来恢复残余的地形信息.本文中用到的TSX/TDX下的bistatic模式的差分干涉相位表达式如下：

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

(5)

主要包含以下几部分：残余地形相位φtop_err、轨道残差相位φorb、大气相位φaps及噪声相位φnoise.由于该双星系统的时间基线几乎为0，因此可认为主从影像成像时刻的大气分布情况一致，因此φaps可忽略不计；TSX/TDX系统有着较常规SAR更精确的轨道信息，同时，我们还可以通过多项式拟合的方式来去除残余轨道相位(Fengetal., 2012;Zhouetal., 2013;Fengetal., 2015a)，因此轨道残差相位φorb也可以忽略不计，从而得到仅含有地形信息的干涉相位.由2.2节可知影响测高精度的主要因素为平台高度H, 水平基线BX，竖直基线BZ，入射角θ，斜距r，由式(4)可知：

(6)

其中：

(7)

由以上精度分析可知，TanDEM-XInSAR干涉测高的精度与基线成反比，基线越长其测高精度越高；而斜距的精度受到硬件设备的影响，如卫星采样时钟抖动及卫星钟的不确定性，该项误差与其产生的高程误差为同一数量级，对高程精度影响较小，可以忽略不计.

3 升降轨数据迭代融合算法提出和DEM技术流程

传统DInSAR技术虽然能够通过干涉相位获取地面的数字高程模型，但是单一轨道的干涉对容易受到高程模糊度、几何畸变、阴影和叠影等因素的影响而无法获取某些区域的地形信息.本文采用升降轨结合的方式，通过不同飞行方式所具有的不同几何关系和参数，如入射角、垂直基线、高程模糊度等，来对因单一轨道产生的阴影、叠影等区域的地形信息进行恢复和补充.目前大部分升降轨融合的方法多考虑与相干性、多视视数相关的干涉相位的概率密度函数，继而得到最后的权重因子(Jiangetal., 2014;Deoetal., 2014).但是TSX/TDX系统的时间基线几乎为0，因此其相干性除了在几何畸变处外，不管是在山区、还是在平坦地区均很好，多视视数越大，其相干性与干涉相位的标准差之间的相干性越小(Toughetal., 1995).因此，本实验中没有将相干性作为定权因子，而是通过提出一种升降轨数据迭代融合的简洁算法来实现大范围DEM数据制作处理.同时，由2.3节可知，不同垂直基线对应着不同的高程敏感性(高程模糊度)，因此，本文除了考虑叠影、阴影在升降轨中的位置信息外，还引入了垂直基线用来定权，从而得到更加精确的DEM模型.

升降轨数据融合的技术流程如图2所示，首先通过对单独的升降轨的CoSSC数据对(CoSSC为DLR通过TanDEM-XProcessor(ITP)数据处理软件进行了一定的辐射校正、精确配准后的产品(Nestoretal., 2010))进行常规的干涉处理，由于CoSSC数据已经进行了精确配准，因此只需要选择一定的多视视数进行干涉得到干涉图，接着引入外部90m分辨率的SRTM数据来削弱区域地形引起的密集条纹.此处，由于SRTM所得到的高程为正高(orthometicheight)，因此需要将高程异常(geoidalheight)去除使其转换至WGS84坐标系下，一般而言在局部区域这个高程异常影响较小可以忽略.然后，对得到的差分干涉图进行滤波和相位解缠，得到残余的地形信息，值得注意的是，由于第一次引入的外部SRTM分辨率较低，从而导致多视视数不能过小(过小的视数会导致外部SRTM无法配准至SAR坐标系下，同时在高楼林立处因外部DEM分辨率不够而导致相位解缠误差)，而最后得到的DEM的分辨率也并不理想.为了获得更为丰富的DEM细节信息，本文通过迭代的思想将第一次处理得到的最终DEM作为第二次的外部DEM，同时减小多视视数，再重复差分干涉处理从而得到更高分辨率的DEM，该操作能够很好地缓解因单视下高楼区的相位解缠误差.此外，本文对TanDEM-X获取的DEM中的水系也进行了一定的处理，通过借鉴PS-InSAR选点的思想，将振幅和相干性双向阈值对水系的轮廓进行一定的判定，记录其位置信息，并对其统一赋值(通过计算其周边像素来获得平均值)达到光滑水面的作用，从而得到更为精确可靠的数字高程模型.另外不同几何条件下获取的SAR影像在地理编码的时候有着一定的偏差，因此，还需要将升降轨得到的高程图统一配准到相同的坐标系统下.本文通过在反编码以后的强度图中选择一定数目(20个)的GCP控制点，再通过最小二乘的思想对其进行一次多项式拟合，得到影像对之间的偏移量关系式.

图2 升降轨数据融合流程图Fig.2 Flowchart of the fusion of ascending and descending data

对配准后的升降轨数据获取的DEM进行一定的融合处理，本文首先分别对升降轨外部DEM进行叠影、阴影区域探测：当坡度角为正的时候，如果坡度角大于等于入射角，此时在雷达影像上将产生叠影，同时也需要考虑坡度角小于入射角时被叠影影响的区域；当坡度角为负的时，如果坡度角的绝对值大于等于入射角时，雷达影像上表现为阴影，同时，当坡度角绝对值小于入射角时，某些区域也受到阴影的影响.接着分别记录探测得到的叠影、阴影区的位置信息为Mi，此处的位置矩阵中仅含有0、1两类元素，以叠影的情况为例，如果升轨中该像素位于叠影区域(记为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 实验数据处理

本文选取珠海、澳门等地作为实验区域，选取了两个不同时间点的升降轨TSX/TDX影像对作为本文的实验数据，详细数据信息如表1所示.干涉对的覆盖范围如图3所示，每景影像的覆盖面积约为33 km×56 km.该实验区域最高海拔约为450 m，高程大于30 m的平均地形坡度为16°，地势总体而言较为平坦，但是地貌和地物特征较为复杂，除了植被茂盛的山区外，还有较多的农田、水系及高楼林立的城市，因而有利于分析在不同地形和地物特征下TanDEM-X的数字高程模型的精度水平.为了验证本文算法获得的DEM的精度和可靠性，还获得了覆盖澳门及部分珠海区域的高精度Lidar数据得到的数字高程模型，其空间分辨率为5 m, 垂直精度高达1 m.

图3 研究区域及影像覆盖黑色实线框分别代表了两个影像时间点的TSX/TDX干涉对，其中一对升轨数据，两对降轨数据；蓝色实线框代表了本次选取的升降轨重复区域.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数据为例，介绍本实验中的计算参数和流程：首先对已经配准好的数据进行常规的差分干涉处理.为了抑制相位噪声，对干涉图进行方位向距离向4×4的多视处理(多视后分辨率约10 m)，再引入外部SRTM(分辨率为90 m)进行差分处理，得到去除基本地形的差分干涉图.接着，根据生成的差分干涉图进行相应的判断，如果信噪比较高的话，可以直接进行相位解缠工作，如果信噪比偏低，建议进行适当滤波处理；然后通过人工屏蔽的方法避免或者减弱低相干区域(海域及孤立的岛屿)对相位解缠的影响，使其不参与解缠计算，对滤波后的差分干涉图中相干性高于阈值(本文选择0.45)的点采用最小费用流(MCF)方法进行解缠，此处，为了后续将相对相位转为绝对相位，本文相位解缠的时候选择的参考对象是整景影像.此外，本文还通过对幅度、相干性设置一定的阈值来获取影像中水系的轮廓.最后，将绝对相位转化为高程信息，并反编码到地球坐标系下，再与重采样后的外部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的空间分辨率，并采用第三节中提出的融合方法进行融合得到最终的高精度的融合DEM产品.

4.2 实验结果及DEM精度分析

图4为我们采用第3节所述的方法将升降轨重复覆盖区域(如图3中蓝色方框所示)进行相应的融合得到的DEM，其中图4a为升轨TSX/TDX得到的DEM，图4b为降轨得到的DEM，图4c为升降轨融合得到的DEM，图4(d,e,f)为黑色框内的放大图，从图中可以看出单一轨道获得的DEM还存在一定的空值(叠影、阴影区及低相干区)，但不同几何条件下的空值位置基本不同，因此通过融合后的DEM能够有效地恢复出几何畸变处的信息，如图5f所示.为了能够对TanDEM-X的DEM进行一定的精度分析，本文选取了含有高精度Lidar数据覆盖的澳门及周边区，Lidar数据的具体情况如图5.从图中，我们可以清晰地看到起伏的山区、密集的澳门城区及新区，图中的两个矩形框分别代表了不同的地物，其中蓝色代表了山区、位于珠海境内，如图5a；红色代表澳门的老城区,如图5b，老城区中分布着密集的居民楼、办公楼，有的楼高达230多米.

图6a为融合后的DEM与Lidar的对比图，从左图中可以看出图中有三处光滑区，均为水域，由于其相干性较低，通过振幅、相干性阈值探测出其位置信息，并赋于周围像素平均值，同时也不参与后面的结果验证.图中的细节信息保留的较好，如山区的棱角较分明，阴影、叠影等地方的细节信息得到较好地恢复，因此，TanDEM-X获得的DEM在城区的信息较为丰富，如右上角处的民房及城镇.此外，由于Lidar数据的获取时间为2010年，在TSX/TDX获取时间之前，因此其地物之间会有一定的差异，且其在除澳门以外的城区均无观测数据.再者，由于升降轨的TanDEM-X数据间相隔较长，约两年，这也会使得地物之间存在较大的变化，如植被的生长、房屋的改建等，因此融合后的DEM也会存在一些因时间跨度较大导致的差异.为了定量的进行对比，本文选取了两个剖面，如图中的黑色实线m，n，其剖面图对应着图6(e, f)，图中升轨、降轨、Lidar及融合的DEM分别用红、绿、蓝及黑色实线表示，从图6e中的小图中可以明显看到，在阴影或者叠影的地方，融合后的DEM明显优于升降轨单独的DEM值，更接近Lidar数据.实线n对应的剖面图如图6f所示，在海拔较高的地方TanDEM-X的DEM与Lidar符合的较好，但是尾部，由于Lidar影像中平地没有进行监测，因而TanDEM-X的DEM较Lidar值偏高，为房屋区.随后，本文也对山区的DEM残差进行了统计分析，如图6i所示，其中红色、绿色及蓝色分别为为升轨、降轨及融合后的DEM与Lidar数据差值的统计，服从正态分布.由图中可以看到融合后的DEM的均值及均方差均有明显改善：均值由原来的1.75 m(升轨)、1.26 m(降轨)降到1.12 m，均方差由最初的3.88(升轨)、4.56(降轨)减少到3.57，具体见表2.

图6(c, d)为澳门老城区的融合DEM及Lidar对比图，图中的绿色平滑区域为水域，通过以上的方法进行识别，并赋予统一数值.从图中可以看出TanDEM-X的DEM结果较粗糙，尤其是在高楼林立的城区，由于雷达的侧视成像而产生很多阴影、叠影等问题在高楼的地方表现的较为明显，因此TanDEM-X获取的DEM在高楼分布较为密集的地方得到的结果较差.然而，在较矮(高程低于30 m)的居民区TanDEM-X获得的DEM的高程精度较高，能够较好地恢复出房屋的高程.此外，本文还对澳门的新城区进行了一定的实验，实验结果表明对于楼房分布较为稀疏的地区，TanDEM的DEM能够较Lidar数据更好地恢复出建筑物的外形及高度，去掉三倍中误差以外的统计值后其标准差基本上在5 m之内.同样，选择两条剖面q，p，从两条剖面图(图6g,h)我们可以看出高程低于30 m的城区监测精度与Lidar数据较符合，但是对于高楼区，大多会产生低估的现象，该现象主要是由于城区高楼之间的阴影、叠影以及相位中心的影响导致相位解缠的时候出现问题，尤其是高楼与平地相接的地方，其相位差可能超过了2π，因而导致解缠错误.通过统计其与Lidar的残差图时(图6j)，无论是单一平台还是融合后的DEM，其标准差均较山区来说大，但是，我们还是可以看到融合后的DEM无论是在残差均值还是标准差方面都有很大的改进，具体统计数值见表2.

表2 高程残差统计对比Table 2 Comparison of height difference

4.3 高精度DEM对于差分干涉的改进

除了对TanDEM获取的DEM进行外部DEM定量比较外，本文还对其进行差分干涉效果的分析.选用一对覆盖澳门、珠海等区域的TerraSAR影像对，主从影像的获取时间分别为2011年12月23日和2012年01月03日，其时间基线为11天，垂直基线为135 m，主影像的入射角为41°.我们选取与升降轨TanDEM所获取DEM重合的区域进行实验，将融合的DEM作为TerraSAR干涉对的外部DEM.为了能够比较分析TanDEM获取的融合DEM对差分干涉的改进，本文还选取了其他两类外部DEM(SRTM和ASTER)进行对比实验.为了使得SRTM能够较好地配准至SAR影像坐标下，本文先选用了方位向距离向为5×5的多视进行常规的DInSAR处理，差分干涉图的方位向距离向分辨率分别为9.47 m×10.39 m，其结果如图7(a—c)所示.从图中的黄色实线圈内可以看到，对于山区来说，TanDEM作为外部DEM时，其在差分干涉图中残存的条纹最少，明显优于另外两种外部DEM.为了能够获取TanDEM对于城区的差分干涉的改进，本文选择了覆盖澳门新城区的一小部分影像在5×5视数的基础上进行方位向距离向2×2的多视迭代，此时的空间分辨率高达3.79 m×4.12 m.同样采用三种不同的外部DEM，其结果如图7(d—f)，从图中可以看出，对于分布较为密集的高楼区，由于TanDEM对于高于30m的建筑物高程恢复能力较弱，因此还残留较多的条纹，如连接老城区的高架桥的周边，但SRTM和ASTER残留着较TanDEM更密的条纹；对于分布较稀疏的高层建筑物，如“威尼斯人度假村”及周边的几个高层建筑，TanDEM能够对其进行改善，从而得到精确去掉地形信息的差分相位，更有利于形变信息的提取(如图7f所示).而SRTM及ASTER作为外部DEM的差分干涉图上仍然残留较多地形误差，其可能原因是由于分辨率过低及配准过程存在的偏移误差.同样，对于2×2视下的山区，如图7(d—f)中黄圈内的山区，TanDEM较其他两者的有着明显的改进，尤其对于阴影、叠影等细节地方的地形也去的较为干净.值得注意的是，本次实验采用的干涉对的垂直基线仅为100多米，如果对于垂直基线较长的干涉对其改进效果将更加显著.

图6 TSX/TDX获取的DEM与Lidar的对比图(a—b)选取的山区的融合DEM及对应的Lidar影像；(c—d)选取的城区的融合DEM及对应的Lidar影像；(e—f)山区内m,n断面图；(g—h)城区内q，p断面图；(i—j)山区及城区的残差统计图.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)中的山区及城区范围，并统计这些范围内相干性大于0.8的点的相位，其结果如图8所示，其中横轴代表这些点的缠绕相位值，单位为rad，纵坐标表示这些参与统计的点的概率密度.由于本次采用的干涉对的时间基线仅为11天，因此可以假设这之间不存在形变信息，因而差分干涉图中所剩下的仅为残余地形、大气及噪声.因为我们所圈取的范围对于三类情况都一样，其大气和噪声的影响都是相同的.因此，我们可以通过其缠绕相位的标准差来判断残余地形相位的多少.图8a为山区范围内的缠绕相位统计结果，其中红色实线为SRTM作为外部 DEM改正后的差分干涉相位的分布，绿色实线和蓝色实线分别代表外部DEM为ASTER及TanDEM的相位统计图，其缠绕相位的标准差依次为：1.4 rad, 1.6 rad及1.4 rad.可以看出外部DEM对差分干涉图中山区的改善程度由高至低依次为：TanDEM>SRTM>ASTER.同样，我们也分别圈取了城区，并做同样的统计，其缠绕相位统计结果如图8b所示.当ASTER作为外部DEM时，其缠绕相位并不稳定，如图8b中的绿色实线的左侧有着较大的突变，而SRTM和TanDEM作为外部DEM时其相位统计分布较为正常，其缠绕相位标准差依次为：1.7 rad，1.3 rad及1.1 rad.其结果与山区范围内的一样，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作用下的差分干涉图的缠绕相位统计图红色实线代表SRTM为外部DEM，绿色和蓝色实线分别代表外部DEM为ASTER和TanDEM的相位统计结果.(a)山区中相干性高于0.8的点的相位统计结果；(b)城区中相干性高于0.8的点的相位统计结果.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能获得更为丰富的细节信息.在高楼分布较为稀疏的新城区TanDEM的DEM精度很高，但是对于高度过高的建筑物(远远超过高程模糊度)仍会存在一定的低估现象.而在高楼林立的老城区，由于高建筑物分布过于密集，几何畸变较为复杂，易产生较为密集的条纹，这也极大地阻碍了后续的相位解缠，从而导致城区的高程反演低估甚至失败.此外，通过TanDEM获取的DEM对DInSAR的改进较SRTM、ASTER等外部DEM更为明显，尤其对于山区，能够较好地去除地形信息，方便后续的形变解译；同时，在分布较为稀疏的新城区或者开发区，TanDEM作为外部DEM去除地形的效果很好，能够利于后续的DInSAR差分相位质量的提高，并为时序InSAR提供更为有利的条件.总的来说，TSX/TDX获得的DEM不管是空间分辨率还是精度都大大优于目前我们常用的SRTM和ASTER，但是城区的细节信息监测仍有待提高.因此高分的TanDEM数据及小视数甚至单视的干涉处理将成为后续的研究重点，通过更高分辨率的数据获取城区精度更高的DEM，同时也为后续的时序InSAR形变监测提供高分辨率外部DEM辅助.

5 结论

本文通过获取的两对覆盖澳门、珠海区域的TSX/TDX干涉对进行了一系列的干涉处理，成功获取了每个干涉对的高精度DEM，同时，通过本文提出的一种结合阴影、叠影、垂直基线的加权融合迭代方法在恢复因侧视成像而导致的几何畸变处的信息的同时得到较HRTI3更高空间分辨率的DEM.实验结果表明几何畸变恢复效果较为突出，通过引入的高精度的Lidar数据与TSX/TDX获得的DEM针对不同地物的精度进行一系列的对比分析，我们发现：针对山区来说，TSX/TDX获取的DEM精度很高，且通过升降轨融合后的DEM精度较单独轨道的精度要高；针对建筑物分布稀疏的城镇，其精度也较高；但对于高楼林立的城区来说，由于建筑物高差大且分布过于密集，从而产生较为严重的几何畸变和解缠误差，导致部分高建筑物的高度产生低估甚至误估，其精度较前两者偏低.与此同时，本文还将TanDEM-X获取的DEM作为TerraSAR的外部DEM来定量分析，通过统计差分干涉图中缠绕相位的标准差来判断其对常规DInSAR的改善效果，结果明显优于SRTM、ASTER.因此本文的研究将为大范围高分辨率、高精度DEM的制作起到参考作用，同时也为今后InSAR的形变监测起到促进作用.

致谢 本研究所用TerraSAR-X/TanDEM-X数据来源于德国宇航局(DLR)项目(XTI_GEOL6191和LAN0820).感谢澳门特别行政区地图绘制地籍局提供的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地震形变监测及参数反演研究. 地球物理学报, 58(5): 1634-1644, doi: 10.6038/cjg20150515.

廖明生, 林辉. 2003. 雷达干涉测量—原理与信号处理基础. 北京: 测绘出版社.

肖金群, 李志伟, 汪长城等. 2012. 相位补偿SAR差分干涉提取山区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监测地表形变.

*通讯作者 冯光财，男，讲师.E-mail: fredgps@csu.edu.cn