Method for topographic correction of optical images of the surface of mars based on shadow typing

Abstract

The present application relates to the technical field of resource exploration, and particularly relates to a Mars surface optical image terrain correction method based on shadow classification, obtaining a Mars optical image of a Mars surface and a digital elevation model matched with the Mars optical image; determining a light factor of sunlight irradiating the Mars surface according to the digital elevation model; dividing regions of the Mars surface into a shadow region, a penumbra region, a fall shadow region and a sunny slope region according to the light factor; performing image segmentation on the Mars optical image which is not complete in atmospheric parameters, and performing semi-empirical link correction, conservative processing or regression compensation on the four regions respectively according to the segmentation result; for the Mars optical image which is complete in atmospheric parameters, determining a slope scattering component and a slope direct component of the Mars surface, performing image segmentation on the Mars optical image, and performing correction on the four regions in each label sub-image according to the segmentation result and in combination with the slope scattering component and the slope direct component.

Description

Technical Field

[0001] This invention belongs to the field of resource exploration technology, and in particular relates to a method for terrain correction of optical images of the Martian surface based on shadow classification. Background Technology

[0002] As one of the core directions of deep space exploration, the accuracy of Mars surface topography data directly determines the success or failure of tasks such as Martian geological structure analysis, mineral resource exploration, and landing site selection. Existing topography correction methods can be divided into semi-empirical model-based and physical model-based methods according to their principles. Both types of methods have their applications in Mars remote sensing image processing, but neither has completely solved the problem of adapting to complex terrain.

[0003] Semi-empirical model-based methods use the illumination coefficient as the core variable and estimate a small number of empirical parameters through statistical regression to normalize and correct the radiance or reflectance of images. The core objective is to eliminate the influence of illumination differences caused by terrain undulations. Typical examples include C-correction, SCS / SCS+C models, and Minnaert models. The core advantages of this type of method are low computational cost, simple engineering implementation, and no reliance on complex auxiliary parameters. However, it has three major limitations in Martian scenes: First, it often uses uniform parameters across the entire map or parameters grouped by large regions, making it difficult to adapt to areas with significant differences in Martian terrain spatial scales, such as canyons and surrounding plains, dunes and impact craters. This can easily lead to local undercorrection at abrupt terrain changes, such as lower brightness on shaded slopes, or overcorrection, such as amplified noise in low-illuminance areas. Second, under low-illuminance conditions, the illumination coefficient IC approaches 0, such as at the bottom of impact craters and on the shaded slopes of canyons. The correction formula is prone to numerical instability, leading to abnormal jumps in reflectance. Third, empirical parameters lack physical interpretability, verifiability, and traceability, making it difficult to reproduce the correction results from different research teams.

[0004] Physical model-based methods are based on radiative transfer theory and emphasize radiative consistency. Through radiative calibration, atmospheric correction and separation of scattering components, they incorporate components such as direct solar radiation, sky scattering radiation and radiation reflected by adjacent terrain into the framework for slope irradiance conversion or reflectance inversion.

[0005] Typical examples include the 6S radiative transfer model coupled with the Perez scattering model and the MODTRAN atmospheric correction combined with the topographic view factor model. The core advantages of this type of method are its clear physical meaning and strong robustness when parameters are complete. However, it faces significant bottlenecks in Mars engineering applications: First, it relies on a large number of auxiliary parameters, such as sensor calibration coefficients, atmospheric scattering parameters, and average surface reflectance. In Mars exploration, parameters are often unavailable, such as atmospheric profiles of unknown regions, forcing model simplification and leading to decreased correction accuracy in complex scenes. Second, it has high computational costs, involving complex steps such as view factor calculation and multi-component radiometric conversion, making it difficult to meet the demands of batch processing massive amounts of Mars imagery. Third, it requires high accuracy of topographic data; when the resolution of the digital elevation model is lower than the image resolution, insufficient topographic detail matching can easily introduce new errors. Summary of the Invention

[0006] In view of this, the present invention aims to provide a terrain correction method for optical images of the Martian surface based on shadow classification. Through a systematic design of "classification-regionalization-linkage", it takes into account stability, information recovery capability, physical consistency and engineering usability in the rugged terrain and complex shadow scenes of Mars, and is superior to the existing single-path terrain correction scheme.

[0007] To achieve the above objectives, the technical solution created by this invention is implemented as follows: A terrain correction method for optical images of the Martian surface based on shadow classification includes: S1: Acquire optical images of the Martian surface and a digital elevation model matching the Martian optical images; determine the illumination coefficient of sunlight illuminating the Martian surface based on the digital elevation model; S2: Based on the illumination coefficient obtained in step S1, the surface of Mars is divided into the umbra, gray shadow, shadow area, and sunny slope area. S3: Determine whether the atmospheric parameters of the Mars optical image in step S1 are complete. If complete, proceed to step S4; otherwise, proceed to step S5. S4: Perform image segmentation on the Mars optical image, and correct the sunny slope area and gray shadow area in each labeled sub-image according to the segmentation results. At the same time, use machine learning methods to correct the shadow area. S5: Determine the slope scattering component and slope direct component of the Martian surface when sunlight shines on it; perform image segmentation on the Martian optical image, and based on the segmentation results, correct the four regions in each labeled sub-image by combining the slope scattering component and the slope direct component. S6: Stitch together the four regions corrected in step S4 or S5 to obtain and output the corrected Mars optical image.

[0008] Furthermore, the process of determining the illumination coefficient of the Martian surface based on the digital elevation model in step S1 includes: Calculate the slope and aspect of the Martian surface for each pixel in the digital elevation model; The preliminary illumination coefficient is calculated using the following formula: IC0=cosσ-cosZ+sinσ×sinZ×cos(β-ω); Where IC0 represents the initial illumination coefficient, σ represents the slope, β represents the aspect, ω represents the azimuth angle when sunlight shines on the surface of Mars, and Z represents the zenith angle when sunlight shines on the surface of Mars. Set an illuminance threshold for low-illuminance pixels in the digital elevation model whose initial illuminance coefficient is less than the lower threshold. Then, apply the following formula to robustly constrain these low-illuminance pixels to obtain their final illuminance coefficients: IC=clip((IC0-T min ) / (T max -T min ),0,1)×IC0; Where IC represents the final illumination coefficient, clip represents the cutoff function, and T min T represents the threshold value under illuminance. max This represents the threshold value for illuminance.

[0009] Further, step S2 includes: the region with an illumination coefficient less than or equal to 0 is the first region; the region with an illumination coefficient greater than 0 but whose brightness in the Mars optical image is significantly lower than the brightness or local average of similar terrain is the second region; both the first and second regions are shaded slope regions, and the remaining regions are sunny slope regions; the set of pixels in the shaded slope region that completely blocks sunlight when it shines directly on the Martian surface is the umbra region; the set of pixels in the shaded slope region with an illumination coefficient greater than 0 but less than or equal to a set intensity threshold is the gray shadow region; ray tracing is performed along the direction of solar incidence in the shaded slope region to obtain the region blocked by adjacent terrain, which is the shadow region.

[0010] Furthermore, step S4 includes: S41: Perform image segmentation on the Mars optical image in step S1 to obtain and output terrain label images; in each label sub-image, use the sunny slope area and gray shadow area obtained in step S2 to establish the regression relationship between each band and the illumination coefficient in the Mars optical image; S42: Correct the gray shadow area and the sunny slope area using the regression relationship obtained in step S41; retain the umbra area obtained in step S2; and correct the umbra area using machine learning based on the sunny slope area adjacent to the shadow area in step S2.

[0011] Furthermore, in step S42: The following rotational correction model is used to correct the gray shadow area and the sunny slope area: I corr (λ)=I(λ)-a k (λ)×(IC-cosZ); Among them, I corr (λ) represents the corrected gray shadow region or sunny slope region, λ represents the band, I(λ) represents the corresponding original gray shadow region or sunny slope region, a k (λ) represents the slope corresponding to the k-th label sub-image in the regression relationship, IC represents the illumination coefficient, and Z represents the zenith angle when sunlight illuminates the surface of Mars.

[0012] Furthermore, the process of correcting the shadow area using machine learning methods in step S42 includes: The SEVI paradigm is used to construct the shadow sensitivity index and obtain the texture information of the shadow area; In the same sub-image of the label, the sunny slope area adjacent to the shadow area is randomly sampled to obtain the sampled image; A random forest regressor was trained using sampled images, shadow sensitivity index, and texture information to obtain a shadow prediction model; The image of the shadow region is input into the shadow prediction model. The obtained predicted image replaces the corresponding image in the shadow region, or the obtained predicted image is added to the corresponding image to obtain the corrected shadow region.

[0013] Furthermore, step S5 includes: S51: Convert the Mars optical image obtained in step S1 into the surface radiance and horizontal apparent reflectance of the Martian surface when sunlight shines on it; determine the cosine of the relative incident angle when sunlight shines on the Martian surface based on the digital elevation model obtained in step S1. S52: Determine the horizontal direct irradiance of sunlight on the surface of Mars based on the solar incident irradiance, and then determine the slope direct irradiance component based on the horizontal direct irradiance and the cosine of the relative incident angle obtained in step S51. S53: Based on the slope and aspect of the Martian surface and the azimuth of the sun, which are reflected in each pixel of the digital elevation model obtained in step S1, the anisotropic model is used to estimate the scattering components of the slope surface. S54: Perform image segmentation on the Mars optical image to obtain a terrain-tagged image; for each pixel in the Mars optical image, calculate the average reflectance and visible sky factor of the pixel neighborhood or the same tag sub-image for each pixel; S55: Using the average reflectance and visible sky factor obtained in step S54, the total irradiance of the slope in the four regions is inverted based on the direct irradiance component obtained in step S52 and the scattering component obtained in step S53. Then, based on the total irradiance of the slope in the four regions and the surface radiance, the slope reflectance of the four regions is calculated.

[0014] Furthermore, step S51 includes: The Mars optical image is converted into an equivalent DN, and then the equivalent DN is converted into top atmospheric radiance according to the sensor gain and bias or radiometric calibration coefficient. Using a radiative transfer model, the top atmospheric radiance obtained in step S51 is converted into surface outgoing radiance and horizontal apparent reflectance. The slope and aspect of the Martian surface are calculated for each pixel in the digital elevation model, and the cosine of the relative incident angle of sunlight on the Martian surface is determined by the following formula: cos(i)=cosσ-cosZ+sinσ×sinZ×cos(β-ω); Where i represents the relative angle of incidence, σ represents the slope, β represents the aspect, ω represents the azimuth angle when sunlight shines on the surface of Mars, and Z represents the zenith angle when sunlight shines on the surface of Mars.

[0015] Furthermore, in step S52: The direct irradiance on the horizontal plane can be obtained using the following formula: E dir,h =E0×cos(Z) Among them, E dir,h E0 represents the direct irradiance on a horizontal surface, while E0 represents the incident irradiance of the sun. The direct sunlight component on the slope is obtained by the following formula: E dir,s =E dir,h ×max(cos(i),0) / cos(Z); Among them, E dir,s This represents the direct sunlight component on the slope.

[0016] Furthermore, in step S55: In calculating the slope reflectance of the gray shadow area and the sunny slope area: The total irradiance of the slope in the gray shadow area or the sunny slope area is calculated by the following inversion formula: E tot_1 =k adj ×F v_1 ×ρ nei_1 +E dir,s_1 +E dif,s_1 ; Among them, E tot_1 E dir,s_1and E dif,s_1 ρ represents the total irradiance, direct irradiance, and diffuse irradiance of the slope in the gray shadow region or the sunny slope region, respectively. nei_1 and F v_1 k represents the average reflectance of the gray shadow area and the visible sky factor of the sunny slope area, respectively. adj Indicates the inversion coefficients; The slope reflectance of the gray shadow area and the sunny slope area is calculated using the following formula: ρ s_1 (λ)=π·L surf_1 (λ) / E tot_1 (λ); Where, ρ s_1 (λ) represents the slope reflectance of the gray shadow area or the sunny slope area, λ represents the band, and L surf_1 (λ) represents the surface radiance of the gray shadow area or the sunny slope area; In calculating the slope reflectance of the umbra and shadow regions: The total irradiance of the slope in the gray shadow area or the sunny slope area is calculated by the following inversion formula: E tot_2 =k adj ×F v_2 ×ρ nei_2 +E dif,s_2 ; Among them, E tot_2 and E dif,s_2 ρ represents the total irradiance and scattering component of the slope in the umbra or shadow region, respectively. nei_2 and F v_2 These represent the average reflectance of the umbra or the shadow region and the visible sky factor, respectively. The slope reflectance of the umbra and shadow regions is calculated using the following formula: ρ s_2 (λ)=π·L surf_2 (λ) / E tot_2 (λ); Where, ρ s_2 (λ) represents the slope reflectance of the umbra or shadow region, L surf_2 (λ) represents the surface radiance of the umbra or shadow region.

[0017] Compared with the prior art, the present invention can achieve the following beneficial effects: (1) The Mars surface optical image topography correction method based on shadow classification created by the present invention further subdivides gray shadow, umbra and shadow from the traditional shadow concept and adopts a differentiated processing strategy, which can avoid parameter deviation caused by misjudging gray shadow as sunny slope or umbra, and at the same time reduce the residual accumulation of shadow boundary and transition zone; the present invention establishes a hybrid correction framework that can switch and combine physical link (step S5) and semi-empirical link (step S4), which has strong physical interpretability when the parameters are complete, and still has good engineering feasibility and batch processing capability when the parameters are incomplete; the Mars surface optical image processed by the present invention can more reliably serve tasks such as mineral identification, landform classification, change detection and landing area engineering interpretation, and reduce the interference of radiation error introduced by topographic effect on subsequent scientific analysis and engineering decision-making; (2) The Mars surface optical image topography correction method based on shadow classification described in this invention, in the semi-empirical link (step S4): by performing parameter adaptive estimation within the topographic label image, the correction parameters can be dynamically adjusted according to changes in surface type, topographic relief and local illumination conditions, significantly reducing the local undercorrection or overcorrection problems caused by uniform parameters across the entire image in the prior art; by implementing constraints such as threshold truncation, weight attenuation and unreliable labeling on low-illuminance pixels, the abnormal reflectance jump, noise amplification and result distortion in low-illuminance areas are effectively suppressed, improving the numerical stability of complex shady slopes, pit bottoms and canyon shaded areas; by adopting conservative strategies such as maintaining the original value, neighborhood uniform smoothing or scatter-dominated small correction for the umbra area, the noise can be mistakenly compensated as false information under extremely low illumination conditions, thereby improving the credibility of the shadow deep area processing results; by establishing a regression prediction compensation link based on shadow sensitivity index, illumination coefficient, local texture features and unit labels, the visible light band reflectance of the shadow area can be effectively restored, making the spectral morphology of the shadow area more consistent with the adjacent similar sunny slope features; (3) The Mars surface optical image terrain correction method based on shadow classification described in this invention comprehensively considers direct radiation, sky scattered radiation and adjacent terrain reflected radiation in the physical link (step S5), which can improve the radiation uniformity of shadow areas, slope break zones and steep slope areas, and reduce the systematic deviation caused by the approximation of a single irradiance component. Attached Figure Description

[0018] The accompanying drawings, which form part of this invention, are used to provide a further understanding of the invention. The illustrative embodiments and descriptions of the invention are used to explain the invention and do not constitute an undue limitation of the invention. In the drawings: Figure 1 A schematic flowchart of the Mars surface optical image terrain correction method based on shadow classification as described in the embodiments of the present invention; Figure 2The flowchart illustrates the method for terrain correction of Mars surface optical images based on shadow classification, as described in an embodiment of the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to the accompanying drawings and specific embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and do not constitute a limitation thereof. It should be noted that, unless otherwise specified, the embodiments and features described herein can be combined with each other.

[0020] The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0021] like Figure 1 and Figure 2 As shown in the embodiment of the present invention, the method for terrain correction of Mars surface optical images based on shadow classification includes: S1: Acquire an optical image of the Martian surface and a digital elevation model (DEM) matching the Martian optical image; determine the illumination coefficient of sunlight illuminating the Martian surface based on the DEM. In this embodiment of the invention, the Martian optical image can be a visible light image or a visible-near-infrared image. Furthermore, step S1 also includes reprojecting the DEM onto the coordinate system of the Martian optical image and resampling it to the resolution of the Martian optical image, and, if necessary, performing slight smoothing on the DEM to suppress noise.

[0022] In some embodiments, the process of determining the illumination coefficient of the Martian surface based on a digital elevation model includes: Calculate the slope and aspect of the Martian surface for each pixel in the digital elevation model; Let the cosine of the relative incident angle of sunlight illuminating the Martian surface be the illumination coefficient, and then calculate the preliminary illumination coefficient using the following formula: IC0=cosi=cosσ-cosZ+sinσ×sinZ×cos(β-ω); Where IC represents the illumination coefficient, σ represents the slope, β represents the aspect, ω represents the azimuth angle when sunlight illuminates the Martian surface, Z represents the zenith angle when sunlight illuminates the Martian surface, and i represents the relative angle of incidence. In this embodiment of the invention, the preliminary illumination coefficient IC0 is further normalized using the following formula to improve consistency across time phases: IC0'=IC0= / cos(Z); Where IC0' represents the normalized initial illumination coefficient; Set an illuminance threshold for low-illuminance pixels in the digital elevation model whose initial illuminance coefficient is less than the lower threshold. Then, apply the following formula to robustly constrain these low-illuminance pixels to obtain their final illuminance coefficients: IC=clip((IC0-T min ) / (T max -T min ),0,1)×IC0; Where IC represents the final illumination coefficient, clip represents the cutoff function, and T min T represents the threshold value under illuminance, preferably taking a value in the range [0.02, 0.08]. max This represents the threshold value for illuminance, preferably within the range of [0.05, 0.15]. The above formula can be understood as weighting and attenuating the initial illuminance IC0 to obtain the final illuminance IC, where the weights are clip((IC0-T)). min ) / (T max -T min ),0,1).

[0023] Existing technologies typically do not set explicit numerical robust constraints for low-light pixels, which can easily lead to the denominator approaching zero, abnormal jumps in reflectivity, and noise amplification when the cosine of the relative incident angle (cosi) is close to zero. This invention introduces mechanisms such as weight attenuation and uncorrected marking to effectively improve the stability of complex shaded areas such as shady slopes, pit bottoms, and steep slopes.

[0024] In some other embodiments, the initial illumination coefficient IC0 can be directly truncated to obtain the final illumination coefficient IC, as shown in the following formula: IC=max(IC0,T min ).

[0025] Understandably, based on the normalized initial illumination coefficient IC0', robust constraints are then applied to the low-light regions in the Mars optical image to obtain the final illumination coefficient.

[0026] In this embodiment of the invention, such low-light pixels are marked as unreliable low-light.

[0027] S2: Based on the illumination coefficient obtained in step S1, the surface of Mars is divided into the umbra, gray shadow, shadow area, and sunny slope area.

[0028] In some embodiments, step S2 includes: a region with an illuminance coefficient less than or equal to 0 is a first region; a region with an illuminance coefficient greater than 0 but whose brightness in the Mars optical image is significantly lower than the brightness or local average of similar terrain is a second region; both the first and second regions are shaded slope regions, and the remaining regions are sunny slope regions; the set of pixels in the shaded slope region that completely blocks sunlight when it directly hits the Martian surface is the umbra region; the set of pixels in the shaded slope region with an illuminance coefficient greater than 0 but less than or equal to a set intensity threshold is the gray shadow region, which is a transitional area where the brightness or texture is close to that of a shadow; ray tracing is performed along the direction of solar incidence in the shaded slope region to obtain the area blocked by adjacent terrain, which is the shadow region, which is the area blocked by adjacent terrain. In this embodiment of the invention: the intensity threshold is the illuminance upper threshold T in step S1. max The umbra region can be obtained using slope geometry discrimination or ray tracing methods. Specifically, this embodiment of the invention adopts the slope geometry discrimination method disclosed in "Luo Qingzhou, Liu Shunxi, Zeng Qihong, et al. Research on the Judgment of Umbra and Shadow Based on DEM [J]. Land Resources Remote Sensing, 2009, (02): 29-31." The main process is to calculate the slope, aspect, and slope normal vector of each grid cell by the digital elevation model, calculate the angle θ between the sun's incident direction and the slope normal vector, and if the cosine value of the angle is less than 0, the slope faces away from the sun and is determined to be the umbra region; the ray tracing method disclosed in "Digital Elevation Model Tutorial" is used to emit rays from the sun position to the ground pixel by pixel, and determine whether it is blocked by higher elevation terrain along the ray path. If it is not blocked and the slope is backlit, it is marked as the umbra.

[0029] The shadow area can also be obtained through morphological expansion and digital elevation model constraints. In this embodiment of the invention, the shadow area is obtained by means of the method described in “Chi Yufeng, Lai Riwen, Yan Qi, et al. Mountain shadow information detection and extraction based on LandSat8 OLI data [J]. Journal of Mountain Science, 2017, 35(04):580-589.DOI:10.16089 / j.cnki.1008-2786.000256.” The general process is as follows: the initial shadow area is obtained by using an initial binary shadow map, morphological expansion (expansion) is performed on it to connect the breaks and fill the holes, and then the digital elevation model is used to constrain only the area with lower elevation (below the projection source) and downstream of the illumination direction to obtain the final shadow area.

[0030] It should be noted that existing technologies typically simplify shadows into a binary division of shadow / non-shadow or simply use a sunny / shaded slope division, which fails to reflect the differences between various shadow formations and radiation states in the complex terrain of Mars. This invention, however, further divides pixels into four categories: sunny slope, gray shadow, umbra, and shadow, and employs differentiated processing and evaluation strategies for each type, thereby reducing errors in shadow boundaries and transition zones.

[0031] S3: Determine whether the atmospheric parameters of the Mars optical image in step S1 are complete. If complete, proceed to step S4; otherwise, proceed to step S5.

[0032] In this embodiment of the invention, the completeness determination in step S3 requires that all of the following three parameters be met simultaneously; failure to meet any one of them indicates that the atmospheric parameters are incomplete: (1) Parameter type completeness judgment (missing items are directly judged as incomplete) To meet the radiative transfer calculation requirements for Martian surface topography correction, it is essential to fully cover four core atmospheric parameters: Martian atmospheric aerosol optical thickness, atmospheric principal component gas column content parameters, atmospheric vertical profile parameters, and atmospheric scattering characteristics parameters. The absence of any one of these parameters will result in an incomplete parameter set. The specific requirements and correction necessity for each parameter type are as follows: For Martian atmospheric aerosol optical depth (AOT): the criteria are that it must contain AOT data that corresponds one-to-one with the imaging band of the image to be corrected, covering the entire imaging band from visible light to near-infrared to short-wave infrared; the necessity of correction: it is used to calculate the transmittance of Martian atmospheric path, correct the attenuation of direct solar radiation in the thin Martian atmosphere, and is the core input for terrain correction radiative transfer calculation, which directly determines the radiative correction accuracy of sunny slopes and shadowed areas. For the atmospheric main component gas column content parameter: the determination requirement is that it must include the concentration of CO2 column, the main component of the Martian atmosphere, and the water vapor column content data of the corresponding imaging band; the necessity of correction: it is used to calculate the radiation attenuation caused by atmospheric gas absorption, correct the difference in incident radiation on slopes with different slopes and aspects, and adapt to the correction needs of complex terrains with large elevation differences such as high-altitude impact craters and mountains on Mars. For atmospheric vertical profile parameters, the criteria are that they must include atmospheric temperature profiles and atmospheric pressure profiles of the area corresponding to the time of image transit, and the number of vertical layers must be greater than or equal to ten. The necessity of correction is to accurately calculate the radiation transfer on the atmospheric vertical path, eliminate the differences in atmospheric attenuation at different altitudes, and solve the correction distortion problem in areas with large elevation differences. For atmospheric scattering characteristic parameters, the determination requirements are that they must include the scattering phase function and single scattering albedo data corresponding to the Martian aerosol type (dust or ice crystals); the necessity of correction is to accurately calculate the atmospheric diffuse radiation component and correct the incident radiation in the umbra, shadow, and gray shadow regions. (2) Parameter value validity determination (invalid parameters are directly determined as incomplete) Assuming the above four essential atmospheric parameters have been covered, the following validity requirements must be met simultaneously; failure to meet any one of these requirements will result in incomplete parameters: Physical range compliance: All parameter values ​​must conform to the objectively reasonable physical range of the Martian atmosphere. Specific thresholds: AOT in the visible-near-infrared band must be in the range of 0.01~2.0 (this can be relaxed to 0.01~4.0 during global dust storms); CO 2 The column concentration needs to be between 5 and 7 mb. Parameters within a km range are invalid; parameters that exceed a reasonable physical range are invalid. Data accuracy compliance: Parameter accuracy must meet the minimum accuracy requirements for Mars terrain correction. Specific judgment thresholds: The absolute error of AOT is less than or equal to 0.05; the temperature measurement error of the atmospheric profile is less than or equal to 2K, and the pressure measurement error is less than or equal to 0.5 mbar; parameters that do not meet the accuracy requirements are invalid. Source traceability: Parameters must originate from authoritative and compliant data sources used in Mars exploration, with a complete traceability path. Compliant data sources include: data from the Mars Reconnaissance Orbiter (MRO) MCS / MARCI instruments, data from the Mars Express SPICAM spectrometer, Mars atmospheric measurement data from the Tianwen-1 orbiter, and officially released data from the Mars Global Reference Atmosphere Model (MGCM); parameters without traceability or from non-authoritative sources are invalid. (3) Spatiotemporal-band matching determination (mismatch is directly determined as incomplete) Provided that the parameter type is complete and the value is valid, the following requirements must also be met. Failure to meet any one of these requirements will result in the parameter being deemed incomplete: Time matching: The difference between the observation time of the parameter and the transit time of the Mars optical image to be corrected is less than or equal to two Martian hours; parameters that are outside the time window cannot represent the true atmospheric state at the time of imaging and are judged as mismatched, i.e. invalid. Spatial matching: The spatial coverage of the parameter must 100% contain the imaging area of ​​the image to be corrected, and the spatial resolution must be compatible with the image resolution; for example, for HiRISE sub-meter high-resolution imagery, the spatial resolution of the corresponding atmospheric parameter must be less than or equal to 1km; parameters with incomplete spatial coverage or mismatched resolution are judged as mismatched. Band matching: The spectral response bands of all atmospheric parameters must correspond one-to-one with the imaging bands of the image to be corrected; for example, if the image contains four imaging bands—blue, green, red, and near-infrared—it must have corresponding AOT parameters for the four bands; parameters that do not correspond to bands are judged as mismatched. The final completeness judgment rule is that atmospheric parameters can be judged as complete only if all of the following conditions are met simultaneously: complete coverage of the four mandatory core atmospheric parameters, with no missing items; all parameters meet the validity requirements of physical range, accuracy, and source; and all parameters meet the matching requirements of time, space, and band. If any one condition is not met, the atmospheric parameters are judged as incomplete.

[0033] Existing technologies typically employ only a single semi-empirical model or a single physical model. The former offers weak physical interpretation, while the latter relies heavily on auxiliary parameters, making it difficult to balance stability and engineering feasibility. This invention constructs a hybrid framework that allows for switching and combination of semi-empirical links (the method in step S4) and physical links (the method in step S5). When parameters are complete, physical consistency processing is performed; when parameters are insufficient, semi-empirical correction or residual correction is applied, thereby improving the overall applicability.

[0034] S4: Perform image segmentation on the Mars optical image, and correct the sunny slope area and gray shadow area in each labeled sub-image based on the segmentation results. At the same time, use machine learning methods to correct the shadow area.

[0035] In some embodiments, step S4 includes: S41: Perform image segmentation on the Mars optical image in step S1 to obtain and output terrain label images; in each label sub-image, use the sunny slope area and gray shadow area obtained in step S2 to establish the regression relationship between each band and the illumination coefficient in the Mars optical image.

[0036] In this embodiment of the invention, the Mars optical image or the principal component feature map after principal component analysis of the Mars optical image is segmented according to the terrain using Mean shift to obtain a terrain label image. It can be understood that each label sub-image in the terrain label image is a segmentation unit. During image segmentation, the terrain in each label sub-image should be as homogeneous as possible. If the number of valid samples in a label sub-image is insufficient, it can be merged with adjacent label sub-images of the same type. For example, if the number of valid samples in a label sub-image is less than a preset minimum sample value (the minimum sample value is adjusted according to the resolution of the Mars optical image; in this embodiment, the minimum sample value is set to be in the range [50, 200]), a strategy of merging segmentation units or borrowing parameters from neighboring regions is used to merge the label sub-images. In other embodiments, Mars terrain semantic segmentation or a lightweight classification model can also be used to obtain the terrain label image. Furthermore, in this embodiment of the invention, robust regression or the Theil-Sen method is used for iterative calculation to obtain the regression relationship between each band and the illumination coefficient, and outliers are removed during the iteration process.

[0037] S42: The gray shadow region and the sunny slope region are corrected using the regression relationship obtained in step S41; the umbra region obtained in step S2 is retained; based on the sunny slope region adjacent to the shadow region in step S2, a machine learning method is used to correct the shadow region. In this embodiment of the invention, only neighborhood uniformity smoothing is performed on the retained umbra region to avoid amplifying noise into false information. Existing technologies often treat the umbra region as an uncorrectable region or directly use uniform stretching and brightness compensation methods, which easily introduces noise amplification and false information. This invention adopts a conservative strategy for the umbra region, including maintaining the original value, neighborhood uniformity smoothing or small correction dominated by scattering, and outputting a quality label to improve the reliability of the results for deep shadow regions.

[0038] In some embodiments, the following rotational correction model is used to correct the gray shadow area and the sunny slope area: I corr (λ)=I(λ)-a k (λ)×(IC-cosZ); Among them, I corr (λ) represents the corrected gray shadow region or sunny slope region, λ represents the band, I(λ) represents the corresponding original gray shadow region or sunny slope region, a k (λ) represents the slope corresponding to the k-th label sub-image in the regression relationship.

[0039] In other embodiments, the existing SCS+C model can be used to correct the gray shadow region and the sunny slope region based on the slope and intercept in the regression relationship. Specifically, for each labeled sub-image of the gray shadow region and the sunny slope region, let the slope α in the regression relationship of that labeled sub-image be... k With intercept b k Dividing by the other yields the adjustment parameter C, i.e., C = b / a. For the slope a... k For labeled sub-images that are close to 0 or have insignificant regression, C is not used for those sub-images. The remaining labeled sub-images are corrected according to the SCS+C formula. Once all labeled sub-images have been corrected, the corrected gray shadow area and sunny slope area are obtained.

[0040] Insignificant regression refers to a univariate linear regression equation established for the band values ​​of Mars optical images and the illumination coefficient (cosine of the solar incidence angle cos(i)) within the tag image, which, after statistical significance testing, cannot prove a true linear correlation between the two at a 95% confidence level. In the embodiments of this invention, the method for determining insignificant regression includes: A regression is considered insignificant if any of the following conditions are met: a. For the regression correlation coefficient (i.e., the slope a) k With intercept b k A two-tailed t-test was performed, and the p-value obtained was >0.05; b. The coefficient of determination R of the regression relationship 2 <0.1; c. The effective sample size for regression is <30.

[0041] For labeled sub-images with slopes close to 0 or insignificant regression, the SCS+C model is not used for correction to avoid correction distortion.

[0042] In some embodiments, the process of correcting the shadow area using machine learning includes: constructing a shadow sensitivity index using the SEVI paradigm and obtaining texture information of the shadow area, where the texture information can be the local mean, variance, or gradient of the shadow area; randomly sampling the sunny slope area adjacent to the shadow area in the same labeled sub-image to obtain sampled images; training a random forest regressor using the sampled images, shadow sensitivity index, and texture information to obtain a shadow prediction model; inputting the image of the shadow area into the shadow prediction model, replacing the corresponding image in the shadow area with the obtained predicted image, or adding the obtained predicted image to the corresponding image to obtain the corrected shadow area. In this embodiment, cross-validation or Bayesian optimization is used to adjust the parameters of the random forest regressor during training. In the prior art, conventional terrain correction of shadow areas has poor effects, often resulting in correction failure, simple brightness stretching, or spectral distortion. This invention establishes a regression prediction compensation link based on shadow sensitivity index, illumination coefficient, local texture, and similar sunny slope training samples to restore the visible light reflectance of the shadow area, thereby enhancing the usability of shadow area information.

[0043] It should be noted that existing technologies mostly use uniform parameters across the entire map or parameters grouped by large regions for regression estimation, which is difficult to adapt to the highly heterogeneous terrain and surface type distribution on the Martian surface. This invention adaptively estimates parameters by performing regression relationships within segmented units or terrain category units, enabling the parameters to respond to local terrain and feature differences, thereby reducing undercorrection and overcorrection.

[0044] S5: Determine the slope scattering component and slope direct radiation component of the Martian surface when sunlight shines on it; perform image segmentation on the Martian optical image, and based on the segmentation results, correct four regions in each labeled sub-image by combining the slope scattering component and the slope direct radiation component.

[0045] In some embodiments, step S5 includes: S51: Convert the Mars optical image obtained in step S1 into the surface radiance and horizontal apparent reflectance of the Martian surface when sunlight shines on it; determine the cosine of the relative incident angle when sunlight shines on the Martian surface based on the digital elevation model obtained in step S1.

[0046] In some embodiments, step S51 includes: Radiometric calibration: Convert the Mars optical image into an equivalent DN, and then convert the equivalent DN into top atmospheric radiance according to the sensor gain and bias or radiometric calibration coefficient; Atmospheric uniformity: Using a radiative transfer model, the obtained top atmospheric radiance is converted into surface outgoing radiance and horizontal apparent reflectance. In other embodiments, the top atmospheric radiance can also be converted into surface outgoing radiance or horizontal apparent reflectance by looking up tables, such as path radiation, up and down transmittance, sky scattering irradiance, etc. The calculation of the slope and aspect of each pixel in the digital elevation model reflecting the Martian surface, and the calculation method for the cosine of the relative incident angle, are described in step S1 and will not be repeated here. In this embodiment of the invention, low-light pixels are determined according to the method in step S1, and such low-light pixels are marked as unreliable low-light pixels.

[0047] In this embodiment of the invention, the work is adapted from "Fernando, J., F. Schmidt, X. Ceamanos, P. Pinet, S. Douté, and Y. Daydou (2013), Surface reflectance of Mars observed by CRISM / MRO: 2. Estimation of surface photometric properties in Gusev Crater and Meridiani Planum, J. Geophys. Res. Planets, 118, 514–539." The radiative transfer model of doi:10.1029 / 2012JE004194. converts the top atmospheric radiance into surface outgoing radiance and horizontal apparent reflectance. The process includes: constructing a radiative transfer equation based on the imaging geometry of the image to be corrected, radiometric calibration coefficients, Martian atmospheric state parameters at the corresponding imaging time, and topographic geometric parameters calculated by the digital elevation model, obtaining four core atmospheric interference parameters: path radiation, up- and down-going atmospheric transmittance, and sky scattering irradiance; removing atmospheric interference from the obtained top atmospheric radiance, inverting to obtain surface outgoing radiance without atmospheric influence, and then normalizing it to convert it into horizontal apparent reflectance; verifying the inversion accuracy using official reference data of flat and uniform areas of Mars, and outputting parameters that meet the requirements of subsequent topographic correction.

[0048] Among them, path radiance, uplink and downlink transmittance, and sky-scattered irradiance are the core intermediate parameters of the Martian atmospheric radiative transfer equation. Path radiance is an atmospheric background interference term in the satellite observation signal that is unrelated to Martian surface reflection, and it is the deduction amount used to invert the surface outgoing irradiance; uplink atmospheric transmittance is the atmospheric attenuation correction coefficient during the upward transmission of surface reflected radiation, and it directly determines the inversion accuracy of surface outgoing irradiance; downlink atmospheric transmittance is the attenuation coefficient during the downward transmission of direct solar radiation, and it determines the total amount of direct radiation received by the surface; sky-scattered irradiance is the total amount of atmospheric scattered radiation received by the surface, and it is the core radiation source for low-illuminance shadowed areas.

[0049] S52: Determine the horizontal direct irradiance of sunlight on the Martian surface based on the solar incident irradiance. Then, determine the slope direct irradiance component based on the horizontal direct irradiance and the cosine of the relative incident angle obtained in step S51.

[0050] In some embodiments, the horizontal direct irradiance is obtained by the following formula: E dir,h =E0×cos(Z) Among them, E dir,h E0 represents the direct irradiance on a horizontal surface, while E0 represents the incident irradiance of the sun. The direct sunlight component on the slope is obtained by the following formula: E dir,s =E dir,h ×max(cos(i),0) / cos(Z); Among them, E dir,s This represents the direct sunlight component on the slope.

[0051] S53: Based on the slope and aspect of the Martian surface and the solar azimuth of each pixel in the digital elevation model obtained in step S1, an anisotropic model is used to estimate the slope scattering component. In this invention, the Perez paradigm is used to estimate the slope scattering component; that is, the input of the Perez paradigm is the slope, aspect, and solar azimuth, and the output is the slope scattering component.

[0052] In some other embodiments, the slope scattering component can also be obtained by isotropic approximation using the following formula: E dif,s ≈E dif,h ×(1+cosσ) / 2.

[0053] S54: Perform image segmentation on the Mars optical image to obtain a terrain-tagged image; for each pixel in the Mars optical image, calculate the average reflectance and visible sky factor of the pixel's neighborhood or the same tagged sub-image. The image segmentation method here is the same as in step S41, and will not be described again.

[0054] In this embodiment of the invention, the average reflectance is the arithmetic mean of the apparent reflectance of the horizontal plane within the specified statistical range in step S51. The statistical range includes two methods: one is a fixed-size pixel neighborhood centered on the target pixel, and the other is the same labeled sub-image obtained by terrain segmentation consistent with the method in step S41. During calculation, invalid values, saturated pixels, and extremely low-light pixels are removed. The visible sky factor, also known as the sky visibility factor, is a standardized parameter that quantitatively describes the impact of terrain occlusion on sky scattering radiation. Its value ranges from 0 to 1 and is calculated based on the digital elevation model obtained in step S1 using the sky hemispherical terrain occlusion integral method. The visible sky factor of the same labeled sub-image is the arithmetic mean of the visible sky factors of all pixels within that unit.

[0055] S55: Using the average reflectance and visible sky factor obtained in step S54, the total irradiance of the slope in the four regions is inverted based on the direct irradiance component obtained in step S52 and the scattering component obtained in step S53. Then, based on the total irradiance of the slope in the four regions and the surface radiance, the slope reflectance of the four regions is calculated.

[0056] In some embodiments, during the calculation of slope reflectance in gray shadow areas and sunny slope areas: The total irradiance of the slope in the gray shadow area or the sunny slope area is calculated by the following inversion formula: E tot_1 =k adj ×F v_1 ×ρ nei_1 +E dir,s_1 +E dif,s_1 ; Among them, E tot_1 E dir,s_1 and E dif,s_1 ρ represents the total irradiance, direct irradiance, and diffuse irradiance of the slope in the gray shadow region or the sunny slope region, respectively. nei_1 and F v_1 k represents the average reflectance of the gray shadow area and the visible sky factor of the sunny slope area, respectively. adj This represents the inversion coefficient, which is obtained empirically and can be determined from a flat region. The slope reflectance of the gray shadow area and the sunny slope area is calculated using the following formula: ρ s_1 (λ)=π·L surf_1 (λ) / E tot_1 (λ); Where, ρ s_1 (λ) represents the slope reflectance of the gray shadow area or the sunny slope area, λ represents the band, and L surf_1(λ) represents the surface radiance of the gray shadow area or the sunny slope area. The slope reflectance is calculated by the above formula to achieve the radiation consistency of the same ground feature under different slope aspects and slopes. In calculating the slope reflectance of the umbra and shadow regions: The total irradiance of the slope in the gray shadow area or the sunny slope area is calculated by the following inversion formula: E tot_2 =k adj ×F v_2 ×ρ nei_2 +E dif,s_2 ; Among them, E tot_2 and E dif,s_2 ρ represents the total irradiance and scattering component of the slope in the umbra or shadow region, respectively. nei_2 and F v_2 These represent the average reflectance of the umbra or the shadow region and the visible sky factor, respectively. The slope reflectance of the umbra and shadow regions is calculated using the following formula: ρ s_2 (λ)=π·L surf_2 (λ) / E tot_2 (λ); Where, ρ s_2 (λ) represents the slope reflectance of the umbra or shadow region, L surf_2 (λ) represents the surface radiance of the umbra or shadow region.

[0057] In this embodiment of the invention, if a correlation residual still exists between slope reflectivity and illumination coefficient, a minor residual correction can be performed based on the slope reflectivity using robust regression of segmentation units according to the semi-empirical link type method (i.e., step S4), correcting only a small trend. Existing semi-empirical models often only make empirical corrections based on the illumination coefficient, failing to adequately consider the contributions of sky scattering and reflection from nearby terrain; single physical models are often difficult to apply stably due to incomplete parameters. This invention explicitly considers direct solar radiation, sky scattering radiation, and radiation reflected from nearby terrain in the physical link (step S5), improving the physical rationality of radiation consistency across complex slopes and shaded areas.

[0058] S6: Stitch together the four regions corrected in step S4 or S5 to obtain and output the corrected Mars optical image. In this embodiment of the invention, step S6 also outputs a corrected quality control (QC) report, which includes correlation analysis, consistency analysis, shadow restoration analysis, and traceability analysis.

[0059] For correlation analysis, the correlation coefficient / regression slope between the statistically corrected Mars optical image and the IC is calculated, with a target result close to 0. For consistency analysis, the mean difference and CV change of similar land features on sunny slopes and in shadows are calculated, with a target result of reduced values. For shadow restoration analysis, the similarity between the predicted shadow area and the spectral morphology of adjacent sunny slopes is calculated, such as using SAM (Spectral Angle Mapper) to quantify the spectral morphological similarity between the predicted shadow area and adjacent sunny slopes. SAM, as a commonly used spectral shape similarity metric in remote sensing, is insensitive to brightness changes caused by illumination and terrain, quantifying only the morphological similarity of the spectral curves, and can accurately verify the restoration effect of the inherent spectral characteristics of land features after shadow correction. The SAM value ranges from 0 to π / 2; the smaller the value, the higher the spectral morphological similarity and the better the shadow restoration effect. For traceability analysis, the sample size and regression determination coefficient R of each tagged sub-image are recorded. 2 Random forest prediction error and quality mask are used to form an evaluable output. In the quality mask: for the corrected Mars optical image obtained in step S4, the quality label Q-umbra=1 is directly output for the umbra region of the Mars optical image, and for other regions, pixels marked as unreliable due to low illumination are labeled as Q-lowIC=1; for the corrected Mars optical image obtained in step S5, the quality label Q-shadow=1 is directly output for the umbra and shadow regions of the Mars optical image, and for other regions, pixels marked as unreliable due to low illumination are labeled as Q-lowIC=1.

[0060] Existing technologies often rely on subjective visual comparison or single statistical measures to evaluate correction effectiveness, making it difficult to distinguish the sources of error in different shadow types and different local scenes. This invention employs methods such as four-category regional grouping statistics, correlation analysis, error indices, spectral consistency testing, and quality labeling records to construct a quantifiable, comparable, and verifiable evaluation system. Furthermore, most existing technologies only output corrected images or provide only a few statistical indicators, lacking parametric maps, quality masks, and classification results. In addition to outputting corrected images, this invention also outputs four types of shadow masks, segmentation unit label maps, unit parameter maps, low-light and shadow quality masks, and QC reports, enhancing result traceability and engineering management capabilities.

[0061] It should be understood that the various forms of processes shown above can be used to reorder, add, or delete steps. For example, the steps described in this invention disclosure can be executed in parallel, sequentially, or in different orders, as long as the desired result of the technical solution disclosed in this invention can be achieved, and this is not limited herein.

[0062] The specific embodiments described above do not constitute a limitation on the scope of protection of this invention. Those skilled in the art should understand that various modifications, combinations, sub-combinations, and substitutions can be made according to design requirements and other factors. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of this invention should be included within the scope of protection of this invention.

Claims

1. A method for terrain correction of Martian surface optical images based on shadow classification, characterized in that, include: S1: Acquire optical images of the Martian surface and a digital elevation model that matches the Martian optical images; The illumination coefficient of sunlight illuminating the surface of Mars was determined based on the digital elevation model. S2: Based on the illumination coefficient obtained in step S1, the surface of Mars is divided into the umbra, gray shadow, shadow area, and sunny slope area. S3: Determine whether the atmospheric parameters of the Mars optical image in step S1 are complete. If complete, proceed to step S4; otherwise, proceed to step S5. S4: Perform image segmentation on the Mars optical image, and correct the sunny slope area and gray shadow area in each labeled sub-image according to the segmentation results. At the same time, use machine learning methods to correct the shadow area. S5: Determine the slope scattering component and slope direct component of the Martian surface when sunlight shines on it; perform image segmentation on the Martian optical image, and based on the segmentation results, correct the four regions in each labeled sub-image by combining the slope scattering component and the slope direct component. S6: Stitch together the four regions corrected in step S4 or S5 to obtain and output the corrected Mars optical image.

2. The method for terrain correction of Martian surface optical images based on shadow classification according to claim 1, characterized in that, Step S1, which involves determining the illumination coefficient of the Martian surface based on the digital elevation model, includes: Calculate the slope and aspect of the Martian surface for each pixel in the digital elevation model; The preliminary illumination coefficient is calculated using the following formula: IC0=cosσ-cosZ+sinσ×sinZ×cos(β-ω); Where IC0 represents the initial illumination coefficient, σ represents the slope, β represents the aspect, ω represents the azimuth angle when sunlight shines on the surface of Mars, and Z represents the zenith angle when sunlight shines on the surface of Mars. Set an illuminance threshold for low-illuminance pixels in the digital elevation model whose initial illuminance coefficient is less than the lower threshold. Then, apply the following formula to robustly constrain these low-illuminance pixels to obtain their final illuminance coefficients: IC=clip((IC0-T min ) / (T max -T min ),0,1)×IC0; Where IC represents the final illumination coefficient, clip represents the cutoff function, and T min T represents the threshold value under illuminance. max This represents the threshold value for illuminance.

3. The method for terrain correction of Martian surface optical images based on shadow classification according to claim 1, characterized in that, Step S2 includes: The region with an illumination coefficient less than or equal to 0 is designated as the first region, and the region with an illumination coefficient greater than 0 but whose brightness in Mars optical images is significantly lower than that of similar terrain or the local average is designated as the second region. Both the first and second areas are shady slopes, while the remaining areas are sunny slopes. The umbra is the set of pixels in the shaded slope region that completely block sunlight from reaching the surface of Mars when it is directly overhead. The set of pixels in the shaded slope area with an illumination coefficient greater than 0 but less than or equal to a set intensity threshold is the gray shadow area. In the shaded area, ray tracing is performed along the direction of the sun's incidence on the outer edge of the umbra. The area blocked by the adjacent terrain is the shadow area.

4. The method for terrain correction of Martian surface optical images based on shadow classification according to claim 1, characterized in that, Step S4 includes: S41: Perform image segmentation on the Mars optical image in step S1 to obtain and output terrain label images; in each label sub-image, use the sunny slope area and gray shadow area obtained in step S2 to establish the regression relationship between each band and the illumination coefficient in the Mars optical image; S42: Correct the gray shadow area and the sunny slope area using the regression relationship obtained in step S41; retain the umbra area obtained in step S2; and correct the umbra area using machine learning based on the sunny slope area adjacent to the shadow area in step S2.

5. The method for terrain correction of Martian surface optical images based on shadow classification according to claim 1, characterized in that, In step S42: The following rotational correction model is used to correct the gray shadow area and the sunny slope area: I corr (λ)=I(λ)-a k (λ)×(IC-cosZ); Among them, I corr (λ) represents the corrected gray shadow region or sunny slope region, λ represents the band, I(λ) represents the corresponding original gray shadow region or sunny slope region, a k (λ) represents the slope corresponding to the k-th label sub-image in the regression relationship, IC represents the illumination coefficient, and Z represents the zenith angle when sunlight illuminates the surface of Mars.

6. The method for terrain correction of Martian surface optical images based on shadow classification according to claim 1, characterized in that, The process of correcting the shadow area using machine learning in step S42 includes: The SEVI paradigm is used to construct the shadow sensitivity index and obtain the texture information of the shadow area; In the same sub-image of the label, the sunny slope area adjacent to the shadow area is randomly sampled to obtain the sampled image; A random forest regressor was trained using sampled images, shadow sensitivity index, and texture information to obtain a shadow prediction model; The image of the shadow region is input into the shadow prediction model. The obtained predicted image replaces the corresponding image in the shadow region, or the obtained predicted image is added to the corresponding image to obtain the corrected shadow region.

7. The method for terrain correction of Martian surface optical images based on shadow classification according to claim 1, characterized in that, Step S5 includes: S51: Convert the Mars optical image obtained in step S1 into the surface radiance and horizontal apparent reflectance of the Martian surface when sunlight shines on it; determine the cosine of the relative incident angle when sunlight shines on the Martian surface based on the digital elevation model obtained in step S1. S52: Determine the horizontal direct irradiance of sunlight on the surface of Mars based on the solar incident irradiance, and then determine the slope direct irradiance component based on the horizontal direct irradiance and the cosine of the relative incident angle obtained in step S51. S53: Based on the slope and aspect of the Martian surface and the azimuth of the sun, which are reflected in each pixel of the digital elevation model obtained in step S1, the anisotropic model is used to estimate the scattering components of the slope surface. S54: Perform image segmentation on the Mars optical image to obtain a terrain-tagged image; for each pixel in the Mars optical image, calculate the average reflectance and visible sky factor of the pixel neighborhood or the same tag sub-image for each pixel; S55: Using the average reflectance and visible sky factor obtained in step S54, the total irradiance of the slope in the four regions is inverted based on the direct irradiance component obtained in step S52 and the scattering component obtained in step S53. Then, based on the total irradiance of the slope in the four regions and the surface radiance, the slope reflectance of the four regions is calculated.

8. The method for terrain correction of Martian surface optical images based on shadow classification according to claim 7, characterized in that, Step S51 includes: The Mars optical image is converted into an equivalent DN, and then the equivalent DN is converted into top atmospheric radiance according to the sensor gain and bias or radiometric calibration coefficient. Using a radiative transfer model, the top atmospheric radiance obtained in step S51 is converted into surface outgoing radiance and horizontal apparent reflectance. The slope and aspect of the Martian surface are calculated for each pixel in the digital elevation model, and the cosine of the relative incident angle of sunlight on the Martian surface is determined by the following formula: cos(i)=cosσ-cosZ+sinσ×sinZ×cos(β-ω); Where i represents the relative angle of incidence, σ represents the slope, β represents the aspect, ω represents the azimuth angle when sunlight shines on the surface of Mars, and Z represents the zenith angle when sunlight shines on the surface of Mars.

9. The method for terrain correction of Martian surface optical images based on shadow classification according to claim 8, characterized in that, In step S52: The direct irradiance on the horizontal plane can be obtained using the following formula: E dir,h =E0×cos(Z) Among them, E dir,h E0 represents the direct irradiance on a horizontal surface, while E0 represents the incident irradiance of the sun. The direct sunlight component on the slope is obtained by the following formula: AND dir,s =E dir,h ×max(cos(i),0) / cos(Z); Among them, E dir,s This represents the direct sunlight component on the slope.

10. The method for terrain correction of Martian surface optical images based on shadow classification according to claim 9, characterized in that, In step S55: In calculating the slope reflectance of the gray shadow area and the sunny slope area: The total irradiance of the slope in the gray shadow area or the sunny slope area is calculated by the following inversion formula: AND tot_1 =k adj ×F v_1 ×ρ nei_1 +E dir,s_1 +E dif,s_1 ; Among them, E tot_1 E dir,s_1 and E dif,s_1 ρ represents the total irradiance, direct irradiance, and diffuse irradiance of the slope in the gray shadow region or the sunny slope region, respectively. nei_1 and F v_1 k represents the average reflectance of the gray shadow area and the visible sky factor of the sunny slope area, respectively. adj Indicates the inversion coefficients; The slope reflectance of the gray shadow area and the sunny slope area is calculated using the following formula: r s_1 (λ)=π·L surf_1 (l) / E tot_1 (l); Where, ρ s_1 (λ) represents the slope reflectance of the gray shadow area or the sunny slope area, λ represents the band, and L surf_1 (λ) represents the surface radiance of the gray shadow area or the sunny slope area; In calculating the slope reflectance of the umbra and shadow regions: The total irradiance of the slope in the gray shadow area or the sunny slope area is calculated by the following inversion formula: E tot_2 =k adj ×F v_2 ×ρ nei_2 +E dif,s_2 ; Among them, E tot_2 and E dif,s_2 ρ represents the total irradiance and scattering component of the slope in the umbra or shadow region, respectively. nei_2 and F v_2 These represent the average reflectance of the umbra or the shadow region and the visible sky factor, respectively. The slope reflectance of the umbra and shadow regions is calculated using the following formula: r s_2 (λ)=π·L surf_2 (l) / E tot_2 (l); Where, ρ s_2 (λ) represents the slope reflectance of the umbra or shadow region, L surf_2 (λ) represents the surface radiance of the umbra or shadow region.

Images