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:
where
and
are the minimum and maximum values of the valid gray values;
is the peak value in histogram with smaller gray values;
is the valley value in the histogram adjacent to
;
is the gray value differences between the peak and valley values;
is a parameter that adjusts
with a value ranging from 0 to 1;
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:
where
and
are the gray values before and after image contrast stretching. Based on our experiences, a value of 0.75 for
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:
where
indicates the 3D ground point’s coordinates on the lunar surface,
are the perspective center coordinates of the best scan line corresponding to the ground point,
is the scale factor,
is the rotation matrix from the LRO NAC frame to the LRO spacecraft body frame,
is the rotation matrix from the LRO spacecraft body frame to the star tracker frame,
is the rotation matrix from the star tracker frame to the J2000 coordinate frame,
is the rotation matrix from the J2000 coordinate frame to the lunar body’s fixed coordinate frame,
are the 2D image point coordinates, and
is the focal length.
and
are constant values,
can be calculated with the planetary constant kernels, and
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:
where
,
, and
and
,
, and
are the measured or initial EO parameters of the camera at time
, respectively;
,
,
,
,
, and
are the error correction terms for the EO parameters of the camera at time
, respectively;
,
,
,
,
, and
are the refined EO parameters of the camera at time
, 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:
where
,
, and
and
,
, and
are the constant error correction terms for the spacecraft positions and camera angles, and
,
, and
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:
where
is the standard deviation of the unit weight,
is the uncertainty of the initial position and attitude observations of the images, and
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
where
is the model’s gray value for the radiometrically calibrated image at the image coordinates
,
is the normal albedo of the object surface at ground point
,
is the emission angle,
is the incidence angle,
is the phase angle, and
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:
where
is the surface normal vector,
is the vector of the viewing direction, and
is the vector of the direction of light. Let
be the height value of a grid point (or pixel) on the DEM, and let
be the sample and line coordinates of the grid point. Then, at each point of the object surface, with normal vector
, the emission angle
and the incidence angle
are a function of
. Specifically, the normal vector
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:
where
is the image gray value in image
(i.e., observation), and
are the residuals of the observed gray value
. The image coordinates
can be determined by the back-projection operation (ground-to-image) of the camera model. Specifically, for each pixel
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
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
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
and the DEM height value
.
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:
where
is the initial terrain,
are the sum of squares of all second-order partial derivatives of
,
is the parameter of the smoothing constraint term, and
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.
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 km2 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 km2 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.