The Generation of High-Resolution Mapping Products for the Lunar South Pole Using Photogrammetry and Photoclinometry

Abstract

High-resolution and high-accuracy mapping products of the Lunar South Pole (LSP) will play a vital role in future lunar exploration missions. Existing lunar global mapping products cannot meet the needs of engineering tasks, such as landing site selection and rover trajectory planning, at the LSP. The Lunar Reconnaissance Orbiter (LRO)’s narrow-angle camera (NAC) can acquire submeter images and has returned a large amount of data covering the LSP. In this study, we combine stereo-photogrammetry and photoclinometry to generate high-resolution digital orthophoto maps (DOMs) and digital elevation models (DEMs) using LRO NAC images for a candidate landing site at the LSP. The special illumination and landscape characteristics of the LSP make the derivation of high-accuracy mapping products from orbiter images extremely difficult. We proposed an easy-to-implement shadow recognition and contrast stretching method based on the histograms of the LRO NAC images, which is beneficial for photogrammetric and photoclinometry processing. In order to automatically generate tie points, we designed an image matching method considering LRO NAC images’ features of long strips and large data volumes. The terrain and smoothness constraints were introduced into the cost function of photoclinometry adjustment, excluding pixels in shadow areas. We used 61 LRO NAC images to generate mapping products covering an area of 400 km². The spatial resolution of the generated DOMs was 1 m/pixel, and the grid spacing of the derived DEMs was 1 m (close to the spatial resolution of the original images). The generated DOMs achieved a relative accuracy of better than 1 pixel. The geometric accuracy of the DEM derived from photoclinometry was consistent with the lunar orbiter laser altimeter (LOLA) DEM with a root mean square error of 0.97 m and an average error of 0.17 m.

1. Introduction

Various studies have shown that there may be water ice at the Lunar South Pole (LSP), especially in the permanently shaded regions (PSRs) [1,2,3]. The water ice at the LSP has important scientific significance and practical value. Specifically, water ice can be used to analyze the source and mechanism of evolution of water on the Moon [4]. It can also be utilized as an important resource for human lunar activities, such as the construction of a future lunar base. Thus, the exploration of the LSP has become the focus of many space agencies. China’s future lunar exploration missions (e.g., Chang’E-7 and Chang’E-8) will explore the LSP in detail [5,6]. The Artemis program of the United States expects to return to the Moon and will conduct specific studies of the LSP [7]. Japan’s Smart Lander for Investigating the Moon (SLIM) probe was landed near the south pole of the Moon. Although SLIM was tilted when it landed, it accomplished its intended mission to some extent. In addition, Russia and India have recently conducted exploration missions in the vicinity of the LSP. High-resolution digital orthophoto maps (DOMs) and digital elevation models (DEMs) are widely used to support landing site selection and rover exploration. However, the existing mapping products of the LSP exhibit a low resolution and insufficient valid coverage, so they cannot meet the engineering requirements.

Currently, the mainstream lunar global mapping products mainly include the near-whole-Moon DEM (GLD100) and whole-Moon DOM produced with the image data from the wide-angle camera (WAC) of the Lunar Reconnaissance Orbiter (LRO), the lunar orbiter laser altimeter (LOLA) DEM produced with the altimeter data from LOLA [8,9], the near-whole-Moon DEM (SLDEM2015) produced with the Terrain Camera (TC) of Japan’s Selenological and Engineering Explorer (SELENE) lunar probe and LOLA data [10], the whole-Moon DEM and DOM produced with image data from China’s Chang’E-2 Stereo Camera (CE2TMap2015), and so on [11]. There are certainly some localized areas where improved products exist, such as the south pole of the Moon or other locations of interest [12]. However, in general, each of these lunar global mapping products has significant weaknesses at the LSP. Alternatively, existing mapping products for the LSP are lacking. For example, SLDEM2015 is partially missing at the LSP, GLD100 only covers latitudes between –79 and 79 degrees, and the image data of Chang’E-2 at the LSP have relatively low accuracy and more image shadows. These problems have resulted in a lack of high-precision and high-quality mapping products for supporting engineering applications and scientific analyses in future lunar exploration missions focused on the LSP.

The LRO NAC can acquire lunar remote sensing images with a submeter spatial resolution based on the pushbroom imaging principle [13,14]. Since its launch in 2009, after years of in-orbit observation, the image data returned from the LRO NAC achieved complete coverage of the entire Moon. In particular, the LSP has been covered many times. In theory, LRO NAC images can be used to generate higher-resolution DOMs and DEMs for the LSP [15]. However, the poor illumination conditions caused by a solar elevation angle of close to zero degrees and the very rugged topographic characteristics of the LSP have caused great difficulties in photogrammetric processing and the generation of mapping products [16,17]. Specifically, most areas of remote sensing images from lunar orbiters are obscured by shadows [18]. Moreover, the appearance of these images also changes greatly with time due to the rapid changes in the solar azimuth angle. Consequently, image matching becomes quite difficult for stereo pairs acquired at different times. Liu et al. (2020) proposed an illumination-invariant matching algorithm to produce a localized DEM of the Chang’E-4 and Chang’E-5 landing areas [19]. Hu et al. (2017) proposed a coupled epipolar rectification method for stereo LRO NAC images, and they introduced a semi-global matching (SGM) algorithm for the matching of LRO NAC images and reduced the discrepancies in LRO NAC images [20]. Chen et al. (2022) proposed a convolutional neural network (CNN)-based topography reconstruction method using LRO NAC images, which generated a DEM mosaic of the Chang’E-4 landing area in the range of 72.8 km × 30.3 km with a resolution of 1.5 m/pixel [21]. Liu et al. (2022) utilized a generative adversarial network (GAN) approach to generate a high-resolution DEM using a low-resolution DEM and high-resolution monocular LRO NAC images [22]. Xin et al. (2018) proposed a high-precision co-registration method to achieve automatic high-precision matching of LRO NAC images with SLDEM2015 and achieved an accuracy of 25.17 m [23].

Stereo-photoclinometry (SPC) is capable of reconstructing the three-dimensional surface of an object through images’ illumination information, and it can achieve pixel-level 3D reconstruction. Beyer et al. (2012) used submeter-scale HiRISE images and SPC to realize fast estimation of the slope of Mars’ surface, and they compared it with a DEM obtained from HiRISE stereo images, which proved that SPC could realize better slope estimation [24]. Liu et al. (2021) used error propagation and simulation analysis to investigate the main reasons affecting the reconstructive effect of SPC when using LRO NAC images; the results showed that the surface albedo had a great influence on the reconstruction quality, and the optimal incidence angle of the sun was about 60 degrees [25]. Liu et al. (2023) combined the atmospheric scattering effect transport model with SPC to effectively reduce the errors caused by aerosols in the Martian atmosphere, and they realized a pixel-level three-dimensional reconstruction of the Martian surface with a geometric accuracy of about two pixels [26].

However, these studies on lunar mapping and SPC reconstruction rarely involved the special terrain and extremely poor illumination conditions present at the LSP. In summary, there are limitations and challenges in the generation of high-quality mapping products for the LSP.

• (1). The poor illumination conditions and the special landscape of the LSP lead to great difficulties in the matching of stereo images. Consequently, the generation of tie points based on sparse feature point matching and the generation of a DEM through dense image matching face significant challenges.

• (2). The LRO NAC was basically perpendicular to the lunar surface when it was photographing. This meant that the intersection angles of stereo pairs and the base-to-height ratio were very small. Thus, the elevation accuracy of DEMs generated with photogrammetric methods is limited.

• (3). Although SPC has been widely used in the planetary mapping community, its effectiveness at the LSP, especially for latitude ranges between −85 degrees and −90 degrees (containing many PSRs), still needs to be further studied.

While the LRO NAC has returned a large number of images of the LSP, there are very few photogrammetrically processed or application-ready cartographic products [27]. In addition, due to the constantly changing positions of the shadows at the LSP, it is necessary to use images of the same position that are photographed at different times over many years to cover the largest possible area [28]. This also increases the difficulty of photogrammetric processing due to the irregular distribution of the tie points and the small overlapping area of stereo pairs. Thus, the photogrammetric processing and terrain reconstruction methods for LRO NAC images at the LSP need further optimization.

In this study, we developed a processing flow that can integrate stereo-photogrammetry and SPC to generate high-resolution mapping products for the LSP based on images from the LRO NAC. The main contributions of this study include the following: an easy-to-implement shadow recognition and image contrast stretching algorithm is proposed for LRO NAC images in the LSP; an image matching method based on an approximate orthophoto space and pyramid matching strategy is developed for the automatic generation of tie points. We made targeted optimization for the characteristics of large data volumes and long strips of LRO NAC images. The terrain and smoothness constraints are introduced into the observation equation of SPC reconstruction. This manuscript is organized as follows: In Section 2, the proposed methodology and the processing flow for LRO NAC images and DEM extraction, respectively, are described. In Section 3, we report the quantitative and qualitative evaluation results of the generated DOMs and DEMs, while in Section 4, we present the conclusion of the work.

2. Methods

The methodology of this study was divided into two parts, encompassing photogrammetric processing and DEM extraction through SPC. The stereo-photogrammetry technique was used to refine the initial position and attitude information of the LRO NAC images. The updated exterior orientation (EO) parameters were used as inputs for SPC, while SPC was used to derive a high-quality DEM with a spatial resolution close to that of the original image. We first describe the photogrammetric processing of LRO NAC images for the LSP.

The high solar incidence angle (close to 90 degrees) and the extremely rugged terrain of the LSP result in a large number of shaded areas. In addition, the illumination in the LSP changes rapidly, which makes stereo images of the same area appear significantly different. As a result, both the generation of tie points based on sparse image matching and DEM extraction through dense image matching always fail at the LSP. In a word, the derivation of mapping products for the LSP is a great challenge. In addition, the spatial resolution of DEMs generated with photogrammetric methods is usually 3–5 times lower than that of the original images. SPC or shape-from-shading (SfS) can be used to derive a DEM at roughly the same resolution as that of the original images [29]. However, SPC can be sensitive to the positional and attitude errors of the camera. Thus, to improve the accuracy of SPC-generated 3D terrain, the initial EO parameters of the images need to be refined using photogrammetric methods. Therefore, the advantages of photogrammetry and SPC are complementary.

Figure 1 presents a flowchart of the geometric processing for LRO NAC images. The details of the processing steps are presented below.

• (1). Data preprocessing (e.g., data format conversion and radiometric calibration) is conducted, and LRO NAC images are selected according to the observation geometry and the spatial coverage, and then the corresponding Spacecraft, Planet, Instrument, Camera-matrix, and Events (SPICE) kernels are attached to the image data.

• (2). Approximate orthophotos are generated using the initial EO parameters and a low-resolution DEM, such as the LOLA DEM.

• (3). The initial geometric accuracy of the LRO NAC images can be automatically evaluated by using image matching on approximate orthophotos. Since the pixel coordinates (samples and lines) of the matched points on approximate orthophotos can be transformed into map projection coordinates, the differences in the map projection coordinates of matching points on approximate orthophotos can reflect the relative accuracy of the LRO NAC images.

• (4). Tie points are generated based on the method of image matching on approximate orthophotos. Then, the initial geometric accuracy is evaluated by analyzing the pixel displacements of conjugate points for adjacent approximate images.

• (5). A pre-adjustment is conducted to remove wrongly matched points in the tie points. Note that a tie point usually contains several measures. To keep sufficient tie points, we do not directly discard the tie points when filtering outliers. Instead, we only remove the measures with large residuals, and other measures with small residuals (e.g., 1 pixel) of the same tie point can still be preserved.

• (6). A combined adjustment is conducted using image observations, non-photogrammetric observations (spacecraft position and camera angles), and LOLA elevation control points. Observation weights are set based on the initial results of the geometric accuracy evaluation and experiences.

• (7). Iterative adjustments are performed until the accuracy requirements of triangulation are met. Then, the refined EO parameters and the triangulated coordinates of tie points are generated. The refined EO parameters can be used for the following SPC processing.

• (8). The refined EO and the LOLA DEM are used to derive DOM and DEM products.

2.1. Automatic Extraction of Shadow Areas and Image Contrast Stretching

The shadows in the LSP significantly influence the photogrammetric processing as well as the SPC reconstruction. However, existing software tools rarely carry out targeted algorithm optimization for this. By analyzing the characteristics of the LRO NAC images in the LSP, we proposed a shadow recognition method that is easy to program and implement. It is based on the statistical analysis results of the image histograms. Figure 2 illustrates a typical histogram of an LRO NAC image in the LSP. As can be seen, the peak values in the histogram with smaller gray values always correspond to the shadow regions in the image. In addition, the gray values of the image are concentrated in the middle region, rather than evenly distributed between the maximum and minimum gray values. Thus, we can make full use of these features for shadow recognition and image enhancement. Specifically, we conduct image contrast stretching on the raw LRO NAC images based on removing image shadows. Firstly, the following equations are used to specify the valid gray values range of raw LRO NAC images:

(1)$\left\{\begin{array}{c}{G}_{left}=G_{peak}+\Delta g\times \lambda \\ \Delta g={G_{valley}-G}_{peak}\\ G_{right}=G_t\end{array}\right.$

where $G_{left}$ and $G_{right}$ are the minimum and maximum values of the valid gray values; $G_{peak}$ is the peak value in histogram with smaller gray values; $G_{valley}$ is the valley value in the histogram adjacent to $G_{peak}$; $\Delta g$ is the gray value differences between the peak and valley values; $\lambda$ is a parameter that adjusts $G_{left}$ with a value ranging from 0 to 1; $G_t$ can be determined by setting a valid pixel percent (e.g., 99.9%). Next, the floating-point gray values in raw LRO NAC images are converted to byte gray values using the following equations:

(2)$G_{new}=\left\{\begin{array}{c}0, { G}_{old}<G_{left} \\ \frac{G_{old}-G_{left}}{G_{right}-G_{left}}\times 255,{ G_{left}\le G}_{old}\le G_{right}\\ 255, { G}_{old}>G_{right}\end{array} \right.$

where $G_{old}$ and $G_{new}$ are the gray values before and after image contrast stretching. Based on our experiences, a value of 0.75 for $\lambda$ yields better results.

It is noted that some image processing software tools such as the ‘isis2std’ application from United States Geological Survey (USGS) Integrated Software for Imagers and Spectrometers (ISIS) version 5.0.1 can perform histogram stretching by directly discarding a portion of the pixels with smaller and larger gray values. But such a method will result in too many shadows on the converted image. As a comparison, the proposed image contrast stretching method can generate a higher-quality contrast stretched image, due to the fact that the pixels in the shadows are taken as invalid pixels. Figure 3 illustrates the shadow extraction and contrast stretching results for an LRO NAC image (i.e., M174930921RE). As can be seen, the proposed method derived a high-quality contrast stretched image and effectively extracted the shadow areas. The pixels in the shadow areas can be marked as invalid values in the following processing flow. This is beneficial for subsequent image matching (tie point generation) and SPC reconstruction.

2.2. Photogrammetric Processing of the LRO NAC Images

Since the LSP is the focus of many lunar exploration missions, the LRO NAC has obtained a large number of images covering it. Thus, given a mapping area at the LSP, the number of the related LRO NAC images is large, from several hundreds to tens of thousands. The first step before practical processing is the selection of valid images according to the requirements of photogrammetry and SPC. This can improve the efficiency of generating mapping products. Conventional photogrammetric processing has certain requirements for observation and illumination conditions to obtain an appropriate base-to-height ratio and improve the success rate of image matching. For example, the photogrammetric processing of OSIRIS images onboard the Rosetta spacecraft involves the selection of images with an incidence angle between 5 and 55 degrees and an emission angle between 0 and 50 degrees [30]. SPC has stricter requirements for the observation and illumination conditions, as it requires that images are acquired from different viewing directions with varying illumination conditions. However, the LRO NAC images of the LSP were acquired in poor illumination conditions. The perpendicular observation mode of the LRO NAC results in a small base-to-height ratio and limited viewing directions. Therefore, the observation and illumination conditions of LRO NAC images of the LSP cannot meet the ideal requirements for both photogrammetry and SPC. Thus, we could not use conventional indicators to select the images.

To maximize the image coverage to generate mapping products for the LSP, we selected images with an incidence angle of less than 88.5 degrees. According to our experience, images with an incidence angle of larger than 88.5 degrees contain too many shadows, which may significantly increase the difficulty of the following processing and contribute less to the final mapping products. In addition, some LRO NAC images have very short strips, and these can be directly removed from the image list through bundle adjustment. To facilitate automatic image selection, we first converted the raw LRO NAC images in the Planetary Data System (PDS) format to cube-format images of the USGS ISIS. Then, the basic information of the LRO NAC images, such as their image resolution, image size, emission angle, incidence angle, and geographic coordinate range, was acquired by using the ‘caminfo’ program in USGS ISIS. Next, the images in the ISIS cube format were converted into common tagged image file format (TIFF) images. The original floating-point TIFF images were converted into 8-bit-per-pixel (BPP)-format grayscale images. The image selection was first conducted by running a shell script, which automatically interpreted the basic information file corresponding to each image and selected the required images. Then, the selected images were further checked by photogrammetric engineers to remove the images with too many shaded areas or poor quality.

2.2.1. Construction of a Camera Model

A camera model was used to conduct a coordinate transformation between 2D image points and 3D ground points, which is the basis for various photogrammetric applications (e.g., DOM and DEM generation). To construct a rigorous sensor model for LRO NAC images, the corresponding SPICE kernels were required. Specifically, the LRO NAC images’ position and attitude information was stored in spacecraft position kernel (SPK) and orientation kernel (CK) files, respectively. In addition, the instrument kernels (IKs) and the planetary constant kernels (PCKs) were also required. The rigorous sensor model for LRO NAC images can be expressed as follows:

(3)$\left[\begin{array}{c}X\\ Y\\ Z\end{array}\right]=\left[\begin{array}{c}{X}_S\\ Y_S\\ Z_S\end{array}\right] +\lambda R_{J2000}^{Moon}{R}_{star}^{J2000}{R}_{body}^{star}{R}_{NAC}^{body}\left[\begin{array}{c}x\\ y\\ -f\end{array}\right]$

where $(X,Y,Z)$ indicates the 3D ground point’s coordinates on the lunar surface, $(X_S,Y_S,Z_S)$ are the perspective center coordinates of the best scan line corresponding to the ground point, $\lambda$ is the scale factor, $R_{NAC}^{body}$ is the rotation matrix from the LRO NAC frame to the LRO spacecraft body frame, $R_{body}^{star}$ is the rotation matrix from the LRO spacecraft body frame to the star tracker frame, $R_{star}^{J2000}$ is the rotation matrix from the star tracker frame to the J2000 coordinate frame, $R_{J2000}^{Moon}$ is the rotation matrix from the J2000 coordinate frame to the lunar body’s fixed coordinate frame, $(x,y)$ are the 2D image point coordinates, and $f$ is the focal length. $R_{NAC}^{body}$ and $R_{body}^{star}$ are constant values, $R_{J2000}^{Moon}$ can be calculated with the planetary constant kernels, and $R_{star}^{J2000}$ is derived from the observation data from the star tracker. Since the LRO NAC is a line-scan sensor, the EO parameters for each scan line can be interpolated using the EO parameters of neighboring sample points.

2.2.2. Automatic Generation of the Tie Points for the Lunar South Pole

High-quality tie points are essential for bundle adjustment. However, for the LSP, an orbiter’s remote sensing images always contain a large number of shaded areas. The illumination conditions of the LSP change rapidly, resulting in significant differences in the appearance of images obtained during different observations. Moreover, the irregular overlapping areas and smaller overlapping ranges also increase the difficulty of image matching. Last but not least, the data volume of LRO NAC images usually reaches the GB level, and many matching algorithms cannot be directly applied.

To further improve the quality and efficiency of the automatic generation of tie points, we developed an image matching method that considered the features of the LRO NAC’s long-strip pushbroom images and the unique terrain of the LSP. Specifically, we first conducted image matching on the low-resolution pyramid level of approximate orthophotos. Then, the matching results could be used to predict the corresponding image chip range when matching at the original resolution. This can provide a good start position for matching conjugate points, which facilitates a decrease in the outliers and improves the computational efficiency. The basic principle of the proposed image matching method is shown in Figure 4. In addition, because the proposed shadow recognition method can mark the shadow areas as invalid pixels, we directly remove feature points adjacent to the shadows. This can further improve the image matching results.

Figure 5 gives the comparison results of image matching between raw image space and approximate orthophoto space. A stereo pair composed of M174984904LE and M174971411RE were used to perform Speeded Up Robust Feature (SURF) matching. The parameters are the same and the number of features to be retained in the matching are 2000, 10,000, and 20,000. As can be seen, the image matching on orthophotos has more tie points. The success rate of image matching was improved from about 20~25% to 35~40%. Since image matching in the LSP faces great difficulties, the proposed image matching strategy can derive more tie points.

The automatic tie point extraction is shown in Figure 6, and the detailed steps are as follows.

• (1). The LRO NAC images were selected based on the method described in Section 2.2 and a small amount of manual inspection. Note that the shaded areas can be marked using the method described in Section 2.1, and the stretched images can be used for further image matching and tie point generation.

• (2). Then, the image overlap information was acquired by invoking USGS ISIS’s ‘findimageoverlaps’ application. Next, an image overlap list file in binary format describing how the images intersected with each other was generated. By parsing this binary-format overlap list file, we generated an image overlap file in text format for the following image matching.

• (3). The original LRO NAC images were firstly geometrically rectified using a coarse-resolution DEM (such as a LOLA DEM) and the initial EO parameters, and corresponding approximate orthophotos were generated. Then, image pyramids were generated for each approximate orthophoto.

• (4). For each set of image pairs to be matched in the overlap list file, image matching based on the SURF algorithm was conducted at the fourth image pyramid level of the approximate orthophotos. The matched points on low-resolution orthophotos were used to construct a coordinate transformation mapping relationship for the current image pair, which could be used to remove outliers and determine the related regions between the reference image and the matched image for image matching at the original resolution. To address the issue of a large data volume, we split the approximate orthophotos into chips (e.g., 2048 × 2048) with some small overlapping area (tens of pixels). Image matching at the original resolution of the approximate orthophoto was conducted on the image chips, and then the matched points on the image chips were merged.

• (5). The matched points on the approximate orthophotos were then converted into the original images. The random sample consensus (RANSAC) algorithm was used to further remove outliers. Considering the large number of shadows at the LSP, we removed feature points very close (e.g., less than two pixels) to the shadows to improve the image matching quality.

• (6). Then, these matched points were used to generate a tie point file for USGS ISIS in the American Standard Code for Information Interchange (ASCII) Parameter Value Language (PVL) format. PVL files are widely used by USGS ISIS for the storage of various data values. The generated tie points in the PVL format were converted into a USGS ISIS tie point file in binary format using USGS ISIS’s ‘cnetpvl2bin’ application. Then, this binary tie point file was directly used by ‘jigsaw’ application of USGS ISIS (version 5.0.1), which is a powerful bundle adjustment software tool that can support a large number of planetary exploration missions.

2.2.3. Absolute Control Using the LOLA DEM

It was noticed that some photogrammetric software, such as the National Aeronautics and Space Administration (NASA) Ames Stereo Pipeline (ASP) version 3.0.0, can use a reference DEM to optimize the results of bundle adjustment by automatically aligning images and the reference DEM [31,32]. However, the matching between LRO NAC images and the LOLA DEM always fails due to the large differences in spatial resolution and the special appearance of LRO NAC images of the LSP. In this study, the LOLA DEM’s elevation control points were introduced into the previously generated tie points file through manual editing by the photogrammetric operators. Specifically, we used USGS ISIS’s ‘qnet’ application to measure the pixel coordinates of a feature point (or measure) on one LRO NAC image according to the same feature point found on a hill-shaded map of the LOLA DEM. Then, the latitude, longitude, and elevation coordinates of the feature point in the LOLA DEM were associated with the measure in the tie point file. Next, the measure was used as a reference to match more measures on other overlapping LRO NAC images using USGS ISIS’s ‘qnet’ application. The hill-shaded LOLA DEM is beneficial for the identification of features (e.g., small craters). It should be noted that there were some errors in the plane direction between the LRO NAC images and the LOLA DEM. Thus, the uncertainties of the elevation control points obtained from the LOLA DEM should have been relatively large.

2.2.4. Bundle Adjustment

For satellite pushbroom images, each scan line has its own EO parameters. Thus, it is impossible and unnecessary to solve all scan lines’ EO parameters directly. In order to conduct the bundle adjustment of LRO NAC images, we used the following mathematical equations to model the trajectory errors of the LRO spacecraft:

(4)$\left\{\begin{array}{c}X\left(t\right)=X_{init}\left(t\right)+∆X\left(t\right)\\ Y\left(t\right)=Y_{init}\left(t\right)+∆Y\left(t\right)\\ Z\left(t\right)=Z_{init}\left(t\right)+∆Z\left(t\right)\\ \omega \left(t\right)={\omega }_{init}\left(t\right)+∆\omega \left(t\right)\\ \phi \left(t\right)={\phi }_{init}\left(t\right)+∆\phi \left(t\right)\\ \kappa \left(t\right)={\kappa }_{init}\left(t\right)+∆\kappa \left(t\right)\end{array}\right.$

where $X_{init}\left(t\right)$, $Y_{init}\left(t\right)$, and $Z_{init}\left(t\right)$ and ${\omega }_{init}\left(t\right)$, ${\phi }_{init}\left(t\right)$, and ${\kappa }_{init}\left(t\right)$ are the measured or initial EO parameters of the camera at time $t$, respectively; $∆X\left(t\right)$, $∆Y\left(t\right)$, $∆Z\left(t\right)$, $∆\mathsf{\omega }\left(t\right)$, $∆\phi \left(t\right)$, and $∆\kappa \left(t\right)$ are the error correction terms for the EO parameters of the camera at time $t$, respectively; $X\left(t\right)$, $Y\left(t\right)$, $Z\left(t\right)$, $\mathsf{\omega }\left(t\right)$, $\phi \left(t\right)$, and $\kappa \left(t\right)$ are the refined EO parameters of the camera at time $t$, respectively. We used both the constant error and velocity error correction terms for the spacecraft positions but only the constant error correction term for the camera angles. Consequently, for each LRO NAC image, there were 9 unknown parameters to be solved. The mathematical models used to correct the image trajectory errors with time t can be expressed as follows:

(5)$\left\{\begin{array}{c}∆X\left(t\right)=X_0+X_1\left(t\right)\\ ∆Y\left(t\right)=Y_0+Y_1\left(t\right)\\ ∆Z\left(t\right)=Z_0+Z_1\left(t\right)\\ ∆\omega \left(t\right)={\omega }_0\\ ∆\phi \left(t\right)={\phi }_0\\ ∆\kappa \left(t\right)={\kappa }_0\end{array}\right.$

where $X_0$, $Y_0$, and $Z_0$ and ${\omega }_0$, ${\phi }_0$, and ${\kappa }_0$ are the constant error correction terms for the spacecraft positions and camera angles, and $X_1\left(t\right)$, $Y_1\left(t\right)$, and $Z_1\left(t\right)$ are the drift (or velocity) error correction terms for the spacecraft positions.

The position and attitude observations of the LRO spacecraft were acquired from Spacecraft, Planet, Instrument, Camera-matrix, and Events (SPICE) kernels. Then, they were converted into the initial EO parameters of the LRO NAC images according to the camera misalignment information. The spacecraft positions and camera angles were used as non-photogrammetric observations for bundle adjustment with image observations. The weight of the image observations (i.e., measures) was always specified as the unit weight (e.g., 1) in the bundle adjustment. The non-photogrammetric observations, such as spacecraft positions and camera angles, needed to be given correct weights as well. The following equation was utilized to specify the weight values of the non-photogrammetric observations:

(6)$P_i=\frac{{\sigma }_0^2}{{\sigma }_i^2}$

where ${\sigma }_0$ is the standard deviation of the unit weight, ${\sigma }_i$ is the uncertainty of the initial position and attitude observations of the images, and $P_i$ is the weight value relative to the corresponding non-photogrammetric observations. Usually, the measurement accuracies of spacecraft positions and camera angles are well known, and they can be used to set the weight values.

2.2.5. DOM Generation

For mapping products near polar areas, polar stereographic projection is always used. Due to the long strips in LRO NAC images and the characteristics of polar stereographic map projection, the generated DOM for a typical LRO NAC image occupied a large amount of disk space (larger than 10 GB). Moreover, the orthorectification process involved a large number of ground-to-image coordinate transformations, which is a very time-consuming operation for pushbroom images [33]. Thus, to improve the computational efficiency of DOM generation for LRO NAC images of the LSP, we adopted an improved back-projection algorithm for pushbroom images, and it was able to significantly decrease the number of iterations needed to determine the best scan line corresponding to a ground point [34]. Based on the generic pushbroom model of USGS ISIS, the interior orientation parameters and EO parameters of LRO NAC images were acquired and stored in a camera file (.cam) and orientation data file (.odf), respectively, according to the file format designed for the airborne linear array camera ADS40 [35]. The orthorectification process was conducted using an inverse method, and the bilinear interpolation method was used to conduct image resampling.

2.3. DEM Extraction

2.3.1. SPC Reconstruction

In short, the SPC terrain reconstruction process can be viewed as a least squares optimization problem, whereas the terrain height or slope and the albedo are considered unknowns. The SPC can use the camera equation to generate a simulated image. Then, for each grid point or pixel on the DEM to be updated, the corresponding gray values on the original LRO NAC image and the simulated image can be determined by ground-to-image coordinates transformation. Then, the SPC is solved by minimizing the square of differences of gray values.

We used the lunar Lambertian model to express the reflectance properties [36]. Since the Moon is airless, there is no need to consider the effects of the atmosphere. Thus, the camera equation can be expressed as

(7)$G_C\left(x^{\prime },y^{\prime }\right)=A_N\left(X,Y,\alpha \right)\left[\Lambda (\alpha )(\frac{2\mathit{cos}\left(i\right)}{\mathit{cos}\left(i\right)+\mathit{cos}\left(e\right)}+\left(1-\Lambda \left(\alpha \right)\right)·\mathit{cos}\left(i\right)\right]$

where $G_C\left(x^{\prime },y^{\prime }\right)$ is the model’s gray value for the radiometrically calibrated image at the image coordinates $p\left(x^{\prime },y^{\prime }\right)$, $A_N\left(X,Y,\alpha \right)$ is the normal albedo of the object surface at ground point $P(X,Y,Z)$, $e$ is the emission angle, $i$ is the incidence angle, $\alpha$ is the phase angle, and

(8)$$\begin{gathered} \Lambda \left(\alpha \right)=1.0+A\alpha +B{\alpha }^2+C{\alpha }^3 \\ A=-0.019, B=2.42\times {10}^{-4},C=-1.46\times {10}^{-6} \end{gathered}$$

The emission angle is the angle between the surface normal vector and the vector of the viewing direction, and the incidence angle is the angle between the surface normal vector and the vector of the direction of light. The emission angle and the incidence angle can be calculated with the following equations:

(9)$\left\{\begin{array}{c}\cos \left(e\right)=\frac{\mathit{n}\cdot \mathit{v}}{\left|\mathit{n}\right|\cdot \left|\mathit{v}\right|}\\ \cos \left(i\right)=\frac{\mathit{n}\cdot \mathit{s}}{\left|\mathit{n}\right|\cdot \left|\mathit{s}\right|}\end{array}\right.$

where $\mathit{n}$ is the surface normal vector, $\mathit{v}$ is the vector of the viewing direction, and $\mathit{s}$ is the vector of the direction of light. Let $Z\left(s, l\right)$ be the height value of a grid point (or pixel) on the DEM, and let $\left(s, l\right)$ be the sample and line coordinates of the grid point. Then, at each point of the object surface, with normal vector $\mathit{n}$, the emission angle $e$ and the incidence angle $i$ are a function of $Z\left(s, l\right)$. Specifically, the normal vector $\mathit{n}$ of a surface point can be determined by the height values of the surrounding grid points.

The basic observation equations for SPC reconstruction are as follows:

(10)$v_j\left(x^{\prime },y^{\prime }\right)=G_C\left(x^{\prime },y^{\prime }\right)-g_j(x^{\prime },y^{\prime })$

where $g_j(x^{\prime },y^{\prime })$ is the image gray value in image $j$ (i.e., observation), and $v_j\left(x^{\prime },y^{\prime }\right)$ are the residuals of the observed gray value $g_j(x^{\prime },y^{\prime })$. The image coordinates $p\left(x^{\prime },y^{\prime }\right)$ can be determined by the back-projection operation (ground-to-image) of the camera model. Specifically, for each pixel $\left(s, l\right)$ on the DEM to be updated by SPC, the map projection coordinates are known. Then, the corresponding longitude and latitude coordinates can be computed through inverse map projection transformation. Next, the spatial Cartesian coordinates $P(X,Y,Z)$ can be computed according to the initial DEM. Then, the image coordinates are computed based on the ground-to-image transformation using Equation (3). Therefore, for each pixel $p\left(s, l\right)$ on the DEM to be updated, an observation equation can be derived for each image. Combining Equations (7)–(10), the unknowns during the least-squares adjustment of SPC are the normal albedo $A_N$ and the DEM height value $Z\left(s, l\right)$.

To prevent the SPC-generated DEM from deviating too far from the initial DEM, a terrain constraint term can be introduced into the cost function of SPC. In addition, considering the good continuity of the terrain on the lunar surface, we also introduce smoothness constraints into the cost function of SPC [37]. Thus, the overall cost function of SPC can be expressed as follows:

(11)$\int \int \sum _j{\left[G_C\left(x^{\prime },y^{\prime }\right)-g_j(x^{\prime },y^{\prime })\right]}^2+\mu {‖{\nabla }^{2}Z\left(s,l\right)‖}^2+\lambda {\left[Z\left(s,l\right)-Z_0\left(s,l\right)\right]}^{2}dsdl$

where $Z_0$ is the initial terrain, ${‖{\nabla }^{2}Z\left(s,l\right)‖}^2$ are the sum of squares of all second-order partial derivatives of $Z\left(s,l\right)$, $\mu$ is the parameter of the smoothing constraint term, and $\lambda$ is the parameter of the terrain constraint term. All input images are used to minimize the cost function and then to reconstruct the DEM. The same weight is used for different images.

We used the original or raw USGS ISIS .cub-format images (float point type) to generate the DEM. The stretched images (byte type) in Section 2.1 were mainly used to conduct image matching to derive tie points for bundle adjustment. When the images are stretched, the threshold value for shadows is determined for each LRO NAC image. This threshold value is then used in the DEM generation stage to exclude the pixels in shadows from the observation equation for SPC reconstruction. Due to the small base-to-height ratios of LRO NAC images and poor illumination in the LSP, the DEM based on dense matching always results in unsatisfactory quality. Thus, we directly use the LOLA DEM as the initial terrain for SPC reconstruction. Because the refined EO parameters are derived from photogrammetric processing and the LRO NAC images are controlled by LOLA DEM, the SPC-generated DEM also benefits from photogrammetric processing. Or rather, the refined DEM is the result of integrating photogrammetry and photoclinometry. Figure 7 presents the SPC reconstruction procedures when using LRO NAC images. The details of the process are given below the figure.

• (1). The image data, initial LOLA DEM, refined EO parameters, and illumination information required for SPC were prepared. The illumination information was acquired from the SPICE kernels.

• (2). Image preprocessing, such as the computation of the incidence angle, emission angle, and phase angle of the images, as well as radiometric calibration, was conducted. Then, the images that could best meet the needs of the observation geometry and illumination conditions for SPC reconstruction were selected.

• (3). SPC reconstruction was performed. Based on the camera equation, simulated images were generated. Terrain and smoothness constraints were introduced into the adjustment procedure.

• (4). Iterative adjustments were performed until the accuracy requirements of SPC were met.

• (5). A DEM and albedo map were generated for a DEM chip. Then, a mosaic of the DEM chips was used to form a large DEM.

To improve the quality of the DEM that was reconstructed through SPC, we automatically specified a gray threshold to indicate the shadow areas by analyzing a histogram of the LRO NAC images. Although the SPC method was able to reconstruct the terrain close to the resolution of the original image, it required relatively more stringent observational conditions for the input images and had higher requirements for computational performance. From the perspective of practical engineering, it is better to conduct DEM extraction using SPC for one small local area at a time.

2.3.2. DEM Mosaicking

We split the entire mapping area into small chips to prevent the SPC processing from exceeding the hardware limitations [38]. If the mapping area covered an area of 10 km × 10 km and the grid spacing of the generated DEM was 1 m, the image size of the DEM to be reconstructed through SPC was 10,000 × 10,000. Thus, the size of each small DEM chip could be 1000 × 1000. The specific size of the DEM chip depended on the available computer hardware resources. To facilitate the mosaicking of small DEM chips, adjacent DEM chips should have a small overlapping area (e.g., 100 m).

In order to derive a large DEM for the entire mapping area, we mosaicked the adjacent DEM chips generated through SPC. Because the adjacent DEM chips were generated with different LRO NAC images, the height value at the same location for adjacent DEM chips could show some small errors. Since the bundle adjustment provided updated EO parameters for LRO NAC images, the errors of geo-positioning in the adjacent strips could largely be removed. Therefore, the subtle differences in the elevation direction between DEM chips could be smoothed by an easy DEM centerline algorithm [39]. Then, a large DEM mosaic with smooth terrain in the boundaries of adjacent DEM chips could be derived. In the stage of mosaicking the DEM chips, we first mosaicked the DEM chips in each row to obtain several strips, and then we mosaicked these strips in columns to obtain a large DEM mosaic.

3. Results

3.1. Test Datasets

The candidate landing area had a latitude range of −85.3 degrees to −89.5 degrees and covered an area of approximately 400 km². According to its geographic coordinates, we downloaded several hundreds of raw LRO NAC images. Using the image selection method described in Section 2.2, 61 valid images were used to generate the mapping products. Figure 8 shows the overall distribution of the raw LRO NAC images in the experimental area. It can be seen that the structure of the distribution of tie points was very irregular, and the overlapping areas of some adjacent observations were very limited. Both open-source and in-house software were used to complete the photogrammetric and SPC processing. USGS ISIS was used to conduct data import and radiometric calibration [40]. The ISIS jigsaw program was used to conduct bundle adjustment. The in-house software included image matching, tie point generation, and orthorectification modules, as well as many shell scripts, to facilitate the complicated photogrammetric and SPC processing. The in-house software was mainly developed using Qt 5.15.2 in the Ubuntu 20.04 operating system. The shell programming language is used in Shell scripts and the C++ programming language is used in Qt. Based on the proposed method, we generated mapping products for a candidate landing site in the LSP covering an area of 400 km². The spatial resolution of the generated DOMs was 1 m/pixel, and the grid spacing of the derived DEMs was 1 m.

Weight setting is also essential for bundle adjustment. Strict weights may limit the range of adjustment of the EO parameters from the initial values, which may result in unsuccessful bundle adjustment results and poor triangulation accuracy. However, weights that are too loose can cause the adjustment results (including both EO parameters and the triangulated points) to deviate too much from the initial values. According to our experiences and the results of the initial evaluation of geometric accuracy derived from image matching on approximate orthophotos, the reconstructed SPICE kernels of the LRO NAC images showed relatively good geometric positioning accuracy. The discrepancies of the geographic coordinates of the matched points on approximate orthophotos were mostly on the orders of several tens of meters, with a small number approaching 200 m. This information facilitated the specification of an appropriate weight. The uncertainties of observations in Table 1 were used to set the weight values in the bundle adjustment.

Figure 9 shows the observation and illumination conditions of the raw LRO NAC images. In Figure 9a, the incidence and emission angles for the images are shown. It can be seen that the incidence angles of most images were between 85° and 90°. This indicates the poor illumination conditions of the LRO NAC images used to generate mapping products for the candidate landing site. Figure 9b shows the spatial resolution of the raw LRO NAC images, which was basically between 0.7 m/pixel and 1.0 m/pixel. Figure 9c shows the distributions of the solar ground azimuth angle and incidence angle, and Figure 9d shows the spacecraft ground azimuth angle and emission angle. These two figures were drawn in polar coordinates, as these can visually represent the observation and illumination conditions. The polar angle in Figure 9c corresponds to the solar ground azimuth, and the polar diameter corresponds to the solar incidence angle (0° means close to the center). It can be seen that the solar azimuth angle was widely distributed in a large range, and the incidence angle was close to 90°, which indicated that the solar illumination conditions varied drastically. This situation made image matching difficult because the appearance of the images could change dramatically. However, this was beneficial for DEM extraction through SPC because images acquired under different illumination conditions can provide complementary information for SPC. The polar angle and polar diameter in Figure 9d correspond to the azimuth angle of the spacecraft and the emission angle, respectively. The emission angle was very close to 0° due to the fact that the camera was approximately perpendicular to the lunar surface when acquiring images.

3.2. Bundle Adjustment Results

Based on the method described in Section 2.2.2, a tie point file was automatically generated from the 61 images. Figure 10 presents the tie points on the local areas of four LRO NAC images. It can be clearly seen that the tie points after bundle adjustment had a high quality and were evenly distributed. The overall distribution of tie points in the images is shown in Figure 11. These results demonstrate that the proposed image matching and automatic tie point generation methods were able to effectively extract enough feature points, even for images with very poor illumination conditions.

The LOLA DEM exhibited high accuracy in the radial direction, and it could be used to improve the height accuracy of the LRO NAC images. Using the method described in Section 2.2.3, we measured a total of 34 elevation control points in the LOLA DEM, as shown with a purple cross in Figure 11. The uncertainty in the height direction of these elevation control points was specified as 1 m, whereas the uncertainty in the plane direction of these elevation control points was specified as 20 m.

In order to remove a small number of outliers in the matched tie points and measures with large residuals, we conducted several iterative adjustments. At the last adjustment, we used 1.5 pixels as the threshold to remove the measures with large residuals. Empirically, this threshold could be specified to be appropriately higher at the initial stage of the bundle adjustment process. Because large errors in the measures had an impact on the surrounding tie points, using too strict of a threshold would remove too many good measures. Finally, the number of ground points was 23,546, and the total number of measures was 49,152. Other information on the results of bundle adjustment is listed in Table 2. The initial EO parameters were used as weighted observations. As described in Section 2.2.4, each image had nine unknowns to be solved in the bundle adjustment, and there were 549 constrained image parameters. The sigma0 value of the bundle adjustment was 0.44, which indicated a successful bundle adjustment result.

To further evaluate the bundle adjustment results, the image residuals and the geometric accuracy of the tie points for one image are shown in Figure 12. As can be seen in Figure 12a, the residuals of most LRO NAC images were between 0.3 and 0.7 pixels for both the sample and line directions. Some images had residuals larger than one pixel. This was mainly because these images had many shadows and showed large differences in appearance compared with the surrounding images. Consequently, the image connection for these images was very weak. Figure 12b shows the geometric accuracy of the tie points in one image. This indicated that the geometric accuracy in the plane direction was less than 1 m, whereas the geometric accuracy of the tie points in the height direction was between 4 and 5 m. The relatively low geometric accuracy in the height direction was mainly due to the small convergence angle and base-to-height ratio of the LRO NAC stereo pairs.

Tie point residuals are important indicators for evaluating the accuracy and results of bundle adjustment. Figure 13 shows two examples of residual vector maps in the image space of the tie points. As can be seen in (a) and (c), the magnitude of the residuals is basically less than one pixel. The corresponding tie point distributions of the original images M112816640RE and M108565000LE can be seen in (b) and (d). Figure 14 shows the vector map of the residuals of the introduced elevation control points. Most of them have residuals within 5 m in the horizontal direction, with some reaching about 10 m. In the elevation direction, the upward red arrows represent having positive residual values, and the downward represent negative. Most of the residuals in the direction of elevation are at the centimeter level, except for very few residuals that reach 0.5 m.

3.3. DOM Accuracy Evaluation

Figure 15 shows the effects of the generated orthophoto mosaic (1 m/pixel) for the candidate landing area. In the right part of Figure 15, it can be seen that the discrepancies in adjacent DOMs are less than one pixel. In particular, the effect of the DOM mosaic can be visualized in the adjacent image regions of the two craters that passed through the impact crater. This indicated that the generated DOMs showed high relative accuracy. It also demonstrated that the bundle adjustment delivered a good result, and the discrepancies between the adjacent LRO NAC images were removed. The updated EO parameters derived from the bundle adjustment provided consistent geometric accuracy for adjacent images. Thus, even if a low-resolution DEM was used in the orthorectification process, the geometric accuracy of the generated DOM would be good. Thus, we were able to generate well-controlled orthophoto mosaics for the candidate landing site.

3.4. DEM Accuracy Evaluation

The LOLA DEM with a grid spacing of 5 m was used as the initial terrain for DEM extraction with SPC. The candidate landing site covered an area of 400 km². Considering that the SPC algorithm is computationally time-consuming and consumes a significant amount of hardware resources, we only processed a small DEM chip of 2000 pixels by 2000 pixels during each iteration of SPC. The overlap between adjacent DEMs was 200 pixels. Thus, in total, there were 121 small DEM chips in 11 rows and 11 columns (see Figure 16). The absence of the white square on the left in Figure 16 is due to the fact that this location is a shaded area in all images. Based on the updated EO parameters derived from the bundle adjustment and the initial terrain, we generated 121 refined DEM chips through SPC. Then, a large DEM was derived by mosaicking these DEM chips (see Figure 17). The results showed that the DEM extracted with SPC was able to reveal more terrain details. It also demonstrated that the SPC technique could be used in the LSP up to a latitude of −88 degrees. Figure 18a shows the solar illumination of the SPC-generated DEM region from Mazarico et al. (2011). It shows the average percentage of illumination at the south pole of the Moon over a period of more than 10 years. The dark blue region in the figure is essentially a shadowed region with no illumination. In order to display the actual illuminated area, we combined the maximum gray value of each DOM pixel as the final gray value of the DOM, and the results are shown in Figure 18b. The unilluminated areas always correspond to low DEM quality. A comparison of Figure 17 and Figure 18 shows that there is also a significant difference in the quality of the SPC-generated DEMs, especially at the edges of the large region of shadows on the left.

Figure 19 compares a generated DEM chip with the corresponding DOM and LOLA DEM, and it provides a qualitative evaluation of the quality of the generated DEM. The grid spacing of the generated DEM was 1 m. The corresponding DOM had a resolution of 1 m/pixel. To facilitate the comparison of the spatial resolution, the LOLA DEM with a grid spacing of 5 m in the same area was resampled to 1 m/pixel. Figure 19a shows that the DEM generated with SPC preserved almost all of the impact craters in the original LRO NAC image. As a comparison, as shown in Figure 19b, even some large craters could not be identified in the LOLA DEM. Figure 19 demonstrates that the generated DEM was able to achieve a spatial resolution close to that of the original image.

Although the DEMs extracted with SPC retained fine-grained details, the quality of the generated DEM was indeed influenced by the input images. As can be seen in the middle of the left side of the shaded DEM map in Figure 19a, there was a small distorted terrain area on the generated DEM. This was largely due to the fact that the available images in this area were limited and did not complement each other well. More specifically, the illumination and observation conditions of the available images were not able to meet the requirements of SPC well, resulting in an inability for the corresponding area to be reconstructed. The distorted terrain could be filtered later.

To quantitatively evaluate the relative height accuracy of adjacent DEM chips, we calculated the height values of some sample points in the overlapping area of two adjacent DEM chips. Then, the height differences were computed. The results are shown in Figure 20. Figure 20a presents the distribution of the sample points. The sampling locations of the overlapping areas of the two adjacent DEM chips are shown at the red vertical lines in Figure 20. The outer boundaries of the two adjacent DEM chips have been marked out in light green and light blue. And the figures are drawn with the hill-shaded map of the generated DEM chips. Figure 20b,c presents a line graph of the height differences at the sample points. The mean value of the height differences was 1.46 m, and the root mean square error (RMSE) was 0.34 m. From the graph, it can be seen that the height difference fluctuates around 1.5 m, and the fluctuation range is approximately 1 m. At the position with the greatest undulation, the difference is relatively large (up to 2 m). It was also noted that there were relatively large discrepancies in the craters. This was caused by the fact that the shadows in the craters varied significantly according to the different sun azimuth angles in the different input images.

To analyze the absolute height accuracy of the SPC-generated DEM, we compared the SPC-generated DEM with the LOLA DEM. Specifically, we computed the height values of the sample points in the SPC-generated DEM and LOLA DEM and then computed the differences in these values. The results are shown in Figure 21. As can be seen, the SPC-generated DEM was very consistent with the LOLA DEM. The height differences had an RMSE of 0.97 m and an average error of 0.17 m. This also indicated that the absolute control points measured in the LOLA DEM played an important role in the bundle adjustment process. Additionally, we derived a pixel-by-pixel height difference map between the SPC-generated DEM and the LOLA DEM, which is shown in Figure 22. It can be seen that the maximum difference between the two is around 50 m over an area of 400 km². The height differences for most of the areas are within 10 m and the larger differences are basically in shaded areas, and the areas with larger slopes (e.g., craters). Note that such differences are mainly due to the fact that the SPC-generated DEM exhibits higher spatial resolution, and can provide more terrain details compared to the LOLA DEM.

4. Conclusions

High-accuracy mapping products for the LSP are essential for landing site selection, rover trajectory planning, and geological mapping for future lunar exploration missions. Due to the poor illumination conditions and the very rugged landscape of the LSP, there are many shadows in remote sensing images from lunar orbiters. Moreover, there are many permanently shadowed regions of the LSP, and the appearance of images from lunar orbiters changes significantly over time. Due to the complicated terrain and special illumination environments of the LSP, conventional photogrammetric methods cannot be used to derive satisfying results. Existing open-source and commercial photogrammetric software faces significant challenges in processing remote sensing images of the LSP.

This study presented a processing flow of generating high-quality mapping products based on the fusion of photogrammetry and SPC for the LSP. Finally, a mosaic of 61 LRO NAC DOM images of the preselected landing area was produced. The discrepancies in adjacent DOMs were less than one pixel. SPC DEMs with an area of up to 20 × 20 km² were produced with a resolution of 1 m/pixel, which is basically the same as the original images. The geometric accuracy of the DEM derived from photoclinometry was consistent with the LOLA DEM with an RMSE of 0.97 m and an average error of 0.17 m. The proposed shadow recognition and image matching method can effectively mark the shadow areas and remove feature points located near the shadow, which greatly improves the quality of the matched tie points. This is one of the main advantages of our method compared to USGS ISIS (version 5.0.1). The experimental results demonstrated that the proposed method was able to deliver a successful bundle adjustment result and that the SPC method can be used at the LSP to derive a high-quality DEM that is close to the resolution of the original image. The detailed accuracy evaluation results demonstrate that the targeted optimization for photogrammetry and photoclinometry works well for LRO NAC images in LSP. The mapping products generated for the candidate landing site can provide good support for engineering tasks. It is also noted that the processing speed of DEM extraction with the SPC method needs to be further optimized. The processing time required to generate a 2 × 2 km² SPC DEM is approximately 4 h on a desktop PC with an Intel i7-10700 processor and 64 GB of RAM. In addition, due to the characteristics of polar stereographic projection and the long-strip LRO NAC images, the data volume of the generated DOM is usually several tens of GB. Thus, the orthorectification algorithm for pushbroom images of the LSP needs to be optimized in the future.

Author Contributions

Conceptualization, X.G.; methodology, X.G. and P.L.; software, X.G. and P.L.; validation, P.L., J.Z. and X.G.; formal analysis, J.Z., Y.W. (Yuying Wang), Z.P., Y.W. (Yinhui Wang), X.M. and Q.W.; investigation, T.L.; resources, P.L., J.Z. and X.G.; data curation, X.M. and Q.W.; writing—original draft, P.L.; writing—review and editing, X.G.; visualization, Z.P.; supervision, X.G.; project administration, X.G.; funding acquisition, X.G. and T.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors on request.

Acknowledgments

We thank the authors of the open-source planetary remote sensing image processing software (i.e., USGS ISIS 5.0.1 and NASA ASP 3.0.0). We also thank the LRO NAC team for providing the planetary images to the public.

Conflicts of Interest

The authors declare no conflict of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Flowchart of the geometric processing for the LRO NAC images.

Figure 2. Contrast stretching based on the statistics of the image histogram for a typical LRO NAC image.

Figure 3. Shadow extraction and contrast stretching results for M174930921RE. (a) Contrast stretching based on the original minimum and maximum gray values. (b) Contrast stretching based on the proposed method. (c) Shadow recognition results based on the proposed method. Red areas in the right figure indicate the shadows.

Figure 4. Illustration of the proposed image matching method for automatic tie point extraction.

Figure 5. Comparison results of image matching between raw images and approximate orthophotos for a stereo pair composed of M174984904LE and M174971411RE. n means the number of features used for matching.

Figure 6. Flowchart of the generation of tie points for LRO NAC images.

Figure 7. Flowchart of SPC reconstruction based on the photogrammetric processing results.

Figure 8. Coverage of LRO NAC images. The red rectangle indicates the candidate landing area.

Figure 9. Observation and illumination conditions of raw LRO NAC images.

Figure 10. Distribution of the tie points on the images of local areas. The green crosses represent matched tie points.

Figure 11. Overall distribution of tie points in the images. The purple crosses represent the absolute elevation control points acquired from the LOLA DEM, and the green crosses represent matched tie points. Overlapping areas of images without tie points are shadowed areas.

Figure 12. Results of the evaluation of geometric accuracy of the bundle adjustment.

Figure 13. Residual vector maps in the image space of the tie points. (a,c) are residual vectors, and (b,d) are the corresponding original image tie points.

Figure 14. Residual vector maps of elevation control points. The explanation of the symbols in the upper left corner is a quantitative description of the magnitude corresponding to the length of the arrow.

Figure 15. DOM mosaic effect. The green rectangle indicates the candidate landing area.

Figure 16. Initial DEM chips for the SPC terrain reconstruction. The spatial resolution of the DEM was 1 m/pixel. The size of each DEM chip was 2 km by 2 km, with an overlap of 200 m between adjacent chips. The blue rectangle indicates the candidate landing area.

Figure 17. Color-shaded map of the generated DEM mosaic (1 m/pixel). The blue rectangle indicates the candidate landing area.

Figure 18. Analysis of the illumination conditions of the available images. (a) Average solar illumination of SPC-generated DEM area. (b) Valid pixels in the candidate landing site. Pixels displayed in black areas have not been illuminated. The blue rectangle indicates the candidate landing area.

Figure 19. Qualitative evaluation of the generated DEM. (a) Comparison of the SPC-generated DEM (left side) and the DOM of an LRO NAC image (right side). The hill-shaded DEM was drawn on top of the DOM of the LRO NAC image (M105925266LE). The gray rectangle on the right side indicates the same area corresponding to the generated DEM on the left side. Red rectangles indicate distorted terrain areas. (b) Comparison between the SPC-generated DEM (1 m/pixel) and resampled LOLA DEM (1 m/pixel). The left side shows the hill-shaded map of the SPC-generated DEM, and the right side shows the hill-shaded map of the LOLA DEM.

Figure 20. Quantitative evaluation results of the height differences in the overlapping areas of adjacent DEM chips. (a) Positions of the sample points in the overlapping area of two adjacent DEM chips. The red line indicates the sample points for calculating the height differences. (b) Elevation profiles of the adjacent DEM chips. (c) Illustration of the height differences between the adjacent DEM chips.

Figure 21. Results of a comparison between the SPC-generated DEM and LOLA DEM. The red line in the top-right subfigure indicates the location of the sample points.

Figure 22. Pixel-by-pixel height difference map between SPC-generated DEM and LOLA DEM.

Table 1

Uncertainties of observations in the bundle adjustment of LRO NAC images.

Observations	Uncertainty
Point latitude	100 (m)
Point longitude	100 (m)
Point radius	10 (m)
LOLA control point latitude	20 (m)
LOLA control point longitude	20 (m)
LOLA control point radius	1 (m)
Spacecraft position	50 (m)
Spacecraft velocity	1 (m/s)
Camera angles	0.1 (°)

Table 2

Bundle adjustment results of the LRO NAC images.

Items	Parameters
Number of images	61
Number of ground points	23,546
Number of total measures	49,152
Constrained image parameters	549
Constrained point parameters	70,638
Iterations	5
Sigma0	0.44

References

1. Li, S.; Lucey, P.G.; Milliken, R.E.; Hayne, P.O.; Fisher, E.; Williams, J.-P.; Hurley, D.M.; Elphic, R.C. Direct evidence of surface exposed water ice in the Lunar Polar Regions. Proc. Natl. Acad. Sci. USA; 2018; 115, pp. 8907-8912. [DOI: https://dx.doi.org/10.1073/pnas.1802345115] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30126996]

2. Deutsch, A.N.; Neumann, G.A.; Head, J.W. New evidence for surface water ice in small-scale cold traps and in three large craters at the north polar region of Mercury from the Mercury Laser Altimeter. Geophys. Res. Lett.; 2017; 44, pp. 9233-9241. [DOI: https://dx.doi.org/10.1002/2017GL074723]

3. Arnold, J.R. Ice in the lunar polar regions. J. Geophys. Res.; 1979; 84, pp. 5659-5668. [DOI: https://dx.doi.org/10.1029/JB084iB10p05659]

4. Biswas, J.; Sheridan, S.; Pitcher, C.; Richter, L.; Reganaz, M.; Barber, S.J.; Reiss, P. Searching for potential ice-rich mining sites on the Moon with the Lunar Volatiles Scout. Planet. Space Sci.; 2020; 181, 104826. [DOI: https://dx.doi.org/10.1016/j.pss.2019.104826]

5. Yu, H.; Rao, W.; Zhang, Y.; Xing, Z. Mission analysis and spacecraft design of Chang’E-7. J. Deep. Space Explor.; 2023; 10, pp. 567-576.

6. Zou, Y.; Liu, Y.; Jia, Y. Overview of China’s upcoming Chang’E series and the scientific objectives and payloads for Chang’E 7 mission. Proceedings of the 51st Annual Lunar and Planetary Science Conference; The Woodlands, TX, USA, 16–20 March 2020; 1755.

7. Evans, M.E.; Graham, L.D. A Flexible Lunar Architecture for Exploration (FLARE) supporting NASA’s Artemis Program. Acta Astronaut.; 2020; 177, pp. 351-372. [DOI: https://dx.doi.org/10.1016/j.actaastro.2020.07.032]

8. Smith, D.E.; Zuber, M.T.; Jackson, G.B.; Cavanaugh, J.F.; Neumann, G.A.; Riris, H.; Sun, X.; Zellar, R.S.; Coltharp, C.; Connelly, J. et al. The Lunar Orbiter Laser Altimeter investigation on the Lunar Reconnaissance Orbiter mission. Space Sci. Rev.; 2010; 150, pp. 209-241. [DOI: https://dx.doi.org/10.1007/s11214-009-9512-y]

9. Barker, M.; Mazarico, E.; Neumann, G.; Smith, D.; Zuber, M.; Head, J.; Sun, X. A new view of the Lunar South Pole from the Lunar Orbiter Laser Altimeter (LOLA). Planet. Space Sci.; 2023; 4, 183. [DOI: https://dx.doi.org/10.3847/PSJ/acf3e1]

10. Barker, M.K.; Mazarico, E.; Neumann, G.A.; Zuber, M.T.; Haruyama, J.; Smith, D.E. A new Lunar Digital Elevation Model from the Lunar Orbiter Laser Altimeter and SELENE Terrain Camera. Icarus; 2016; 273, pp. 346-355. [DOI: https://dx.doi.org/10.1016/j.icarus.2015.07.039]

11. Kirk, R.L.; Archinal, B.A.; Gaddis, L.R.; Rosiek, M.R. Lunar cartography: Progress in the 2000s and prospects for the 2010s. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci.; 2012; XXXIX-B4, pp. 489-494.

12. Barker, M.K.; Mazarico, E.; Neumann, G.A.; Smith, D.E.; Zuber, M.T.; Head, J.W. Improved LOLA elevation maps for south pole landing sites: Error estimates and their impact on illumination conditions. Planet Space Sci.; 2021; 203, 105119. [DOI: https://dx.doi.org/10.1016/j.pss.2020.105119]

13. Robinson, M.S.; Brylow, S.M.; Tschimmel, M.; Humm, D.; Lawrence, S.J.; Thomas, P.C.; Denevi, B.W.; Bowman-Cisneros, E.; Zerr, J.; Ravine, M.A. et al. Lunar Reconnaissance Orbiter Camera (LROC) instrument overview. Space Sci. Rev.; 2010; 150, pp. 81-124. [DOI: https://dx.doi.org/10.1007/s11214-010-9634-2]

14. Speyerer, E.J.; Wagner, R.V.; Robinson, M.S.; Licht, A.; Thomas, P.C.; Becker, K.; Anderson, J.; Brylow, S.M.; Humm, D.C.; Tschimmel, M. Pre-flight and on-orbit geometric calibration of the Lunar Reconnaissance Orbiter Camera. Space Sci. Rev.; 2016; 200, pp. 357-392. [DOI: https://dx.doi.org/10.1007/s11214-014-0073-3]

15. Di, K.; Jia, M.; Xin, X.; Liu, B.; Liu, Z.; Peng, M.; Yue, Z. High resolution seamless DOM generation over Chang’E-5 landing area using LROC NAC images. Int. Arch. Photogramm. Remote Sens. Spatial Inf. Sci.; 2018; XLII-3, pp. 271-276. [DOI: https://dx.doi.org/10.5194/isprs-archives-XLII-3-271-2018]

16. Mazarico, E.; Neumann, G.A.; Smith, D.E.; Zuber, M.T.; Torrence, M.H. Illumination conditions of the Lunar Polar Regions using LOLA topography. Icarus; 2011; 211, pp. 1066-1081. [DOI: https://dx.doi.org/10.1016/j.icarus.2010.10.030]

17. Cohen, B.A.; Hayne, P.O.; Greenhagen, B.; Paige, D.A.; Seybold, C.; Baker, J. Lunar flashlight: Illuminating the Lunar South Pole. IEEE Aerosp. Electron. Syst. Mag.; 2020; 35, pp. 46-52. [DOI: https://dx.doi.org/10.1109/MAES.2019.2950746]

18. Zhang, Y.; Liu, B.; Di, K.; Liu, S.; Yue, Z.; Han, S.; Wang, J.; Wan, W.; Xie, B. Analysis of illumination conditions in the Lunar South Polar Region using multi-temporal high-resolution orbital images. Remote Sens.; 2023; 15, 5691. [DOI: https://dx.doi.org/10.3390/rs15245691]

19. Liu, W.C.; Wu, B. An integrated photogrammetric and photoclinometric approach for illumination-invariant pixel-resolution 3D mapping of the Lunar surface. ISPRS J. Photogramm. Remote Sens.; 2020; 159, pp. 153-168. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2019.11.017]

20. Hu, H.; Wu, B. Precision 3D surface reconstruction from LRO NAC images using semi-global matching with couple epipolar rectification. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci.; 2017; XLII-3/W1, pp. 55-61. [DOI: https://dx.doi.org/10.5194/isprs-archives-XLII-3-W1-55-2017]

21. Chen, H.; Hu, X.; Glaser, P.; Xiao, H.; Ye, Z.; Zhang, H.; Tong, X.; Oberst, J. CNN-based large area pixel-resolution topography retrieval from single-view LROC NAC images constrained with SLDEM. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2022; 15, pp. 9398-9416. [DOI: https://dx.doi.org/10.1109/JSTARS.2022.3214926]

22. Liu, Y.; Wang, Y.; Di, K.; Peng, M.; Wan, W.; Liu, Z. A generative adversarial network for pixel-scale Lunar DEM generation from high-resolution monocular imagery and low-resolution DEM. Remote Sens.; 2022; 14, 5420. [DOI: https://dx.doi.org/10.3390/rs14215420]

23. Xin, X.; Liu, B.; Di, K.; Jia, M.; Oberst, J. High-precision co-registration of orbiter imagery and Digital Elevation Model constrained by both geometric and photometric information. ISPRS J. Photogramm. Remote Sens.; 2018; 144, pp. 28-37. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2018.06.016]

24. Beyer, R.A.; Kirk, R.L. HiRISE photoclinometry of final MSL landing sites. Proceedings of the 43rd Lunar and Planetary Science Conference; The Woodlands, TX, USA, 19–23 March 2012; 2694.

25. Liu, W.C.; Wu, B. Influence of solar incidence angle on single-image photoclinometry for precision Lunar topographic mapping. ISPRS J. Photogramm. Remote Sens.; 2021; 182, pp. 208-227. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2021.10.011]

26. Liu, W.C.; Wu, B. Atmosphere-aware photoclinometry for pixel-wise 3D topographic mapping of Mars. ISPRS J. Photogramm. Remote Sens.; 2023; 204, pp. 237-256. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2023.09.017]

27. Chen, C.; Ye, Z.; Xu, Y.; Liu, D.; Huang, R.; Zhou, M.; Xie, H.; Tong, X. Large-scale block bundle adjustment of LROC NAC images for Lunar South Pole mapping based on topographic constraint. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2024; 17, pp. 2731-2746. [DOI: https://dx.doi.org/10.1109/JSTARS.2023.3346199]

28. Henriksen, M.R.; Manheim, M.R.; Burns, K.N.; Seymour, P.; Speyerer, E.J.; Deran, A.; Boyd, A.K.; Howington-Kraus, E.; Rosiek, M.R.; Archinal, B.A. et al. Extracting accurate and precise topography from LROC Narrow Angle Camera stereo observations. Icarus; 2017; 283, pp. 122-137. [DOI: https://dx.doi.org/10.1016/j.icarus.2016.05.012]

29. Jiang, C.; Douté, S.; Luo, B.; Zhang, L. Fusion of photogrammetric and photoclinometric information for high-resolution DEMs from Mars in-orbit imagery. ISPRS J. Photogramm. Remote Sens.; 2017; 130, pp. 418-430. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2017.06.010]

30. Preusker, F.; Scholten, F.; Matz, K.-D.; Roatsch, T.; Willner, K.; Hviid, S.; Knollenberg, J.; Jorda, L.; Gutiérrez, P.; Kührt, E. et al. Shape model, reference system definition, and cartographic mapping standards for comet 67P/Churyumov-Gerasimenko. Stereo-photogrammetric analysis of Rosetta/OSIRIS image data. Astron. Astrophys.; 2015; 583, A33. [DOI: https://dx.doi.org/10.1051/0004-6361/201526349]

31. Beyer, R.A.; Alexandrov, O.; McMichael, S. The Ames Stereo Pipeline: NASA’s open source software for deriving and processing terrain data. Earth Space Sci.; 2018; 5, pp. 537-548. [DOI: https://dx.doi.org/10.1029/2018EA000409]

32. Shean, D.E.; Alexandrov, O.; Moratto, Z.M.; Smith, B.E.; Joughin, I.R.; Porter, C.; Morin, P. An automated, open-source pipeline for mass production of Digital Elevation Models (DEMs) from very-high-resolution commercial stereo satellite imagery. ISPRS J. Photogramm. Remote Sens.; 2016; 116, pp. 101-117. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2016.03.012]

33. Geng, X.; Xu, Q.; Wang, J.Y.; Lan, C.; Qin, F.; Xing, S. Generation of large-scale orthophoto mosaics using MEX HRSC images for the candidate landing regions of China’s first Mars mission. IEEE Trans. Geosci. Remote Sens.; 2021; 60, pp. 1-20. [DOI: https://dx.doi.org/10.1109/TGRS.2021.3128518]

34. Geng, X.; Xu, Q.; Lan, C.; Xing, S.; Hou, Y.; Lyu, L. Orthorectification of planetary linear pushbroom images based on an improved back-projection algorithm. IEEE Geosci. Remote Sens. Lett.; 2019; 16, pp. 854-858. [DOI: https://dx.doi.org/10.1109/LGRS.2018.2885086]

35. Geng, X.; Xu, Q.; Xing, S.; Lan, C. A Generic pushbroom sensor model for planetary photogrammetry. Earth Space Sci.; 2020; 7, e2019EA001014. [DOI: https://dx.doi.org/10.1029/2019EA001014]

36. Lohse, V.; Heipke, C.; Kirk, R.L. Derivation of planetary topography using multi-image shape-from-shading. Planet. Space Sci.; 2006; 54, pp. 661-674. [DOI: https://dx.doi.org/10.1016/j.pss.2006.03.002]

37. Min, S.; Rama, C.; Tal, S. Reconstructing a 3-D depth map from one or more images. Comput. Vis. Image Underst.; 1991; 53, pp. 219-226.

38. Kirk, R.L.; Mayer, D.P.; Dundas, C.M.; Wheeler, B.H.; Beyer, R.A.; Alexandrov, O. Comparision of Digital Terrain Models from two photoclinometry methods. Int. Arch. Photogramm. Remote Sens. Spat. Inf. Sci.; 2022; XLIII-B3-2022, pp. 1059-1067. [DOI: https://dx.doi.org/10.5194/isprs-archives-XLIII-B3-2022-1059-2022]

39. Korokhin, V.; Velikodsky, Y.; Shkuratov, Y.; Kaydash, V.; Mall, U.; Videen, G. Using LROC WAC data for Lunar surface photoclinometry. Planet. Space Sci.; 2018; 160, pp. 120-135. [DOI: https://dx.doi.org/10.1016/j.pss.2018.05.020]

40. Edmundson, K.L.; Cook, D.A.; Thomas, O.H.; Archinal, B.A.; Kirk, R.L. Jigsaw: The ISIS3 bundle adjustment for extraterrestrial photogrammetry. ISPRS Ann. Photogramm. Remote Sens. Spat. Inf. Sci.; 2012; I–4, pp. 203-208. [DOI: https://dx.doi.org/10.5194/isprsannals-I-4-203-2012]