利用相干系数改进的数字高程模型反演方法

乐 颖，夏元平，钱文龙

(1.东华理工大学测绘工程学院，江西 南昌 330013；2.福州市勘测院，福建 福州 350108)

0 引言

数字高程模型(digital elevation model，DEM)蕴含了丰富的地表高程信息，在测绘、地质地貌、工程建设、环境治理等诸多领域具有极其重要的作用[1-4]。由于我国南北方气候差异明显，尤其是在我国西南部区域，经常伴有多云多雨天气，通过摄影测量、光学遥感等常规手段很难获取这些区域的DEM[5]。此外，相较于光学遥感，合成孔径雷达干涉测量(interferometric synthetic aperture radar，InSAR)反演DEM时，不受天气状况的限制，可全天候全天时进行数据采集，并且具有高精度、高分辨率的优势，从而成为了获取 DEM 的常用手段之一[6]。针对不同地区、不同应用领域，分析其影响因素，并尝试各种改进方法获取DEM，从而改善DEM精度仍是目前研究的热点[7]。

InSAR技术是目前获取DEM的常用手段，通过获取同一区域两幅SAR影像，并利用影像上的相位信息进行差分干涉，从而反演出DEM。然而，InSAR技术由于轨道误差、大气影响、空间与时间失相关等固有缺陷，在实际应用中受到一定的影响。一些学者基于InSAR技术进行了改进研究，从而提高 InSAR 数据处理的效率和精度，并将其有效地用来反演DEM[8]。文献[9]以DEM数据为基础，结合强局部加权回归算法，利用反距离加权方法进行空间插值，从而获得了可靠性更好的河床重构数据。文献[10]提出了一种基于DEM辅助后向投影模型的InSAR高程反演方法，显著降低了干涉条纹的密度，有效地提高了DEM反演精度。文献[11]采用InSAR技术获取矿区DEM，并对定南县岭北离子型稀土矿区数据进行实验，得到了该矿区的地形地貌特征。但是，上述文献在InSAR反演DEM过程中，主要依靠人工凭借其经验选取轨道精炼控制点，这对操作者的专业性要求比较严格，否则难以选取误差较小的控制点。鉴于此，本文提出利用相干系数改进的数字高程模型反演方法。

1 传统InSAR技术反演DEM原理

SAR影像经由信号采集与数据处理之后，影像中每一像素都将蕴含与斜距相关的相位信息以及地面分辨元的雷达后向散射强度信息。InSAR技术主要从干涉相位与干涉相关数据中获取感兴趣的可用信息。而干涉相位一般与传感器到目标的距离有直接关系，是InSAR数据处理与信号提取的关键部分。在实际应用中，InSAR干涉相位ψ主要包括参考椭球面相位φref、地形相位φtop、形变相位φdef、大气相位φatm和噪声相位φnoi等多个分量。其中φref、φtop和φdef为主要组成部分，通常采用较短的时间基线和较长的空间基线干涉对；φnoi可采用低通相位滤波的方法有效抑制，φatm可根据相位的高通滤波进行一定的削弱。因此，干涉相位ψ可表示为[12]：

ψ=φref+φtop+φdef+φatm+φnoi。

(1)

φref是由参考椭球面本身引起的干涉相位，计算公式为：

(2)

式(2)中，λ表示波长，r1和r2分别为影像中任一像元等斜距投影到参考椭球面后对应于主、副影像的斜距，R1为任一像元对应的主影像斜距。

由于B∥=δR=r1-r2=Bsin(θ0-α)，将其代入式(2)，φref则可表示为：

(3)

式(3)中，θ0为主影像卫星雷达至地面点目标的侧视角，B和α分别是经由卫星轨道数据和配准结果计算所得到的对应基线长度和基线倾角。

φtop是由目标大地高h所引起的相位分量，计算公式为：

(4)

由于斜距r和卫星高度H可远大于观测目标大地高h，即由目标高程引起的侧视角变化量δθ将很小，因此，δθ可表示为：

(5)

将式(5)代入式(4)中，得到φtop分量的表达式为：

(6)

式(6)中，B⊥为干涉垂直基线长度。

φdef计算的是沿雷达视线方向的一维地表形变量，计算公式为：

(7)

式(7)中，Δr是由于观测目标位移造成的沿雷达视线方向的斜距变化量。

最后，通过获取同一区域的两景SAR影像形成一个干涉对，而对于地面上的每一分辨元来说，ψ直接对应着地面分辨元与两传感器间的距离之差的精确信息。因此，可利用传感器位置数据、SAR系统参数和基线参数等，借助于ψ计算每一地面分辨元的三维位置，即为DEM的反演。

2 利用相干系数改进的DEM反演方法

2.1 利用相干系数的DEM反演方法原理

将同一区域的主、副两幅SAR影像所对应的像素相位值进行差分，便可得获得一个一次差分相位图，即干涉相位图。在传统InSAR技术反演DEM的数据处理过程中，主要包括以下几个步骤：影像配准、干涉图生成、自适应滤波和相干性计算、最小费用流相位解缠、控制点选取、轨道精炼和重去平、相位转形变及地理编码等，处理流程图如图1所示。

图1 实验流程图Fig.1 Experimental flow chart

上述每个数据处理步骤过程中产生的误差都将对干涉结果的精度和可靠性产生巨大的影响，而利用相干系数改进的数字高程模型反演方法的具体操作步骤为：

1) 影像配准。利用轨道数据和参考DEM获取局部非参数偏移估计，并计算每个窗口的互相关函数，取最大值作为所选位置的参考位移；然后依据方位角和距离像元位置信息，经由多项式计算残差参数偏移，将其与原始局部非参数偏移估计相加，从而获得细化的残差参数偏移；最后，将局部非参数偏移和改进的基于互相关拟合进行配准，得到亚像素级精度。

2) 干涉图生成。根据两景配准后的SAR影像之间的相位差(Φ)，可计算得到地球上某地面目标与传感器位置之间的距离差。干涉相位的表达式

(8)

式(8)中，Image(I)和Real(I)分别表示干涉图的虚部和实部。

3) 自适应滤波和相干性计算。每个稳定区域中所包含的像素选择，主要来源于该稳定区域中所有像素的均值(Mall)与新像素值(Mnew)之间的差值，该差值对Mall进行了归一化处理，是可供选择的候选值。采用区域增长方法识别新的候选像素，相似均值因子是由一个数字在线性尺度上的表示，公式为：

(9)

一般情况下，需要多次过程迭代，才可获取最佳相似平均因子。

4) 最小费用流相位解缠。SAR影像中所有像素所对应的相位皆有整周模糊度问题，干涉图的相位只能是2π的倍数。在干涉处理中，本文采用最小费用流(minimum cost flow，MCF)方法为所有像素确定干涉相位的整周未知数。此类方法采用正方形的格网，综合了图像上所有的像元，并将相干性小于阈值的像元进行掩膜处理。将解缠最小相干性阈值设置为 0.15，若像元的相干系数大于0.15，则将该像元进行解缠。相较于区域增长法，当存在大范围的低相干或其他限制增长的因素导致解缠困难时，最小费用流算法可以取得更优的结果。

5) 控制点选取。在缺少地面控制点信息并对研究区形变没有先验认识的情况，盲目选择地面控制点将大概率引入误差，从而导致最终结果产生偏差。因此，为了减少轨道误差对DEM反演的影响，本文通过借助相干系数来剔除相干性低的像元点，保留相干性高的像元点作为控制点进行轨道精炼[13]。相干图可描绘两幅影像的相关程度，通过计算干涉条纹图的相干系数，能够定量评价SAR影像精确配准后所得到干涉条纹图的质量。相干系数于1993年最初由普拉蒂提出，相干系数的取值范围为[0，1]，值越接近1说明相干性越高，其理论模型公式为[14]：

(10)

式(10)中，E[·]表示数学期望，u*表示共轭复数，u1与u2为主副图像的信号。在实际的InSAR处理过程中，在同一像素中无法获取计算所需数量的采样值。因此，根据式(10)的理论模型，基于SAR影像复数数据计算获得相干系数的标准表达式为[15]：

(11)

式(11)中，u1(n，m)、u2(n，m)分别为主、副影像数据块内某个坐标(n，m)处的复数值；|·|2为对应的二阶范数；M与N分别为计算相干性的数据块尺寸大小；m与n为数据块内对应的行列号。

6) 轨道精炼和重去平。该步骤主要是通过输入控制点对轨道形态进行估算，而对于未解缠的相位信息，是否能够正确转换为高度(或位移)值也取决于该步骤是否成功。它可用于改进轨道(即纠正可能的误差)，也可用于计算相位偏移(即获得绝对相位值)，以及消除可能存在的相位偏移。

7)控制点选择精度评定。为了更系统地验证本文方法选取控制点的可靠性，首先统计干涉相干性系数图中符合要求的像元数量，然后根据相干系数阈值选取轨道精炼控制点，最后计算各点绝对残余误差，计算公式为：

(12)

式(12)中，γ是干涉相干性测度，H是高程模糊度。

8) 相位转形变及地理编码。绝对校准和解缠相位与合成相位重新组合，并将其转换为高度并编码成地图投影。通过对比距离-多普勒方法、相关的大地测量与制图转换方法，发现距离-多普勒方程可同时应用于两个天线，不仅能够提供各像素的高度，而且可得到其在给定的地图和大地参考系统中的具体位置。

2.2 评价指标

本文以标准差和均方根偏差作为精度评价指标，对DEM的精度进行定量评定[16]。标准差(standard deviation，SD)是方差的算术平方根，用σ表示，在概率统计中最常作为统计分布程度上的测量依据。标准差主要反映数据的离散程度，值越小，则说明数据越稳定。其计算公式如下：

(13)

均方根偏差(root mean square error，RMSE)是一种常用的测量数值之间差异的量度，主要体现观测值与真值之间的误差，值越小，则说明方法所得结果越好。其计算公式如下：

(14)

3 实例验证

3.1 研究区域概况

赣州市地处江西省南部，占据着省内最大的地域面积，总面积大约为整个江西省总面积的1/4，并且在江西省内拥有最多的区县数量。定南县作为赣州的区县之一，总面积达到了1 321平方公里，坐落于江西省最南部。定南县离子型稀土资源储量非常丰富，涵盖了20多种矿藏，并且县区内河流分布广泛，主要分属于定南水和桃江两水系[17]。

基于此，本文将定南县岭北镇作为研究区域，地理坐标位置为东经114°47′49″E～115°22′48″E，北纬 24°33′37″N～25°03′21″N。图2是赣州市定南县的土壤类型分布图，数据来源于世界土壤数据库[18]。尽管定南县的稀土种类丰富、品位高，但受限于落后的开采技术，单一的监测技术，乱采乱挖、釆富弃贫的现象仍然大规模的存在。长期以往导致大量矿产资源浪费，矿山地质环境愈发严重。因此，利用InSAR技术对DEM进行反演，从而获取最新的数字高程模型，对于地表形变监测具有十分重要的意义。

图2 赣州市定南县土壤类型分布Fig.2 Distribution of soil types in Dingnan County, Ganzhou City

3.2 数据源

3.2.1雷达影像数据

本文使用的雷达影像数据为C波段的Sentinel-1A卫星的SLC数据，选取了2020年两景影像数据，雷达影像数据的基本参数如表1所示。哨兵1号(Sentinel-1)卫星是欧洲航天局哥白尼计划中的地球观测卫星，由两颗卫星组成，载有C波段合成孔径雷达，可提供连续图像。轨道信息是InSAR数据处理中非常重要的信息，从最初的图像配准到最后的形变图像生成都有着重要的作用。含有误差的轨道信息造成基线误差以残差条纹的形式存在于干涉图中。因此，本文使用了哨兵1号相应的卫星精密轨道数据对轨道信息进行修正，可去除因轨道误差引起的系统性误差。

表1 雷达影像数据的基本参数Tab.1 Basic parameters of radar image data

3.2.2外部高程数据

在雷达影像数据处理过程中需要借助外部DEM，用来提供参考地形或者参考地理坐标系。本文选取的是SRTM 30 m数据，SRTM(shuttle radar topography mission)即航天飞机雷达地形测绘使命，是美国太空总署(NASA)和国防部国家测绘局(NIMA)以及德国与意大利航天机构共同合作完成联合测量，由美国发射的“奋进”号航天飞机上搭载SRTM系统完成[19]。SRTM系统共有约9.8万亿字节雷达影像的数据量，经过数据处理，制成了数字地形高程模型；于2003年开始公开发布，目前最新的版本为V4.1版本；分为30 m分辨率的SRTM1和90 m分辨率的SRTM3[20]。

3.3 实验结果分析

3.3.1DEM反演结果

本文根据相干系数阈值选取了35个轨道精炼控制点，然后分别计算传统InSAR技术反演DEM方法与本文方法选取轨道精炼控制点的绝对残余误差，结果对比如图3所示。从图3中可以明显看出，本文方法用于选取轨道精炼控制点整体的绝对残余误差较小，只有个别点的绝对残余误差偏大，误差符合要求；而传统InSAR技术反演DEM方法选取控制点的绝对残余误差整体偏大，后续需要筛选绝对残余误差较大的点。此外，文中为了对比分析InSAR反演DEM方法与本文方法用于选取轨道精炼控制点的准确性，引入了标准差指标，并计算得到两种方法绝对残余误差的标准差分别为0.71和0.34，足以定量说明本文方法准确性更高。因此，在缺少地面控制点信息或研究区形变先验认识的情况下，若需要选择可靠的地面控制点，可借助相干性系数进行剔除低相干性点。最终，通过轨道精炼和重去平步骤得到利用相干性系数的InSAR技术反演DEM与InSAR技术反演DEM这两种方法的轨道精炼多项式分别为-4.67+0.000 017R+0.002 5A和-10.37-0.000 041R+0.000 88A，其中，R代表距离向，A代表方位向。

图3 轨道精炼控制点精度对比Fig.3 Precision comparison of track refining control points

利用InSAR技术获取的DEM分辨率为15 m，数据整体表现优良，对地形的细节刻画明显优于SRTM DEM[5]。图4(a)、(b)和(c)分别为SRTM DEM、InSAR反演的DEM和本文方法反演的DEM；(d)、(e)和(f)分别是SRTM DEM山体阴影图、InSAR反演DEM山体阴影图和本文方法反演DEM山体阴影图。由于时间跨度不同具有部分差异，但是这三种DEM结果在空间分布上整体一致，说明利用InSAR技术反演得到最新的DEM是可行的。在相同分辨率情况下，本文方法比InSAR反演的DEM方法结果在纹理信息和高程变化分布上与SRTM DEM更相符。

3.3.2结果分析

为了对比分析InSAR技术反演DEM结果和本文方法反演DEM结果与SRTM DEM之间的差异，分别从定性和定量两方面对这两种方法反演的DEM结果进行对比分析。

1) 定性分析

相干系数图对DEM反演结果的质量有很大的影响，倘若整体的相干系数都处于一个很低的情况下，最终很难反演出准确的DEM结果。为了使得本文研究区域的相干性系数尽量高，本文选取了秋冬季节的两景影像，避免了夏季受植被大范围覆盖的干扰，图5为本文方法反演DEM的相干性系数图。相干系数的取值范围为在0和1之间，且值越接近1，代表着两幅影像的相干性越高。从图5中可以看出，整个研究区域整体的相干性都较高，从而说明，相干系数图的质量较优，所得DEM反演结果准确性较高。与此同时，将本文方法反演DEM的等高线叠加在SRTM DEM上，结果见图6。图6表明这两者的DEM结果高程起伏形态相同，且与SRTM DEM相比，本文方法反演DEM的分辨率更高，纹理信息更丰富，利用InSAR技术进行DEM反演结果是可靠的。利用相干系数的InSAR技术反演DEM可在减少轨道误差的情况下，更准确且精细地反演出赣州市定南县岭北镇最新的数字高程模型。

图5 本文方法反演的DEM相干性系数图Fig.5 DEM coherence coefficient diagram retrieved by this method

图6 本文方法反演的DEM等高线叠加SRTM DEMFig.6 DEM contour line inversion by this method is superimposed with SRTM DEM

2) 拟合分析

为了方便分析两组数据之间的内在联系，常常通过数学模型进行曲线拟合。首先将InSAR技术反演DEM结果、本文方法反演DEM结果和SRTM DEM裁剪至同一区域范围，并重采样使这三种结果的行列数相同，最后利用线性模型、二次模型、指数模型、几何模型、双曲模型和对数平方模型这六种常见数学模型对三种DEM结果分别进行两两拟合，并比较选取最优拟合模型，其结果见图7。本文方法反演DEM结果和SRTM DEM最优的拟合模型为对数平方模型，其表达式为Y=7 982.60-6 711.10×alog(x)+1 452.20×alog2(x)，拟合结果见图7(a)；InSAR技术反演DEM结果和SRTM DEM最优的拟合模型为二次模型，其表达式为Y=238.92+0.35x+0.000 4x2，拟合结果见图7(b)；本文方法反演DEM结果和InSAR技术反演DEM结果最优的拟合模型为二次模型，其表达式为Y=-159.20+1.12x-0.000 1x2，拟合结果见图7(c)。结果表明，本文方法反演DEM结果比InSAR技术反演DEM结果与SRTM DEM更相似，相关性更高。

图7 三种DEM间的拟合分析结果Fig.7 Fitting analysis results of three kinds of DEM

3.4 结果评估

利用InSAR技术对DEM进行反演可以获取最新的数字高程模型，对于地表形变监测具有十分重要的意义。由于赣州市定南县岭北镇的矿产资源丰富，因此本文结合已有矿权数据选取33个特征点代表InSAR技术反演稀土矿区DEM的精度对结果进一步定量评估，特征点的分布情况见图8。以SRTM DEM作为参考DEM，分别对比InSAR技术反演DEM结果、本文方法反演DEM结果和参考DEM之间的标准差和均方根差，其结果对比见图9。从图中可以看出，这些特征点的高程趋势高度一致，本文方法反演DEM结果比InSAR技术反演DEM结果更接近参考DEM的高程。SRTM DEM、InSAR技术反演DEM方法和本文方法的标准差分别为86.92、83.28和76.17，赣州市定南县多为山区和居民地，因此地势起伏差异会很大。InSAR技术反演DEM方法相较于SRTM DEM的均方根偏差和标准差分别是152.42 m和24.27 m；本文方法相较于SRTM DEM的均方根偏差和标准差分别是36.92 m和24.07 m。结果表明，相较于InSAR技术反演DEM方法，本文方法的均方根偏差和标准差分别降低了75.78 %和8.53%。由于InSAR反演DEM方法和本文方法所用数据和处理方法一致，仅仅是选取的控制点不同，因此，可以将这部分归因于轨道精炼误差导致的。利用相干系数改进的数字高程模型反演方法可以减少轨道精炼引起的误差，使得最终反演的DEM结果更准确。

图8 特征点选取分布图Fig.8 Distribution diagram of feature points selection

图9 SRTM、InSAR反演DEM和本文方法反演DEM结果对比Fig.9 Comparison of DEM inversion results from SRTM and InSAR and DEM inversion results from this method

4 结论

本文提出利用相干系数改进的数字高程模型反演方法。该方法通过设置一定的相干系数阈值，筛选出适宜的高相干性像元作为控制点进行轨道精炼处理，最后反演出DEM结果。以赣州市定南县岭北镇作为研究区域，借助 2020 年两景的Sentinel-1A卫星数据，以SRTM DEM作为参考DEM进行对比验证，得到以下结论：1) 利用InSAR技术进行DEM反演是可行的，其结果在空间分布和纹理信息上与SRTM DEM结果一致；2) 利用相干系数改进的数字高程模型反演方法可以减少轨道精炼引起的误差，其绝对残余误差的标准差为0.34，与InSAR技术反演DEM方法相比降低了52.11%，引入本文方法选取的控制点进行数据处理可以使最终反演出的DEM结果更准确；3) 以SRTM DEM为参考DEM，通过选取33个特征点对InSAR技术反演DEM方法和本文方法分别进行精度验证，结果表明，相较于InSAR技术反演DEM方法，本文方法的均方根偏差和标准差分别降低了75.78%和8.53%。本实验仅用免费的雷达数据进行了InSAR技术的反演DEM研究，对于其他雷达数据反演DEM的适用性还有待验证。后续可收集其他雷达卫星的影像数据进行DEM反演研究，或者将利用相干系数反演的DEM作为参考DEM，用于DInSAR技术中进行地表形变监测研究。