A method for generating planetary surface DEMs by integrating photogrammetry and stereophotometry

Abstract

This invention relates to the field of 3D reconstruction technology, specifically to a method for generating a planetary surface DEM by integrating photogrammetry and stereophotometry. The method involves acquiring reconstruction preparation data, including a remote sensing image dataset, various parameter values ​​corresponding to each remote sensing image in the dataset, and an illumination model and its parameters. A target region is constructed, and the remote sensing images in the dataset are filtered based on the parameter values ​​corresponding to each image to obtain the optimal image dataset and an initial low-resolution DTM for the target region. For a single ground point within the target region, a stereophotometric error equation is constructed based on the optimal image dataset, the initial low-resolution DTM, and the illumination model and its parameters. The reconstruction result for each ground point is determined iteratively, ultimately resulting in a refined reconstructed image. This invention integrates photogrammetry and stereophotometry, effectively improving the 3D reconstruction effect of planetary surfaces.

Description

Technical Field

[0001] This invention relates to the field of planetary surface three-dimensional construction technology, specifically to a method for generating a planetary surface DEM by integrating photogrammetry and stereophotometry. Background Technology

[0002] High-resolution, high-precision digital elevation models (DEMs) of the surfaces of extraterrestrial bodies such as the Moon and Mars are crucial supporting data for engineering missions such as landing site selection, flight navigation, and rover path planning. Existing DEM products for extraterrestrial surfaces are typically created using classical photogrammetric methods based on data from earlier exploration missions, and their resolution remains in the tens to hundreds of meters range, which is no longer sufficient to meet the needs of deep space exploration engineering applications.

[0003] The main technical approaches for 3D reconstruction of planetary surfaces include image-based 3D reconstruction methods and lidar-based remote sensing methods. Image-based 3D reconstruction methods have the advantage of a wide operating range for the mapping camera and relatively high resolution for the reconstructed terrain; however, the quality of the reconstructed model is limited by lighting conditions. Since lidar is an active remote sensing device, lidar-based remote sensing methods offer the advantages of high elevation accuracy and are not limited by lighting conditions; however, their operating range and reconstruction resolution are limited.

[0004] Image-based methods for planetary 3D reconstruction mainly include photogrammetry, Structure from Motion (SfM), and Stereophotoclinometry (SPC). Photogrammetry and SfM are similar in principle, both using image matching to reconstruct 3D terrain information. Their main advantages are their mature engineering applications, high 3D reconstruction efficiency, and the ability to directly calculate absolute elevation. However, their disadvantages include relatively low resolution (typically 3-5 times the original image resolution) and poor reconstruction results in areas with weak texture. Furthermore, planetary remote sensing images often have sparse surface textures, making feature point matching difficult. Stereophotometry uses image grayscale information to calculate slope and albedo, thereby reconstructing 3D terrain. It can generate detailed surface terrain with the same resolution as the original image and is not limited by the sparseness of image textures, giving it advantages in planetary surface 3D reconstruction. However, stereophotometry calculates terrain slope rather than absolute elevation, and accuracy control in large-scale processing requires reference to the terrain as a constraint.

[0005] In recent years, stereophotometry has been increasingly applied to 3D reconstruction in planetary exploration missions. Its main advantage is the ability to obtain more detailed 3D models, better supporting engineering tasks such as flight navigation and landing site selection. Based on stereophotometric 3D reconstruction technology, at least 17 small celestial bodies have been created, and researchers are gradually applying stereophotometry to larger celestial bodies such as the Moon and Mars. Robert Gaskell, a researcher at the Planetary Science Institute (PSI), developed the SPC software toolkit based on the principles of stereophotometric reconstruction. Initially used primarily for optical navigation in deep space exploration missions, it has since been gradually applied to 3D reconstruction of planetary surfaces. PSI used stereophotometry for the 3D reconstruction of the asteroid Bennu, creating a low-resolution global DTM with a grid spacing of 3m based on images from the probe's approach phase. Before landing, a global DTM with a grid spacing of 5cm and a locally detailed DTM with a grid spacing of 0.5cm were created using a large amount of image data. The Jet Propulsion Laboratory (JPL) in the United States created a 3D model of Ceres with a grid spacing of 100m using 38,000 Dawn images based on stereophotometry.

[0006] Although stereophotometric-based 3D reconstruction methods have certain advantages in the field of planetary 3D reconstruction, the main drawbacks of existing stereophotometric algorithms are that they are computationally expensive. When modeling large areas, external constraint information is required to ensure the overall geometric accuracy of the reconstructed model. Moreover, the selection of lighting models and model parameters, as well as image shadows, are also prominent problems that restrict the application of stereophotometry to large-scale planetary surface 3D reconstruction, making it difficult to guarantee the resolution and accuracy of planetary 3D reconstruction. Summary of the Invention

[0007] The purpose of this invention is to provide a method for generating a planetary surface DEM that integrates photogrammetry and stereophotometry, in order to solve the problem of low resolution and accuracy of existing stereophotometric methods for three-dimensional reconstruction of planetary surfaces.

[0008] To address the aforementioned technical problems, this invention provides a method for generating a planetary surface DEM that integrates photogrammetry and stereophotometry, comprising the following steps:

[0009] Obtain reconstruction preparation data, which includes at least: remote sensing image dataset, various parameter values ​​corresponding to each remote sensing image in the remote sensing image dataset, and illumination model and its model parameters. The various parameter values ​​include at least: incident angle, observation angle, camera azimuth angle and image intersection angle.

[0010] Construct a target region, and based on the parameter values ​​corresponding to each remote sensing image in the remote sensing image dataset, filter the remote sensing images in the remote sensing image dataset to obtain the optimal image dataset for the target region, and use the optimal image dataset to construct an initial low-resolution DTM.

[0011] For a single ground point in the target area, a stereophotometric error equation is constructed based on the best image dataset of the target area, the initial low-resolution DTM, the illumination model and its model parameters, and the reconstruction result of the single ground point is determined iteratively based on the stereophotometric error equation.

[0012] The reconstruction results of individual ground points in the target area are merged to generate a refined reconstruction map.

[0013] Furthermore, the remote sensing images in the aforementioned remote sensing image dataset are filtered to obtain the optimal image dataset for the target area, including:

[0014] Based on the various parameter values ​​corresponding to each remote sensing image in the remote sensing image dataset, at least two first-type remote sensing images and at least three second-type remote sensing images for stereo positioning are determined from the remote sensing image dataset. The first-type remote sensing images have an incident angle ranging from 0 to 20°, an observation angle ranging from 35° to 50°, and an image intersection angle ranging from 70° to 110°. The second-type remote sensing images have an incident angle ranging from 30° to 50°, an observation angle ranging from 0° to 20°, and camera azimuth angles including the four directions of east, south, west, and north.

[0015] The set consisting of all the first-class and second-class remote sensing images is determined as the optimal image dataset for the target area.

[0016] Furthermore, the error equation for stereophotometry is constructed, including:

[0017] Based on the optimal image dataset of the target area, the initial low-resolution DTM, and the illumination model and its parameters, the image coordinates, remote sensing image grayscale values, and model grayscale values ​​corresponding to the individual ground point are determined, and then the stereophotometric error equation is constructed.

[0018] Furthermore, when the illumination model is the Lunar-Lambert illumination model, the calculation formula corresponding to the stereophotometric error equation is:

[0019]

[0020]

[0021] Among them, G C(x′,y′) represents the model grayscale value corresponding to the ground point in the j-th image; (x′,y′) represents the image coordinates corresponding to the ground point; A N denoted by α, i represents the incident angle, e represents the observation angle; Λ(α) represents the influence factor of the illumination model; T represents the multiplicative factor of atmospheric transport effects; H represents the additive factor of atmospheric transport effects; g j (x′, y′) represents the grayscale value of the ground point in the j-th image; v j (x′, y′) represents g j The residual of (x′,y′); Indicates the elevation value; The normalized albedo is represented by k and l, which represent the row and column indices of the grid points corresponding to the ground points, respectively.

[0022] Furthermore, the method also includes:

[0023] Based on each candidate lighting model and its model parameters, at least two different methods are used to perform stereophotometric 3D reconstruction to obtain the 3D reconstruction results corresponding to each candidate lighting model and its model parameters.

[0024] Based on the 3D reconstruction results corresponding to each candidate lighting model and its parameters, an applicability analysis of the lighting model and its parameters is conducted to determine the optimal lighting model and its parameters.

[0025] Furthermore, based on the selected illumination models and their parameters, stereophotometric 3D reconstruction was performed in the experimental area using the original resolution images, and global stereophotometric 3D reconstruction was performed using the sampled images obtained by downsampling the original resolution images.

[0026] Furthermore, the formula for calculating the overall cost function of stereophotometric 3D reconstruction is as follows:

[0027]

[0028] Where φ(X,Y) represents the terrain; I k (φ)(X,Y) represents the image grayscale value corresponding to the ground point when back-projected onto the k-th image, T k Represents the exposure parameters on the k-th image; A(X,Y) represents the albedo of the pixel on the k-th image when the ground point is back-projected onto it; R k (φ)(X,Y) represents the reflection coefficient; φ(X,Y) represents the sum of squares of all second-order partial differential terms; μ represents the smoothing term parameter; λ represents the adjustment parameter for the degree of deviation between the reconstructed terrain and the reference terrain; φ0(X,Y) represents the reference terrain.

[0029] This invention offers the following advantages: Photogrammetry boasts high efficiency in 3D reconstruction and the ability to directly calculate absolute elevation, but the resolution of the generated terrain is relatively low. Stereophotometry, on the other hand, can reconstruct more detailed 3D terrain, but it calculates terrain slope rather than absolute elevation, and accuracy control during large-scale processing requires reference to the terrain as a constraint. Therefore, this invention effectively improves the accuracy and quality of planetary surface digital elevation model reconstruction by fusing photogrammetry and stereophotometry to construct a high-resolution planetary surface DEM. Attached Figure Description

[0030] To more clearly illustrate the technical solutions and advantages in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0031] Figure 1 This is a flowchart of a planetary surface DEM generation method that integrates photogrammetry and stereophotometry according to an embodiment of the present invention;

[0032] Figure 2 This is a schematic diagram of the image incident angle, observation angle, and phase angle according to an embodiment of the present invention;

[0033] Figure 3 This is a flowchart of the stereophotogrammetry method according to an embodiment of the present invention;

[0034] Figure 4 This is a flowchart of a stereophotometric three-dimensional reconstruction method according to an embodiment of the present invention;

[0035] Figure 5 This is a flowchart illustrating the determination of the lighting model and model parameters in an embodiment of the present invention.

[0036] Figure 6 This is a schematic diagram of a geometrically corrected remote sensing image of Vesta, according to an embodiment of the present invention.

[0037] Figure 7 This is a reconstructed local DTM of Vesta from an embodiment of the present invention.

[0038] Figure 8 This is a high-resolution DTM of the lunar south polar region according to an embodiment of the present invention. Detailed Implementation

[0039] To further illustrate the technical means and effects adopted by the present invention to achieve the intended purpose, the specific implementation methods, structures, features, and effects of the technical solution proposed according to the present invention are described in detail below with reference to the accompanying drawings and preferred embodiments. In the following description, different "one embodiment" or "another embodiment" do not necessarily refer to the same embodiment. Furthermore, specific features, structures, or characteristics in one or more embodiments can be combined in any suitable form.

[0040] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art to which this invention pertains. Furthermore, all parameters or indices in the formulas discussed herein are normalized values ​​to eliminate the influence of dimensions.

[0041] To address the low resolution and accuracy of existing stereophotometric methods for 3D reconstruction of planetary surfaces, this embodiment provides a method for generating a planetary surface DEM that integrates photogrammetry and stereophotometry. This method combines photogrammetry and stereophotometry for 3D planetary reconstruction, and the corresponding flowchart is shown below. Figure 1 As shown, it includes the following steps:

[0042] Obtain reconstruction preparation data, which includes at least: remote sensing image dataset, various parameter values ​​corresponding to each remote sensing image in the remote sensing image dataset, initial low-resolution DTM, and illumination model and its model parameters. The various parameter values ​​include at least: incident angle, observation angle, camera azimuth angle, and image intersection angle.

[0043] Construct a target region, and based on the parameter values ​​corresponding to each remote sensing image in the remote sensing image dataset, filter the remote sensing images in the remote sensing image dataset to obtain the optimal image dataset for the target region.

[0044] For a single ground point in the target area, a stereophotometric error equation is constructed based on the best image dataset of the target area, the initial low-resolution DTM, the illumination model and its model parameters, and the reconstruction result of the single ground point is determined iteratively based on the stereophotometric error equation.

[0045] The reconstruction results of individual ground points in the target area are merged to generate a refined reconstruction map.

[0046] The following, in conjunction with the accompanying drawings, provides a detailed description of the specific implementation steps for each of the above steps.

[0047] First, in the aforementioned method for generating a planetary surface DEM by fusing photogrammetry and stereophotometry, in order to obtain a high-quality 3D shape model, the images need to be screened according to the requirements of stereophotometry and photogrammetry techniques. That is, at least five images are used for 3D reconstruction of the target area. Four of these images require the detector to take pictures from the east, south, west, and north directions relative to the target area at an observation angle of about 30°. The ideal value for the azimuth angle variation of these four images is about 90°, which is beneficial for extracting terrain slope information in the longitude and latitude directions and can supplement shadow areas. The other image requires a low incidence angle (i.e., close to local noon), but the phase angle should be avoided so that the grayscale information of the image is mainly affected by the ground albedo rather than the terrain undulation.

[0048] In terms of illumination models, commonly used models in stereophotometric planetary 3D modeling include Lambert, Lunar–Lambert, Hapke, and Chron models. Due to significant differences in surface physical properties and albedo among different types of asteroids, the applicability of existing illumination models and parameters to planets often requires extensive experimentation. Currently, the illumination models and parameters used in stereophotometric planetary 3D reconstruction are mainly determined based on experience and remote sensing data. This embodiment adaptively determines the optimal illumination model and parameters by performing 3D reconstruction on a small number of experimental areas and quantitatively analyzing the effect of the 3D reconstructed terrain.

[0049] Specifically, the above-mentioned method for generating a planetary surface DEM by fusing photogrammetry and stereophotometry mainly includes the following three parts:

[0050] Part 1: Stereo Pair Optimization

[0051] Assuming the light source is a point source at infinite distance and its direction is known (the light source for planetary surface remote sensing images is usually the sun, thus satisfying this condition), the grayscale values ​​in the image space can be calculated according to the imaging equation, that is, the grayscale values ​​can be expressed as a function of terrain slope and surface illumination properties. Therefore, the model grayscale values ​​of ground points can be calculated from albedo and terrain information based on the imaging equation. Stereophotometry uses the grayscale values ​​of the remote sensing image itself as observation values, combines them with the model grayscale values ​​calculated by the imaging equation to establish an error equation (or observation equation), and then performs iterative adjustment to reconstruct the three-dimensional terrain.

[0052] Selecting the optimal image dataset from a large number of remote sensing images is the first problem to be solved in asteroid 3D reconstruction processing. Image screening can reduce the processing time of subsequent steps and improve processing efficiency; on the other hand, it can also remove invalid or redundant images, improving the accuracy of 3D reconstruction. During image screening, image metadata information can be used, combined with the requirements of photogrammetry and stereophotometry for lighting and observation conditions, to select the best image dataset for the target area.

[0053] Figure 2 This diagram illustrates the incident angle, observation angle, and phase angle of an image during remote sensing. Where i is the incident angle, e is the observation angle, and g is the phase angle. The incident angle is the angle between the line connecting the light source (sun) to a point on the ground and the normal to the ground point. The observation angle is the angle between the line connecting the camera ray to a point on the ground and the normal to the ground point. The phase angle is the angle between the line connecting the ground point to the light source (sun) and the line connecting the ground point to the camera ray. The solar azimuth angle is the angle between the projection of sunlight onto the ground plane and the local meridian, typically measured clockwise from the north direction of the target point to the direction of sunlight incidence.

[0054] Given the detector's position and attitude, as well as the light source's position, the incident angle and observation angle can be calculated from the DTM's height value:

[0055]

[0056] Where i represents the incident angle; e represents the observation angle; n represents the normal vector of the ground point P(X,Y,Z); s represents the light source direction vector of the ground point P(X,Y,Z); and v represents the observation direction vector of the ground point P(X,Y,Z).

[0057] Photogrammetric processing requires sufficient illumination during remote sensing image acquisition to ensure the images contain enough information for matching. A large angle of incidence results in poor illumination and numerous shadowed areas; conversely, a small angle of incidence allows direct sunlight to distort the surface radiation. The observation angle indicates the camera's tilt. Photogrammetric processing requires a certain tilt in the stereo image to create a suitable intersection angle, but this tilt cannot be too large to limit image distortion and resolution differences. Based on experience from Rosetta, OSIRIS-REx, and other small body exploration projects, as well as photogrammetric processing in Earth observation, the illumination and observation conditions required for asteroid remote sensing image photogrammetric processing are shown in Table 1. Furthermore, to facilitate image matching, the difference in illumination conditions between stereo images should generally not be too large; for example, the difference in angle of incidence is typically less than 10°.

[0058] Table 1 Optimal lighting and observation conditions for photogrammetric processing.

[0059] Angle of incidence 5–70° Observation angle 0–55° Phase angle 5–125° Stereoscopic image intersection angle 20–45°

[0060] Stereophotometry typically requires at least five images to cover the target area, with two for stereo geometric localization and the other three for acquiring ground radiometric information from different light source locations. Stereophotometry requires remote sensing images of the target area to be acquired from different photographic directions and illumination angles. Therefore, to achieve the best 3D reconstruction results using stereophotometry, it is necessary to analyze not only the incident and observation angles of the images but also information such as the detector azimuth angle (camera azimuth angle). Based on existing literature and preliminary practical experience, the requirements for remote sensing images in this embodiment of stereophotometry are shown in Table 2.

[0061] Table 2 Optimal illumination and observation conditions for stereophotometric processing.

[0062]

[0063] Part Two: Preliminary DEM Generation via Photogrammetry

[0064] Stereophotogrammetry (SPG) is theoretically mature and reliable in engineering applications, but it has strict requirements for input data, requiring information such as camera geometry, probe position, and attitude. The processing flow for planetary 3D reconstruction based on stereophotogrammetry is as follows: Figure 3 As shown, the specific steps are as follows:

[0065] (1) Preprocessing: Preprocess the planetary remote sensing images (radio calibration, etc.), and use the pre-calibrated sensor placement error parameters and small celestial body rotation parameters to convert the position and attitude information of the detector into the position and attitude information of the camera (i.e., exterior orientation elements); and use the pre-calibrated camera parameters (i.e., interior orientation elements) to build a camera model.

[0066] (2) Stereo image pair selection: Combining the metadata information of the image (incident angle, emission angle, phase angle, resolution, etc.) and the initial camera position and attitude information, the overlap and intersection angle of the stereo images are calculated, and image pairs that meet the stereo photogrammetry conditions are automatically selected, such as: the incident angle is 5 to 70°, the stereo image intersection angle is 10 to 50°, the resolution difference is less than 2 times, and the image overlap is greater than 20%.

[0067] (3) Establishing a tie point control network: The construction of the tie point control network is the basis of photogrammetric bundle adjustment. In order to improve the accuracy and reliability of tie point matching, a matching strategy based on the space of the approximate orthophoto is adopted. First, the original image is corrected to an approximate orthophoto using the initial camera position and attitude data. Then, matching is performed on the approximate orthophoto. The matching results are then converted to the original image for photogrammetric processing. The scale variation problem in the matching of small celestial bodies is solved by comprehensively using methods such as normalized cross correlation (NCC), image pyramid matching, and SIFT / SURF feature matching. Gross errors in tie points are eliminated by using the RANSAC algorithm, terrain continuity conditions, and epipolar geometric constraints.

[0068] (4) Bundle Adjustment: An adjustment mathematical model is established based on the collinearity condition equation. The initial position and attitude information of the camera are adjusted together as weighted observations. The second-order polynomial error correction model is used to correct the orbit position and attitude measurement errors. A large number of remote sensing images of small celestial bodies are pre-adjusted to identify individual images with abnormal accuracy. Then, virtual control point technology is used to convert reliable connection points in the overall control network to individual images with abnormal accuracy, thereby improving the overall adjustment result quality and forming refined exterior orientation elements.

[0069] (5) Generate 3D digital elevation model products: Based on the refined exterior orientation element information, the stereo image is densely matched to generate digital elevation model products, and then the DEM is spliced ​​and mosaicked.

[0070] Third step: Stereophotometric 3D reconstruction

[0071] The calculation of the stereophotometric imaging equation requires the use of a lighting model. Taking the Lunar-Lambert lighting model, commonly used in the field of 3D reconstruction of planetary surfaces, as an example, its mathematical formula is as follows:

[0072]

[0073] Where i represents the incident angle; e represents the observation angle; α is the phase angle; R(i,e,α) is the irradiance corresponding to the ground point P(X,Y,Z); A N (X,Y,α) represents the normalized or standardized albedo; the Lunar-Lambert illumination model is essentially a combination of the Lambert model and the Lommel-Seeliger model, and Λ(α) represents the factor that adjusts the influence of the Lambert and Lommel-Seeliger illumination models, with an empirical value of Λ(α) = 1.0 + A·α + B·α 2 +C·α 3, A=-0.019; B=0.242E-3; C=-1.46E-6.

[0074] The imaging equation is further extended here, taking into account the camera offset and gain effects, as well as the two spatially invariant parameters T and H describing atmospheric effects, resulting in the following imaging equation:

[0075]

[0076] Where G(x′,y′) represents the model gray value of image point p calculated by the imaging equation; (x′,y′) represents the image coordinates of image point p, which can be obtained by back projection from ground point P(X,Y,Z) according to the collinearity condition equation; O(x′,y′) represents the camera offset; k(x′,y′) represents the camera gain factor; n represents the light attenuation index; γ represents the angle between the optical axis and the line connecting image point p and the corresponding ground point P(X,Y,Z); d represents the camera aperture; f represents the lens focal length; E s The irradiance of the scene is represented by T; T and H represent the multiplicative and additive factors of the atmospheric transport effect, respectively.

[0077] Let the grayscale value of the model of the radiometric calibration image be G. C Given (x′, y′), the mathematical form of the above imaging equation can be expressed as:

[0078]

[0079] For the Moon and other celestial bodies without atmospheric effects, the above equation can be further simplified to:

[0080] At this point, the model grayscale value of the inspection image depends on the standardized albedo A. N The factors include the influence factor Λ(α), the incident angle, and the observation angle, which can be calculated from the ground elevation value.

[0081] Let DTM grid point (X) k ,Y l The elevation value of ) is Z(X) k ,Y l ), where k and l represent the row and column indices of the grid points, respectively. The error equation for stereophotometry is:

[0082]

[0083] Among them, g j (x′,y′) represents the gray value of image point p(x′,y′) on image j, which can be obtained by interpolation of the gray values ​​of the four adjacent pixels; v j (x′, y′) is defined by g j The residuals of the observed values; It is the elevation value of the terrain surface (unknown); It is the normalized albedo (unknown).

[0084] The above formula relates to elevation values. The nonlinear function needs to be linearized and solved iteratively, therefore elevation values ​​are required. The initial values ​​can be obtained by constructing an initial terrain using photogrammetry in a dense matching manner. This initial terrain is also called the initial low-resolution DTM, and the elevation values ​​are determined based on the initial low-resolution DTM. The initial value, and It is in linear form in the formula and does not require initial values.

[0085] The process of stereophotometric 3D reconstruction of asteroid surface is as follows: Figure 4 As shown, the specific processing steps are as follows:

[0086] ① Prepare input data, including radiometrically calibrated remote sensing image datasets, image orientation elements inside and outside the image, photography time, light source location, and various information such as illumination model and model parameters.

[0087] ② Construct a Maplet (i.e., a small rectangular target area) and filter images. Based on the sampling interval of the DTM to be generated, Maplets are typically constructed at 100×100 pixels, with adjacent Maplets overlapping by 30% in both longitude and latitude. Observation conditions (observation angle, camera azimuth angle, etc.) and illumination conditions (incident angle, solar azimuth angle, etc.) are calculated by combining information such as photography time, light source position, and exterior orientation elements of the images. The intersection angle of the stereo images is determined by matching corresponding points and calculating the angles between the photographic rays from those points. Images that meet the requirements of stereophotometric processing are then selected. Specifically, based on the optimal illumination and observation conditions for stereophotometric processing in Table 2 above, the radiometrically calibrated remote sensing image dataset is filtered to select the best image dataset that meets the conditions in the table. The above stereo photogrammetry method is used to construct an initial low-resolution DTM using the best image dataset, also known as the initial DTM.

[0088] ③ For each ground point on the Maplet, the image coordinates and remote sensing image grayscale values ​​corresponding to the ground point are calculated using the existing collinearity condition equation and the initial low-resolution DTM back projection. The model grayscale value of the ground point is calculated using the imaging equation, and regularization terms such as smoothness constraints and photogrammetric constraints are added. Then, the stereophotometric error equation is constructed.

[0089] ④ Based on the stereophotometric error equation, a single Maplet is subjected to stereophotometric adjustment according to the least squares principle. The adjustment output can refine the digital elevation model, thereby obtaining the slope map and albedo map. The slope map is integrated to obtain the elevation map. When the various indicators in the iteration process meet the accuracy requirements or reach the maximum number of iterations, the iteration calculation is terminated.

[0090] ⑤ Establish the transformation relationship between adjacent regions by using feature points in overlapping regions, merge the 3D reconstruction results of multiple Maplets, and generate a high-resolution large-scale DTM.

[0091] The stereophotometric processing described above demonstrates that the determination of the illumination model and its parameters directly affects the 3D reconstruction results. For asteroids with unknown surface physical properties to be explored, this embodiment first constructs an illumination model and its parameter combination based on experience. Then, two methods are used to quantitatively analyze the results of the stereophotometric 3D reconstruction. The process is as follows: Figure 5 As shown. Method 1 utilizes original resolution imagery in several small test areas to perform coarse-resolution 3D reconstruction using photogrammetry (approximately 3-5 times the horizontal resolution of the original imagery). Method 2 uses low-resolution imagery obtained by downsampling the original imagery (5-10 times the resolution of the original imagery) for global stereophotometric 3D reconstruction. Finally, using the topographic information reconstructed by photogrammetry as a reference, the optimal lighting model and model parameters for stereophotometric reconstruction are determined by comparing and analyzing the number of iterations, the magnitude of the stereophotometric adjustment residuals, and the consistency of geometric accuracy of the reconstructed topography in adjacent areas. The geometric accuracy of the reconstructed topography in adjacent areas is determined by rendering the DTM (Digital Transformer Model) with lighting, converting it into a simulated image, and then matching the DTM-rendered images. By analyzing the differences in the coordinates of the matching points, the geometric accuracy is automatically calculated.

[0092] Photometric measurements essentially solve for slope information. To ensure the absolute accuracy of stereophotometric 3D reconstruction over large areas, a regularization term constraining the deviation between the reconstructed terrain and the reference terrain needs to be added to the cost function. This involves fusing a low-resolution 3D terrain model generated by photogrammetry. This low-resolution 3D terrain model refers to the 3D model determined using photogrammetric reconstruction methods based on the best image dataset. Additionally, to ensure a certain degree of continuity in the reconstructed terrain, a regularization term controlling the smoothness of the terrain can also be added to the cost function. Therefore, the overall cost function of stereophotometric processing can be expressed as follows:

[0093]

[0094] Where φ(X,Y) represents the terrain; I k(φ)(X,Y) represents the image grayscale value corresponding to the ground point backprojected onto the k-th image; T k Represents the exposure parameters on the k-th image; A(X,Y) represents the albedo of the pixels projected onto the k-th image from the ground points; R k (φ)(X,Y) represents the reflection coefficient; φ(X,Y) represents the sum of squares of all second-order partial differential terms; μ represents the smoothing term parameter; λ represents the adjustment parameter for the degree of deviation between the reconstructed terrain and the reference terrain; φ0(X,Y) represents the reference terrain.

[0095] Experimental verification of three-dimensional reconstruction of planets based on stereophotometry. Figure 6 and Figure 7 This is an experiment using stereophotometric processing of Vesta images obtained from the Dawn spacecraft. Among them, Figure 6 This is a geometrically corrected remote sensing image of Vesta; Figure 7 The reconstructed local DTM of Vesta (rendered using the Hillshade method).

[0096] To address the need for high-precision navigation maps of the lunar south polar landing area during the fourth phase of the lunar exploration program, 75 US LRO NAC lunar remote sensing images (resolution 0.8–1.0 m / pixel) were processed using stereophotometry to create a high-resolution DTM for the lunar south polar region with a grid spacing of 1 m. Figure 8 As shown. In Figure 8 In the image, the background is an LRO NAC orthophoto map with a resolution of 1.0 m / pixel. The color display area is a rendered DTM with a range of 1000 m × 1000 m.

[0097] This invention establishes an adaptive stereophotometric illumination model of the planetary surface, integrating the advantages of photogrammetry and stereophotometric 3D reconstruction. Image position and attitude data refined through photogrammetric processing can be used to calculate the imaging equations for stereophotometry. The image position and attitude include the X, Y, and Z coordinates of the image in the world coordinate system and the attitude of the image coordinate system relative to the world coordinate system, typically represented using rotation matrices or quaternions. Furthermore, the low-resolution digital elevation model generated by photogrammetry can serve as initial values ​​and external terrain constraint information for stereophotometric 3D modeling, while stereophotometry theoretically can reconstruct high-resolution, detailed terrain. Therefore, this invention, integrating photogrammetry and stereophotometry, can acquire more refined 3D terrain information of the planetary surface, providing technical support for probe landing site selection, rover path planning, and scientific research on planetary geology and geomorphology.

[0098] It should be noted that the above embodiments are only used to illustrate the technical solutions of this application, and are not intended to limit them. Although this application has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of this application, and should all be included within the protection scope of this application.

Claims

1. A method for generating a planetary surface DEM by integrating photogrammetry and stereophotometry, characterized in that, Includes the following steps: Obtain reconstruction preparation data, which includes at least: remote sensing image dataset, various parameter values ​​corresponding to each remote sensing image in the remote sensing image dataset, and illumination model and its model parameters. The various parameter values ​​include at least: incident angle, observation angle, camera azimuth angle and image intersection angle. Construct a target region, and based on the parameter values ​​corresponding to each remote sensing image in the remote sensing image dataset, filter the remote sensing images in the remote sensing image dataset to obtain the optimal image dataset for the target region, and use the optimal image dataset to construct an initial low-resolution DTM. For a single ground point in the target area, a stereophotometric error equation is constructed based on the best image dataset of the target area, the initial low-resolution DTM, the illumination model and its model parameters, and the reconstruction result of the single ground point is determined iteratively based on the stereophotometric error equation. The reconstruction results of individual ground points in the target area are merged to generate a refined reconstruction map; Specifically, based on the various parameter values ​​corresponding to each remote sensing image in the remote sensing image dataset, at least two first-type remote sensing images and at least three second-type remote sensing images for stereo positioning are determined from the dataset. The first-type remote sensing images have an incident angle ranging from 0 to 20°, an observation angle ranging from 35° to 50°, and an image intersection angle ranging from 70° to 110°. The second-type remote sensing images have an incident angle ranging from 30° to 50°, an observation angle ranging from 0° to 20°, and camera azimuth angles including east, south, west, and north. The set consisting of all first-type and second-type remote sensing images is determined as the optimal image dataset for the target area. Specifically, based on each candidate lighting model and its parameters, at least two different methods are used to perform stereophotometric 3D reconstruction to obtain the 3D reconstruction results corresponding to each candidate lighting model and its parameters; based on the 3D reconstruction results corresponding to each candidate lighting model and its parameters, an applicability analysis of the lighting model and its parameters is performed to determine the optimal lighting model and its parameters. Specifically, based on the selected illumination models and their parameters, stereophotometric 3D reconstruction is performed in the experimental area using the original resolution images, and global stereophotometric 3D reconstruction is performed using the sampled images obtained by downsampling the original resolution images.

2. The method for generating a planetary surface DEM by fusing photogrammetry and stereophotometry according to claim 1, characterized in that, Constructing the error equation for stereophotometry, including: Based on the optimal image dataset of the target area, the initial low-resolution DTM, and the illumination model and its parameters, the image coordinates, remote sensing image grayscale values, and model grayscale values ​​corresponding to the individual ground point are determined, and then the stereophotometric error equation is constructed.

3. The method for generating a planetary surface DEM by fusing photogrammetry and stereophotometry according to claim 2, characterized in that, When the illumination model is the Lunar-Lambert illumination model, the calculation formula corresponding to the stereophotometric error equation is: ； ； in, This represents the grayscale value of the model corresponding to the ground point in the j-th image; This represents the image coordinates corresponding to the ground point; Represents standardized albedo; Indicates the angle of incidence. Indicates the observation angle; Indicates the influence factors of the illumination model; The multiplicative factor representing atmospheric transport effects; An additive factor representing atmospheric transport effects; This represents the grayscale value of the ground point in the j-th image; express The residual; Indicates the elevation value; Represents normalized albedo; These represent the row and column indices of the grid point corresponding to the ground point.

4. The method for generating a planetary surface DEM by fusing photogrammetry and stereophotometry according to claim 1, characterized in that, The formula for calculating the overall cost function of stereophotometric 3D reconstruction is as follows: ； in, Indicates terrain; This indicates projecting the ground point back onto the first... The corresponding grayscale value on the image. Indicates the first Exposure parameters on the image; This indicates projecting the ground point back onto the first... Albedo of pixels in an image; Indicates the reflection coefficient; express The sum of squares of all second-order partial differential terms; Indicates the smoothing term parameter; An adjustment parameter indicating the degree of deviation between the reconstructed terrain and the reference terrain; Indicates the reference terrain.

Images