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。另外NL-InSAR可以利用较快的并行计算方案，在较快时间内完成高质量滤波处理。

在山区和地表起伏较大的复杂地形区域，雷达的侧视成像造成叠掩和阴影^[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$ 和 $z'$ 是两幅配准图像上对应的像素点，根据Goodman模型^[20]， $z$ 与 $z'$ 服从如下分布：

其中， ${Σ} $ 是 $2 \times 2$ 的相干系数矩阵。

其中， $R$ 与 $R'$ 是散射系数， $\beta $ 是真实相位， $D$ 是相干系数。

非局部相位滤波^[12]对3个参数 $R$ , $\beta $ , $D$ 联合估计，其结果一般会优于仅利用相位信息的滤波方法，SAR成像过程表明相位突变的地方通常伴随着幅度或者相干性的变化，因此不同的参数相互支撑彼此的估计。其核心是基于图像分块间的相似性质去除噪声，采用较大的搜索窗搜索相似的像素点，相似的像素不再局限于某个局部区域，当前像素的估计值由图像中与它具有相似邻域结构的像素加权得到。实际处理过程中通常分为3个步骤^[21]，如图1(a)所示，首先根据相似性准则搜索出相似块，块之间相似性通常定义为在每个对应像素的相似度的总和 $\Delta \left( {s,t} \right) = \displaystyle\sum\nolimits_\tau {\delta \left( {v\left( {s + \tau } \right),v\left( {t + \tau } \right)} \right)} $ , $\delta $ 为所选取的像素的相似性；选出相似块之后，这些块在第2步中结合起来估计雷达特征，如相干系数等；最后将第2步中的估计值投影到雷达图像空间。

假设两幅图像的散射特性相同，即 $R = R'$ ，这种假设在相干性较高的区域是成立的，同时令 $A = \left| z \right|$ , $A' = \left| {z'} \right|$ ; $\phi = \arg \left( {z{{z'}^*}} \right)$ 为去平地后的相位，则有：

NL-InSAR使用权重最大似然估计(Weighted Maximum Likelihood Estimation, WMLE)来估计权重，在分块s中，将参数估计记作：

在文献[22]中Seymour和Cumming区别于传统的最大似然估计假设，令 $R = R'$ ，将公式扩展到WMLE中，令 ${\varTheta _s} = \left( {{R_s},{\beta _s},{D_s}} \right)$ 以及 ${\mathcal{O}_t} = \left( {{A_t},{{A'}\!\!_t},{\phi _t}} \right)$ 。

其中， $a \!=\!\!\! \displaystyle\sum\nolimits_t\!\!\! {w\left( {s,t} \right)\!\frac{{{{\left| {{z_t}} \right|}^2} \!+\! {{\left| {{{z'}\!_t}} \right|}^2}}}{2}} $ , $x \!=\!\!\! \displaystyle\sum\nolimits_t \!\!\!{w\left(\! {s,t} \right)\!{z_t}{{z'}^*}{\!\!\!_t}}\; $ , $N = \displaystyle\sum\nolimits_t {w\left( {s,t} \right)} $ 。 $w\left( {s,t} \right)$ 是基于分块搜索计算的权重，如图1(b)所示。

由于噪声的存在，基于局部小窗口的相干性估计不仅会受到噪声的影响，而且会产生一定的偏差。而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所示。

其中， ${\gamma _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滤波的结果，从中可以明显看出，NL-InSAR滤波效果较好，在有效抑制噪声的同时对于细节特征的保持更为丰富；而基于窗口的Goldstein滤波含噪更多，在人造建筑区域如道路和建筑的细节损失较为明显。

同时，NL-InSAR滤波对于相干性的估计更加准确，选取相干系数图部分结果如图9所示，图9(a)采用NL-InSAR估计得到的相干系数图，图9(b)是Goldstein滤波估计的相干系数图。图像所在区域中存在相干性比较高的建筑物、失相干区域的水域、相干系数较低的草地，在非局部滤波的相干系数图可以明显看出这些差别，而Goldstein滤波之后的相干系数整体偏高，同时从相干系数统计直方图可以看出，Goldstein滤波器得到的相干系数较集中分布在高相干性区域，NL-InSAR得到的相干系数更加符合所在区域的相干性分布特征，Goldstein滤波器相干系数存在过估计的现象，非局部滤波产生无偏估计的相干系数的同时保持了图像的边缘细节。

为了降低相位解缠的难度，减少相位解缠误差，利用外部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%。

为了对所获取的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。