CN114964176B - Topography mapping method for moon permanent shadow area - Google Patents

Abstract

The invention relates to a topographic mapping method for a moon permanent shadow area, which comprises the following steps: step S1, performing fixed step adjustment on each laser section of a lunar orbit laser altimeter LOLA within the range of an accepted view field, taking the rest laser observation values in a monorail laser data neighborhood as constraints, minimizing the weighted Root Mean Square Error (RMSE) of interpolation elevation and observation elevation difference values, and realizing track-by-track self-constraint adjustment of the laser sections; s2, updating the adjusted laser profile into initial state data for iteration, and continuously converging the plane position adjustment value; s3, removing abnormal values from the laser data after iterative adjustment to obtain moon laser altimetry data after error correction; and S4, obtaining a digital elevation model representing the three-dimensional terrain of the permanent shadow area of the moon based on the moon laser altimetry data after error correction. Compared with the prior art, the method eliminates a large amount of topographic artifacts in the original image, and improves the topographic mapping quality of the moon permanent shadow area.

Description

Topography mapping method for moon permanent shadow area

Technical Field

The invention relates to the technical field of lunar surveying and mapping, in particular to a lunar permanent shadow area topographic surveying and mapping method.

Background

Due to the slight inclination of the axis of rotation of the moon, it is difficult for sunlight to reach directly to certain areas of the lunar polar region, which are known as permanent shadow areas (PSRs). The lunar laser altimeter can realize fine mapping of lunar topography, and particularly has the advantage of active mapping on a lunar Permanent Shadow Region (PSRs), which provides basic topography data for planetary detection tasks such as Liu Xuanzhi, lander line and water ice detection. However, various factors such as spacecraft orbit control, laser pointing and hardware faults can lead to uncertainty of the geographical position of a laser spot, which can cause errors in laser data and a great deal of artifacts in three-dimensional terrain products. The above-mentioned problems also exist with lunar orbit machine laser altimeter (LOLA) data issued by the national astronaut agency (NASA) as the highest precision laser altimeter data on the moon by now.

The Digital Elevation Model (DEM) can reflect the three-dimensional topography of the lunar surface, is the basis and key of lunar field exploration, can help human beings and inspection devices to clearly see the surrounding environment and identify important features, and provides references for exploration task planning and decision making. Unlike optical remote sensing, which uses the sun as a light source, laser imaging is independent of external light and can map PSRs topography with its own light source.

Therefore, the lunar orbit device is almost all equipped with laser load for high-precision topographic mapping and lunar surface control. Although the laser observation has the advantage of high precision, errors exist in the laser measurement due to uncertainty of factors such as track reconstruction quality, laser pointing deviation and the like, and the errors can be directly mapped into high-resolution terrain products to present terrain artifacts distributed in accordance with the laser footprint tracks. Thus, geolocation calibration of laser footprints is particularly important for laser high precision mapping.

The method for correcting the geographic positioning error of the laser footprint is a cross-track method. The crossed rail is mainly used for adjusting the laser footprint position by minimizing the distance between laser observation data of different rails at the crossing point; the direct height measurement method is based on priori high-precision terrain data, and each laser observation is optimally matched with the priori terrain, so that the purpose of correcting errors is achieved. However, for lunar laser altimetry data, the LOLA is taken as the current highest-precision lunar apparent survey data, and priori topographic data with higher precision than the lunar apparent survey data cannot be obtained; and the method of cross-track is computationally intensive for dense observer zone intersections and the geometric constraints are weak.

Barker et al modified the LOLA DEM product in a self-constraining manner using the reduced LOLA DEM as reference terrain data, but this method uses the generation of all remaining points to reduce the DEM when adjusting the laser profile in batches, which makes the area away from the target profile time consuming and redundant when constructing the DEM with laser points. In addition, if the spatial distribution of the laser profiles of the same batch is very dense, the points used for generating the reference DEM are very sparse, so that the generated reduced DEM error is larger.

In view of the above, it is desirable to devise a new method for mapping the topography of a permanently shaded area of the moon to correct errors in the laser spot, thereby accurately mapping the three-dimensional topography of the permanently shaded area to eliminate a significant amount of topographical artifacts in the original product.

Disclosure of Invention

The invention aims to overcome the defects of the prior art and provide a method for mapping the topography of the permanently shaded area of the moon, which eliminates the geographic positioning error to the greatest extent and eliminates the abnormal points at the same time, thereby eliminating a great amount of topography artifacts.

The aim of the invention can be achieved by the following technical scheme:

the invention provides a topographic mapping method of a moon permanent shadow area, which comprises the following steps:

step S1, performing fixed step adjustment on each laser section of a lunar orbit laser altimeter LOLA within the range of an accepted view field, taking the rest laser observation values in a monorail laser data neighborhood as constraints, minimizing the weighted Root Mean Square Error (RMSE) of interpolation elevation and observation elevation difference values, and realizing track-by-track self-constraint adjustment of the laser sections;

s2, updating the adjusted laser profile into initial state data for iteration, and continuously converging the plane position adjustment value;

s3, removing abnormal values from the laser data after iterative adjustment to obtain moon laser altimetry data after error correction;

and S4, obtaining a digital elevation model DEM representing the three-dimensional topography of the permanent shadow area of the moon based on the moon laser altimetry data after error correction.

Preferably, the step S1 comprises the following sub-steps:

step S101, data preprocessing is carried out on LOLA laser points;

step S102, selecting monorail RDR data as a target section T, and expanding k outwards in all directions according to the distribution range ₁ A step of cutting out a reference laser point set R serving as reference topographic data from the rest laser points;

step S103, generating a KD tree for nearest neighbor searching by using the reference laser point set R;

step S104, establishing a radius of k ₂ A preset adjustment area of rice, and adjusting a target profile;

step S105, adopting inverse distance weight interpolation according to k ₃ Nearest neighbor k in meter range ₄ A plurality of points, and the elevation value z after each position adjustment is interpolated _interp ；

Step S106, according to the actual observed elevation value and the interpolated elevation value z of the target profile T _interp And calculating a weighted Root Mean Square Error (RMSE) with the expression:

wherein v is the elevation difference, and p is a Huber weight function;

step S107, finishing each step of adjustment in the preset adjustment area, and recording the plane position adjustment value (Deltax, deltay) and the elevation difference v corresponding to the minimum RMSE value _T Updating plane coordinates (x _T ,y _T ,z _T ) The expression is:

wherein, (x) ₀ ,y ₀ ,z ₀ ) For the original plane coordinates, Δz is v _T Is a weighted average of (2);

step S108, after the current single target profile T is adjusted, the self-constrained iterative adjustment of the next laser profile is performed.

Preferably, the data preprocessing in step S101 includes noise rejection, projective transformation and track-by-track storage.

Preferably, the step S104 specifically includes: based on the conversion relation between polar projection coordinate system and local coordinate system in the directions of vertical and vertical tracks, the step length k is respectively used in the vertical and vertical track directions ₅ Adjusting a target section T, wherein the coordinate conversion expression is as follows:

wherein θ is the angle between the laser profile T and the x-axis of the polar projection coordinate system, x _aC 、y _aC Is the coordinates in the local coordinate system, x _SP 、y _SP Is the coordinates of the laser spot in the polar projection system.

Preferably, said k ₁ 、k ₂ 、k ₃ 、k ₄ And k ₅ Set to 150, 50, 100, 10 and 2.5, respectively.

Preferably, the root mean square error RMSE expression weighted in step S106 is:

wherein v is the elevation difference, p is Huber weight function, and n is the number of laser points.

Preferably, the step S3 specifically includes: and constructing a standardized trending slope map and residual statistics, and eliminating abnormal values of pseudo-topography observation.

Preferably, said step S3 comprises the following sub-steps:

step 301, generating a gradient map with preset resolution by using the adjusted laser points, and calculating a trending gradient S corresponding to each point after median filtering processing is performed by adopting the size of a preset pixel window:

S＝S ₁ -S ₂

wherein S is ₁ ，S ₂ Local to each point respectivelyGradient and filtering gradient;

step S302, adopting roughness R to normalize the trending slope S, and constructing statistics delta ₁ ：

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

Step S303, calculating the absolute median difference MAD of the elevation difference v obtained after the last round of adjustment in the self-constrained iterative adjustment process, and selecting the value of MAD with the absolute value larger than K times to form a new statistic delta ₂ ；

Step S304, according to statistics delta ₁ And delta ₂ And eliminating the laser points distributed at the two tails according to a preset quantile value.

Preferably, the roughness R in the step S302 is a median gradient.

Preferably, K in step S303 is 3.

Compared with the prior art, the invention has the following advantages:

1) According to the invention, a self-constraint iterative adjustment method is adopted to adjust all LOLA laser points in a research area track by track, external high-precision reference topographic data is not needed in the adjustment process, LOLA height measurement profile data are integrally adjusted through geometric constraints among LOLA adjacent distribution data, geographic positioning errors of the laser points are eliminated, a large amount of topographic artifacts in DEM topographic data are eliminated, and the three-dimensional topographic map quality of a moon permanent shadow area is improved;

2) The invention uses the normalized trend-removing gradient and the last iteration elevation residual error to construct statistics to filter out abnormal values which can cause abnormal protrusions or pits in the DEM.

Drawings

FIG. 1 is a flow chart of a track-by-track self-constrained iterative adjustment of the present invention;

FIG. 2 is a diagram showing experimental data of the present invention; wherein, (a) is LOLA RDR monorail laser profile data, (b) is LOLA DEM with 5m resolution, (c) is DEM with 5m resolution improvement issued by NASA, (d) is LROC NAC image with 1.5 m resolution;

FIG. 3 is a schematic illustration of monorail self-constraining adjustment;

FIG. 4 is a graph of the laser spot plane position adjustment values for each track in each iteration;

FIG. 5 is a DEM results display;

FIG. 6 is a graph of comparative analysis of results;

FIG. 7 is a topographical cross-sectional view; the left image is a global terrain section, and the right image is a local enlarged image;

FIG. 8 is a graphical representation of the results of NASA-improved DEM and DEM differential characterization of the present invention;

fig. 9 is a (partial) three-dimensional terrain rendering of a sarton impingement pit.

Detailed Description

The following description of the embodiments of the present invention will be made clearly and fully with reference to the accompanying drawings, in which it is evident that the embodiments described are some, but not all embodiments of the invention. All other embodiments, which can be made by those skilled in the art based on the embodiments of the present invention without making any inventive effort, shall fall within the scope of the present invention.

Examples

The embodiment provides a method for surveying and mapping the topography of a moon permanent shadow area, which comprises the following steps:

1. obtaining RDR data products and LDEM data distributed track by track;

the research object of the invention is LOLA laser altimetry data. The LOLA data science team performs hierarchical processing and release on the data products, wherein EDR (Experiment Data Records) is original scientific observation data, and RDR (Reduced Data Records) data products are generated after calculation and calibration.

The laser observation values in the RDR data product are arranged according to time to form monorail laser altimetry data containing information such as time stamps, laser footprint three-dimensional positions, laser pulse characteristic parameters and the like. After median filtering and spatial interpolation processing are carried out on laser points in RDR, a Lunar Digital Elevation Model (LDEM) data product is generated, and a large amount of topographic artifacts exist in the product.

According to the invention, RDR data products distributed track by track are used as experimental data, LDEM data are used as comparison analysis data, meanwhile, improved DEM data newly issued by NASA are used as precision evaluation data, and LROC (Lunar Reconnaissance Orbiter Camera) NAC (Narrow Angle Camera) images with higher resolution are used as auxiliary data for visually judging the authenticity of individual topographic features.

2. Track-by-track self-constrained iterative adjustment of monorail RDR laser profiles

And for all laser points in the research area, taking a single-rail RDR laser profile as an object, and realizing self-restraint adjustment one by one. The main idea is as follows: the geometrical space constraint exists between the monorail laser profile and the laser observation points distributed nearby, the position of the target profile which is best fit with the terrain in the neighborhood range can be found by utilizing the relation, and the criterion of the fitting degree is measured by using weighted Root Mean Square Error (RMSE). As shown in fig. 3 (a), gray points represent laser points (projected on a two-dimensional plane) distributed in the neighborhood as reference terrain, other points represent spatial distribution of the target profile corresponding to three different plane positions, corresponding error bars represent residuals of the observed elevation values and the interpolated elevations, and fig. 3 (b) is a distribution diagram corresponding to the elevation residuals. Obviously, the B section is more consistent with the terrain, so the corresponding plane position is closer to the true position of the target section.

The main implementation flow of the method is shown in figure 1. The detailed steps are as follows:

1) All laser spots are preprocessed. Noise rejection and projection transformation are included, and then track-by-track storage is performed;

2) Selecting single-rail RDR data as a target section T, and expanding 150 meters outwards in each direction according to the distribution range of the single-rail RDR data to cut out a laser point set R serving as reference topographic data from the rest laser points;

3) Generating a kd-tree by using the reference point set R for nearest neighbor searching;

4) And establishing a preset adjustment area with the radius of 50 m. And establishing a conversion relation between the polar position projection coordinate system and a local coordinate system in the directions of the vertical track and the vertical track, and respectively adjusting the target section T in the directions of the vertical track and the vertical track by 2.5m steps. The coordinate conversion relationship is as follows:

wherein θ is the angle between the laser section T and the x-axis of the polar projection coordinate system, x _AC ，y _AC Is the coordinates in the local coordinate system, x _SP ，y _SP Coordinates of the laser point under a polar projection system;

5) Interpolation of inverse distance weight is adopted, and the elevation value z after each position adjustment is interpolated according to the nearest ten points in the range of 100 meters _interp ；

6) According to the actual observed elevation value and the interpolated elevation value z of T _interp Calculating the weighting

Where v is the elevation difference and p employs a Huber weight function.

7) Each step of adjustment in the preset adjustment area is completed, and the plane position adjustment value (delta x, delta y) corresponding to the minimum RMSE value and the elevation difference v at the moment are recorded _T . The original plane coordinates (x ₀ ,y ₀ ,z ₀ ) Update to (x) _T ,y _T ,z _T ) The expression is:

wherein Δz is v _T Is a weighted average of (c).

8) The target profile T is adjusted. The adjustment process described above is performed for the next laser profile.

In consideration of errors in the laser spot data referenced in the self-constrained adjustment process, the laser spot data after the adjustment is completed is updated, so that the iterative process of the above procedure is needed. The purpose of the iteration is to bring the adjustment value of the position of the data plane of the laser profile per track towards 0, i.e. eventually to the optimal position. In general, all data can reach convergence after 4-6 iterations.

3. Outlier rejection

Some abnormal values exist in the laser data after the iterative adjustment is completed. These outliers are highly influential in the topographical data products because they can lead to some unusual micro-pits or spike pseudo-features in the DEM, and thus, it is necessary to eliminate them. For a single distribution of outliers, abrupt changes in the slope of the terrain at these points are typically caused by significant differences in elevation values from the surrounding terrain. According to this feature, a gradient map of 5m resolution is first generated using the adjusted laser spot, and then median filtering processing is performed using a 25 pixel window size. Thus, the detrack gradient S corresponding to each point can be calculated:

S＝S ₁ -S ₂

wherein S is ₁ ，S ₂ The local gradient and the filtered gradient for each point, respectively.

Trending slope S may characterize a point of great slope change but does not represent an abnormal slope change, as normal fluctuations in terrain may also result in a change in slope, which is normalized by roughness R. There are many methods for calculating the roughness, and the median of the gradient is used here. Statistics can thus be constructed:

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

in addition, in the self-constrained iterative adjustment process, larger points in the elevation difference v obtained after the last round of adjustment are eliminated, because the points obviously have great differences from the surrounding terrain. The method of rejection is to first calculate the absolute median difference MAD of the elevation difference v. Selecting values with absolute values greater than 3 MAD to form new statistic delta ₂ 。

Finally, according to delta ₁ ，δ ₂ And the two statistics are used for eliminating the laser points distributed at the two tails according to the quantile value of 0.001.

4. Experimental results and discussion

Data around the periphery of the sarton impingement pit is processed using the methods presented herein and the results are displayed and analyzed. Fig. 5 shows a partial close-up of the shakerton impingement pit edge. Wherein, (a) is a DEM product issued by NASA official, (b) is a DEM obtained by self-constrained iterative adjustment of laser points, and (c) is a DEM obtained by removing abnormal points from laser data in (b).

It can be seen that the laser points processed by the method provided herein eliminate geographical positioning errors to the greatest extent, and meanwhile eliminate abnormal points, so that it is obvious from DEM that a great deal of original topographic artifacts are eliminated. These topographical artefacts can lead to tens of meters of error in elevation.

To verify the correctness of the results, the results herein were quantitatively compared with the latest improved DEM issued by NASA (https:// pgda. Gsfc. NASA gov/products/78.). Taking a site01 area with the size of 16km which is published by the site01 area as an object, wherein the range of the area under a moon south pole stereoscopic projection coordinate system is X-19000-3000 m, Y-4000-20000 m. Firstly, projecting the laser point processed by the method onto an improved DEM, and calculating an average absolute deviation (MAE) of the elevation value of the laser point and the elevation value of the DEM to be 0.25 m and a Root Mean Square Error (RMSE) to be 0.46 m (left in fig. 6); meanwhile, the DEM generated by the laser spot processed herein is subtracted therefrom, and the distribution of the difference DEM pixel values is counted (fig. 6, right).

As can be seen from the topographical profile of fig. 7, the DEM and improved DEM herein also have higher consistency and better topographical compliance than the original DEM product (LDEM).

Although the consistency with NASA to improve DEM is higher, there are also very individual pixels with a large difference in elevation. Therefore, we analyzed the differences in combination with LROC NAC imaging. As shown in fig. 8, (a) is a NASA-improved DEM, (b) is a DEM generated herein, A, B, C is three topographical features that differ significantly. Specifically described as a NASA-modified DEM, a pit is present at both A, B and a spur is present at C. This results in elevation differences in these local areas of up to tens of meters. In the upper right-hand black dashed box of fig. 8, the first column DEM is an enlarged view at A, B, C and the second column is a geo-registered LROC NAC image. The principal cause of this discrepancy is found by projecting the removed outliers onto the map, as the laser points corresponding to these features are considered outliers and not true topographical features when processed using the methods herein, so these pits and spikes are not present on the DEM herein. By visual comparison with LROC NAC images, corresponding topographical features of such large scale are not found at the locations corresponding to these target points.

At the same time, it can be seen in the generation of topographical profiles (FIG. 8) across these features that the topographical trends of the DEM and LDEM of the present invention are more consistent where there are significant differences in both A and B. As for the spike in C, it appears to be a pseudo feature because it is distributed over the smooth inner wall of the merle pit, which appears to be a general trend out of topography.

Finally, taking a partial area of the antarctic sarton impingement pit as an example, the laser data of the area is adjusted by using the method disclosed herein, and then a topographic map is regenerated. As a result, as shown in fig. 9, the left graph is a rendered topographical map generated by the LDEM product, and the right graph is a rendered topographical map generated by the DEM herein.

In summary, the invention provides a permanent shadow region topography mapping method based on lunar satellite-borne laser altimetry data aiming at the characteristics of lunar satellite-borne laser altimetry data. Firstly, a self-constrained iterative adjustment method is adopted to adjust all LOLA laser points in a research area track by track, the LOLA height measurement profile data is integrally adjusted by geometric constraints among LOLA adjacent distribution data without external high-precision reference topographic data in the adjustment process, the geographic positioning error of the laser points is eliminated, and a large amount of topographic artifacts are eliminated on DEM topographic data as a result. At the same time, the normalized detritus slope and the last iteration elevation residual are used to construct statistics to filter out outliers that may lead to outlier protrusions or pits in the DEM. In the results, taking the example of a shackleton impingement pit, a local topography of the moon's south pole permanent shadow region generated from the RDR dataset processed by the method herein is shown.

While the invention has been described with reference to certain preferred embodiments, it will be understood by those skilled in the art that various changes and substitutions of equivalents may be made and equivalents will be apparent to those skilled in the art without departing from the scope of the invention. Therefore, the protection scope of the invention is subject to the protection scope of the claims.

Claims (9)

1. A method of topographic mapping a permanently shaded area of the moon, the method comprising:

step S1, performing fixed step adjustment on each laser section of a lunar orbit laser altimeter LOLA within the range of an accepted view field, taking the rest laser observation values in a monorail laser data neighborhood as constraints, minimizing the weighted Root Mean Square Error (RMSE) of interpolation elevation and observation elevation difference values, and realizing track-by-track self-constraint adjustment of the laser sections;

s2, updating the adjusted laser profile into initial state data for iteration, and continuously converging the plane position adjustment value;

s3, removing abnormal values from the laser data after iterative adjustment to obtain moon laser altimetry data after error correction;

and S4, obtaining a digital elevation model representing the three-dimensional terrain of the permanent shadow area of the moon based on the moon laser altimetry data after error correction.

2. A method of topographic mapping a permanently shaded area of the moon according to claim 1, wherein said step S1 comprises the sub-steps of:

step S101, data preprocessing is carried out on LOLA laser points;

step S102, selecting monorail RDR data as a target section T, and expanding k outwards in all directions according to the distribution range ₁ A step of cutting out a reference laser point set R serving as reference topographic data from the rest laser points;

step S103, generating a KD tree for nearest neighbor searching by using the reference laser point set R;

step S104, establishing a radius of k ₂ A preset adjustment area of rice, and adjusting a target profile;

step S105, interpolation by inverse distance weightAccording to k ₃ Nearest neighbor k in meter range ₄ A plurality of points, and the elevation value z after each position adjustment is interpolated _interp ；

Step S106, according to the actual observed elevation value and the interpolated elevation value z of the target profile T _interp And calculating a weighted Root Mean Square Error (RMSE) with the expression:

wherein v is the elevation difference, and p is a Huber weight function;

step S107, finishing each step of adjustment in the preset adjustment area, recording the plane position adjustment value (Deltax, deltay) and the elevation difference v corresponding to the minimum RMSE value, and updating the plane coordinate (x) of the target profile T _T ,y _T ,z _T ) The expression is:

wherein, (x) ₀ ,y ₀ ,z ₀ ) For the original plane coordinates, Δz is a weighted average of v;

step S108, after the current single target profile T is adjusted, the self-constrained iterative adjustment of the next laser profile is performed.

3. The method according to claim 2, wherein the data preprocessing in step S101 includes noise rejection, projective transformation, and track-by-track storage.

4. The method of claim 2, wherein the step S104 is specifically: based on the conversion relation between polar projection coordinate system and local coordinate system in the directions of vertical and vertical tracks, the step length k is respectively used in the vertical and vertical track directions ₅ Adjusting a target section T, wherein the coordinate conversion expression is as follows:

/>

wherein θ is the angle between the laser profile T and the x-axis of the polar projection coordinate system, x _AC 、y _AC Is the coordinates in the local coordinate system, x _SP 、y _SP Is the coordinates of the laser spot in the polar projection system.

5. The method of claim 4, wherein k is ₁ 、k ₂ 、k ₃ 、k ₄ And k ₅ Set to 150, 50, 100, 10 and 2.5, respectively.

6. The method for topographic mapping of permanently shaded areas of the moon according to claim 2, wherein said step S3 is specifically: and constructing a standardized trending slope map and residual statistics, and eliminating abnormal values of pseudo-topography observation.

7. The method of claim 6, wherein said step S3 comprises the sub-steps of:

step 301, generating a gradient map with preset resolution by using the adjusted laser points, and calculating a trending gradient S corresponding to each point after median filtering processing is performed by adopting the size of a preset pixel window:

S＝S ₁ -S ₂

wherein S is ₁ ，S ₂ Local gradient and filtering gradient of each point respectively;

step S302, adopting roughness R to normalize the trending slope S, and constructing statistics delta ₁ ：

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

Step S303, calculating the absolute median difference MAD of the elevation difference v obtained after the last round of adjustment in the self-constrained iterative adjustment process, and performing absolute adjustmentSelecting the value of MAD with the value larger than K times to form new statistic delta ₂ ；

Step S304, according to statistics delta ₁ And delta ₂ And eliminating the laser points distributed at the two tails according to a preset quantile value.

8. The method according to claim 7, wherein the roughness R in the step S302 is a median gradient.

9. The method according to claim 7, wherein k in the step S303 is 3.

Images