CN114964176A - 一种月球永久阴影区域地形测绘方法 - Google Patents

Abstract

本发明涉及一种月球永久阴影区域地形测绘方法，包括：步骤S1、对月球轨道器激光高度计LOLA的每个激光剖面在其接受视场范围内进行定步长调整，以单轨激光数据邻域内其余的激光观测值作为约束，最小化插值高程和观测高程差值的加权均方根误差RMSE，实现激光剖面的逐轨自约束调整；步骤S2、将调整后的激光剖面更新为初始状态数据进行迭代，不断收敛平面位置调整值；步骤S3、对迭代调整后的激光数据进行异常值剔除，得到误差修正后的月球激光测高数据；步骤S4、基于误差修正后的月球激光测高数据，得到表征月球永久阴影区域三维地形的数字高程模型。与现有技术相比，本发明消除了以消除原始图像中大量的地形伪影，提高了月球永久阴影区域地形测绘的质量。

Description

一种月球永久阴影区域地形测绘方法

技术领域

本发明涉及月球测绘技术领域，尤其是涉及一种月球永久阴影区域地形测绘方法。

背景技术

由于月球自转轴的微小倾斜，阳光很难直接到达月球极地的某些区域，这些区域被称为永久阴影区(PSRs)。月球激光高度计能够实现月球地形精细测绘，尤其对月球永久阴影区域(PSRs)具备主动测绘的优势，这为着陆选址、着陆器当行和水冰探测等行星探测任务提供基础地形数据。然而，航天器轨道控制、激光指向和硬件故障等多种因素会导致激光光斑地理位置的不确定性，这会使激光数据存在误差，三维地形产品出现大量伪影。美国国家宇航局(NASA)发布的月球轨道器激光高度计(LOLA)数据作为截止目前月球上最高精度的激光测高数据，也存在上述问题。

数字高程模型(DEM)能够反映月表三维地形，是月球实地勘探的基础和关键，它可以帮助人类和巡视器清楚地看到周围环境并识别重要特征，为勘探任务规划和决策提供参考。与使用太阳作为光源的光学遥感不同，激光成像不依赖于外部光，可以利用自身的光源测绘PSRs地形。

因此，月球轨道器几乎都配备了用于高精度地形测绘和月表控制的激光载荷。虽然激光观测本身具有精度高的优势，但由于轨道重建质量和激光指向偏差等因素的不确定性导致激光测量存在误差，并且这些误差可以直接映射到高分辨率地形产品中，呈现出与激光足印轨迹一致分布的地形伪影。因此，激光足印的地理定位校准对于激光高精度测绘尤为重要。

修正激光足印的地理定位误差的方法有交叉轨法。交叉轨主要是通过最小化不同轨激光观测数据间在交叉点的距离对激光足印位置的调整；直接测高法是基于先验的高精度地形数据，将每一次的激光观测与先验地形最佳匹配，从而达到修正误差的目的。但对于月球激光测高数据而言，LOLA作为目前最高精度的月表观测数据，无法获取到比其更高精度的先验地形数据；并且交叉轨的方法对于密集观测极区交叉点的计算量大，并且几何约束弱。

Barker等人利用缩减LOLA DEM作为参考地形数据，以自约束的方式改进了LOLADEM产品，但是，该方法在分批调整激光剖面时，会用所有剩余点的生成缩减DEM，这使得远离目标剖面的区域在利用激光点构建DEM时耗时且冗余。此外，如果同一批激光剖面空间分布十分密集，则剩余作为生成参考DEM的点非常稀疏，因此会导致生成的缩减DEM误差较大。

针对上述问题，亟需设计一种新的月球永久阴影区域地形测绘方法来修正激光点的误差，从而精确绘制永久阴影区域三维地形，以消除原始产品中大量的地形伪影。

发明内容

本发明的目的就是为了克服上述现有技术存在的缺陷而提供了一种月球永久阴影区域地形测绘方法，该方法最大程度上消除了地理定位误差，同时剔除了异常点，从而消除大量地形伪影。

本发明的目的可以通过以下技术方案来实现：

本发明提供了一种月球永久阴影区域地形测绘方法，该方法包括：

步骤S1、对月球轨道器激光高度计LOLA的每个激光剖面在其接受视场范围内进行定步长调整，以单轨激光数据邻域内其余的激光观测值作为约束，最小化插值高程和观测高程差值的加权均方根误差RMSE，实现激光剖面的逐轨自约束调整；

步骤S2、将调整后的激光剖面更新为初始状态数据进行迭代，不断收敛平面位置调整值；

步骤S3、对迭代调整后的激光数据进行异常值剔除，得到误差修正后的月球激光测高数据；

步骤S4、基于误差修正后的月球激光测高数据，得到表征月球永久阴影区域三维地形的数字高程模型DEM。

优选地，所述步骤S1包括以下子步骤：

步骤S101、对LOLA激光点进行数据预处理；

步骤S102、选取单轨RDR数据作为目标剖面T，根据其分布范围，各方向向外扩展k₁米，从剩余激光点中裁剪出作为参考地形数据的参考激光点集R；

步骤S103、以参考激光点集R生成用于进行最近邻搜索的KD树；

步骤S104、建立半径为k₂米的预设调整区域，调整目标剖面；

步骤S105、采用反距离权重插值，根据k₃米范围内的最近邻k₄个点，插值出每次调整位置后的高程值z_interp；

步骤S106、根据目标剖面T的实际观测高程值与插值高程值z_interp计算加权的均方根误差RMSE，表达式为：

其中，v为高程差，p采用Huber权函数；

步骤S107、完成预设调整区域内的每一步调整，记录最小的RMSE值对应的平面位置调整值(Δx,Δy)以及高程差v_T，更新目标剖面T的平面坐标(x_T,y_T,z_T)，表达式为：

其中，(x₀,y₀,z₀)为原始平面坐标，Δz为v_T的加权平均值；

步骤S108、完成当前单个目标剖面T调整后，进行下一个激光剖面的自约束迭代调整。

优选地，所述步骤S101中数据预处理包括噪声剔除、投影变换以及逐轨存储。

优选地，所述步骤S104具体为：基于极地投影坐标系与以延轨和垂轨方向的局部坐标系之间的转换关系，分别在垂轨和延轨方向以步长k₅调整目标剖面T，坐标转换表达式为：

式中，θ为激光剖面T与极地投影坐标系x轴的夹角，x_aC、y_aC是在局部坐标系下的坐标，x_SP、y_SP为激光点在极地投影系下的坐标。

优选地，所述k₁、k₂、k₃、k₄和k₅分别设置为150、50、100、10和2.5。

优选地，所述步骤S106中加权的均方根误差RMSE表达式为：

其中，v为高程差，p采用Huber权函数，n为激光点数量。

优选地，所述步骤S3具体为：构建标准化去趋势坡度图和残差统计量，剔除伪地形观测的异常值。

优选地，所述步骤S3包括以下子步骤：

步骤S301、利用调整后的激光点生成预设分辨率的坡度图，采用预设像素窗口大小进行中值滤波处理后，计算每个点对应的去趋势坡度S：

S＝S₁-S₂

式中，S₁，S₂分别为每个点的局部坡度和滤波坡度；

步骤S302、采用粗糙度R对去趋势坡度S标准化，构建统计量δ₁：

δ₁＝S/R＝(S₁-S₂)/R

步骤S303、计算自约束迭代调整过程中最后一轮调整后得到的高程差v的绝对中位差MAD，将绝对值大于K倍MAD的值挑选出来构成新的统计量δ₂；

步骤S304、根据统计量δ₁和δ₂，按照预设分位值对分布在双尾的激光点进行剔除。

优选地，所述步骤S302中的粗糙度R为坡度中位数。

优选地，所述步骤S303中的K为3。

与现有技术相比，本发明具有以下优点：

1)本发明采用自约束迭代调整的方法对研究区域内的所有LOLA激光点逐轨进行调整，在调整过程中无需借助外部高精度参考地形数据，通过LOLA邻近分布数据之间的几何约束，整体调整LOLA测高剖面数据，消除了激光点的地理定位误差，消除了在DEM地形数据中的大量地形伪影，提高了月球永久阴影区域三维地形图质量；

2)本发明通过标准化去趋势坡度和最后一次迭代的高程残差用于构建统计量，以过滤掉一些会导致DEM中异常突起或坑洞的异常值。

附图说明

图1为本发明的逐轨自约束迭代调整流程图；

图2为本发明的实验数据展示图；其中，(a)为LOLA RDR单轨激光剖面数据，(b)为5米分辨率的LOLA DEM，(c)为NASA发布的5米分辨率改善的DEM，(d)为1.5米分辨率的LROCNAC影像；

图3为单轨自约束调整示意图；

图4为每轮迭代中各轨激光点平面位置调整值；

图5为DEM结果展示；

图6为结果比较分析图；

图7为地形剖面图；其中，左图为地形剖面全局，右二为局部放大图；

图8为NASA改善的DEM和本发明的DEM差异特征分析结果示意图；

图9为沙克尔顿撞击坑(局部)三维地形渲染图。

具体实施方式

下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例是本发明的一部分实施例，而不是全部实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动的前提下所获得的所有其他实施例，都应属于本发明保护的范围。

实施例

本实施例给出了一种月球永久阴影区域地形测绘方法，包括：

1、获取逐轨分布的RDR数据产品以及LDEM数据；

本发明的研究对象为LOLA激光测高数据。LOLA数据科学团队将数据产品进行分级处理并发布，其中EDR(Experiment Data Records)产品为原始的科学观测数据，经过解算和校准后生成RDR(Reduced Data Records)数据产品。

RDR数据产品中激光观测值按照时间排列，形成包含时间戳、激光足印三维位置、激光脉冲特征参数等信息的单轨激光测高数据。对RDR中的激光点经过中值滤波及空间插值处理后，生成了月球数字高程模型(LDEM)数据产品，该产品中存在大量的地形伪影。

本发明以逐轨分布的RDR数据产品作为实验数据，LDEM数据作为对比分析数据，同时，以NASA新发布的改善后的DEM数据作为精度评价数据，以更高分辨率的LROC(LunarReconnaissance Orbiter Camera)NAC(Narrow Angle Camera)影像辅助用作目视判别个别地形特征的真伪。

2、对单轨RDR激光剖面进行逐轨自约束迭代调整

对研究区域内的所有激光点，以单轨RDR激光剖面作为对象，逐个实现自约束调整。主要思想为：利用单轨激光剖面与其邻近分布的激光观测点存在几何空间上的约束，利用这种关系必定能找到目标剖面在邻域范围内与地形最佳拟合的位置，拟合程度的准则用加权均方根误差(RMSE)衡量。如图3(a)中灰色点表示邻域内分布的作为参考地形的激光点(投影在二维平面)，其他点表示目标剖面在三个不同的平面位置对应的空间分布，对应的误差棒表示观测高程值与插值高程的残差，图3(b)是高程残差对应的分布图。显然，B剖面与地形更加吻合，所以其对应的平面位置更接近于目标剖面的真实位置。

该方法主要的实现流程如图1所示。详细步骤如下：

1)对所有激光点预处理。包括噪声剔除、投影变换，然后逐轨存储；

2)选取单轨RDR数据作为目标剖面T，根据其分布范围，各方向向外扩展150米从剩余激光点中裁剪出作为参考地形数据的激光点集R；

3)以参考点集R生成kd-tree，用于最近邻搜索；

4)建立50m为半径的预设调整区域。建立极地方位投影坐标系与以延轨和垂轨方向的局部坐标系之间的转换关系，分别在垂轨和延轨方向以2.5m步长调整目标剖面T。坐标转换关系如下：

其中，θ为激光剖面T与极地投影坐标系x轴的夹角，x_AC，y_AC是在局部坐标系下的坐标，x_SP，y_SP为激光点在极地投影系下的坐标；

5)采用反距离权重插值，根据100米范围内的最近邻十个点，插值出每次调整位置后的高程值z_interp；

6)根据T的实际观测高程值与插值高程值z_interp计算加权的

其中v为高程差，p采用Huber权函数。

7)完成预设调整区域内的每一步调整，记录最小的RMSE值对应的平面位置调整值(Δx,Δy)及此时的高程差v_T。此时目标剖面T的原始平面坐标(x₀,y₀,z₀)更新为(x_T,y_T,z_T)，表达式为：

其中，Δz为v_T的加权平均值。

8)目标剖面T完成调整。对下一个激光剖面进行上述调整过程。

考虑到由于在自约束调整过程中参考的激光点数据也存在误差，完成调整后激光点数据会被更新，因此需要再进行上述流程的迭代过程。迭代的目的是为了使每轨激光剖面数据平面位置的调整值趋于0，即最终调整到了最佳的位置。一般情况下，所有数据在经过4-6轮迭代即可达到收敛。

3、异常点剔除

完成迭代调整后的激光数据中存在一些异常值。这些异常点对于地形数据产品影响很大，因为它们会导致DEM中出现一些异常的微小坑洞或突刺伪特征，因此剔除这些异常点是十分有必要的。对于单个分布的异常点，通常由于其高程值明显与周围地形不符，所以会导致在这些点地形坡度会发生突变。根据这一特征，首先利用调整后的激光点生成了5m分辨率的坡度图，然后用25像素窗口大小进行中值滤波处理。由此，可以计算每个点对应的去趋势坡度S：

S＝S₁-S₂

其中，S₁，S₂分别为每个点的局部坡度和滤波坡度。

去趋势坡度S可以表征出坡度变化大的点，但并不能代表坡度异常变化，因为地形的正常起伏也会导致坡度的变化，所以用粗糙度R对其进行标准化。粗糙度的计算方法有很多，这里采用的是坡度的中位数。因此可以构建统计量：

δ₁＝S/R＝(S₁-S₂)/R

除此之外，在自约束迭代调整过程中，对于最后一轮调整后得到的高程差v中较大的点进行剔除，因为这类点明显与周围地形存在很大的差异。剔除的方法是首先计算高程差v的绝对中位差MAD。将绝对值大于3*MAD的值挑选出来构成新的统计量δ₂。

最后，根据δ₁，δ₂两个统计量，按照分位值为0.001对分布在双尾的激光点进行剔除。

4、实验结果与讨论

利用本文提出的方法，对沙克尔顿撞击坑周边的数据进行了处理，并对结果进行展示和分析。图5展示了沙克尔顿撞击坑边缘的局部特写。其中，(a)是NASA官方发布的DEM产品，(b)是自约束迭代调整后的激光点得到的DEM，(c)是对(b)中的激光数据进行了异常点剔除后得到的DEM。

可以看出，经过本文提出的方法处理后的激光点，最大程度上消除了地理定位误差，同时剔除了异常点，因此，从DEM中明显看出原有的大量地形伪影被消除了。这些地形伪影会导致高程上存在数十米的误差。

为了验证结果的正确性，将本文的结果与NASA发布最新的改善后的DEM(https://pgda.gsfc.nasa.gov/products/78.)进行了定量比较分析。以其发布的16km*16km大小的site01区域为对象，该区域在月球南极立体投影坐标系下的范围为X:-19000～-3000米，Y:-4000～-20000米。首先，将使用本方法处理后的激光点投影至改善后的DEM上，计算了激光点高程值和DEM的高程值的平均绝对偏差(MAE)为0.25米和均方根误差(RMSE)为0.46米(图6，左)；同时，用本文处理后的激光点生成的DEM与其相减，统计了差值DEM像素值的分布(图6，右)。

通过图7地形剖面图可以看出，与原始的DEM产品(LDEM)相比，本文的DEM和改善的DEM一致性也更高，地形符合程度更好。

虽然与NASA改善DEM的一致性较高，但是也存在极个别高程相差较大的像素。所以，结合LROC NAC影像，我们对差异进行了分析。如图8所示，(a)为NASA改善的DEM，(b)为本文生成的DEM，A、B、C为明显差异的三个地形特征。具体描述为在NASA改善的DEM中，A、B两处都存在一个坑洞，C处存在一个突刺。这导致了在这些局部区域的高程差达到了数十米。图8中右上角的黑色虚线框中，第一列DEM在A、B、C处的放大图，第二列是地理配准后的LROCNAC影像。通过将剔除掉的异常点投影在图上可以发现，导致这种差异的主要原因是在利用本文方法对激光点进行处理时，这些特征对应的激光点被视为异常点，并非是真实的地形特征点，所以本文的DEM上没有这些坑洞和突刺。通过与LROC NAC影像进行目视比较，在这些目标点对应的位置，并没有发现相对应的如此大尺度的地形特征。

同时，生成跨越这些特征的地形剖面(图8)中可以看到，在存在显着差异的A和B两处，本发明的DEM和LDEM的地形趋势更加一致。至于C中的突刺似乎更像是一种伪特征，因为它分布在陨石坑光滑的内壁上，似乎不符合地形的总体趋势。

最后，以南极沙克尔顿撞击坑的部分区域为例，使用本文所提出的方法对该区域的激光数据进行调整后，再生成地形图。结果如图9所示，左图为由LDEM产品生成渲染地形图，右图为本文DEM生成的渲染地形图。

综上，针对月球星载激光测高数据的特点，本发明提出了一种基于月球星载激光测高数据的永久阴影区域地形测绘方法。首先采用自约束迭代调整的方法对研究区域内的所有LOLA激光点逐轨进行调整，在调整过程中不借助外部高精度参考地形数据，而是以LOLA邻近分布数据之间的几何约束，整体调整LOLA测高剖面数据，消除激光点的地理定位误差，结果在DEM地形数据上表现为消除了大量地形伪影。同时，标准化去趋势坡度和最后一次迭代的高程残差用于构建统计量，以过滤掉一些会导致DEM中异常突起或坑洞的异常值。在结果中，以沙克尔顿撞击坑为例，展示了本文方法处理后的RDR数据集生成的月球南极永久阴影区域的局部地形图。

以上所述，仅为本发明的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明揭露的技术范围内，可轻易想到各种等效的修改或替换，这些修改或替换都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应以权利要求的保护范围为准。

Claims (10)

1.一种月球永久阴影区域地形测绘方法，其特征在于，该方法包括：

步骤S1、对月球轨道器激光高度计LOLA的每个激光剖面在其接受视场范围内进行定步长调整，以单轨激光数据邻域内其余的激光观测值作为约束，最小化插值高程和观测高程差值的加权均方根误差RMSE，实现激光剖面的逐轨自约束调整；

步骤S2、将调整后的激光剖面更新为初始状态数据进行迭代，不断收敛平面位置调整值；

步骤S3、对迭代调整后的激光数据进行异常值剔除，得到误差修正后的月球激光测高数据；

步骤S4、基于误差修正后的月球激光测高数据，得到表征月球永久阴影区域三维地形的数字高程模型。

2.根据权利要求1所述的一种月球永久阴影区域地形测绘方法，其特征在于，所述步骤S1包括以下子步骤：

步骤S101、对LOLA激光点进行数据预处理；

步骤S102、选取单轨RDR数据作为目标剖面T，根据其分布范围，各方向向外扩展k₁米，从剩余激光点中裁剪出作为参考地形数据的参考激光点集R；

步骤S103、以参考激光点集R生成用于进行最近邻搜索的KD树；

步骤S104、建立半径为k₂米的预设调整区域，调整目标剖面；

步骤S105、采用反距离权重插值，根据k₃米范围内的最近邻k₄个点，插值出每次调整位置后的高程值z_interp；

步骤S106、根据目标剖面T的实际观测高程值与插值高程值z_interp计算加权的均方根误差RMSE，表达式为：

其中，v为高程差，p采用Huber权函数；

步骤S107、完成预设调整区域内的每一步调整，记录最小的RMSE值对应的平面位置调整值(Δx，Δy)以及高程差v_T，更新目标剖面T的平面坐标(x_T，y_T，z_T)，表达式为：

其中，(x₀，y₀，z₀)为原始平面坐标，Δz为v_T的加权平均值；

步骤S108、完成当前单个目标剖面T调整后，进行下一个激光剖面的自约束迭代调整。

3.根据权利要求2所述的一种月球永久阴影区域地形测绘方法，其特征在于，所述步骤S101中数据预处理包括噪声剔除、投影变换以及逐轨存储。

4.根据权利要求2所述的一种月球永久阴影区域地形测绘方法，其特征在于，所述步骤S104具体为：基于极地投影坐标系与以延轨和垂轨方向的局部坐标系之间的转换关系，分别在垂轨和延轨方向以步长k₅调整目标剖面T，坐标转换表达式为：

式中，θ为激光剖面T与极地投影坐标系x轴的夹角，x_AC、y_AC是在局部坐标系下的坐标，x_SP、y_SP为激光点在极地投影系下的坐标。

5.根据权利要求4所述的一种月球永久阴影区域地形测绘方法，其特征在于，所述k₁、k₂、k₃、k₄和k₅分别设置为150、50、100、10和2.5。

6.根据权利要求2所述的一种月球永久阴影区域地形测绘方法，其特征在于，所述步骤S106中加权的均方根误差RMSE表达式为：

其中，v为高程差，p采用Huber权函数，n为激光点数量。

7.根据权利要求2所述的一种月球永久阴影区域地形测绘方法，其特征在于，所述步骤S3具体为：构建标准化去趋势坡度图和残差统计量，剔除伪地形观测的异常值。

8.根据权利要求7所述的一种月球永久阴影区域地形测绘方法，其特征在于，所述步骤S3包括以下子步骤：

步骤S301、利用调整后的激光点生成预设分辨率的坡度图，采用预设像素窗口大小进行中值滤波处理后，计算每个点对应的去趋势坡度S：

S＝S₁-S₂

式中，S₁，S₂分别为每个点的局部坡度和滤波坡度；

步骤S302、采用粗糙度R对去趋势坡度S标准化，构建统计量δ₁：

δ₁＝S/R＝(S₁-S₂)/R

步骤S303、计算自约束迭代调整过程中最后一轮调整后得到的高程差v的绝对中位差MAD，将绝对值大于K倍MAD的值挑选出来构成新的统计量δ₂；

步骤S304、根据统计量δ₁和δ₂，按照预设分位值对分布在双尾的激光点进行剔除。

9.根据权利要求8所述的一种月球永久阴影区域地形测绘方法，其特征在于，所述步骤S302中的粗糙度R为坡度中位数。

10.根据权利要求8所述的一种月球永久阴影区域地形测绘方法，其特征在于，所述步骤S303中的K为3。

Images