LOLA DEM Assisted Photogrammetric Processing of Lunar Reconnaissance Orbiter NAC Images for the Lunar South Pole

Because of its special significance, exploration missions and research at the Lunar South Pole (LSP) are of considerable importance. The mapping products of the LSP such as Digital Ortho Maps (DOMs) and Digital Elevation Models (DEMs) are essential to support exploration missions. However, in the LSP region, the available data are relatively limited. The large amounts of shadows and the special environment of the LSP has an important impact on both image matching and generation of mapping products. In this paper, we propose a photogrammetric processing method to improve the accuracy of the LRO NAC images by using LOLA DEM as auxiliary control data. The LOLA DEM was first rendered using the same illumination conditions as the corresponding LRO NAC raw images, and then matched to the LRO NAC DOMs. Next, the corresponding elevations of the tie points between LOLA DEMs and LRO NAC images were extracted and automatically added to the control network for bundle adjustment. In order to verify the feasibility of the method, we selected 17 LRO NAC images near the LSP at 85 ° south latitude. The results show that the proposed method can efficiently acquire uniformly distributed elevation control points for bundle adjustment and the generated DEM has a good consistency with the LOLA DEM.


Introduction
According to research, there may be water ice in the cold traps in the permanently shadowed regions (PSRs) of the Lunar South Pole (LSP) (Biswas et al., 2020;Li et al., 2018b).Water is an important substance for the origin of life and contributes to the study of the evolutionary processes of the Moon.In addition, water is a direct aid to in situ exploration and the establishment of scientific research stations at the lunar south pole.In recent years, research and exploration missions about the LSP have also become increasingly hot.The fourth phase of China's Lunar Exploration Program will carry out focused research on the South Pole, including in-situ research and utilization of South Pole resources, as well as research on the distribution of water and volatile components (Yu et al., 2023;Zou et al., 2020).The U.S. Artemis program also plans to land near the LSP and conduct related research (Evans et al., 2020).Other countries, such as India and Japan, also have corresponding plans to carry out exploration missions to the LSP (Ohtake et al., 2021).High precision positioning and navigation maps can be used to select the landing area and ensure the smooth landing of the lander, which plays an important supporting role in the exploration mission of the LSP.
Photogrammetry is a mature method widely used in the production of digital orthophotos map (DOM) and digital elevation model (DEM) of remote sensing images for deep space exploration missions.However, the LSP has complex terrain and harsh illumination conditions (Cohen et al., 2020), which have a significant impact on the production of key mapping products such as DOM and DEM.Specifically, the numerous impact craters at the LSP have resulted in rugged terrain and drastic changes in the slope of the lunar surface.The location in the polar region leads to an extremely low solar altitude angle, resulting in a nearly 90 degrees angle of solar incidence (Zhang et al., 2023).

* Corresponding author
The excessively large incidence angle and rugged terrain ultimately result in a wide and discontinuous distribution of shadow that can be seen everywhere.These discontinuous shadows severely reduce the range of available image and make the image discontinuous, which can have a critical impact on image matching.The result of image matching directly affects the quality of bundle adjustment, ultimately leading to poor positioning accuracy.At present, there are different advantages and disadvantages in various types of lunar mapping products (e.g., spatial resolution, lunar coverage, positioning accuracy and three-dimensional coverage) due to the specific scientific objectives of their respective exploration missions.Among the more commonly used lunar source data, the Chang'E-2 image with 7 m/pixel resolution has a stereoscopic viewpoint, which enables three-dimensional coverage of the photographed area, and has achieved full lunar surface coverage.However, due to its short in-orbit shooting time, and due to a large number of shadows in the lunar polar regions, it is actually difficult to utilize most of the images at the location of the LSP.The Japanese Kaguya Terrain Camera (TC) instrument does not have complete coverage at the LSP and some low and middle latitudes.
Lunar Reconnaissance Orbiter (LRO) Narrow Angle Camera (NAC) images have a resolution of 0.5~1.5 m/pixel, which is the highest resolution remote sensing image of the Moon at present (Robinson et al., 2010).The LRO was launched in 2009, and after more than ten years of in-orbit photography, the images returned by LRO NAC has basically covered the entire moon, with over two million images.Designed for focused observations of the LSP, LRO achieved multiple coverage of the LSP region over a long time series.However, there are some difficulties accompanying LRO, which belongs to the long strip linear array pushbroom remote sensing images, and the high-resolution leads to its huge data volume.The linear array remote sensing image and large data volume lead to the complexity of its camera model, which is more difficult and time-consuming in photogrammetric processing.The Lunar Orbiter Laser Altimeter (LOLA) DEM is generated from the elevation points recorded by the laser altimeter carried by the LRO (Smith et al., 2011).The LOLA achieves a nominal measurement accuracy of 10 cm (Smith et al., 2010a).It has good positioning accuracy and a large amount of dense coverage near the LSP (Smith et al., 2010b).A key reason for the low positioning accuracy of lunar remote sensing mapping products is the lack of absolute control points in the lunar surface.The only absolute control data currently available are the five lunar laser ranging retroreflector (LRRR) points placed on the lunar surface by exploration missions of the former Soviet Union and the United States.However, the locations of these reflectors are all distributed in the lunar sea region on the front side of the Moon, and there is no available absolute control data in the LSP region.LOLA DEM is regarded as a better choice for the current control data because of its good ranging and positioning accuracy.Therefore, LOLA DEM is also widely used as control data in photogrammetric processing of lunar images (Barker et al., 2015).(Xie et al., 2022;Barker et al., 2023) focused on the PSRs of the LSP, and they targeted and optimized the accuracy of the LOLA DEM to attenuate the streaking problem in these regions.To solve the weak convergence problem of LRO NAC images, (Chen et al., 2024) added a reference DEM to the bundle adjustment model as a terrain constraint, and outliers are eliminated using iterative reweighting, resulting in a large-scale bundle adjustment with a relative positioning accuracy of less than 1 m, and the inter-image inconsistency is reduced to 0.5 pixels.(Barker et al., 2021) optimized LOLA laser altimetry trajectory data by iterative bundle adjustment method, generated new DEMs of four landing zones, reduced geo-positioning error by more than ten times after optimization, and investigated height inconsistency and its error propagation law.(Xin et al., 2018) used affine scale-invariant feature transformations, rational function models, and least-squares to co-registrate rendered SLDEM2015 with LRO NAC images.(Wu et al., 2017) improved the aperture offset of LRO NAC images based on the geometric model by least-squares adjustment using the triplematching tie points between NAC images, which reduced the target-space inconsistency to the meter level.(Hu and Wu, 2018) proposed a method of block adjustment and coupled epipolar rectification for NAC stereo images, which merges the LRO NAC left and right images into one stereo model in the disparity space, eliminating the internal inconsistency in the disparity.The internal inconsistency in the bundle adjustment is eliminated.
In the case of photogrammetric processing or bundle adjustment of control network for large scale or even global imagery, there are more problems leading to worse accuracy (Di et al., 2017).Using other data as a reference for global bundle adjustment is one way (Di et al., 2018;Archinal et al., 2023).Global remote sensing images typically consist of years of data, and it's difficult to correct their errors to the same level due to the complicated error sources between images (Geng et al., 2021).Alternatively, error control can be achieved by analysing the internal and external orientation elements (Ostrach et al., 2023), and building the corresponding fitting model (Ren et al., 2019).In summary, the lack of control data is currently a major issue in the LSP mapping.The complex terrain and poor illumination also make matching a challenge in photogrammetric processing.In this paper, we propose a LOLA DEM assisted photogrammetric processing method of LRO NAC images for the LSP, which helps to improve the positioning accuracy of mapping products.
The manuscript is organized in the following way: the second section is the proposed method, which describes the construction of the initial control network, the extraction of control points, the introduction of control points and the generation of the final DEM.The third section is the experimental part to validate the proposed method, including the accuracy assessment of the DOM and DEM products.The fourth section shows the summary of our work.

LOLA DEM Assisted Photogrammetric Processing Method
This method proposed in this paper combines the high positioning accuracy of LOLA DEM and the high-resolution characteristics of LRO NAC.Firstly, match the rendered LOLA DEM with the DOM of LRO NAC.Secondly, convert the matching tie points into 3D coordinates and automatically add them to the initial control network.Finally, higher positioning accuracy is achieved through further adjustment.Figure 1 is a flowchart of LOLA DEM assisted photogrammetric processing method.The proposed method is specifically divided into the following steps: (1) The LRO NAC images are processed with conventional photogrammetric processing flow (i.e., image preprocessing, spice initialization, image matching for tie points), and initial control network is generated.
(2) Then, the illumination information corresponding to the exposure time of each LRO NAC image is obtained, such as the incidence angle and the azimuth angle of the sun.
(3) Use the illumination conditions of each LRO NAC image to render LOLA DEM and generate corresponding hillshading maps.
(4) Match DOMs of LRO NAC images with LOLA DEM simulated images based on a thumbnail to obtain the tie points files.
(5) Extract the elevation information corresponding to the LOLA DEM from the pixel coordinates of the tie points, and convert the matched tie points into three-dimensional coordinates.
(6) Add control points with latitude, longitude, and elevation attributes to the previously constructed control network, and we specify some uncertainties (e.g., 20m, 20m, 1m) to the control points as weights.
(7) Perform iterative bundle adjustment again on the control network with added elevation control points, and generate the final DEMs and DOMs.
Matching the LOLA DEM shade maps with the geometrically corrected LRO NAC images reduces the matching error to some extent, since the orthorectified images have eliminated geometric errors and can be resampled as the same spatial resolution of the corresponding LOLA DEM.Moreover, the hill shaded LOLA DEM maps are rendered by the same illumination information as the LRO NAC images to be matched.This facilitates to improve the image matching result, and can obtain more control points.

Generation of the initial control network
The original LRO NAC images are stored in PDS format, and they need to be converted to other formats that can be interpreted by photogrammetric software.The Integrated Software for Imagers and Spectrometers (ISIS) software of U.S. Geological Survey (USGS) can be used to convert LRO NAC raw images into internally usable Cube image files.Afterwards, preprocessing such as radiometrically calibration and echo removing is required.After completing these steps, SPICE information can be bound to LRO NAC image data.This step involves matching the metadata information such as external orientation elements recorded during image exposure.After binding SPICE information to LRO NAC images, the images have initial positional information.
With this positional coordinate information, it is possible to calculate the outline of each image's external boundary.The image outline can determine the overlapping relationship between different images, generating an overlapping relationship file that records each overlapping area.The initial control network was built based on this overlapping relationship file.Firstly, establish tie points in all overlapping areas, which can establish connections between different images.Afterwards, some tie points are removed according to certain parameters, such as matching error, density, etc.By completing these, a preliminary control network with accurate matching and overall uniformity is obtained.Further use Jigsaw to iteratively adjust multiple times, gradually removing outliers and optimizing the external orientation elements of the images (Edmundson et al., 2012).Finally, a satisfactory control network is achieved, which is the initial control network in the proposed method.

Automatic Acquisition of Elevation Control Points
Illumination information was extracted from the pre-processed LRO NAC image, including the incident angle of light and the azimuth angle of the sun.Illumination information play a crucial role in the rendering of the LOLA DEM as well as the subsequent image matching between LRO NAC images and hill shaded LOLA DEM images.Based on the image matching results, the tie points are obtained, and the three-dimensional coordinate information of the tie points is extracted through the corresponding pixel coordinates of the tie points in the DEM.

Illumination Information Extraction for LRO NAC Image:
The pre-processed LRO NAC image contains various image metadata, which is a supplementary explanation of the image, including detailed image exposure time, sample and line numbers of pixels, center point longitude and latitude, image range, sample and line direction resolution, and illumination information.The illumination information includes the incidence angle, sun azimuth, etc.In this method, accurate illumination information is used to render LOLA DEM, where the incidence angle of the light and the sun azimuth are the two most important parameters.
The angle of incidence of light is closely related to the altitude angle of the sun, usually adding up to 90°.The solar altitude angle determines the length of the lunar shadow and also affects the area of the shadow in the image.The azimuth of the sun determines the direction of the source of sunlight, which also determines the direction of shadows.The shape and distribution of lunar shadows are jointly shaped by the two key indicators of solar altitude and azimuth.At the LSP, it has a consistently low solar altitude angle and rapidly changing solar azimuth, as well as drastically changing terrain.This illumination environment inevitably leads to rapid changes in the distribution of lunar shadows, and correspondingly, the radiation information of LRO NAC images also undergoes relatively significant changes at different exposure times.
In this case, a slight change in the solar azimuth will lead to a more noticeable change in the appearance of the LSP.So, LOLA DEM rendered by accurate illumination information, to some extent, can ensure the similarity between the rendered image and the captured image as much as possible, thereby greatly improving the likelihood of successful matching between the rendered image and the DOM image of LRO NAC.

Generation of LOLA DEM Hillshaded Images:
With the accurate illumination information extracted from LRO NAC images, simulation rendering of LOLA DEMs can be completed through illumination models.In addition to precise solar altitude and azimuth, elevation information from DEM is also required to calculate the slope and aspect of each pixel.Formula (1) is the key calculation process for rendering LOLA DEM.
The  ℎℎ here is the grayscale value of the pixels calculated in the rendered image, used to reflect the simulated radiance;  ℎ is the zenith angle, which is related to the altitude angle;   is the slope of the current pixel;  _ℎ is the mathematical form of the solar azimuth, starting from a different position than the solar azimuth;   is the slope direction of the current pixel. (3) The relationship between the solar zenith angle and altitude angle can be obtained from formula (2).Formula (3) converts the solar azimuth from geographical units to mathematical units, where   is the solar azimuth in geographical units represented by the compass direction.
Formula ( 4) can calculate the slope value of each pixel, where respectively, which can be obtained from the elevation values of the 8 pixels around the pixel.The direction of slope can be obtained from formula (5).

Image Matching and Elevation Control Points
Acquisition: LRO NAC is a linear array pushbroom sensor with dual-view non-stereoscopic coverage, and its images are characterized by large amount of data and relatively complex camera model (Geng et al., 2020).In addition, the image matching is extremely difficult due to the rugged terrain and drastic changes in illumination at the LSP.Moreover, matching LRO NAC images with simulated LOLA DEM hill shaded images poses great difficulties due to significant differences in scale and texture.
First of all, the resolution of LRO NAC image is around 1m/pixel, while the resolution of LOLA DEM at the region of interest is 20m/pixel.The large difference in resolution will significantly reduce the success rate of image matching.Second, there is a certain difference in appearance between the simulated LOLA DEM hill shade images and the LRO NAC images.Moreover, the LSP is characterized by many shadows, leading to a further reduction in the valid image data.These factors reduce the probability of a successful match.
To automatically acquire more elevation control points from LOLA DEM, we optimize for the image matching between LRO NAC images and LOLA DEM hill shaded images.First, the LRO NAC DOM image and the LOLA DEM are uniformly resampled to a resolution of 10m/pixel using cubic convolution interpolation resampling.This resolution lies in between the two, and the cubic convolution interpolation resampling can better ensure the continuity with pixels.Afterwards, perform LOLA DEM illumination rendering and generate the corresponding hill shade images.A three-layer pyramid is built for the NAC DOMs and the rendered images.Then, the low-resolution thumbnail image matching is performed, and the positional relationship of the corresponding blocks in the thumbnails is determined by the Hough transform.Next, the feature point matching on the original scale is performed based on the determined positional relationship of blocks.This can ensure the accuracy of matched points to a certain extent, and improve the efficiency.
For the problem of excessive shadows, we analysed the shadow threshold of each image through statistical histograms.Then use this threshold to round off a certain proportion of extreme values and re-stretch the remaining pixel values to a pixel depth of 8 bytes.This substantially improves the contrast of the image in non-shaded regions, thereby increasing the number of matching points.

Introduction of Control Points
The tie point has corresponding pixel coordinates in the LOLA DEM, and elevation information can be automatically extracted through this position.Then the pixel coordinates of image matching can be transformed to obtain three-dimensional control points with longitude, latitude, and elevation.These control points have elevation information from the LOLA DEM.
By adding measures, it is possible to bind the elevation control points to the initial control network generated by previous photogrammetric processing.The introduction of elevation control points can play a constraining role in the subsequent bundle adjustment process, thereby improving the accuracy of the entire control network.Specifically, the elevation control points require an appropriate weight to determine the range of adjustment, i.e. uncertainty, during bundle adjustment.Setting the uncertainty requires a deep understanding of the source of the control points.Based on the accuracy level of the LOLA DEM in plane and height directions, 20m, 20m, and 1m are suitable uncertainties for the longitude, latitude, and elevation of the introduced control points.

Generation of DEM
After adjusting to an acceptable level using the introduced control points, the adjusted images can be used to generate DEM and DOM.The emission angle of the LSP LRO NAC image is small (less than two degrees), and the intersection angle is also too small for stereo images, which can lead to poor reconstruction results.Here we use the Shape from Shading (SfS) method to derive DEM.Since SfS does not have very strict requirements for the intersection angle, and can generate accurate terrain close to the original image resolution.
The image with adjusted exterior orientation elements is the main input of SfS, and LOLA DEM is also required as a reference DEM.An accurate shadow threshold is a necessary prerequisite for generating high-quality DEM.In addition, we also introduce some smoothing constraints and initial terrain constrains to improve the quality of the SfS generated DEM.

Test Datasets
To verify the effectiveness of the method, in the vicinity of the LSP 85° south latitude, we used 17 LRO NAC images and LOLA DEM data for photogrammetric processing.Table 1 shows the metadata information of the LRO NAC images, including the resolution of the image, the solar incidence angle, the emission angle, and other information.Figure 2 shows the distribution of LRO NAC Images and the tie points of initial control network.

Image Matching Results
Figure 3 shows the image matching results of a pair of LRO NAC DOM and LOLA DEM rendered image.In the thumbnail image matching stage, the number of feature points matched for this pair of images is 8.The number of feature points matched in the original image is 48.Usually there are fewer points matched in the thumbnail images.But the matching results in the thumbnail images can be used for subsequent matching on the original scale, and helps to determine the corresponding image blocks.The right side of Figure 3 shows the detailed information of the rendered LOLA DEM and LRO NAC real terrain.Despite using the same illumination conditions for rendering, there are still differences in the appearance, especially for areas near the shadows.The shadow edges in LRO NAC images are sharper, while the rendered images' shadow edges appear relatively smooth.In addition, the illumination model is also a factor that may lead to this phenomenon.
In the end, a total of 192 roughly uniformly distributed feature points were matched.And these points will be converted into elevation control points and then introduced into the initial control network.The automatic extraction of LOLA DEM control points is completed through self-developed software using Qt 5.15.2 under Ubuntu 20.04.Open source software such as OpenCV and GDAL are used to facilitate the image processing and feature matching.

Bundle Adjustment Results of Control Network
The main indicators for the adjustment results of the experimental data are shown in Table 2. Figure 4 shows the root mean square (RMS) of the images in the adjustment process, which is basically around 1 pixels.

DOM Accuracy Evaluation
The stitching effect of adjacent DOMs is an important and direct means to reflect the bundle adjustment results and accuracy of the control network.Figure 5 shows the stitching effect of DOM after adding LOLA DEM control points to the control network.
For the mosaicked DOMs, the stitching accuracy is basically within one pixel, with some areas near the shadows slightly higher than one pixel.As shown in Figure 5, there is a display of the stitching details in the small subfigure above Figure 5.It can be seen that the crater here happens to be at the junction of the adjacent images, but the stitching of this crater is clearly complete at the pixel level.This also indicates that a satisfactory bundle adjustment result was achieved after the introduction of LOLA DEM control points.

DEM Accuracy Evaluation
Based on the refined exterior orientation elements, we conducted SfS reconstruction for two test regions, each covers an area of 1km×1km with a grid spacing of 1 m/pixel.As shown in Figure 6, the SfS generated DEMs' terrain is relatively fine and complete.It is obvious that the reconstructed DEM has more details.Many small impact craters that do not exist in the LOLA DEM have also been reconstructed, and the terrain has good integrity and connectivity.In order to further evaluate the effectiveness of introducing LOLA DEM elevation control points, we extracted two profile lines from the reconstructed SfS DEM and extracted the same profile lines from the corresponding LOLA DEM.We evaluated the effectiveness of introducing control points by comparing the differences between SfS generated DEM and LOLA DEM.A more detailed analysis is shown in Figure 9.It shows the difference between the reconstructed SfS DEM and LOLA DEM profiles, and the fluctuation of the line reflects the magnitude of the elevation difference between the two.From the graph, it can be seen that fluctuation range of both lines is basically within 2 meters, indicating that the difference between the reconstructed DEM and LOLA DEM is very small.The average elevation difference between region 1 and region 2 is 0.7 and 0.4 meters, respectively, with standard deviations of 1.0 and 0.7 meters, both of which do not exceed 1 meter or one pixel.

Conclusions
Compared to the current manual editing method of introducing control points in USGS ISIS, the proposed method can automatically add elevation control points with higher efficiency.And the experimental results demonstrate that the method can improve the positioning accuracy of mapping products.However, there is still room for improvement, such as image matching with poor illumination conditions, weakly textured areas and large differences in scale.These are currently the main difficulties in mapping the LSP and will also be the direction of our future research.
of elevation change in the x and y directions,

Figure 2 .
Figure 2. Distribution of LRO NAC Images and tie points.Yellow crosses are ignored tie points located in shaded areas.

Figure 3 .
Figure 3. Matching results of image M1166617859RE.The area in the green line on the right is the corresponding area.

Figure 4 .
Figure 4. RMS of images after bundle adjustment.

Figure 7 .
Figure 7. DEM profile chart of reconstruction area 1.

Figure 8 .
Figure 8. DEM profile chart of reconstruction area 2.

Figures 7
Figures7 and 8show the comparison between the SfS generated DEM and LOLA DEM profiles at two locations.The red line in the small square colored image at the corner shows the corresponding position of the profile line.Firstly, it can be seen that the trends of SfS DEM and LOLA DEM in the profile lines of the two positions are basically the same, indicating that the introduced elevation control points have played a good role.In addition, SfS generated DEM is more detailed in the profile compared to LOLA DEM, while LOLA DEM is smoother.This is because SfS DEM is capable of reconstructing terrain close to the original image.

Figure 9 .
Figure 9. Differences in profiles between SfS DEM and LOLA DEM.

Table 2 .
Main indicators of bundle adjustment results.