基于区域网平差的干涉合成孔径雷达数字高程模型校正

章 鑫, 刘世杰, 李彬彬, 童小华

(1.同济大学 测绘与地理信息学院，上海 200092; 2.同济大学 土木工程防灾国家重点实验室，上海 200092)

合成孔径雷达干涉测量(Interferometric Synthetic Aperture Radar，InSAR)是获取数字高程模型(Digital Elevation Model，DEM)数据的重要手段，能够描述地面的起伏形态特征，广泛地应用于地形分析、地质、水文监测、生态、农业以及军事应用等领域[1-2]，但其精度受到轨道定位、SAR影像的配准、相位解缠等精度的影响[3-4]。为了提高DEM的利用价值，需要对DEM进行精度提升。

大量学者对于InSAR 得到的DEM产品进行了质量提升研究。Crosetto[5]利用控制点及升降轨数据重叠区的公共点对ERS-1(Earth Remote Sensing Satellite-1)数据进行优化，分析并提出了高程误差模型，DEM精度由18.1 m提升至11.4 m，但是没提及DEM平面位置精度。Wessel等[6]总结了InSAR生成的DEM的校准方法，给出了TanDEM-X任务的DEM误差校正方案，在模拟与实际DEM中进行了DEM的区域网平差处理，也只实现了高程精度的提升。Gruber等[7-8]以激光测高(Ice、Cloud and Land Elevation Satellite，ICESat)数据的高程作为地面控制点，使用区域网平差的方式对TanDEM-X DEM进行系统误差校正，结果证明了通过DEM区域网平差可显著提高相邻DEM间的相对高程精度，但是没考虑DEM初始定位误差对ICESat约束的影响。和会等[9]利用ICESat 对南极ERS数据的 InSAR DEM做了纠正，但采用的是简单的线性高程纠正模型，需要密集的控制点分布，且没有利用DEM重叠区的一致性约束。Hueso等[3]结合E-SAR系统对TanDEM-X的校正方法进行实验验证，并给出TanDEM-X的校正方法过程，分析了DEM系统误差的成分特性。Hossein等[10]用ICESat数据作为控制点，利用移动平均插值算法得到DEM间的差值面来提高先进星载热发射和反射辐射仪传感器生成的全球DEM的高程精度，只利用了控制数据来进行高程上的误差拟合。张永俊[11]利用了多源控制点基于定标的方法逐景对星载分布式InSAR DEM做了矫正。已有研究采用了简单线性高程纠正模型，依赖密集的控制点，且没有考虑重叠区的约束。而在引入区域网平差方法时，直接将控制点的平面位置绑定到DEM对应位置的栅格进行高程校正，没有考虑DEM初始平面误差造成的影响，这会使控制点产生错误的高程约束。

ICESat-2[12-13]是2018年发射的搭载了光子计数激光测高仪的传感器，相比2003年发射的ICESat所搭载的全波形激光雷达系统,具有更高的精度和更小的激光足印。笔者同时考虑InSAR DEM的平面和高程改正，基于星载InSAR DEM误差特性分析，构建InSAR DEM区域网平差方法，利用ICESat-2激光数据产品进行高程约束，对DEM的平面、高程精度进行校正，输出大范围的高精度的DEM产品。并对比了有无平面控制情况下，区域网平差优化的DEM结果。

1 算法基础

DEM平差流程如图1所示。首先收集Sentinel A IW影像数据对，进行干涉处理生成对应的DEM数据；再对收集的ICESat-2激光点筛选生成区域网平差所需的高程控制点和检核点；最后针对生成的DEM，先后进行连接点匹配、区域网平差、DEM纠正处理过程。本文将区域网平差分解为平面定位六参数的优化和高程误差模型的平差求解两个独立的平差过程，联同连接点的地理坐标约束和高程约束来提高DEM的平面定位和高程精度。

图1 DEM平差流程

1.1 激光点筛选

针对目标DEM范围内的ICESat-2点，首先根据ATL08产品的属性标签进行筛选，基于Li等[14]提出的筛选流程来提取初始的ATL08高精度数据，提取的ATL08点分成两部分，分别用于控制和检核。此外，本文引入两个新的约束：① 由于DEM可能存在初始平面定位误差，保留的高程控制数据应该位于平坦的地形上，且激光点对应的DEM栅格应该处于平坦区域[15]，如式(1)所示，通过设定阈值对DEM栅格周围像素高程的方差进行限制。② 过于密集的点对平差结果没有明显改善，因此本文通过格网均匀地选取高精度的激光高程控制数据。对于部分格网内没有控制点分布的情况，本文引入虚拟控制点作为辅助约束。

(1)

式中：R和C为目标栅格周边邻域的行列号；μ为目标栅格的高程值。由此来统计其邻域栅格的高程差异变化。

1.2 DEM 连接点匹配

DEM在重叠区域存在高程误差和定位误差，提取连接点需综合考虑两者误差的影响。DEM特征包括高程、坡度和纹理等，其中高程信息最为直接地显示DEM地形信息。利用滑动匹配过程对DEM重叠区公共连接点进行初始提取，对于窗口内获得的整像素，进一步利用双线性插值方法，使定位误差校正于亚像素级，其中匹配测度函数采用比较稳健且效率较高的高程差标准差[3]，实现连接点高精度提取。

1.3 区域网平差

DEM平面平差是利用连接点的平面约束进行平面区域网平差。对于DEM，每个栅格的地理坐标由角点的大地坐标以及像元分辨率确定。平差前后的DEM的分辨率不会改变，因此对DEM整体做平面刚体变换纠正，如式(2)所示。

(2)

式中：X和Y分别为DEM某个栅格的地理坐标；x和y为对应的行列号；dx和dy为DEM的角点经纬度，即平面区域网平差过程将优化DEM的定位六参数中角点的平面坐标;矩阵A为旋转参数和比例参数(分辨率)形成的矩阵。此时，通过引入DEM连接点的平面坐标一致性约束，建立约束方程为

(3)

式中:X和Y分别为地理坐标，包含同名点的DEMJ和K对应的地理坐标在优化后应该达成一致。通过至少两对连接点建立的约束方程，即可通过最小二乘求解对应的DEM平面定位参数纠正模型的系数。

DEM高程区域网平差是为了校正不同数据集之间的相对误差，各DEM的高程误差改正模型参数通过高程控制点进行约束优化，同时利用重叠区的连接点高程值一致条件建立约束。通过Gruber等[7]的研究，针对InSAR生成 DEM，其误差分布模型符合：

(4)

式中：x为DEM的列号；y为行号。即由于SAR的成像模式，导致其方位向和距离向的误差分布规律不一致，在任意方位时刻干涉的各项参数误差恒定，因此距离方向高程纠正采用一阶模型，而方位向的阶数由于卫星的位置变化通常采用多项式模型[6]。对于DEMJ和K上某个确定的连接点对(xJ,yJ)和(xK,yK)，平差约束则为它们的原始高程加上对应的高程改正项达成一致，如式(5)所示。

hDEM_J(xJ,yJ)-gJ(xJ,yJ)=hDEM_K(xK,yK)-gK(xK,yK)

(5)

式中：gJ(xJ,yJ)和gK(xK,yK)分别为DEMJ和K的高程改正模型；hDEM_J(xJ,yJ)为DEMJ在(xJ,yJ)处的初始高程；hDEM_K(xK,yK)为DEMK在(xK,yK)处的初始高程。结合对控制点的约束，观测方程可写为

V=AX-L

(6)

式中：X为高程校正模型的系数向量。矩阵A为

(7)

式中：A的前两行表示 ICESat-2 测高数据，分别用于改进DEMJ和K。第3行表示相邻 DEM 连接点之间的高程限制。然后利用最小二乘原理消除各DEM的高程系统误差，得到误差校正模型系数。针对高程控制点缺失的部分，为了解决缺少控制点部分导致的局部精度低的问题且使平差稳定，引入高程虚拟控制点，降低对应的权重，来达到虚拟约束。

2 实验结果

2.1 实验数据

本文收集了9对间隔12天的Sentinel A IW影像数据，通过干涉处理生成了对应的DEM数据，如图2所示。初始的DEM由于初始系统和干涉过程引入的相位误差，导致DEM重叠部分存在明显的纹理不一致现象。

2.2 实验结果

首先对收集的ICESat-2激光点进行筛选，筛选前后DEM范围内的激光点分布对比如图3所示。筛选前6866个点，筛选后设置200栅格大小的格网约束。考虑到平坦地区的坡度要求为小于2°，且DEM分辨率为20 m，本文设置DEM栅格的邻域方差阈值为(20×tan2°)2。由此保留了1230个激光点，其中按照三分之一的比例从每景DEM上分布的激光点中随机选择检核点，剩余的作为控制点。随后匹配出508对同名连接点。平面和高程区域网平差依次进行后，得到了平差纠正后的DEM数据，平差前后的DEM对比如图4所示。DEM重叠区的高程不一致现象被明显纠正，通过连接点的高程一致性约束，平差的DEM在重叠区高程达到了很好的吻合。同时，为了对比本文方法和传统不考虑平面的平差方法，利用同一数据只进行了高程区域网平差，平差后DEM不一致性也能被良好地纠正。但在不考虑平面误差时，由于初始的平面误差，激光检核点分布的位置也会存在一定的错位。图5为DEM平差后局部一致性对比，选取两个DEM重叠区纹理进行了放大。图5(a)中初始的DEM平面和高程不一致性十分明显。不考虑平面偏移时，经过高程平差后的图5(b)结果表明高程一致性得到一定提升，但是平面上仍存在着初始偏移。而由图5(c)可知，考虑平面偏移的平差后DEM重叠区的高程和平面纹理的一致性都十分吻合。区域二的初始高程误差更为明显，整体对比也更为直观。图6为DEM平差前后检核点高程误差直方图对比，不考虑平面误差时，高程检核点精度也有提升，误差趋于0均值的正态分布，但是残差接近于0的部分相较本文方法相对较少。表1为DEM平差结果的精度评估，统计了检核点误差和连接点定位坐标的误差RMSE。经过平面平差，连接点的平面定位相对精度也有明显提升，从初始的182.462 m提升至14.887 m，优于一个像素。而只进行高程平差时，DEM的高程误差RMSE从平差前的19.372 m提升到了7.865 m，这是由于检核点点位分布仍存在错误的定位，经过本文提出的考虑平面平差的一体化平差方法，最后的高程精度提升至3.459 m。

图2 初始9景Sentinel DEM数据

图3 区域内所有ICESat-2数据

图4 平差前后DEM局部提升效果对比

表1 DEM平差结果的精度评估

图5 DEM平差后局部一致性对比

图6 DEM平差前后检核点高程误差直方图对比

4 结束语

提出了一种针对InSAR DEM数据的区域网平差纠正方法，同时考虑了DEM的平面位置偏移和高程误差，并通过ICESat-2激光测高数据进行高程约束和检核。对哨兵1A数据干涉生成的DEM数据进行了验证，通过和传统的不考虑平面误差的DEM纠正方法进行对比，证明本方法显著地提升了DEM的高程精度(提升至82%以上)以及平面的纹理一致性(提升至92%)。后续考虑将试验推广到更多传感器生成的DEM数据中，研究多传感器数据间的区域网平差方法。