Mars multi-source orbiter image joint adjustment method based on generalized control constraint

Abstract

The invention relates to a Mars multi-source orbiter image joint adjustment method based on generalized control constraint, which comprises the following steps of: extracting plane control points for a low-resolution reference DOM (Document Object Model) and a Mars orbiter optical image, and matching to obtain connection points; based on the fitting residual error of the image rational function model of the Mars orbiter, establishing an image space compensation model of a superimposed polynomial and a cubic spline function, and obtaining an initial compensation parameter; and according to the plane control points, the connection points and the initial compensation parameters, combining plane and elevation control constraints to construct an adjustment model, and solving the adjustment model by adopting a self-adaptive weighting strategy to obtain geometric positioning parameters of the image. Compared with the prior art, the method has the advantages that multi-source image united bundle adjustment is carried out through cooperation of the plane weak control constraint and the elevation control constraint, the absolute and relative positioning precision of the Mars multi-source orbiter image can be effectively improved, and accurate registration of the multi-source orbiter image is achieved.

Description

Technical Field

[0001] The present invention relates to the technical field of image registration, and in particular to a Mars multi-source orbiter image joint adjustment method based on generalized control constraints. Background Art

[0002] High-resolution Mars topography is an important means to understand the topography and structure of Mars, carry out Mars exploration projects and scientific research. With the development of Mars mapping technology and the development of exploration activities, a large number of multi-task and multi-coverage Mars remote sensing images have been accumulated.

[0003] Currently, among the orbiter optical images used to generate Mars planetary mapping products, the medium-resolution images that can cover the entire moon can only support the acquisition of meter-scale surface features at best, while the sub-meter-level high-resolution images used to generate high-precision local mapping products have limited coverage and uneven distribution. Therefore, the coordinated use of multi-source images is conducive to capturing more comprehensive and accurate information on the morphological features of Mars, and repeated observations of multi-source images can eliminate random errors, which is conducive to further improving the geometric positioning accuracy.

[0004] However, due to the limited accuracy of orbiter attitude and orbit determination, there are widespread geometric positioning inconsistencies between images. Bundle adjustment is the main method to eliminate inconsistencies between images and realize the coordinated use of multi-source images. In the absence of measured ground control points on Mars, it is necessary to introduce control information in the adjustment process to eliminate orbit and attitude errors, improve solution stability and image geometric positioning accuracy. Using only the Mars Orbiter Laser Altimeter (MOLA) and its Digital Elevation Model (DEM) as the Mars reference for control adjustment will be affected by their limited resolution and planar accuracy. The planar inconsistency between the two can be minimized by selecting identifiable surface features in optical images and aligning them with the prior DEM. However, this method has a low degree of automation and is affected by factors such as image distortion, jitter, imaging model differences, and illumination changes, making it difficult to directly obtain reliable feature points between heterogeneous data for accurate alignment. Summary of the invention

[0005] The purpose of the present invention is to overcome the defects of the above-mentioned prior art and to provide a Mars multi-source orbiter image joint adjustment method based on generalized control constraints, which can improve the geometric positioning accuracy and consistency of multi-source orbiter images and quickly and automatically align multi-source images to a unified reference benchmark.

[0006] The object of the present invention can be achieved by the following technical solution: A Mars multi-source orbiter image joint adjustment method based on generalized control constraints, comprising the following steps:

[0007] S1. Extract plane control points and match them to obtain connection points based on low-resolution reference DOM and Mars orbiter optical images;

[0008] S2. Based on the fitting residual of the rational function model of the Mars orbiter image, an image-space compensation model of superimposed polynomial and cubic spline function is established and initial compensation parameters are obtained;

[0009] S3. According to the plane control points, connection points and initial compensation parameters, combined with the plane and elevation control constraints, an adjustment model is constructed, and an adaptive weighted strategy is used to solve the adjustment model to obtain the geometric positioning parameters of the image.

[0010] Furthermore, the step S1 includes the following steps:

[0011] S11, preprocessing the optical image of the Mars orbiter to obtain a preprocessed optical image;

[0012] S12, extracting connection points between images through least squares template matching processing for the pre-processed optical image;

[0013] S13, combining the pre-processed optical image and the low-resolution reference DOM, outputting the crater center through image segmentation, mask segmentation, and crater fitting processing;

[0014] The crater centers are then screened, and the matched crater centers are obtained through point set registration to obtain the plane control points.

[0015] Furthermore, in step S12, for optical images with different illumination, a Histogram of Oriented Phase Consistency (HOPC) is used as a similarity measure to perform least squares template matching;

[0016] For optical images with consistent lighting conditions, the normalized cross correlation (NCC) is used as the similarity measure for least squares template matching.

[0017] Furthermore, the specific process of step S13 includes:

[0018] Contrast-limited adaptive histogram equalization and gamma correction are used to perform radiometric enhancement on the preprocessed optical image and low-resolution reference DOM to strengthen the crater rim.

[0019] The enhanced image is extracted and segmented using the SAM model (Segment Anything Model).

[0020] In the process of crater extraction, the segmentation masks that do not meet the shape requirements are screened out, and the roundness and ellipticity of the mask are evaluated by calculating the area and circumference of the mask. When the ratio of the circumference to the circumference and the ratio of the major and minor axes exceed the threshold of 0.9, the mask is judged to be a circle; when the ratio of its major and minor axis area to the actual area exceeds the threshold of 0.8, the mask is judged to be an ellipse; the fitting accuracy of the circle or ellipse is evaluated by edge detection, and only the masks with high fitting accuracy and accurate center position are retained;

[0021] Then, an improved point set registration method that retains global and local structures is used to realize the registration of crater positions between images. First, based on the prior geographic information of the craters, the target centers that lack candidate targets within the distance range are screened out. At the same time, local features are enhanced by constructing a feature descriptor that integrates the aspect ratio, fitting area ratio, and shape context features to ensure that the global and local structures are retained during the matching process to complete the registration.

[0022] Furthermore, the specific process of step S2 is as follows:

[0023] The elevation of the image area is divided into several elevation surfaces, and a three-dimensional virtual control grid of the image is established in the image space with a set interval as the virtual control point of the image space. On this basis, the image coordinates of the virtual control grid are projected onto each elevation surface according to the strict imaging model to obtain the virtual control point of the object space, and the Rational Polynomial Coefficient (RPC) is calculated by fitting, and the geometric relationship between the image space and the object space is expressed by the rational function model.

[0024] The image space compensation model of linear and nonlinear errors of images is constructed by using superimposed polynomials and cubic spline functions to obtain initial compensation parameters for compensating nonlinear image errors.

[0025] Furthermore, the expression of the rational function model is:

[0026]

[0027] Among them, (X o ,Y o ) and (X s ,Y s ) represent the offset and scale value of the normalized pixel coordinates (X, Y), (P, L, H) are the latitude, longitude and elevation of the normalized object point, Num s , Num l ,Den s and Den l It is a third-order polynomial with a total of 20 coefficients, namely the rational function model parameter RPC.

[0028] Furthermore, the expression of the image space compensation model is:

[0029]

[0030] Among them, the parameters (e0, e1, e2, f0, f1, f2) represent the affine transformation used to compensate for the linear deviation. According to the time variation characteristics caused by posture jitter, the average residual of the image points in the same row is used as the optimal estimate of the RPC nonlinear deviation. The corrected pixel coordinates (X ′ ,Y ′ ) is calculated as follows:

[0031]

[0032] in, are the unknown parameters of the cubic spline function, which are used to compensate for nonlinear image errors, and i=1,2…m-1 is the number of connection points.

[0033] Further, the process of constructing the adjustment model in step S3 includes:

[0034] The error equation is established and the regional block adjustment is performed. According to the tie points, virtual control points, plane weak control constraints based on the center of the crater, terrain elevation constraints and rational function model parameters RPC obtained by multi-source image matching, the error equation of the joint adjustment model is constructed, which includes the following four observation equations:

[0035]

[0036] Among them, TPs are connection points, ECPs are elevation control points, VCPs are virtual control points, PCPs are plane control points, V TPs is the observation equation of the connection point, V ECPs is the observation equation of the elevation control constraint, V PCPs is the observation equation of the plane control constraint, V VCPs is the observation equation of the virtual control point, t is the correction number of the adjustment parameter, x is the unknown ground coordinate of the connection point, matrices A and B are the error equation coefficient matrices, L is the constant term, P i,i＝1,2,3,4 For the corresponding weight, for the error equation, construct the normal equation, and use the least squares method to fit and obtain the object coordinates of the most approximate connection points;

[0037] Then, terrain constraints based on the prior DEM are introduced, and elevation control is achieved by minimizing the sum of squares of the distances from the connection points to the DEM surface, which corresponds to the second observation equation V ECPs The elevation error equation V h The expression is as follows:

[0038]

[0039] Among them, h0 is the elevation of the connection point, (lat0, lon0) is the geographic plane coordinate, and DEM (lat0, lon0) is the elevation value obtained by interpolation at (lat0, lon0). and are the slopes of DEM in latitude and longitude respectively, (Δlat, Δlon, Δh) are the corrections of the ground coordinates of the connection points;

[0040] A plane weak control constraint is constructed to introduce the terrain change information near the plane control point into the plane control constraint, and the object space 3D coordinates of the control point are adaptively adjusted to achieve automatic registration of the image with the center of the crater in the DOM. The difference between the pixel coordinates of the plane control point and the corresponding 3D coordinates in the reference map is minimized by the following formula:

[0041]

[0042] Where (X, Y) is the pixel coordinate of the center of the crater in the orbiter image, RPC X (lat,lon,DEM) and RPC Y (lat, lon, DEM) is the pixel coordinate of the center of the crater obtained by reprojecting the object space coordinates in the reference DOM and DEM onto the image, which is consistent with the third observation equation V PCPs The corresponding plane error equation v X and v0 is expressed as follows:

[0043]

[0044] in, as well as and The derivatives of the image coordinates with respect to latitude, longitude and elevation are obtained for reprojection. (X0, Y0) are the pixel coordinates of the plane control point. (ΔX′, ΔY′) are the differences between the image coordinates calculated by reprojecting the object coordinates of the plane control point and the corresponding actual image coordinates.

[0045] Virtual control points (VCPs) are introduced to effectively reduce the degrees of freedom of the adjustment model. VCPs are strategically selected by predefined thresholds and grid sizes and introduced as weighted observations into the fourth observation equation V VCPs middle.

[0046] Further, the process of solving the adjustment model using the adaptive weighting strategy in step S3 includes:

[0047] According to the prior matching accuracy σ RPs , define the weight P1 of the connection point as:

[0048]

[0049] The prior terrain is used as control information, and the error and point inconsistency caused by the terrain are used as the "gross error" in the control information. In order to quantitatively reflect the impact of the terrain on the accuracy, the weight function of the elevation control constraint is defined as follows:

[0050]

[0051] Among them, img res is the image resolution, d h is the elevation difference, h mean is the average elevation difference of the data set, σ h is the standard deviation of the elevation difference. When the elevation difference is less than σ h When the height difference exceeds three times σ h When , the weight gradually approaches zero;

[0052] Calculate the scale factor for each adjustment area separately and Where Num tie is the number of connection points, Num pcp is the number of plane control points, Num vcp is the number of virtual control points, considering the crater radius r and the fitting standard deviation f, the weight function of the plane control constraint is defined as follows:

[0053]

[0054]

[0055] Among them, d p is the plane difference, p mean is the average value of the plane difference, σ p is the standard deviation of the plane difference, R max is the normalized radius factor, F max is the normalized fitting factor, the smoothing function σ f As f approaches the set threshold f threshold Gradually reduce the weight to zero, the steepness of the transition curve is controlled by the constant k, when the standard deviation is lower than f threshold When the value is higher than f, the weight decreases rapidly. threshold When , the weight decreases more slowly;

[0056] A dynamic weighting function is established to alleviate the imbalance between the number of virtual control points, and then the unknown ground coordinates of the connection points are calculated according to the least squares method to solve the refined geometric positioning parameters of each image.

[0057] Furthermore, the expression of the dynamic weighting function is:

[0058]

[0059] Where N represents the number of iterations.

[0060] Compared with the prior art, the present invention has the following advantages:

[0061] The present invention first extracts plane control points and matches connection points for low-resolution reference DOM and Mars orbiter optical images; then, based on the fitting residual of the rational function model of the Mars orbiter image, an image square compensation model of superimposed polynomials and cubic spline functions is established and initial compensation parameters are obtained; then, according to the plane control points, connection points and initial compensation parameters, combined with plane and elevation control constraints, an adjustment model is constructed, and an adaptive weighted strategy is used to solve the adjustment model to obtain the geometric positioning parameters of the image. Thus, by coordinating plane weak control constraints and elevation control constraints to carry out multi-source image joint bundle method adjustment, the absolute and relative positioning accuracy of the Mars multi-source orbiter image can be effectively improved, and the accurate registration of the multi-source orbiter image can be achieved.

[0062] When extracting plane control points and matching connection points on the surface of Mars, the present invention adopts contrast-limited adaptive histogram equalization and gamma correction to perform radiation enhancement on sub-satellite images and low-resolution reference DOM to strengthen the edge of the crater, thereby facilitating the subsequent extraction and segmentation of the image using the SAM model. On the other hand, the area and perimeter of the mask are calculated to evaluate its roundness and ellipticity, and then the fitting accuracy of the circle or ellipse is evaluated through edge detection, and only the mask with high fitting accuracy and accurate center position is retained. In addition, in view of the problem that it is difficult to directly match feature points due to differences in imaging between heterogeneous images, an improved point set registration method that retains global and local structures is adopted, which can achieve accurate registration of crater positions between images.

[0063] The present invention adopts a rational function model to express the geometric relationship between the image space and the object space, which is conducive to the joint processing of multi-source images. In view of the problem that the traditional RPC model and low-order polynomials cannot effectively compensate for nonlinear errors, the present invention adopts superimposed polynomials and cubic spline functions to construct a unified compensation model for image linear and nonlinear errors, effectively considering the linear and nonlinear distortion of orbiter images, and ensuring the robustness of the model during adjustment calculation.

[0064] The present invention constructs an adjustment model combining plane weak control constraints, elevation control constraints, image square compensation model, image connection points and virtual control points. Among them, by minimizing the elevation difference between the connection point triangulation result and the prior DEM, an elevation error equation is established, the plane distance is minimized by the terrain change information of the center of the crater, the plane error equation is established, and the reference data is used to construct an adaptive weight function to improve the adjustment accuracy and convergence speed, so as to realize the correction of image sensor and track errors.

[0065] The invention takes into account the different intersection conditions of the connection points between multi-source images, and the non-convergence or instability may occur during the adjustment process. Therefore, an adaptive weighting strategy combining connection points, elevation constraints, plane weak control constraints and virtual control points is designed to solve the adjustment model, which can ensure the robustness of the adjustment, and dynamically coordinate various constraints in combination with the introduced adaptive weighting function, so as to effectively improve the relative and absolute geometric positioning accuracy of multi-source images. BRIEF DESCRIPTION OF THE DRAWINGS

[0066] Figure 1 It is a schematic diagram of the method flow of the present invention;

[0067] Figure 2 It is a schematic diagram of the application process of the embodiment;

[0068] Figure 3 Schematic diagram of crater extraction, screening and registration in the embodiment;

[0069] Figure 4a Schematic diagram of elevation control constraints in the embodiment;

[0070] Figure 4b Schematic diagram of plane weak control constraint in the embodiment;

[0071] Figure 5a The comparison of geometric consistency between DOM before and after adjustment of Gale crater in the embodiment;

[0072] Figure 5b The geometric consistency comparison between the DOM before and after the adjustment of the Jezero crater in the embodiment;

[0073] Figure 5c The geometric consistency comparison between DOM before and after adjustment of the Utopia Plain region in the embodiment;

[0074] Figure 5d This is a 6 m / pixel CTX image mosaic of the Gale crater area in the embodiment;

[0075] Figure 5e This is a 6 m / pixel CTX image mosaic of the Jezero crater area in the embodiment;

[0076] Figure 5f This is a mosaic of 6m / pixel CTX images of the Utopia Plain area in the embodiment;

[0077] Figure 5g The HiRISE mosaic of Gale crater at 1 m / pixel in the embodiment;

[0078] Figure 5h This is a HiRISE mosaic of the Jezero crater area at 1 meter / pixel in the embodiment;

[0079] Figure 5i This is a HiRISE mosaic image of the Utopia Planitia region at 1 meter / pixel in the embodiment. DETAILED DESCRIPTION

[0080] The present invention is described in detail below with reference to the accompanying drawings and specific embodiments.

[0081] Example

[0082] like Figure 1 As shown, a Mars multi-source orbiter image joint adjustment method based on generalized control constraints includes the following steps:

[0083] S1. Extract plane control points and match them to obtain connection points based on low-resolution reference DOM and Mars orbiter optical images;

[0084] S2. Based on the fitting residual of the rational function model of the Mars orbiter image, an image-space compensation model of superimposed polynomial and cubic spline function is established and initial compensation parameters are obtained;

[0085] S3. According to the plane control points, connection points and initial compensation parameters, combined with the plane and elevation control constraints, an adjustment model is constructed, and an adaptive weighted strategy is used to solve the adjustment model to obtain the geometric positioning parameters of the image.

[0086] This embodiment applies the above solution, such as Figure 2 As shown in the figure, it mainly includes three aspects: plane control point extraction and connection point matching, image space compensation model of superposition polynomial and cubic spline function, and adjustment method of comprehensive geometric constraints and adaptive weighting strategy, which can improve the geometric positioning accuracy and consistency of multi-source orbiter images, and quickly and automatically align multi-source images to a unified reference benchmark.

[0087] Specific:

[0088] 1. Plane control point extraction and connection point matching

[0089] In view of the significant difference in imaging mode between the low-resolution reference DOM and the Mars orbiter optical image, which makes it difficult to effectively match them, this scheme proposes an efficient Mars surface plane control point extraction and tie point matching strategy to obtain reliable image relationships and control points. Contrast-limited adaptive histogram equalization and gamma correction are used to perform radiation enhancement on the sub-satellite image (i.e., the image taken vertically, included in the pre-processed optical image) and the low-resolution reference DOM to strengthen the crater edge, so as to facilitate the subsequent extraction and segmentation of the image using the SAM model (Segment Anything Model).

[0090] In the process of crater extraction, considering that most craters are circular or elliptical, it is necessary to screen out segmentation masks that do not conform to these shapes. The specific strategy is to evaluate the roundness and ellipticity of the mask by calculating the area and circumference of the mask. When the ratio of the circumference to the circumference and the ratio of the major and minor axes exceed the threshold of 0.9, the mask is judged to be a circle, and when the ratio of its major and minor axis area to the actual area exceeds the threshold of 0.8, the mask is judged to be an ellipse. The fitting accuracy of the circle or ellipse is evaluated by edge detection, and only masks with high fitting accuracy and accurate center position are retained.

[0091] In order to solve the problem that it is difficult to directly match feature points due to differences in imaging between heterogeneous images, this solution adopts an improved point set registration method that retains global and local structures to achieve accurate registration of crater positions between images. First, based on the prior geographic information of the crater, the target center that lacks candidate targets within the distance range is screened out; at the same time, by constructing a feature descriptor that integrates the major-minor axis ratio, the fitting area ratio, and the shape context features, the local features are enhanced to ensure that the global and local structures are retained during the matching process to complete the registration. Figure 3 The results of crater extraction, screening and registration are shown. The Histogram of Oriented Phase Consistency (HOPC) is used as the similarity measure for least squares template matching between CTX and HiRISE images with different illumination conditions. When the illumination conditions are relatively consistent, the Normalized Cross Correlation (NCC) is used as the similarity measure for least squares template matching to achieve reliable connection point extraction between images.

[0092] 2. Image Square Compensation Model of Superimposed Polynomial and Cubic Spline Function

[0093] This scheme uses a rational function model to express the geometric relationship between the image space and the object space, which is conducive to the joint processing of multi-source images. Its expression is:

[0094]

[0095] Among them, (X o ,Y o ) and (X s ,Y s ) represent the offset and scale value of the normalized pixel coordinates (X, Y), (P, L, H) are the latitude, longitude and elevation of the normalized object point, Num s , Num l ,Den s and Den lIt is a third-order polynomial with 20 coefficients, also known as the rational function model parameter (RPC). By dividing the elevation of the image area into several elevation surfaces, a three-dimensional virtual control grid of the image is established at a certain interval in the image plane as the image plane virtual control point. On this basis, the virtual control grid image coordinates are projected onto each elevation surface according to the strict imaging model to obtain the object space virtual control point, and the RPC can be calculated by fitting.

[0096] In view of the problem that traditional RPC models and low-order polynomials cannot effectively compensate for nonlinear errors, this solution uses superimposed polynomials and cubic spline functions to construct a unified compensation model for image linear and nonlinear errors to ensure the robustness of the model during adjustment calculation. The specific form of the compensation model is as follows:

[0097]

[0098] The parameters (e0, e1, e2, f0, f1, f2) represent the affine transformation used to compensate for the linear deviation. In view of the temporal variation characteristics caused by posture jitter, the average residual of the image points in the same row is taken as the optimal estimate of the RPC nonlinear deviation. The corrected pixel coordinates (X ′ ,Y ′ ) is calculated as follows:

[0099]

[0100] in, are the unknown parameters of the cubic spline function, which are used to compensate for nonlinear image errors, and i=1,2…m-1 is the number of connection points.

[0101] 3. Adjustment method based on geometric constraints and adaptive weighting strategy

[0102] (1) Construction of adjustment model combining plane and elevation control constraints

[0103] The error equation is established and the regional block adjustment is performed. The error equation of the joint adjustment model is constructed based on the tie points, virtual control points, plane weak control constraints based on the center of the crater, terrain elevation constraints, and rational function model parameters obtained by multi-source image matching. It mainly consists of the following four observation equations:

[0104]

[0105] Among them, TPs are connection points, ECPs are elevation control points, VCPs are virtual control points, PCPs are plane control points, V TPs is the observation equation of the connection point, V ECPs is the observation equation of the elevation control constraint, V PCPsis the observation equation of the plane control constraint, V VCPs is the observation equation of the virtual control point. t is the correction number of the adjustment parameter, and x is the unknown ground coordinate of the connection point. The matrices A and B are the error equation coefficient matrices, L is the constant term, and P i,i＝1,2,3,4 For the error equation, construct the normal equation, and use the least squares method to fit and obtain the object coordinates of the most approximate connection points.

[0106] In order to solve the problem of difficulty in obtaining measured control points on Mars and limited absolute positioning accuracy, this scheme introduces terrain constraints based on prior DEM and realizes elevation control by minimizing the sum of squares of distances from connection points to the DEM surface, such as Figure 4a As shown. Corresponding to the second observation equation V ECPs The elevation error equation V h The expression is as follows:

[0107]

[0108] Among them, h0 is the elevation of the connection point, (lat0, lon0) is the geographic plane coordinate, and DEM (lat0, lon0) is the elevation value obtained by interpolation at (lat0, lon0). and are the slopes of DEM in latitude and longitude respectively, and (Δlat, Δlon, Δh) are the corrections of the ground coordinates of the connection points.

[0109] In order to solve the problem of inconsistent fitting positions due to imaging differences of multi-source images in the process of extracting the center of the crater, a plane weak control constraint is constructed. By introducing the terrain change information near the plane control point into the plane control constraint, the object space 3D coordinates of the control point are adaptively adjusted, thereby realizing automatic registration of the image with the crater center in DOM and improving the absolute plane positioning accuracy. Figure 4b The difference between the pixel coordinates of the plane control points and the corresponding three-dimensional coordinates in the reference map can be minimized by the following formula:

[0110]

[0111] Where (X, Y) is the pixel coordinate of the center of the crater in the orbiter image, RPC X (lat,lon,DEM) and RPC Y (lat, lon, DEM) is the pixel coordinate of the center of the crater obtained by reprojecting the object space coordinates in the reference DOM and DEM onto the image. Therefore, it is consistent with the third observation equation V PCPs The corresponding plane error equation v X and v Y The expression is as follows:

[0112]

[0113] in, as well as and The derivatives of the image coordinates with respect to latitude, longitude and elevation are obtained for reprojection. (X0, Y0) are the pixel coordinates of the plane control point, and (ΔX′, ΔY′) are the differences between the image coordinates calculated by reprojecting the object coordinates of the plane control point and the corresponding actual image coordinates.

[0114] In addition, virtual control points (VCPs) are introduced to effectively reduce the degrees of freedom of the adjustment model, thereby solving the rank deficiency problem of the normal equation coefficient matrix and improving the robustness and convergence speed of the adjustment. VCPs are strategically selected through predefined thresholds and grid sizes and introduced as weighted observations into the fourth observation equation V VCPs middle.

[0115] (2) Adaptive weighted strategy method

[0116] Due to the different intersection conditions of the tie points between multi-source images, non-convergence or instability may occur during the adjustment process. To ensure the robustness of the adjustment, it is necessary to develop an adaptive weighting strategy that combines tie points, elevation constraints, plane weak control constraints, and virtual control points. The weight P1 of the tie point can be calculated based on the prior matching accuracy σ TPs The definition is as follows:

[0117]

[0118] Among them, σ TPs is the prior matching accuracy.

[0119] The control effect of the prior terrain as control information varies with the terrain undulation. In flat or gently sloping areas, due to the small terrain undulation, the plane coordinate correction number is small, and the control effect on the plane displacement is weak; while in areas with obvious terrain undulation, the plane coordinate correction number increases and the control effect is enhanced. The impact of local terrain changes on the adjustment results is reflected by formulas (5) and (7). The errors and point inconsistencies caused by the terrain can be regarded as "gross errors" in the control information and need to be effectively eliminated during the adjustment process. In order to quantitatively reflect the impact of terrain on accuracy, the weight function of the elevation control constraint is defined as follows:

[0120]

[0121] Among them, img res is the image resolution, d h is the elevation difference, h mean is the average elevation difference of the data set, σ h is the standard deviation of the elevation difference. hWhen the height difference exceeds three times σ h , the weight gradually approaches zero.

[0122] Due to the imbalance in the number of virtual control points, plane control points, and tie points, it is possible that the elevation control is strong and the plane control is weak. Therefore, the scale factor is calculated for each adjustment area separately. and Where Num tie is the number of connection points, Num pcp is the number of plane control points, Num vcp is the number of virtual control points. Considering the crater radius r and the fitting standard deviation f, the weight function of the plane control constraint is defined as follows:

[0123]

[0124] Among them, d p is the plane difference, p mean is the average value of the plane difference, σ p is the standard deviation of the plane difference. R max is the normalized radius factor, F max is the normalized fitting factor. Smoothing function σ f As f approaches the set threshold f threshold (Empirically set to 0.35) Gradually reduce the weight to zero. The steepness of the transition curve is controlled by the constant k, which is set to 9. When the standard deviation is lower than f threshold When the value is higher than f, the weight decreases rapidly. theeshold When , the weight decreases more gently, which effectively alleviates the influence of outliers and enhances the robustness of the adjustment process.

[0125] In order to speed up the convergence and avoid the problem of weakening the plane control effect of the crater center due to the fixed weights given to the virtual control points, a dynamic weighting function is established to alleviate the imbalance between the number of virtual control points. The specific expression is as follows:

[0126]

[0127] Where N is the number of iterations. The unknown ground coordinates of the connection points are calculated using the least squares method, thereby solving the refined geometric positioning parameters of each image.

[0128] In order to verify the effectiveness of this solution, this embodiment selects HiRISE (High-Resolution Imaging Science Experiment) and CTX (Context Camera) images of three areas as experimental data. Area 1 is near the Gale crater, area 2 is near the Jezero crater, and area 3 is near the Utopia Planitia. Eight images are selected for each area. The experiment uses the THEMIS mosaic with a resolution of 100 meters per pixel and the HRSC / MOLA Blended DEM with a resolution of 200 meters per pixel as plane and elevation reference data to obtain reliable control information.

[0129] The rational function model before and after the generalized control constraint adjustment is used to project the triangulated object space 3D coordinates of the matching tie points into the image space, and the residual between the back-projected image space coordinates and the actual image space coordinates in the row and column directions of the image is calculated. Table 1 summarizes the image space compensation and the image space back-projection root mean square error of the tie points before and after the adjustment. The results show that due to the influence of linear and nonlinear deviations on the initial RPCs, the maximum root mean square of the reprojection residual is about 13 pixels in the row direction and about 7 pixels in the column direction. After applying image space compensation, the geometric consistency of the image coordinates has improved, but it still does not meet the requirements of high-precision applications. After the generalized control constraint adjustment, the root mean square error is reduced to less than 0.5 pixels, which improves the geometric consistency of the image and effectively compensates for the linear and nonlinear deviations in the RPC.

[0130] Table 1. Root mean square error of back-projection adjustment with generalized control constraints (unit: pixel)

[0131]

[0132] The relative positioning root mean square error before and after the generalized control constraint adjustment is calculated, as shown in Table 2. The relative positioning error is the difference between the object space three-dimensional coordinates corresponding to the triangulation calculation activities of the connection points between multiple pairs of images. The plane root mean square error of the multi-source image before adjustment is about 80 meters, and some areas exceed 200 meters. The elevation root mean square error ranges from 15 meters to 110 meters. Through the generalized control constraint adjustment, the relative plane and elevation root mean square errors are reduced, among which the plane root mean square error is less than 7 meters, and the elevation root mean square error is less than 3 meters.

[0133] Table 2. Root mean square error of relative positioning before and after generalized control constraint adjustment (unit: meter)

[0134]

[0135] 20 checkpoints were manually measured from the reference DEM and DOM as the true values ​​of their object space 3D coordinates. The absolute positioning accuracy of the images before and after the adjustment was statistically analyzed by calculating the difference between the object space 3D coordinates obtained by triangulation of the checkpoints before and after the adjustment and the true values ​​of the object space 3D coordinates. Table 3 shows the root mean square error of absolute positioning before and after the generalized control constraint adjustment. The deviation range of the initial RPC in the X and Y directions relative to the reference DOM is more than 300 meters, while the deviation range in the Z direction relative to the reference terrain is between 50 and 240 meters. The generalized control constraint adjustment method reduces the deviation in the X and Y directions to within 23 meters and the deviation in the Z direction to within 5 meters by introducing the terrain change information of the prior terrain. However, it should be noted that the manually selected checkpoints themselves have certain errors.

[0136] Table 3. Root mean square error of absolute positioning before and after generalized control constraint adjustment (unit: meter)

[0137]

[0138] Using RPC optimized based on adjustment results, HRSC / MOLABlended DEM was used as the reference DEM, and the orthophoto map was generated using an equidistant projection coordinate system. Among them, CTX DOM was resampled with a uniform 6-meter resolution, while HiRISE DOM was resampled to a 1-meter resolution. Figure 5a to Figure 5i The figure shows the consistency comparison between DOM before and after adjustment and the image mosaic result of the experimental area. Registration deviation often occurs at high-frequency terrain features (such as crater edges and curved structures), resulting in pixel offsets between images. The pixel offset of the experimental area image before adjustment was on the order of hundreds of meters, but it has been corrected to the sub-pixel level after generalized control constraint adjustment, indicating that the method proposed in this scheme has a significant effect in solving image inconsistency.

[0139] In summary, this scheme aims to solve the problem of lack of measured control points on the surface of Mars and the widespread inconsistency between multi-source data. It proposes a joint adjustment method for Mars multi-source orbiter data under generalized control constraints. It uses a priori MOLADEM and low-resolution digital orthophoto map (DOM) as control benchmarks, and focuses on breakthroughs in key technologies such as Mars surface point feature extraction and image matching, image square compensation model of superimposed polynomials and cubic spline functions, adjustment model construction based on geometric constraints and adaptive weighting strategy, so as to improve the geometric positioning accuracy of multi-source orbiter images.

[0140] On the one hand, the large segmentation model is used to extract craters on the Martian surface, and the major-minor axis ratio and area ratio of the crater ellipse fitting are combined to optimize the point set registration method that retains the global and local structure, so as to achieve accurate matching of point features on the Martian surface; on the other hand, an image square deviation compensation model based on cubic splines and affine functions is constructed to effectively consider the linear and nonlinear distortion of the orbiter image; in addition, an adjustment model combining plane weak control constraints, elevation control constraints, image square compensation model, image connection points and virtual control points is constructed.

[0141] Therefore, by introducing the image space compensation model, using affine transformation and cubic spline function orbiter images to correct linear and nonlinear errors, and integrating connection points, virtual control points, and plane and elevation constraints, the geometric optimization of images is achieved: the SAM model is used to extract and segment craters, and the improved point set registration algorithm is used to match the center of the crater to generate reliable plane control points, thereby realizing plane weak constraint control; combining the slope information of the prior terrain, the elevation constraint control is established; and the adaptive weighting function is further introduced to dynamically coordinate various constraints to improve the relative and absolute geometric positioning accuracy of multi-source images. The method proposed in this scheme can improve the absolute and relative positioning accuracy of Mars multi-source orbiter images, and realize the accurate registration of multi-source images under the unified Mars reference benchmark.

Claims

1. A joint adjustment method for Mars multi-source orbiter images based on generalized control constraints, characterized in that: The following steps are involved: S1. Extract plane control points and match them to obtain connection points based on low-resolution reference DOM and Mars orbiter optical images; S2. Based on the fitting residual of the rational function model of the Mars orbiter image, an image-space compensation model of superimposed polynomial and cubic spline function is established and initial compensation parameters are obtained; S3. According to the plane control points, connection points and initial compensation parameters, combined with the plane and elevation control constraints, an adjustment model is constructed, and an adaptive weighted strategy is used to solve the adjustment model to obtain the geometric positioning parameters of the image.

2. A Mars multi-source orbiter image joint adjustment method based on generalized control constraints according to claim 1, characterized in that: The step S1 comprises the following steps: S11, preprocessing the optical image of the Mars orbiter to obtain a preprocessed optical image; S12, extracting connection points between images through least squares template matching processing for the pre-processed optical image; S13, combining the pre-processed optical image and the low-resolution reference DOM, outputting the crater center through image segmentation, mask segmentation, and crater fitting processing; The crater centers are then screened, and the matched crater centers are obtained through point set registration to obtain the plane control points.

3. A Mars multi-source orbiter image joint adjustment method based on generalized control constraints according to claim 2, characterized in that: In the step S12, for the optical images with different illumination, the directional phase consistency histogram HOPC is used as a similarity measure to perform least squares template matching; For optical images with consistent lighting conditions, the normalized cross correlation (NCC) is used as the similarity measure for least squares template matching.

4. A Mars multi-source orbiter image joint adjustment method based on generalized control constraints according to claim 2, characterized in that: The specific process of step S13 includes: Contrast-limited adaptive histogram equalization and gamma correction are used to perform radiometric enhancement on the preprocessed optical image and low-resolution reference DOM to strengthen the crater rim. The SAM model is used to extract and segment the enhanced image; In the process of crater extraction, the segmentation masks that do not meet the shape requirements are screened out, and the roundness and ellipticity of the mask are evaluated by calculating the area and circumference of the mask. When the ratio of the circumference to the circumference and the ratio of the major and minor axes exceed the threshold of 0.9, the mask is judged to be a circle; when the ratio of its major and minor axis area to the actual area exceeds the threshold of 0.8, the mask is judged to be an ellipse; the fitting accuracy of the circle or ellipse is evaluated by edge detection, and only the masks with high fitting accuracy and accurate center position are retained; Then, an improved point set registration method that retains global and local structures is used to realize the registration of crater positions between images. First, based on the prior geographic information of the craters, the target centers that lack candidate targets within the distance range are screened out. At the same time, local features are enhanced by constructing a feature descriptor that integrates the aspect ratio, fitting area ratio, and shape context features to ensure that the global and local structures are retained during the matching process to complete the registration.

5. The Mars multi-source orbiter image joint adjustment method based on generalized control constraints according to claim 1, characterized in that: The specific process of step S2 is: The elevation of the image area is divided into several elevation surfaces, and a three-dimensional virtual control grid of the image is established in the image space with a set interval as the virtual control point of the image space. On this basis, the image coordinates of the virtual control grid are projected onto each elevation surface according to the strict imaging model to obtain the virtual control point of the object space, and the rational function model parameters RPC are calculated by fitting, and the geometric relationship between the image space and the object space is expressed by the rational function model. The image space compensation model of linear and nonlinear errors of images is constructed by using superimposed polynomials and cubic spline functions to obtain initial compensation parameters for compensating nonlinear image errors.

6. The Mars multi-source orbiter image joint adjustment method based on generalized control constraints according to claim 5, characterized in that: The expression of the rational function model is: Among them, (X o ,Y o ) and (X s ,Y s ) represent the offset and scale value of the normalized pixel coordinates (X, Y), (P, L, H) are the latitude, longitude and elevation of the normalized object point, Num s , Num l ,Den s and Den l It is a third-order polynomial with a total of 20 coefficients, namely the rational function model parameter RPC.

7. The method for joint adjustment of Mars multi-source orbiter images based on generalized control constraints according to claim 6, characterized in that: The expression of the image space compensation model is: Among them, the parameters (e0, e1, e2, f0, f1, f2) represent the affine transformation used to compensate for the linear deviation. According to the time variation characteristics caused by posture jitter, the average residual of the image points in the same row is used as the optimal estimate of the RPC nonlinear deviation. The corrected pixel coordinates (X ′ ,Y ′ ) is calculated as follows: in, are the unknown parameters of the cubic spline function, which are used to compensate for nonlinear image errors, and i=1,2…m-1 is the number of connection points.

8. The Mars multi-source orbiter image joint adjustment method based on generalized control constraints according to claim 7, characterized in that: The process of constructing the adjustment model in step S3 includes: The error equation is established and the regional block adjustment is performed. According to the tie points, virtual control points, plane weak control constraints based on the center of the crater, terrain elevation constraints and rational function model parameters RPC obtained by multi-source image matching, the error equation of the joint adjustment model is constructed, which includes the following four observation equations: Among them, TPs are connection points, ECPs are elevation control points, VCPs are virtual control points, PCPs are plane control points, V TPs is the observation equation of the connection point, V ECPs is the observation equation of the elevation control constraint, V PCPs is the observation equation of the plane control constraint, V VCps is the observation equation of the virtual control point, t is the correction number of the adjustment parameter, x is the unknown ground coordinate of the connection point, matrices A and B are the error equation coefficient matrices, L is the constant term, P i,i＝1,2,3,4 For the corresponding weight, for the error equation, construct the normal equation, and use the least squares method to fit and obtain the object coordinates of the most approximate connection points; Then, terrain constraints based on the prior DEM are introduced, and elevation control is achieved by minimizing the sum of squares of the distances from the connection points to the DEM surface, which corresponds to the second observation equation V ECPs The height error equation V h The expression is as follows: Among them, h0 is the elevation of the connection point, (lat0, lon0) is the geographic plane coordinate, and DEM (lat0, lon0) is the elevation value obtained by interpolation at (lat0, lon0). and are the slopes of DEM in latitude and longitude respectively, (Δlat, Δlon, Δh) are the corrections of the ground coordinates of the connection points; A plane weak control constraint is constructed to introduce the terrain change information near the plane control point into the plane control constraint, and the object space 3D coordinates of the control point are adaptively adjusted to achieve automatic registration of the image with the center of the crater in the DOM. The difference between the pixel coordinates of the plane control point and the corresponding 3D coordinates in the reference map is minimized by the following formula: Where (X, Y) is the pixel coordinate of the center of the crater in the orbiter image, RPC X (lat,lom,DEM) and RPC Y (lat, lon, DEM) is the pixel coordinate of the center of the crater obtained by reprojecting the object space coordinates in the reference DOM and DEM onto the image, which is consistent with the third observation equation V PCPs The corresponding plane error equation v X and v Y The expression is as follows: in, as well as and The derivatives of the image coordinates with respect to latitude, longitude and elevation are obtained for reprojection. (X0, Y0) are the pixel coordinates of the plane control point. (ΔX′, ΔY′) are the differences between the image coordinates calculated by reprojecting the object coordinates of the plane control point and the corresponding actual image coordinates. Virtual control points (VCPs) are introduced to effectively reduce the degrees of freedom of the adjustment model. VCPs are strategically selected by predefined thresholds and grid sizes and introduced as weighted observations into the fourth observation equation V VCPs middle.

9. The Mars multi-source orbiter image joint adjustment method based on generalized control constraints according to claim 8, characterized in that: The process of using the adaptive weighting strategy to solve the adjustment model in step S3 includes: According to the prior matching accuracy σ TPs , define the weight P1 of the connection point as: Taking the prior terrain as control information, the error caused by the terrain and the point inconsistency as the "gross error" in the control information, in order to quantitatively reflect the impact of the terrain on the accuracy, the weight function of the elevation control constraint is defined as follows: Among them, img res is the image resolution, d h is the elevation difference, h mean is the average elevation difference of the data set, σ h is the standard deviation of the elevation difference. When the elevation difference is less than σ h When the height difference exceeds three times σ h When , the weight gradually approaches zero; Calculate the scale factor for each adjustment area separately and Where Num tie is the number of connection points, Num pcp is the number of plane control points, Num vcp is the number of virtual control points, considering the crater radius r and the fitting standard deviation f, the weight function of the plane control constraint is defined as follows: Among them, d p is the plane difference, p mean is the average value of the plane difference, σ p is the standard deviation of the plane difference, R max is the normalized radius factor, F max is the normalized fitting factor, the smoothing function σ f As f approaches the set threshold f threshold Gradually reduce the weight to zero, the steepness of the transition curve is controlled by the constant k, when the standard deviation is lower than f threshold When the value is higher than f, the weight decreases rapidly. threshold When , the weight decreases more slowly; A dynamic weighting function is established to alleviate the imbalance between the number of virtual control points, and then the unknown ground coordinates of the connection points are calculated according to the least squares method to solve the refined geometric positioning parameters of each image.

10. The Mars multi-source orbiter image joint adjustment method based on generalized control constraints according to claim 9, characterized in that: The expression of the dynamic weighting function is: Where N represents the number of iterations.

Images