CN114689015A - 一种提高光学卫星立体影像dsm高程精度的方法 - Google Patents

Abstract

本发明涉及光学卫星立体影像技术。本发明针对缺少地面实测控制点影响了光学卫星立体影像生成DSM的高程精度，进一步限制了高分辨率DSM在测绘、水文、地质等领域的应用问题，提供一种提高光学卫星立体影像DSM高程精度的方法,利用卫星立体影像创建DSM，获取同一区域的卫星激光测高数据；对二者进行空间坐标变换后，进行配准，以二者的高程差的标准差最小时的配准状态作为配准结果，提高DSM的高程精度。利用ICESat‑2卫星测高数据高程精度高，在地面沿着直线轨迹密集分布，且轨迹间距离较大的特点，将多条ICESat‑2轨迹组合在一起可以近似为一个高程精度优越、沿激光轨迹方向分辨率高的DSM，从而利用DSM配准技术提高光学卫星立体影像生成DSM的几何精度。

Description

一种提高光学卫星立体影像DSM高程精度的方法

技术领域

本发明涉及光学卫星立体影像技术，特别涉及光学卫星立体影像高程精度技术。

背景技术

随着遥感和摄影测量技术的不断发展，高分辨率立体影像成为获取大范围内高分辨率DSM(Digital Surface Model)的主要途径。这得益于两方面技术的进步，第一，亚米级立体遥感影像成像技术，并且传感器类型趋于多样化；第二，高分辨率遥感影像处理技术，极大促进了DSM的获取效率。同时，随着研究的深入，在测绘、地质，水文等科学研究领域，要求获取更高精度的DSM。根据卫星线阵推扫式成像传感器的工作原理，在卫星设计中，影响卫星影像定位精度的主要因素有定轨精度、定姿精度、传感器精度及各时间信息精度等。摄影测量中提高卫星影像定位精度的常规技术有：a.传感器几何定标；b.区域网平差计算；前者需要对成像传感器和星敏感器等设备系统误差特性进行长期的分析来确保遥感影像在全球范围内定位精度的一致性；后者需要利用额外的信息来改正目标影像的相关参数中的误差,难点是并不总是可以获取地面实测控制点，在高原山区等交通不便的区域，一般无法获取数量足够的控制点，这就导致区域网平差技术的使用受到地域的限制；c.使用三维点云配准的方法，来抵偿成像条件对DSM造成的误差，同样能够提高DSM的几何精度。由于高精度姿态测定系统的发展，在无控制点条件下卫星影像平面定位精度可以得到理想的结果，但是高程精度的提升效果弱于平面精度，比如：Worldview-2影像在无控制点条件下，平面定位精度可以达到2.3m，高程精度为5m。

为获取全球范围内的高程数据，量化冰川消融对海平面上升的影响，估计全球生物量及监测地形变化，美国宇航局(NASA)于2018年9月15发射了第二颗冰、云和陆地高程卫星ICESat-2(Ice,Cloud,and land Elevation Satellite-2),其上搭载的先进地形激光测高系统(Advanced Topographic Laser Altimeter System，ATLAS)采用了微脉冲光子计数激光雷达技术，测量精度达到10cm。沿着轨迹方向相邻激光点相距70cm，激光足印直径约为17m，随着能量的消耗，最终直径能达到20m，可同时发射6束激光，每两束强、弱激光为一组，在垂直于地面轨迹方向上，相邻两组之间的距离为3.3km，每组内的激光相距90m。强弱激光束组合的模式便于ATLAS在低反射率，与高反射率地表获取足够的返回光子的同时节约能量，3组强弱光束形成交叉测量，能够有效提高对坡度和地表高程变化的检测能力。ATLAS数据产品共分为5个等级，所有数据均可在National Snow&Ice Data Center(NSIDC)获取。ATL08为3A级产品，记录了陆地及植被冠层的高程信息，还包括每个激光点的经纬度坐标及相关的质量评价参数，沿轨迹方向的空间分辨率为100m，时间分辨率为91天，数据覆盖全球。ATL08是在ATL03数据的基础上，采用DRAGANN滤波算法处理得到。ATL03为L2级全球地理定位光子数据，空间分辨率为0.7m，记录的信息包含每个光子位于WGS84椭球上的椭球高、纬度、经度和光子飞行时间，ATL03产品是低级别、针对仪器的产品(ATL01/02)和高级别、针对特定表面、以科学研究为中心的产品(ATL06及以上)之间的桥梁，提供了更高层次数据产品所需的所有信息。此外，ATL03数据中将每个光子分为的信号光子或背景光子，信号光子又被标识为陆冰、海冰、海洋、陆地和内陆水域五类，并提供了这些分类的置信度评估。ATL03数据应用多种地球物理参数校正，优化每个光子的高程，允许用户根据需要对这些地球物理参数进行编辑。

自ICESat-2数据发布以来，已经在近海岸水深测量，冰川消融及高度变化，森林火灾区域检测，森林生物量估计，森林冠层高度测量等领域表现出重要作用。ICESat-2精确指向定位(Precision Pointing Determination,PPD)技术使ICESat-2激光点平面地理定位精度优于6.5m。ATL08产品的表达的高程位置位于3DEP(3-D elevation program)DEM(Digital Elevation Model,简称DEM)和实际地表高程之间。实验表明，选择适当的参数进行筛选可以从ATL08中获取高精度的控制点。然而，对于ICESat-2激光测高数据，因为没有记录激光点足印影像，所以在区域网平差中的应用受到限制。另外，如果通过空间插值的方法采用离散点云构建DEM，结果容易受到离散点云空间分布状态的影响，不适合于ICESat-2数据。ATLAS的平面定位精度优于6.5m，高程精度在平原地区为0.2m，山区为2.0m，优于遥感影像直接对地定位时的高程精度。可见，利用ICESat-2卫星激光测高数据提高光学立体影像高程精度，是一种可行的思路，但是必须对两种数据进行精确匹配。

发明内容

本发明针对缺少地面实测控制点影响了光学卫星立体影像生成DSM的高程精度，进一步限制了高分辨率DSM在测绘、水文、地质等领域的应用问题，提供一种提高光学卫星立体影像高程精度的方法，利用ICESat-2卫星测高数据高程精度高，在地面沿着直线轨迹密集分布，且轨迹间距离较大的特点，将多条ICESat-2轨迹组合在一起可以近似为一个高程精度优越、分辨率高的DSM，从而利用DSM配准技术提高光学卫星立体影像生成DSM的几何精度。

本发明解决上述技术问题采用的技术方案是，一种提高光学卫星立体影像DSM高程精度的方法,包括以下步骤：利用卫星立体影像创建DSM，获取同一区域的卫星激光测高数据；对二者进行空间坐标变换后，进行配准，以二者的高程差的标准差最小时的配准状态作为配准结果，提高DSM的高程精度。

具体的，利用卫星立体像创建DSM时，影像几何模型采用RFM模型,输出分辨率为影像不超过辨率的2倍，并以规则格网形式进行存储。

具体的，包括以下步骤：获取同一区域的卫星激光测高数据后进行预处理，获得高质量激光测高数据,再进行空间坐标变换；所述预处理，包括以下步骤：

剔除卫星激光测高数据中的弱能量激光束数据，在剩余的数据中选取地表覆盖类型为陆地时的中、高置信度光子；所述中、高置信度光子沿着激光轨迹线性聚集分布，中、高置信度光子分布在垂直于激光轨迹方向线段上，取同一线段内光子的平均高程值作为此处的最终高程；剔除高程异常的光子。

具体的，剔除高程异常的光子，包括以下步骤：

设置光子与DSM高程差的阈值，将超过阈值的光子作为高程异常的光子剔除。

进一步的，对二者进行空间坐标变换后，进行配准，以二者的高程差的标准差最小时的配准状态作为配准结果，提高DSM的高程精度，包括以下步骤：

对DSM及高质量激光测高数据进行空间坐标变换后，设置最大平移旋转量，建立对应点关系，以高质量激光测高数据为源数据，DSM作为参考数据，绕X、Y、Z轴执行空间旋转变换；且在X-O-Y平面内进行平移变换，计算不同变换状态下对应点高程差的标准差；

选取高程差的标准差最小时的旋转平移量作为变换参数；

将DSM作为源数据，根据变换参数进行反向变换，从而提高DSM的高程精度。

进一步的，还包括以下步骤：在不同变换状态下，剔除该状态下高质量激光测高数据中的高程异常的光子后再重新计算该变换状态下对应点高程差的标准差。

进一步的，所述不同变换状态下对应点高程差的标准差，计算公式如下：

ΔH＝{Δh(i,j)/Δh(i,j)＝H_ATL(i,j)-H_DSM(i,j),i∈1,2…n,j∈1,2…m}；

其中，H_ATL(i,j)、H_DSM(i,j)分别为高质量激光测高数据的激光点与DSM构成的对应点的高程，Δh(i,j)为每一次变换后对应点的高程差，ΔH为对应点在所有变换状态下的高程差集合；i为对应点数量，j为高质量激光测高数据执行空间变换的次数；

其中，Δh_sd(j)为高质量激光测高数据在第j次空间变换下对应点的高程差的标准差，ΔH_SD为所有变换状态下对应点的高程差的标准差集合。

进一步的，所述变换参数的计算公式如下：

其中，T为平移矩阵；R为旋转矩阵；α、β、γ、Δx、Δy、Δz分别为x、y、z轴的旋转平移量；z轴方向的平移参数Δz为高程差的标准差最小时的高程差均值。

再进一步的，剔除高程异常的光子的计算公式如下：

其中，H_{u_ATL}(i,j)为第j次变换状态下剔除高程异常的高质量激光测高数据，H_{ut_ICESAT2}为所有变换状态下剔除高程异常点后的高质量激光测高数据集合；

则有：

ΔH_u＝{Δh_u(i,j)/Δh_u(i,j)＝H_{u_ATL}(i,j)-H_DSM(i,j),i∈1,2…n,j∈1,2…m}；

其中，ΔH_u为所有变换状态下重新计算的高程差集合，Δh_u(i,j)为第j次变换下重新计算的高程差。

再进一步的，所述最小高程差的标准差的计算公式如下：

其中，Δh_usd(j)为同一变换状态下重新计算的高程差的标准差；Δh_{usd_min}为所有变换状态下最小的高程差的标准差。

本发明的有益效果是，在本发明方案中，利用ICESat-2卫星测高数据高程精度高，在地面沿着直线轨迹密集分布，且轨迹间距离较大的特点，将多条ICESat-2轨迹组合在一起可以近似为一个高程精度优越、沿激光轨迹方向分辨率高的DSM，并提出一种基于地形相似性的DSM配准策略(Based on Terrain Similarity，BOTS)，在两个DSM间建立目标函数，通过移动、旋转DSM使目标函数最小化，将两个DSM在三维空间中对齐。从而在无地面实测控制点的条件下，实现ICESat-2激光测高数据辅助提高遥感卫星光学立体影像生成DSM的高程精度。

附图说明

图1是本发明实施例中实验区域分布图。

图2是本发明实施例中各实验区中的地形及ATL03数据与立体影像范围。

图3为本发明实施例中基于地形相似性的DSM配准流程图。

图4为本发明实施例中ATL03数据预处理流程。

图5为本发明实施例中同一轨迹中的ATL03光子点云图。

图6为本发明实施例中基于空间坐标变换的DSM配准示意图。

图7为本发明实施例中基于地形相似性DSM配准方法流程图。

图8为本发明实施例中利用SGM算法由5对立体像提取的DSM；a、b、c、d、e分别为：WV2_DSM、ZY3_DSM_N32E91、GF7_DSM、ZY3_DSM_N39E115、SV_DSM。

图9为本发明实施例中配准后DSM中心平面平移量；左图为每个DSM中心在东西方向的移动量，右图为每个DSM中心在南北方向的移动量。

图10为本发明实施例中DSM高程残差分布；a、b、c、d、e分别为各个DSM在ICP与BOTS算法配准前后高程残差绝对值的四分位数与均值。

图11为本发明实施例中DSM高程精度评价。

图12为本发明实施例中BOTS算法配准后DSM高程残差空间分布；a、b、c、d、e分别表示WV2_DSM、ZY3_DSM_N32E91、GF7_DSM、ZY3_DSM_N39E115、SV_DSM经过BOTS算法配准后高程残差绝对值的空间分布；ATL03DSM高程残差检查点符号的大小取绝于高程残差绝对值的大小，也就是符号越大表明高程残差越大。

具体实施方式

下面结合实施例及附图，详细描述本发明的技术方案。

本发明所述的一种提高光学卫星立体影像高程精度的方法,包括以下步骤：利用卫星立体影像创建DSM，获取同一区域的卫星激光测高数据；对二者进行空间坐标变换后，进行配准，以二者的高程差的标准差最小时的配准状态作为配准结果，提高DSM的高程精度。利用ICESat-2卫星测高数据高程精度高，在地面沿着直线轨迹密集分布，且轨迹间距离较大的特点，将多条ICESat-2轨迹组合在一起可以近似为一个高程精度优越、沿激光轨迹方向分辨率高的DSM，并提出一种基于地形相似性的DSM配准策略(Based on TerrainSimilarity，BOTS)，在两个DSM间建立目标函数，通过移动、旋转DSM使目标函数最小化，将两个DSM在三维空间中对齐。从而在无地面实测控制点的条件下，实现ICESat-2激光测高数据辅助提高遥感卫星光学立体影像生成DSM的高程精度。

实施例

为了进一步充分说明本发明的可行性，本例选择不同的地形地貌条件，采用Worldview-2、SV-1、GF-7、ZY-3立体影像验证了本文所提方法的可行性和有效性，并与目前最常用的迭代最近点算法(Iterative Closest Point，ICP)进行了比较。

为了充分说明本发明的有效性，选择位于不同地点、地形及地表覆盖类型不同的区域作为研究区。如图1所示，研究区A域位于柴达木盆地南缘，昆仑山东段北坡，北纬36°12'～36°20'，东经95°42'～96°3'的范围内。因为海拔在2700m-4300m之间，空气稀薄，导致区域内植被稀少，大部分山体上基岩裸露，除少数河床谷地有少量红柳、芦苇等灌木生长外，其他地域无任何植被，属高原荒漠山区。区域内地形复杂，山势陡峭，海拔落差大，是研究高原山区光学遥感影像对地定位精度的典型区域，实验结果具有代表性。研究区B位于青藏高原，北纬32°45'～33°21'东经91°40'～92°21'的范围内。海拔在4700m-6100m之间，区域内植被以高原草地为主，无高大植被，山顶有积雪覆盖，区域内地形复杂，山势陡峭，海拔落差大。研究区C位于长江中下游，北纬27°9'～27°20'，东经75°46'～75°57'的范围内。海拔在0m-1400m之间，区域内覆盖高大植被，属于山地地形，山势陡峭，海拔落差大。研究区D位于华北平原，北纬39°7'～39°40'，东经115°19'～116°1'的范围内。海拔在0m-1400m之间，区域内植被以农作物为主，地形以平原为主，西北区域为山地地形。研究区E位于印度，北纬27°10'～27°20'，东经75°46'～75°56'的范围内。海拔在380m-750m之间，区域内植被以农作物为主，地形以平原为主，东部区域为丘陵。表1描述了各研究区内的立体影像及主要参数。

表1：研究区立体影像及主要参数

基于以下原因，本例选取ATL03作为实验用的参考点云数据。a.虽然ATL08与ATL03都提供了地形表面高程数据，但是ATL08沿着轨迹方向的空间分辨率偏低，激光点数量少，不能详细的描述沿着轨迹方向的地形起伏；b.ATL08中记录的高程为地形高程与植被冠层的高度，所用DSM通过光学遥感影像获取，它并不能穿过地物获取地形高程，而ATL03中记录的高程包含了地表地物的高程，这一点与DSM相同；c.ATL03沿着轨迹方向空间分辨率高，能详细的表达地形起伏，点云数量庞大。所需ATL03数据通过所描述的工具获取，时间覆盖范围为2018年到2021年。但应当说明的是，无论采用哪一种数据进行配准，均应当落入本发明的保护范围。各实验区中的地形及ATL03数据与立体影像范围如图2所示。

本例以卫星全色立体像对为基础数据提取DSM，然后获取同一区域的ATL03，利用基于地形相似性的DSM配准策略，对二者进行空间坐标变换后，进行配准，以二者的高程差的标准差最小时的配准状态作为配准结果，提高DSM的高程精度。

具体DSM配准策略分为以下4个步骤，如图3所示，(1)(2)数据预处理。剔除置信度低的ATL03激光点，统一多源数据之间的空间参考系统；(3)基于地形相似性的DSM配准。以DSM为参考数据，ATL03为源数据，计算变换参数；(4)DSM逆变换；将DSM作为源数据，根据变换参数进行反方变换，使DSM与ATL03配准，提高DSM的几何精度。

其中，利用卫星立体像创建DSM，采用半全局匹配(SGM)作为密集匹配方法，最大可能的减少高程异常值的出现。根据实验，密集数据集作为参考数据时的配准结果比稀疏数据集作为参考数据具有更高的准确度，所以本文以立体影像提取的DSM作为参考数据，ATL03激光测高数据作为待配准数据，待求得配准参数后，再把转换参数逆应用到DSM上。影像几何模型采用RFM模型通过优化有理函数模型(rational function model，RFM)，输出分辨率不超过影像分辨率的2倍，DSM存储形式为规则格网，这是为了方便处理计算。

虽然ICESat-2激光测高的精度较高，但是受地形、采集时间、云层及气溶胶等因素的影响，激光点的垂直精度并不稳定，存在误差较大的数据点。所以配准开始前有必要对ATL03进行预处理，提取高质量ATL03激光高程数据，可以有效提高配准的准确度。本例中ATL03数据预处理过程如图4所示具体包括以下几个步骤：

首先，对ATL03筛选数据：

剔除卫星激光测高数据中的弱能量激光束数据，在剩余的数据中选取地表覆盖类型为陆地时的中、高置信度光子。研究显示，在ATL08数据中h_te_uncertainty(100m段内拟合高程总的不确定度)参数值较小的激光点质量较高，而在实验区域内强激光束的h_te_uncertainty参数值普遍低于弱激光束，且弱激光束在山体区域有大范围缺失。此外，如图5所示，在ATL03数据集中将所有光子分为信号光子与背景噪声，并通过signal_conf_ph参数给出了每个信号光子在陆地、海洋、海冰、陆地冰和内陆水中的置信度。当signal_conf_ph的取值为-2、-1、0、1时分别表示该光子属于TEP(transmitter echo pulse)事件光子，与特定表面类型无关的光子，噪声光子，缓冲光子，当signal_conf_ph的取值为2、3和4时分别表示该光子属于低置信度，中置信度与高置信度光子。综上所述，首先从ATL03数据中剔除弱能量激光束数据，然后在剩余的数据中选取地表覆盖类型为陆地时的中、高置信度ATL03光子作为实验数据。

signal_conf_ph参数的取值与ATL03光子的类别存在以下关系，signal_conf_ph＝-2、-1、0、1分别表示该光子属于TEP事件光子、与特定表面类型无关的光子、噪声光子、缓冲光子，signal_conf_ph＝2、3、4分别表示该光子属于低置信度光子，中置信度光子与高置信度光子；这里将，signal_conf_ph<2的光子统称为噪声光子；中、高置信度光子沿着激光轨迹线性聚集分布。

其次，对激光点稀释：

如图5所示，根据参数筛选出的中、高置信度光子沿着激光轨迹线性聚集分布。即聚在一起的中、高置信度光子分布在垂直于激光轨迹方向的近似于直线的直线段上，线段长度约为0.5m，线段与线段之间的距离为0.7m，同一线段上的激光点个数约为3-5个，且光子最大高程与最小高程最大相差5m。所以这里取同一线段内光子的平均高程值作为此处的最终高程。

最后，剔除高程异常数据：

经过参数筛选ATL03数据质量有了极大提高，但是此时ATL03数据中包含多种类型的光子，比如：植被冠层光子，冰面光子，地形光子及水体下的光子，而DSM表达了地形之上地物的高程，所以，在执行空间变换之前还需要剔除高程异常的光子。通常为位于地物以下的光子作为异常点去除。核心思想是通过设置ATL03激光点与DSM高程差的阈值，将超过阈值的激光点作为异常数据从ATL03中剔除。对于阈值，取地形相对平缓且地表裸露区域的ATL03激光点与DSM高程差的平均值。在本例中将采取反复剔除ATL03高程异常数据点的策略来最低程度的减小ATL03高程异常数据的影响，这一点将在下文中详细描述。

如图6所示，传感器在采集数据时，姿态、轨道测量值存在不可避免的误差，导致在相同坐标系下，同一区域的DSM数据集之间坐标值不一致，两者不能重合。对于同一数据集，各数据点采集条件相同，所受系统误差相等，此时两DSM数据集间的坐标不一致可等效为，使两个DSM数据集重合时，其中一个DSM数据集在坐标系中的平移旋转量。基于此思想，对数据处理完成后，对二者进行基于地形相似性的DSM配准，包括以下步骤：

对DSM及高质量ATL03数据进行空间坐标变换，设置最大平移旋转量，建立对应点关系，以高质量ATL03数据为源数据，DSM作为参考数据，绕X、Y、Z轴执行空间旋转变换；且在X-O-Y平面内进行平移变换，计算不同变换状态下对应点高程差的标准差。选取高程差的标准差最小时的旋转平移量作为变换参数；将DSM作为源数据，根据变换参数进行反方变换，从而提高DSM的高程精度。具体如图7所示。

首先，对DSM及高质量激光测高数据进行空间坐标变换，变换后的坐标可依据式(1)所示的变换矩阵计算,因为在本文中DSM与ATL03都表示WGS84坐标下的绝对坐标，所以在坐标变换中不考虑比例因素。

式中：T为平移矩阵；R为旋转矩阵；αβγ为x、y、z轴的旋转量；ΔxΔyΔz分别为x、y、z轴的平移量。z轴方向的平移参数Δz为高程差的标准差最小时的高程差均值。

其次，设置最大平移旋转量，目的是确定平面移动范围，减少不必要的计算量，提高配准效率。平面最大移动距离由立体遥感影像直接对地定位精度与ICESat-2的平面精度共同决定，依据影像类型而有所不同。

建立对应点关系，具体是用ATL03在DSM上内插，获取DSM高程值，并与ATL03构成对应点。为了削弱ATL03数据与DEM(Digital Elevation Model,简称DEM)之间的相关性，用对应点高差表达两个数据之间的关系。

对于同一范围内的地形，两个DSM数据所表达的地形起伏变化应当是相似的，当两者在二维平面方向对齐的时候，其不同平面坐标处的高差应当相等，也就是高程差的标准差最小。基于此，本文DSM配准以高程差的标准差为目标函数，坐标变换参数的计算分为两步，第一步，计算α、β、γ、Δx、Δy参数，第二步，计算z轴的平移参数Δz。

以ATL03作为源数据集，DSM作为参考数据集，绕x、y、z轴执行空间旋转变换且在X-O-Y平面内进行平移变换，依据(1)式计算变换后的坐标值，利用(2)(3)式计算不同变换状态下对应点高程差的标准差。

所述不同变换状态下对应点高程差的标准差，计算公式如下：

ΔH＝{Δh(i,j)/Δh(i,j)＝H_ATL03(i,j)-H_DSM(i,j),i∈1,2…n,j∈1,2…m}；(2)

其中，H_ATL03(i,j)、H_DSM(i,j)分别为ATL03的激光点与DSM构成的对应点的高程，Δh(i,j)为每一次变换后对应点的高程差，ΔH为对应点在所有变换状态下的高程差集合；i为对应点数量，j为高质量激光测高数据执行空间变换的次数；

其中，Δh_sd(j)为ATL03在第j次空间变换下对应点的高程差的标准差，ΔH_SD为所有变换状态下对应点的高程差的标准差集合。

为了减小ATL03高程异常数据的影响，本例中在不同变换状态下反复剔除ATL03高程异常数据点。具体为，在不同变换状态下，依据公式(4)剔除该状态下ATL03中的高程异常的光子后，再依据式(5)重新计算该变换状态下对应点高程差的标准差。

其中，H_{u_ATL03}(i,j)为第j次变换状态下剔除高程异常的高质量激光测高数据，H_{ut_ICESAT2}为所有变换状态下剔除高程异常点后的高质量激光测高数据集合；

则有：

ΔH_u＝{Δh_u(i,j)/Δh_u(i,j)＝H_{u_ATL03}(i,j)-H_DSM(i,j),i∈1,2…n,j∈1,2…m}；(5)

其中，ΔH_u为所有变换状态下重新计算的高程差集合，Δh_u(i,j)为第j次变换下重新计算的高程差。最后，依据式(6)及式(1)选取高差标准差最小时的旋转平移量作为α βγ Δx Δy参数。此时，计算高程差的标准差最小时的高差均值，作为z轴方向的平移参数Δz，以获取完整的6个变换参数。将DSM作为源数据，根据6个变换参数进行反方变换，从而提高DSM的高程精度。其中，最小高程差的标准差的计算公式如下：

其中，Δh_usd(j)为同一变换状态下重新计算的高程差的标准差；Δh_{usd_min}为所有变换状态下最小的高程差的标准差。

本例从统计学角度出发，以高差标准差最小为准则建立目标函数；从配准参数的解算方法来讲，在规定的平移旋转范围内以固定的步长遍历算有的坐标变换状态，相比于最小二乘原理迭代计算配准参数，本例方法避免了局部最优解的缺点。同时，与常用的ICP相比复杂度低，易实现，在有效配准两个点云的同时能够抵抗大量的异常点影响。

由研究区内的5对立体影像，采用SGM匹配方法提取了5幅DSM，如图8所示，a、b、c、d、e分别为：WV2_DSM，空间分辨率1.0m、ZY3_DSM_N32E91，空间分辨率5.0m、GF7_DSM，空间分辨率1.5m、ZY3_DSM_N39E115，空间分辨率5.0m、SV_DSM,空间分辨率1.5m。选用WGS84模型作为DSM的坐标系统，高程基准面为WGS84椭球高，与ATL03的空间坐标系统保持一致。为了全面、详细的评估本文所用配准方法的效果，将DEM的优化结果从平面、垂直两个方向与点到平面的ICP配准结果进行对比分析。

在二维平面方向，通过对比分析用ICP算法与本例方方法配准优化后的DSM中心在平面方向的移动向量，来评估本例基于地形相似性DSM配准算法的效果。表2统计了5幅DSM应用ICP算法与本例算法配准优化后，DSM中心在正东，正北方向的平移量，并统计了平移量差值。从中可以看出对于WV2_DSM、SV_DSM、GF7_DSM经过两种配准算法优化后，DSM中心在东、北两个方向的平移量差值达到了亚米级，小于DSM的空间分辨率，对于ZY3_DSM_N32E91及ZY3_DSM_N39E115DSM中心在东、北两个方向的平移量差值大于1m，但是也是小于DSM的空间分辨率。由此可以看出，在允许误差小于一个像素时，在平面方向本例配准算法与ICP点云配准算法效果相当，且DSM空间分辨率越小，本例配准效果越接近于ICP算法。图9显示了应用本例算法与ICP算法配准优化后，DSM中心平面平移量的空间分布，相比于表格的形式可以更加直观的看出DSM经过两种算法优化后，DSM中心平面移动向量的异同。

表2 DSM配准优化后中心平面平移量 单位：m

在垂直方向，采用交叉验证、对比分析的方法评估本例基于地形相似性(BOTS)DSM配准算法的效果。交叉验证的具体方案为：在最终筛选出的ATL03数据中，选取均匀分布于整个DSM的激光点作为高程检查点，用于统计DSM的高程残差，评估高程精度，剩余的点则用于优化DSM几何精度。图12显示了各个DSM中高程检查点的空间分布。图10中a、b、c、d、e分别统计了WV2_DSM、ZY3_DSM_N32E91、GF7_DSM、ZY3_DSM_N39E115、SV_DSM经过两种配准算法优化前后高程残差绝对值的百分位数及平均值。其中，各个DSM经过两种配准算法优化后高程残差有了较大改善；对于应用两种配种算法优化后的DSM，在0％～75％的DSM高程残差之间，残差变化量小，表明高程残差分布集中，但是DSM高程残差在75％～100％百分位之间时，残差变化量突然发生较大变化，表明高程残差分布离散，所以这里认为位于75％～100％百分位之间的高程残差点属于异常数据。此外，在各个百分位处，每个DSM应用两种配种算法优化后的高程残差绝对值及高程残差绝对平均值接近一致，表明在垂直方向基于地形相似性(BOTS)DSM配准算法可以达到ICP算法的效果。

根据图10的分析结果，在剔除了75％～100％百分位之间的异常高程残差数据后，表3统计了应用两种配准算法优化前后各个DSM中高程残差的均方根误差(RMSE)，以此作为指标，评价各个DSM的高程精度，通过对比RMSE来说明基于本例地形相似性(BOTS)DSM配准算法的效果。表3表明，在剔除异常检查点后，各个未配准优化的DSM高程残差RMSE值的相对大小与实际相符；在经过ICP算法与BOTS算法优化后，WV2_DEM、SV_DEM的高程残差RMSE值达到了亚米级，理论上GF7_DSM的高程残差RMSE也可以达到亚米级，但是由于GF7_DSM区域内覆盖茂密植被，进一步损失了GF7_DSM与ATL03的高程质量，导致GF7_DSM的高程残差RMSE大于1m；对于ZY3_DSM_N32E91与ZY3_DSM_N39E115，认为是空间分辨率小的原因导致高程残差RMSE值大于1m。

表3 DEM高程残差RMSE统计单位：m

图11直观的对比了各DSM应用两种配准算法优化前后DSM高程残差的RMSE值，结合表3的统计分析，表明在垂直方向上经过基于地形相似性(BOTS)DSM配准算法后的DSM的高程残差RMSE比未配准优化的DSM的高程残差RMSE表现出明显的优势，且效果与ICP算法相当。只是由于各立体像对的成像条件不同，导致优化的幅度有所差异。

图12中a、b、c、d、e分别显示了WV2_DSM、ZY3_DSM_N32E91、GF7_DSM、ZY3_DSM_N39E115、SV_DSM经过基于地形相似性(BOTS)地形三维点云配准算法优化后，高程残差绝对值的空间分布。可以看出检查点在各DSM上基本均匀分布，保证了DSM高程精度分析的可靠性。DSM高程优化结果受地形、地表覆盖类型及DSM空间分辨率影响较大。各个DSM中高程残差大的点位多分布于地形陡峭、复杂的山区，而平缓的区域高程残差较小。

综上所述，得益于地球轨道空间传感器位姿测量技术的发展，目前多数高分辨率卫星遥感影像直接对地定位精度优于10m，但是在许多高精度大地测量应用、地理环境监测研究中还达不到要求。卫星摄影测量中通常采用在区域网平差中加入少量的控制点来提高卫星遥感影像的定位精度，但是这项工作流程中需要准确识别控制点在影像中的位置，而在复杂的地形区域无法实现控制点实地测量。因此，本发明介绍了一种基于地形相似性(BOTS)的DSM与ICESat-2激光测高数据配准算法来优化了卫星DSM的几何精度。通过实验发现，在二维平面方向，BOTS算法与ICP算法配准精度的差距小于一个像素且DSM分辨率越高差距越小，在误差允许的范围内可达到与ICP方法同样的配准精度；在垂直方向BOTS算法与ICP算法配准精度接近一致。从算法复杂性来讲，BOTS算法比ICP算法更加简单易实现。由于BOTS算法通过遍历所有的可能的坐标进行配准参数获取，避免了ICP算法局部最优解的限制。在实验中发现，在ATL03数据中存在大量异常值或ATL03与DSM的初始空间姿态相距大时ICP配准会失败，显示本例的配准方法更具有鲁棒性。融合ICESat-2优化后的DSM高程RMSE(Root mean square error)可以达到亚米级，高程精度提高了73％～92％，其中Worldview-2 DSM优化后的高程RMSE为0.71m。

Claims (10)

1.一种提高光学卫星立体影像DSM高程精度的方法,其特征在于，包括以下步骤：利用卫星立体影像创建DSM，获取同一区域的卫星激光测高数据；对二者进行空间坐标变换后，进行配准，以二者的高程差的标准差最小时的配准状态作为配准结果，提高DSM的高程精度。

2.根据权利要求1所述的一种提高光学卫星立体影像DSM高程精度的方法,其特征在于，利用卫星立体像创建DSM时，影像几何模型采用RFM模型,输出分辨率不超过影像分辨率的2倍，并以规则格网形式进行存储。

3.根据权利要求1或2所述的一种提高光学卫星立体影像DSM高程精度的方法,其特征在于，包括以下步骤：获取同一区域的卫星激光测高数据后进行预处理，获得高质量激光测高数据,再进行空间坐标变换；所述预处理，包括以下步骤：

剔除卫星激光测高数据中的弱能量激光束数据，在剩余的数据中选取地表覆盖类型为陆地时的中、高置信度光子；所述中、高置信度光子沿着激光轨迹线性聚集分布，中、高置信度光子分布在垂直于激光轨迹方向线段上，取同一线段内光子的平均高程值作为此处的最终高程；剔除高程异常的光子。

4.根据权利要求3所述的一种提高光学卫星立体影像DSM高程精度的方法,其特征在于，剔除高程异常的光子，包括以下步骤：

设置光子与DSM高程差的阈值，将超过阈值的光子作为高程异常的光子剔除。

5.根据权利要求4所述的一种提高光学卫星立体影像DSM高程精度的方法,其特征在于，对二者进行空间坐标变换后，进行配准，以二者的高程差的标准差最小时的配准状态作为配准结果，提高DSM的高程精度，包括以下步骤：

对DSM及高质量激光测高数据进行空间坐标变换后，设置最大平移旋转量，建立对应点关系，以高质量激光测高数据为源数据，DSM作为参考数据，绕X、Y、Z轴执行空间旋转变换；且在X-O-Y平面内进行平移变换，计算不同变换状态下对应点高程差的标准差；

选取高程差的标准差最小时的旋转平移量作为变换参数；

将DSM作为源数据，根据变换参数进行反向变换，从而提高DSM的高程精度。

6.根据权利要求5所述的一种提高光学卫星立体影像DSM高程精度的方法,其特征在于，还包括以下步骤：在不同变换状态下，剔除该状态下高质量激光测高数据中的高程异常的光子后再重新计算该变换状态下对应点高程差的标准差。

7.根据权利要求6所述的一种提高光学卫星立体影像DSM高程精度的方法,其特征在于，所述不同变换状态下对应点高程差的标准差，计算公式如下：

ΔH＝{Δh(i,j)/Δh(i,j)＝H_ATL(i,j)-H_DSM(i,j),i∈1,2…n,j∈1,2…m}；

其中，H_ATL(i,j)、H_DSM(i,j)分别为高质量激光测高数据的激光点与DSM构成的对应点的高程，Δh(i，j)为每一次变换后对应点的高程差，ΔH为对应点在所有变换状态下的高程差集合；i为对应点数量，j为高质量激光测高数据执行空间变换的次数；

其中，Δh_sd(j)为高质量激光测高数据在第j次空间变换下对应点的高程差的标准差，ΔH_SD为所有变换状态下对应点的高程差的标准差集合。

8.根据权利要求7所述的一种提高光学卫星立体影像DSM高程精度的方法，其特征在于，所述变换参数的计算公式如下：

其中，T为平移矩阵；R为旋转矩阵；α、β、γ、Δx、Δy、Δz分别为x、y、z轴的旋转平移量；z轴方向的平移参数Δz为高程差的标准差最小时的高程差均值。

9.根据权利要求8所述的一种提高光学卫星立体影像DSM高程精度的方法，其特征在于，剔除高程异常的光子的计算公式如下：

其中，H_{u_ATL}(i，j)为第j次变换状态下剔除高程异常的高质量激光测高数据，H_{ut_ICESAT2}为所有变换状态下剔除高程异常点后的高质量激光测高数据集合；

则有：

ΔH_u＝{Δh_u(i，j)/Δh_u(i，j)＝H_{u_ATL}(i，j)-H_DSM(i，j)，i∈1，2…n，j∈1，2…m}；

其中，ΔH_u为所有变换状态下重新计算的高程差集合，Δh_u(i，j)为第j次变换下重新计算的高程差。

10.根据权利要求9所述的一种提高光学卫星立体影像DSM高程精度的方法，其特征在于，所述最小高程差的标准差的计算公式如下：

其中，Δh_usd(j)为同一变换状态下重新计算的高程差的标准差；Δh_{usd_min}为所有变换状态下最小的高程差的标准差。

Images