Remote sensing image terrain radiation correction method, device and equipment

Abstract

The invention relates to the field of remote sensing image processing, and provides a remote sensing image terrain radiation correction method, device and equipment, and the method comprises the steps: obtaining a to-be-corrected remote sensing image, calculating the terrain parameter of each pixel, dividing the to-be-corrected remote sensing image into a shadow region and a non-shadow region, and determining the direction weighting function of each pixel; calculating a topographic modulation spherical harmonic primary function corresponding to each pixel; training samples are selected in the shadow area and the non-shadow area respectively, a linear model is constructed, and spherical harmonic illumination coefficients representing overall three-dimensional illumination distribution of the shadow area and the non-shadow area are solved; and based on the topographic modulation spherical harmonic primary function, combining the spherical harmonic illumination coefficient corresponding to the area to which each pixel belongs, reconstructing the average incident radiance of each pixel, and then performing inversion to obtain the real reflectivity of the earth surface. According to the method, the problem of inaccurate radiation distortion correction of the remote sensing image under the complex terrain in the prior art is solved, and accurate inversion of the real reflectivity of the earth surface is realized.

Description

Technical Field

[0001] This invention relates to the field of remote sensing image processing, and in particular to a method, apparatus and equipment for topographic radiometric correction of remote sensing images. Background Technology

[0002] Remote sensing technology has become a core pillar of modern Earth observation systems, providing an irreplaceable technical means for acquiring large-scale, multi-temporal surface information. Especially in mountainous and hilly areas with complex terrain, remote sensing images can provide crucial data such as land cover, vegetation status, and environmental changes, supporting quantitative applications in many fields, including resource exploration, ecological monitoring, and disaster assessment. However, due to the influence of topographic relief, the surface radiation signals received by remote sensing sensors undergo significant geometric and radiometric distortions. These distortions directly limit the accuracy of quantitative applications of remote sensing data in complex terrain areas.

[0003] To eliminate the influence of topography, various topographic radiation correction methods have been developed. Early solutions, such as the band ratio method and empirical statistical models, often rely on specific assumptions (e.g., proportional changes in reflectivity across different bands) or require flat area references, resulting in poor universality under complex terrain and easy loss of spectral information. Classical models based on the Lambertian assumption are widely used, but they often lead to undercorrection or overcorrection in shaded areas because they ignore the non-Lambertian reflectance characteristics of the surface and the incident-observation geometry. The later-developed Minnaert model introduces an empirical exponent to simulate non-Lambertian effects, but the physical meaning of this parameter is unclear, making it difficult to accurately characterize the biaxial reflectance characteristics under complex lighting conditions. Some more advanced models begin to consider radiation component decomposition or nearby topographic reflection, but they still use simplified or uniform treatments for the contributions of sky scattering or nearby reflection, failing to accurately quantify the differences in radiation under different incident directions.

[0004] Therefore, how to construct a topographic radiation correction method that can strictly follow the physical radiation transfer mechanism and efficiently handle the multi-directional illumination effects under complex terrain, so as to accurately restore the true spectral characteristics of the land surface, has become a key technical bottleneck for improving the quantitative application level of remote sensing data in mountainous areas, and is also an important issue that urgently needs to be solved in the field of remote sensing preprocessing. Summary of the Invention

[0005] This invention provides a method, apparatus, and device for topographic radiometric correction of remote sensing images, which solves the problem of inaccurate radiometric distortion correction of remote sensing images under complex terrain in the prior art, and realizes accurate inversion of the true reflectance of the land surface.

[0006] This invention provides a method for topographic radiometric correction of remote sensing images, comprising the following steps: Acquire the remote sensing image to be corrected and its corresponding digital elevation model, and calculate the terrain parameters of each pixel in the remote sensing image to be corrected based on the digital elevation model; The remote sensing image to be corrected is divided into shadowed and unshadowed areas based on the terrain parameters. The orientation weight function of each pixel in each discrete incident direction is determined based on the terrain parameters; Based on the directional weight function of each pixel itself, calculate its corresponding directional spherical harmonic basis function; based on the directional spherical harmonic basis function of each pixel in all discrete incident directions, calculate its corresponding terrain modulation spherical harmonic basis function. Representative pixels were selected as training samples from both the shaded and unshaded areas. Based on the reflection characteristics of the training samples and the terrain-modulated spherical harmonic basis function, a linear model is constructed, and the spherical harmonic illumination coefficients characterizing the overall three-dimensional illumination distribution of the shadowed and unshadowed areas are solved. Based on the topographically modulated spherical harmonic basis function of each pixel, and combined with the spherical harmonic illumination coefficient corresponding to its region, the average incident radiance of each pixel is reconstructed. Based on the reflected radiance observed by the sensor and the reconstructed average incident radiance, the true surface reflectance after eliminating the influence of topography is obtained by inversion.

[0007] According to a method for topographic radiometric correction of remote sensing images provided by the present invention, the step of calculating the topographic parameters of each pixel in the remote sensing image to be corrected based on the digital elevation model includes: determining the surface normal vector of each pixel in the remote sensing image to be corrected based on the digital elevation model; calculating the topographic parameters of each pixel according to the surface normal vector, wherein the topographic parameters include at least slope, aspect, solar incidence angle, and maximum shading angle in multiple discrete azimuth directions.

[0008] According to a remote sensing image topographic radiometric correction method provided by the present invention, the method determines the directional weight function for each discrete incident direction of each pixel based on the topographic parameters, including: selecting a corresponding weight calculation mode based on whether each discrete incident direction is the direction of direct sunlight and whether the current pixel is located in a shadow area; the weight calculation mode includes at least: direct sunlight and non-shadow mode, direct sunlight and shadow mode, non-direct sunlight and non-shadow mode, and non-direct sunlight and shadow mode; in the non-shadow mode, the weight function is jointly determined by the proportion of direct sunlight component and the sky scattering visibility factor; in the shadow mode, the weight function is jointly determined by the proportion of environmental reflectance irradiance, the neighboring reflectance visibility enhancement factor, and the neighboring reflectance ratio.

[0009] According to a remote sensing image topographic radiometric correction method provided by the present invention, representative pixels are selected as training samples in the shaded area and the non-shaded area respectively. The method includes: in the shaded area, generating a cumulative histogram based on the solar incidence angle of all pixels, removing extreme values, and dividing the remaining solar incidence angle range into several intervals; in each interval, selecting the pixel whose brightness is closest to the average brightness of the pixels in the current interval as the training sample of the current interval; and selecting training samples in the non-shaded area using the same principle as selecting training samples in the shaded area.

[0010] According to the present invention, a method for topographic radiometric correction of remote sensing images includes the following steps: constructing a linear model based on the reflectance characteristics of the training samples and the topographic modulated spherical harmonic basis functions, and solving for the spherical harmonic illumination coefficients characterizing the overall three-dimensional illumination distribution of the shadowed and unshadowed regions. The method comprises: performing radiometric calibration and atmospheric correction on the remote sensing image to be corrected to obtain a preliminary surface reflectance; calculating the reflected radiance in the sensor observation direction based on the preliminary surface reflectance using a radiative transfer model; constructing and calculating a set of real spherical harmonic basis function vectors; for each training sample, calculating the incident radiance of each training sample based on the reflected radiance, and calculating the corresponding topographic modulated spherical harmonic basis function based on the real spherical harmonic basis function vectors; establishing an overdetermined linear equation set with incident radiance as the dependent variable and the topographic modulated spherical harmonic basis functions as the independent variable, based on the training sample sets of the shadowed and unshadowed regions respectively; and solving the overdetermined linear equation set using the least squares method to obtain the spherical harmonic illumination coefficient vectors of the shadowed and unshadowed regions respectively.

[0011] According to a remote sensing image topographic radiometric correction method provided by the present invention, the step of calculating the incident radiance of each training sample based on the reflected radiance and calculating the topographic modulation spherical harmonic basis function corresponding to each training sample based on the real spherical harmonic basis function vector includes: using the real spherical harmonic basis function vector and the direction weight function of each training sample in each discrete incident direction; calculating the directional spherical harmonic basis function of each training sample in each discrete direction based on the direction weight function; and performing a weighted average of the directional spherical harmonic basis functions of all incident directions for each training sample to obtain the topographic modulation spherical harmonic basis function corresponding to the training sample.

[0012] According to the method for topographic radiometric correction of remote sensing images provided by the present invention, the empirical reflectance of each training sample in the sensor observation direction is calculated based on the bidirectional reflectance distribution function parameters. The method includes: calculating the bidirectional reflectance of the training sample in the sensor observation direction and any discrete incident direction in the upper hemisphere based on the kernel-driven model weight parameters and the pre-calculated volume scattering kernel function matrix and geometric optics kernel function matrix; and taking a weighted average of the bidirectional reflectances in all discrete incident directions according to the solid angle to obtain the empirical reflectance of the training sample.

[0013] According to a remote sensing image topographic radiometric correction method provided by the present invention, the method reconstructs the average incident radiance of each pixel based on the topographic modulation spherical harmonic basis function of each pixel and the spherical harmonic illumination coefficient corresponding to its region. This includes: for pixels in shaded areas, weighting the real spherical harmonic basis function vector using its directional weight function in each discrete incident direction to obtain the directional spherical harmonic basis function in each incident direction, and then averaging the directional weight function by solid angle to obtain the topographic modulation spherical harmonic basis function of the pixel in the shaded area; for pixels in non-shaded areas, using its directional weight function in each discrete incident direction to obtain the directional spherical harmonic basis function in each incident direction. The directional weighting function for each discrete incident direction is used to weight the real spherical harmonic basis function vector to obtain the directional spherical harmonic basis function for each incident direction. The topographic modulation spherical harmonic basis function for the non-shaded area is obtained by weighting the solid angle. The topographic modulation spherical harmonic basis function for the shaded area is then multiplied by the spherical harmonic illumination coefficient vector of the shaded area to obtain the average incident radiance of each pixel in the shaded area. The topographic modulation spherical harmonic basis function for the non-shaded area is then multiplied by the spherical harmonic illumination coefficient vector of the non-shaded area to obtain the average incident radiance of each pixel in the non-shaded area.

[0014] According to a method for topographic radiometric correction of remote sensing images provided by the present invention, the step of inverting the surface reflectance after eliminating the influence of topography based on the reflected radiance observed by the sensor and the reconstructed average incident radiance includes: taking the reflected radiance observed by the sensor for each band of the remote sensing image to be corrected as the numerator; taking the reconstructed average incident radiance of the corresponding band as the denominator; performing ratio calculation pixel by pixel and band by band to obtain the topographically corrected surface reflectance of each band; and combining the topographically corrected surface reflectance of each band to form the surface reflectance after eliminating the influence of topographic undulation.

[0015] The present invention also provides a remote sensing image topographic radiometric correction device, comprising the following modules: The terrain parameter calculation module is used to acquire the remote sensing image to be corrected and its corresponding digital elevation model, and to calculate the terrain parameters of each pixel in the remote sensing image to be corrected based on the digital elevation model. The shadow region segmentation module is used to divide the remote sensing image to be corrected into shadow regions and non-shadow regions according to the terrain parameters. The weighting function determination module is used to determine the directional weighting function of each pixel in each discrete incident direction based on the terrain parameters. The terrain-modulated spherical harmonic basis function determination module is used to calculate the corresponding directional spherical harmonic basis function based on the directional weight function of each pixel itself; and to calculate the corresponding terrain-modulated spherical harmonic basis function based on the directional spherical harmonic basis function of each pixel in all discrete incident directions. The training sample selection module is used to select representative pixels as training samples in the shaded area and the non-shaded area respectively. The illumination coefficient solution module is used to construct a linear model based on the reflection characteristics of the training samples and the terrain modulation spherical harmonic basis function, and solve for the spherical harmonic illumination coefficients that characterize the overall three-dimensional illumination distribution of the shadow area and the non-shadow area. The average incident radiance reconstruction module is used to reconstruct the average incident radiance of each pixel based on the terrain-modulated spherical harmonic basis function of each pixel and the spherical harmonic illumination coefficient corresponding to its region. The true reflectance inversion module is used to invert the true surface reflectance after eliminating the influence of topography based on the reflected radiance observed by the sensor and the reconstructed average incident radiance.

[0016] The present invention also provides an electronic device, including a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the remote sensing image topographic radiometric correction method as described above.

[0017] The present invention also provides a non-transitory computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the remote sensing image topographic radiometric correction method as described above.

[0018] The present invention also provides a computer program product, including a computer program that, when executed by a processor, implements the remote sensing image topographic radiometric correction method as described above.

[0019] This invention provides a method, apparatus, and device for topographic radiometric correction of remote sensing images, which offers the following advantages: It precisely calculates the topographic parameters of each pixel based on a digital elevation model, accurately dividing shadowed and unshadowed areas accordingly, thus laying the foundation for subsequent differential processing. Furthermore, by establishing a directional weight function for each pixel, it precisely quantifies the combined contributions of direct sunlight, sky scattering, and reflection from adjacent terrain from different incident directions, overcoming the limitations of traditional methods in simplifying complex lighting conditions. It models the three-dimensional illumination distribution by introducing spherical harmonic functions. Representative training samples are selected from both shadowed and unshadowed areas, and linear models are constructed using their reflection characteristics and directional weight functions to solve for the spherical harmonic irradiance coefficients characterizing the overall illumination pattern of each region. This process achieves a quantitative mapping from physical radiative transfer to mathematical function expression, making the characterization of irradiance distribution under complex terrain both physically rigorous and computationally feasible. By calculating the topographic modulated spherical harmonic basis function of each pixel using its own directional weight function, and combining it with the spherical harmonic irradiance coefficients of the region, the incident irradiance that each pixel should receive is accurately reconstructed. By calculating the ratio between the reflected radiance observed by the sensor and the reconstructed average incident radiance, the irradiance variation caused by topographic relief is effectively eliminated, and the true surface reflectance, which is closer to the inherent spectral characteristics of the ground objects, is retrieved. Therefore, this invention significantly improves the radiometric correction accuracy of remote sensing images in complex terrain areas and enhances the spectral consistency and spatial uniformity of the images. Attached Figure Description

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

[0021] Figure 1 This is a flowchart illustrating the remote sensing image topographic radiometric correction method provided by the present invention.

[0022] Figure 2 This is a schematic diagram of the overall process of remote sensing image topographic radiometric correction provided by the present invention.

[0023] Figure 3 The present invention provides the Landsat8 OLI remote sensing image to be corrected, the corresponding ALOS DSM, and the shadow recognition results.

[0024] Figure 4 This is a schematic diagram of the distribution of weight parameters of MCD43A1 and BRDF after interpolation to the OLI scale provided by the present invention (taking the blue band as an example).

[0025] Figure 5a , Figure 5b This is a band-by-band display diagram of the original reflectance, the reconstructed average incident radiance, and the resulting reflectance image provided by the present invention.

[0026] Figure 6 This is a comparison image of various color composite effects of the original reflectance image and the resulting reflectance image provided by the present invention.

[0027] Figure 7 This is a schematic diagram of the structure of the remote sensing image topographic radiation correction device provided by the present invention.

[0028] Figure 8 This is a schematic diagram of the structure of the electronic device provided by the present invention. Detailed Implementation

[0029] To make the objectives, technical solutions, and advantages of this invention clearer, the technical solutions of this invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of this invention. All other embodiments obtained by those skilled in the art based on the embodiments of this invention without creative effort are within the scope of protection of this invention.

[0030] The following is combined with Figures 1-8 The embodiments of the present invention are described in detail.

[0031] The remote sensing image topographic radiometric correction method provided in this embodiment of the invention is executed by a remote sensing image topographic radiometric correction device, which can be configured in a computer. The computer can be a local computer or a cloud computer. The local computer can be a computer, tablet, etc., and no specific limitation is made here.

[0032] Figure 1 This is a flowchart illustrating the remote sensing image topographic radiometric correction method provided by the present invention, as shown below. Figure 1 As shown, the method includes the following steps: S110. Obtain the remote sensing image to be corrected and its corresponding digital elevation model, and calculate the terrain parameters of each pixel in the remote sensing image to be corrected based on the digital elevation model.

[0033] According to the method for topographic radiometric correction of remote sensing images provided by the present invention, radiometric calibration and atmospheric correction are performed on the remote sensing image to be corrected to obtain a preliminary surface reflectance; based on the radiative transfer model, the reflected radiance in the sensor observation direction is calculated according to the preliminary surface reflectance.

[0034] Specifically, step 1: data preprocessing.

[0035] (1) Radiometric calibration and FLAASH atmospheric correction were performed on the remote sensing images to obtain the surface reflectance based on the Lambertian assumption and ignoring the radiation effect of topographic relief. The 6S (Second Simulation of the Satellite Signal in the Solar Spectrum) radiative transfer model was used to simulate the atmospheric conditions and imaging geometry corresponding to the remote sensing images, and to calculate the direct solar irradiance received by the Earth's surface under the assumption of a horizontally homogeneous Lambertian body. Sky diffuse irradiance and ambient irradiance The total incident irradiance is thus calculated. and the percentage of each component , and .

[0036] (1) in, Under the assumption of a horizontally uniform Lambertian body, the total incident irradiance received by the Earth's surface is typically expressed in watts per square meter per micrometer. The subscript "flat" indicates "horizontal," and λ indicates that the value is related to wavelength (band). Under the same assumptions, this represents the direct solar irradiance received by the Earth's surface. This represents the amount of sky diffuse irradiance received by the Earth's surface under the same assumptions. "dif" is an abbreviation for diffuse. env is an abbreviation for environment, representing the ambient reflected irradiance (reflected light from the surrounding terrain) received by the Earth's surface under the same assumptions.

[0037] And further calculate the corresponding reflected radiance in the sensor's observation direction. : (2) in, This represents the reflected radiance received by the sensor in the direction of observation. The unit is typically watts per square meter per steradian per micrometer. 'r' indicates "reflected". The raw value represents the initial surface reflectance obtained after radiometric calibration and atmospheric correction.

[0038] (2) Based on the multi-band remote sensing image, it is first converted into a single-band broadband albedo image, and the ratio between the MCD43A1 pixel size and the remote sensing image pixel size is calculated. .

[0039] For a given pixel center in MCD43A1, calculate its corresponding position in the albedo image, and use this position as the center. A sub-image of albedo is cropped based on the radius. Pixel pairs are randomly selected from this sub-image, and the Euclidean distance and albedo semivariance between each pixel pair are calculated. Then, the different distances are statistically analyzed. Corresponding mean semivariance By fitting a spherical variogram model using nonlinear least squares: (3) Here, γ(d) is the semivariance, a core indicator for measuring spatial correlation in geostatistics. d is the Euclidean distance between two pixels. C0 is the nugget value, representing the semivariance when the distance is 0. Theoretically, it should be 0, but in actual data, it exists due to measurement errors or variations smaller than the sampling scale, reflecting random noise or microstructural variations. C is the structural variance (a part of the sill value), representing the variation caused by spatial autocorrelation structure. h is the range, the effective distance of spatial correlation. When d > h, there is no longer spatial correlation between the two points, and the semivariance reaches a stable value (sill value). C0 + C is the sill value, representing the stable upper limit of the semivariance when the distance is sufficiently large.

[0040] Get Gold Value , base value and variable range This method is used to calculate the range of all pixels in the MCD43A1 image.

[0041] For a certain pixel in a remote sensing image, in the context of... Search for the center of an MCD43A1 pixel within a neighborhood of the radius; if the distance between the remote sensing image pixel and the center of an MCD43A1 pixel is less than the range of that MCD43A1 pixel, calculate the distance-to-range ratio; assign the BRDF weight parameters of the MCD43A1 pixel with the smallest ratio to the remote sensing image pixel. This interpolation method ensures that all pixels in the remote sensing image receive the three weight parameters of the kernel-driven BRDF. , and .

[0042] (3) Within the upper hemisphere, the zenith angle is... Discretize the interval as Each sampling point, with azimuth angle as Discretize the interval as 100 sampling points, thus constructing a collection of 100 sampling points. A discrete sample set of independent directions is generated. Then, the Ross-Thick volume scattering kernel function between each pair of incident and reflected directions is calculated. With Li-Sparse Reciprocal geometric optics kernel function : (4) in, , The volume scattering nucleus is used to simulate the reflection characteristics caused by multiple scattering of photons within the vegetation canopy. The zenith angle is the angle between the direction of the sun and the normal to the Earth's surface. The zenith angle is the angle between the sensor's observation direction and the Earth's surface normal, representing the reflected (observed) zenith angle. ε is the relative azimuth angle, the angle between the projections of the solar incident direction and the sensor observation direction onto the horizontal plane. ε is the scattering phase angle, the angle between the solar incident direction and the sensor observation direction.

[0043] (5) It is a geometric optics kernel that simulates the reflection characteristics caused by the shadows and shading effects produced by the three-dimensional structure of the Earth's surface. "geo" is an abbreviation for geometric.

[0044] The formulas involved in this calculation process are summarized as follows: (6) in, ′, ε′ represents the corrected incident and observed zenith angles, which were transformed from the original angles based on the canopy shape parameters to accommodate the "sparse canopy" assumption in the Li-Sparse model. ε′ represents the angle based on the corrected angles. ′, The calculated corrected scattering phase angle. sec(x) is the secant function, sec(x) = 1 / cos(x). O is a measure representing the overlapping region of the two projections (the projections of the sun and the sensor direction onto the horizontal plane). t is an intermediate angle variable related to the geometry of the overlapping region, used to calculate the overlapping area. D is a distance parameter describing the distance between the projection points of the sun direction and the sensor direction on the horizontal projection plane. b / r is the shape parameter ratio, where b is the vertical semi-axis of the ellipsoid (used to simulate the canopy), and r is the horizontal semi-axis. A ratio of 1 indicates that the canopy shape is simplified to a sphere. h / b is the height parameter ratio, where h is the height from the center of the ellipsoid to the ground. A ratio of 2 is the standard parameterization of the Li-SparseReciprocal kernel, indicating that the canopy height is twice its vertical semi-axis.

[0045] By radiometrically calibrating and performing FLAASH atmospheric correction on the original remote sensing images, the interference of atmospheric scattering and absorption on the surface signal was effectively eliminated, resulting in a relatively pure initial value of surface reflectance. This laid an accurate radiometric benchmark for all subsequent quantitative analyses. The 6S radiative transfer model was used to accurately simulate the components and proportions of direct solar radiation, sky scattering, and environmental reflected irradiance during imaging. This not only provided a standard radiometric reference under horizontal surface conditions but, more importantly, quantified the contribution weights of different radiation sources to surface illumination, providing key physical parameters for subsequent detailed modeling of the differentiated impacts of terrain on various types of illumination. An innovative spatial interpolation method was used to adapt MODIS BRDF product parameters to high-resolution remote sensing images, enabling each pixel to obtain kernel-driven model weight parameters characterizing its non-Lambertian reflectance properties. This operation successfully combined macroscopic BRDF prior knowledge with local pixel characteristics, maintaining the rigor of the physical mechanism while overcoming the limitations of insufficient direct observation data. By pre-calculating the discrete directional sample set and corresponding kernel function matrix in the upper hemisphere space, the complex BRDF directionality calculation is transformed into a data structure that can be efficiently called, significantly improving the computational efficiency of the entire correction process. In summary, step 1, through systematic data preparation and physical parameterization, constructs a data processing framework that combines radiometric accuracy, physical completeness, and computational feasibility, providing reliable data support and theoretical assurance for subsequent terrain lighting modeling and correction.

[0046] According to the present invention, a method for topographic radiometric correction of remote sensing images is provided, which calculates the topographic parameters of each pixel in the remote sensing image to be corrected based on a digital elevation model, including: determining the surface normal vector of each pixel in the remote sensing image to be corrected based on the digital elevation model; and calculating the topographic parameters of each pixel based on the surface normal vector, wherein the topographic parameters include at least slope, aspect, solar incidence angle, and maximum shading angle in multiple discrete azimuth directions.

[0047] Specifically, a three-dimensional spatial coordinate system is established with the center of the top-left pixel in the remote sensing image as the origin: the X-axis points due east, the Y-axis points due north, and the Z-axis points towards the zenith. Based on the pixel spatial resolution and digital surface model of the remote sensing image, the three-dimensional spatial coordinates corresponding to the center of each pixel in the image are determined. Using a triangular mesh model, the corresponding surface normal vector is calculated based on the three-dimensional coordinates of each pixel center, thereby obtaining the slope, aspect, and solar incidence angle relative to the sun for each pixel.

[0048] Define a terrain occlusion analysis search radius. For any pixel in the image, search along each azimuth direction in the discrete direction sample set constructed in step 1. Within this search radius, calculate and record the maximum terrain shading angle that can shade the pixel. Simultaneously, the positions of neighboring pixels that generate the maximum occlusion angle are recorded. This process yields the maximum occlusion angle of all pixels in the remote sensing image in each azimuth direction, as well as the corresponding occluded pixel location.

[0049] This step provides a precise spatial geometry and illumination shading analysis foundation for topographic radiometric correction. By establishing a three-dimensional coordinate system and calculating the surface normal vector for each pixel, the two-dimensional image information is transformed into topographic parameters (slope, aspect, solar incidence angle) suitable for three-dimensional spatial analysis, achieving precise quantification of the microscopic geometry of the land surface. By introducing a search radius and systematically calculating the maximum shading angle in each azimuth direction, not only is the complete spatial pattern of each pixel's shading by surrounding terrain accurately depicted, but the location of the shading source is also recorded simultaneously, providing crucial spatial correlation information for subsequent quantification of the reflection contribution of neighboring terrain.

[0050] Step 2: Terrain Parameter Calculation and Shadow Determination A three-dimensional spatial coordinate system is established with the center of the top-left pixel in the remote sensing image as the origin, the X-axis pointing due east, the Y-axis pointing due north, and the Z-axis pointing towards the zenith. Based on the pixel size and DSM of the remote sensing image, the three-dimensional coordinates of each pixel center are determined, and the normal vector is calculated using a triangular mesh model to obtain the slope, aspect, and solar incidence angle of each pixel.

[0051] Set the search radius. For a specific remote sensing image pixel, search along a certain azimuth direction. (The discrete direction sample set from step 1) Within the search radius, find the maximum occlusion angle. Simultaneously record the pixel position that generates the maximum occlusion angle. This allows us to obtain the maximum occlusion angle and the location of the occluded pixel in each azimuth direction of the remote sensing image.

[0052] S120. Divide the remote sensing image to be corrected into shaded and unshaded areas according to the terrain parameters.

[0053] Specifically, for any pixel in the image, shadow identification is performed. If the maximum shading angle of the pixel in the direction of the solar azimuth is greater than the solar altitude angle at the time of image formation, the pixel is determined to be in the terrain shadow area; if the solar incidence angle at the pixel is greater than 90°, the pixel is designated as the umbra. Finally, the shadow area and the umbra are merged, thereby dividing the entire remote sensing image into shadowed and unshadowed areas.

[0054] The automatic segmentation of shadow areas is based on strict geometric criteria (comparison of maximum shading angle and solar altitude angle, solar incidence angle threshold), which overcomes the shortcomings of subjective and inaccurate shadow identification in traditional methods and achieves objective and accurate separation of shadow and non-shadow areas.

[0055] Step 3: Construct and compute a set of real spherical harmonic basis function vectors.

[0056] Set the number of zones in the spherical harmonic expansion Then it has a number Sum of rank The range of values ​​for is: , Number The corresponding total number of spherical harmonic functions for: No. The indices of the spherical harmonic functions are: For a specific incident direction in the discrete direction sample set of step 1 Calculate the first spherical harmonic function : (7) in, The zenith angle is the angle between the direction of the sun and the normal to the Earth's surface. The zenith angle is the angle between the sensor's observation direction and the Earth's surface normal, representing the reflected (observed) zenith angle. For each azimuth direction in the discrete direction sample set. As the normalization factor, In the formula For the accompanying Legendre polynomial ( The following step-by-step recursive algorithm is used to calculate the result: (1) Initial conditions (2) Calculate the diagonal polynomial The double factorial is defined as: (3) Recursive generation (4) General recursive relations generate higher-order ones. Thus, the incident direction is obtained. The real spherical harmonic function vector below: Further calculate the real spherical harmonic function vectors for all incident directions in the discrete direction sample set.

[0057] S130. Determine the directional weight function of each pixel in each discrete incident direction based on the terrain parameters.

[0058] According to the method for terrain radiometric correction of remote sensing images provided by the present invention, the method determines the directional weight function of each pixel in each discrete incident direction based on terrain parameters, including: selecting the corresponding weight calculation mode according to whether each discrete incident direction is the direction of direct sunlight and whether the current pixel is located in a shadow area; the weight calculation mode includes at least: direct sunlight and non-shadow mode, direct sunlight and shadow mode, non-direct sunlight and non-shadow mode, and non-direct sunlight and shadow mode; in the non-shadow mode, the weight function is jointly determined by the proportion of direct sunlight component and the sky scattering visibility factor; in the shadow mode, the weight function is jointly determined by the proportion of environmental reflectance irradiance, the neighboring reflectance visibility enhancement factor, and the neighboring reflectance ratio.

[0059] Specifically, step 4: Calculate the direction weight function.

[0060] For a given pixel in a remote sensing image, calculate a certain azimuth angle in a discrete direction sample set. Directional weighting function It comprehensively characterizes the accessibility of direct sunlight, the visibility of sky scattering, and the enhanced reflection effect of nearby terrain: (8) In the formula, The angle of solar incidence during imaging. Azimuth Maximum occlusion angle in the direction (obtained in step 1); , , The ratio of direct solar irradiance, sky scattering, and ambient reflected irradiance under the assumption of a horizontal and uniform Earth surface (obtained in step 1); and Azimuth angles Sky scattering visibility factor and proximity reflection visibility enhancement factor in the direction of reflection; Azimuth The reflectance ratio of adjacent directions.

[0061] Sky scattering visibility factor for: (9) Mapping the "level surface" assumption of 6S to actual undulating terrain Mapping the "uniform environment" assumption of 6S to the non-uniform reflection scenario of actual terrain: (10) Pick , Azimuth Within the visible area (from the target pixel to the occluded pixel) The average reflectance of ground features between (between) The uniform surface reflectance was set for the 6S simulation.

[0062] Further calculations yielded the directional weight function values ​​for each pixel in the remote sensing image in each azimuth direction.

[0063] This step addresses the challenge of quantifying the differential impact of terrain on illumination from different sources by establishing a refined directional weighting function calculation model. This function adaptively selects four weighting calculation modes based on whether the incident direction is direct sunlight and whether the pixel is in shadow, accurately characterizing the accessibility of direct sunlight, the visibility of sky scattering, and the reflection enhancement effect of adjacent terrain in each direction. Thus, the complex three-dimensional shading and reflection effects caused by terrain are transformed into a series of calculable physical parameters, laying a precise physical foundation for subsequently constructing a quantitative mapping relationship between radiative transfer and spherical harmonic functions.

[0064] S140. Calculate the corresponding directional spherical harmonic basis function based on the directional weight function of each pixel itself; calculate the corresponding terrain modulation spherical harmonic basis function based on the directional spherical harmonic basis function of each pixel in all discrete incident directions.

[0065] S150. Select representative pixels from both the shaded and unshaded areas as training samples.

[0066] According to the present invention, a method for topographic radiometric correction of remote sensing images is provided, which selects representative pixels as training samples in both shadowed and unshadowed areas. The method includes: in the shadowed area, generating a cumulative histogram based on the solar incidence angle of all pixels, removing extreme values, and dividing the remaining solar incidence angle range into several intervals; in each interval, selecting the pixel whose brightness is closest to the average brightness of the pixels in the current interval as the training sample of the current interval; and selecting training samples in the unshadowed area using the same principle as selecting training samples in the shadowed area.

[0067] Specifically, step 5: Select the training sample set.

[0068] Training sample sets are selected from the shadow and non-shadow regions of the remote sensing image to serve the subsequent calculation of the illumination coefficients of the shadow and non-shadow regions.

[0069] For the shadowed areas of the remote sensing image, a cumulative histogram is first generated for the solar incidence angles of all pixels within the area (calculated in step 2). Extreme values ​​below 1% and below 99% are removed, and the remaining solar incidence angle range is divided into several intervals at fixed intervals. If the number of pixels in a certain interval is lower than a set threshold, it is merged with neighboring intervals. For a specific band in the remote sensing image, the mean pixel brightness within each interval is calculated, and the pixel closest to the mean is selected as the training sample for that band in that interval. To ensure sufficient sampling of the illumination direction and the stability of the regression model, a certain number of training samples is required. .

[0070] Similarly, a training sample set of non-shaded regions is constructed, and the number of samples is denoted as . .

[0071] This step employs a training sample selection method based on the statistical distribution of solar incidence angles to ensure the representativeness and stability of the model training. Samples are selected within both shaded and unshaded regions by removing extreme values ​​and dividing the data into intervals, ensuring that the samples adequately cover the main differences in illumination conditions within each region. This sampling strategy not only guarantees the data quality for subsequent linear regression analysis but also avoids overfitting or underfitting problems caused by uneven sample distribution or insufficient representativeness, thereby improving the robustness and generalization ability of the obtained illumination coefficient vector.

[0072] S160. Based on the reflection characteristics of the training samples and the terrain-modulated spherical harmonic basis function, a linear model is constructed, and the spherical harmonic illumination coefficients characterizing the overall three-dimensional illumination distribution of the shadow and non-shadow areas are solved.

[0073] According to the present invention, a method for correcting the topographic radiometric properties of remote sensing images is provided. Based on the reflectance characteristics of training samples and the topographic modulation spherical harmonic basis function, a linear model is constructed, and the spherical harmonic illumination coefficients characterizing the overall three-dimensional illumination distribution of shadowed and unshadowed areas are solved. The method includes: for each training sample, calculating the incident radiance of each training sample based on the reflected radiance, and calculating the topographic modulation spherical harmonic basis function corresponding to each training sample based on the real spherical harmonic basis function vector; establishing an overdetermined linear equation system with incident radiance as the dependent variable and the topographic modulation spherical harmonic basis function as the independent variable, based on the training sample sets of shadowed and unshadowed areas respectively; solving the overdetermined linear equation system using the least squares method to obtain the spherical harmonic illumination coefficient vectors of the shadowed and unshadowed areas respectively.

[0074] According to the method for terrain radiometric correction of remote sensing images provided by the present invention, the terrain modulation spherical harmonic basis function corresponding to each training sample is calculated based on the real spherical harmonic basis function vector. The method includes: using the real spherical harmonic basis function vector and the directional weight function of each training sample in each discrete incident direction, calculating the directional spherical harmonic basis function of each training sample in each discrete direction; and taking a weighted average of the directional spherical harmonic basis functions of all incident directions for each training sample to obtain the terrain modulation spherical harmonic basis function corresponding to the training sample.

[0075] According to the method for topographic radiometric correction of remote sensing images provided by the present invention, the empirical reflectance of each training sample in the sensor observation direction is calculated based on the bidirectional reflectance distribution function parameters. The method includes: calculating the bidirectional reflectance of the training sample in the sensor observation direction and any discrete incident direction in the upper hemisphere based on the kernel-driven model weight parameters and the pre-calculated volume scattering kernel function matrix and geometric optics kernel function matrix; and performing a weighted average of the bidirectional reflectances in all discrete incident directions according to the solid angle to obtain the empirical reflectance of the training sample.

[0076] Step 6: Calculate the empirical reflectance of the training samples.

[0077] For a given training sample, based on its corresponding BRDF weight parameters , and (Obtained in Step 1), and the Ross Thick volume scattering kernel and Li-Sparse Reciprocal geometric optics scattering kernel matrix (Obtained in Step 1), calculate their incidence direction in a certain incident direction in the upper hemisphere. Sensor observation direction Diaxial reflectivity: (11) Where brdf is the bidirectional reflectance distribution function, and λ is the wavelength. The angle of incidence is the zenith angle. This refers to the reflected (observed) zenith angle. Relative azimuth angle. . For isotropic scattering weight parameters. For volume scattering weight parameters. This is the geometric optics scattering weighting parameter. For volume scattering kernel function. For geometric optics kernel functions.

[0078] Then, the empirical reflectance of the training samples is obtained by taking the solid angle-weighted average of the dihedral reflectances of all incident directions in the discrete direction sample set. : (12) in, This is the empirical reflectance. Here, n represents the zenith angle and azimuth angle in the direction of sensor observation. represents the number of discrete sampling points, indicating the total number of sampling points obtained after discretizing the zenith angle and azimuth angle of the upper hemisphere, respectively. a and b are summation indices; a iterates through the zenith angle sampling points, and b iterates through the azimuth angle sampling points. Let be the a-th discrete incident zenith angle. Solid angle element. , , The discretization interval for the zenith angle and azimuth angle is denoted as .

[0079] The empirical reflectance of each sample in the sensor observation direction is further obtained by using the above method to train the shadowed area training sample set and the unshaded area training sample set.

[0080] By combining a kernel-driven BRDF model with a kernel function matrix, the empirical reflectance of each training sample in the sensor observation direction is calculated. This calculation effectively integrates prior knowledge of non-Lambertian reflectance of ground features, extending the reflectance characteristics of a pixel in a specific observation direction to represent the average reflectance capability across all possible incident directions in the upper hemisphere. The resulting "empirical reflectance" links the reflected radiance observed by the sensor with the total incident radiance received by the pixel, providing a crucial conversion bridge for accurately retrieving the true incident radiance of the training samples.

[0081] Step 7: Calculate the incident radiance of the training samples.

[0082] For a given training sample, based on its reflected radiance in the sensor observation direction... (Step 1) and empirical reflectance Calculate the incident radiance of the training samples : (13) In the formula, It represents the weighted average of the radiance incident on this pixel from all directions in the upper hemisphere space; The reflectance radiance of the training samples.

[0083] The incident radiance of each sample in the training sample set of the shaded region and the training sample set of the unshaded region is calculated using the method described above.

[0084] This step uses the empirical reflectivity obtained in step 6 to invert the observed reflected radiance, thus obtaining the average incident radiance received by the training samples. This value characterizes the weighted sum of radiant energy actually received by a pixel from all directions in the upper hemisphere under the influence of complex terrain. It provides the most crucial "true value" label for subsequent training of the spherical harmonic illumination model, making it possible to fit the three-dimensional illumination distribution using mathematical methods. It is the core data node connecting physical observation and mathematical model.

[0085] Step 8: Calculate the terrain-modulated spherical harmonic basis functions of the training samples.

[0086] For a given training sample, calculate a certain incident direction. Directional spherical harmonic basis functions: (14) in, These are the directional spherical harmonic basis functions. This is the direction weighting function. is the vector of real spherical harmonic basis functions.

[0087] By taking a weighted average of the directional spherical harmonic basis functions for all incident directions, the terrain-modulated spherical harmonic basis function corresponding to this sample is obtained: (15) in, Modulate the spherical harmonic basis function vector for terrain. This is the tilt factor.

[0088] The terrain-modulated spherical harmonic basis function values ​​for each sample in the training sample set of the shaded region and the training sample set of the unshaded region are calculated using this method.

[0089] This step constructs a unique terrain-modulated spherical harmonic basis function for each training sample by combining the directional weight function with the real spherical harmonic basis function. This operation encodes the physical weights characterizing terrain lighting conditions into the mathematical basis function characterizing spatial lighting distribution, generating a "fingerprint" or "signature" of the terrain lighting characteristics for each sample point. This "signature" not only reflects the local terrain attributes of the sample point but is also expressed in the standard mathematical language of spherical harmonic functions, laying a formally unified and physically meaningful input foundation for subsequent use of linear models to uniformly solve for the overall regional lighting pattern.

[0090] Step 9: Calculate the spherical harmonic illumination coefficient.

[0091] For the training sample set in the shaded area A series of samples combine to form a linear system: (16) in, for 3D sample incident radiance vector, for 3D terrain-modulated spherical harmonic basis function matrix, for Vector of spherical harmonic illumination coefficients in the shaded region.

[0092] Next, the least squares method is used to solve for the spherical harmonic illumination coefficient vector in the shaded region: (17) Similarly, calculate the spherical harmonic illumination coefficient vector for the non-shaded area. .

[0093] This step uses the least squares method to fit the training sample sets for both shaded and unshaded areas, solving for their respective spherical harmonic illumination coefficient vectors. This process completes the extraction of a mathematical model representing the three-dimensional spatial distribution of illumination across the entire region (shaded or unshaded area) from the individual illumination characteristics of a large number of discrete sample points. The resulting coefficient vector is essentially a compact and powerful illumination distribution code, capable of efficiently and accurately describing the complex irradiance variation patterns caused by terrain in the region with a small number of parameters. This represents a key leap from point-to-area illumination modeling.

[0094] S170. Based on the topographic modulation spherical harmonic basis function of each pixel, combined with the spherical harmonic illumination coefficient of its corresponding region, the average incident radiance of each pixel is reconstructed.

[0095] According to a remote sensing image topographic radiometric correction method provided by the present invention, based on the topographic modulation spherical harmonic basis function of each pixel and combined with the spherical harmonic illumination coefficient corresponding to its region, the average incident radiance of each pixel is reconstructed. The method includes: for pixels in shaded areas, using the directional weight function of the pixel in each discrete incident direction to weight the real spherical harmonic basis function vector to obtain the directional spherical harmonic basis function in each incident direction, and then averaging the directional spherical harmonic basis function of the pixel in the shaded area by solid angle weighting; for pixels in non-shaded areas, using the directional weight function of the pixel in each discrete incident direction to weight the real spherical harmonic basis function vector to obtain the directional spherical harmonic basis function in each incident direction, and then averaging the directional spherical harmonic basis function of the pixel in the shaded area by solid angle weighting; The directional weighting function in the incident direction is used to weight the real spherical harmonic basis function vector to obtain the directional spherical harmonic basis function in each incident direction. The topographic modulation spherical harmonic basis function of the non-shaded area pixel is obtained by weighting the solid angle. The topographic modulation spherical harmonic basis function of the shaded pixel is multiplied by the spherical harmonic illumination coefficient vector of the shaded area to obtain the average incident radiance of each pixel in the shaded area. The topographic modulation spherical harmonic basis function of the non-shaded pixel is multiplied by the spherical harmonic illumination coefficient vector of the non-shaded area to obtain the average incident radiance of each pixel in the non-shaded area.

[0096] Specifically, step 10: reconstruct the average incident radiance image.

[0097] For a pixel in the shadow region of the remote sensing image, calculate a certain incident direction according to formula (14) in step 8. Directional spherical harmonic basis functions Then, according to formula (15) in step 8, the directional spherical harmonic basis functions of all incident directions are integrated and averaged to obtain the terrain modulation spherical harmonic basis function corresponding to the pixel. The incident radiance of the pixel in the shadow area for: (18) Similarly, the pixels in the non-shaded areas of a remotely sensed image are calculated. and incident radiance: (19) The incident radiance value of each pixel in the shadow and non-shadow regions of the remote sensing image is calculated separately, and the incident radiance image corresponding to each band of the remote sensing image is reconstructed.

[0098] This step applies the regional illumination coefficient model obtained in step 9 to every pixel in the entire image. By calculating the topographic modulation spherical harmonic basis function of each pixel and linearly combining it with the illumination coefficient of the corresponding region, the average incident radiance distribution map of the entire image is reconstructed. This process enables the learned regional illumination rules to be extended to a massive number of pixels, thereby efficiently and consistently estimating the total incident light energy received by each pixel under real terrain conditions. The generated average incident radiance image is the core intermediate product for eliminating the influence of terrain, and it intuitively shows the radiance differences caused by terrain in a visual form.

[0099] S180. Based on the reflected radiance observed by the sensor and the reconstructed average incident radiance, the true surface reflectance after eliminating the influence of topography is obtained by inversion.

[0100] According to the method for topographic radiometric correction of remote sensing images provided by the present invention, the true surface reflectance after eliminating the influence of topography is obtained by inverting the reflected radiance observed by the sensor and the reconstructed average incident radiance. The method includes: taking the reflected radiance observed by the sensor for each band of the remote sensing image to be corrected as the numerator; taking the reconstructed average incident radiance of the corresponding band as the denominator; performing ratio calculation pixel by pixel and band by band to obtain the topographically corrected surface reflectance of each band; and combining the topographically corrected surface reflectance of each band to form the true surface reflectance after eliminating the influence of topographic undulation.

[0101] Step 11: Invert the true reflectance image.

[0102] For a single pixel in the remote sensing image, the true surface reflectance after terrain correction is calculated based on the reflected radiance (step 1) from the sensor observation direction and the reconstructed average incident radiance (step 10): (20) In the formula, the true surface reflectance after terrain correction is... It effectively eliminates the changes in irradiance caused by topography and can more realistically reflect the inherent reflective properties of the earth's surface.

[0103] This step involves dividing the original reflected radiance observed by the sensor by the reconstructed average incident radiance, ultimately retrieving the true surface reflectance after eliminating the influence of topographic irradiance. This ratio calculation is physically equivalent to normalizing the reflected signal observed on a sloping surface to standard horizontal illumination conditions, thereby eliminating topographic-induced illumination variations and allowing the pixel values ​​in the corrected image to more purely reflect the inherent spectral characteristics of the land cover. The resulting reflectance product significantly improves the spectral consistency of similar land cover features across different slope aspects and gradients, enhancing the accuracy and reliability of the image for quantitative applications such as land cover classification and parameter inversion.

[0104] To verify the topographic radiometric correction effect of this invention, Landsat 8 OLI Level 1 remote sensing image (cropped to a 300 × 300 pixel mountain sub-region) was used as the data to be corrected. Its path number and row number were 123 / 32, and the imaging date was October 30, 2017. Six bands were selected from this image: blue (B), green (G), red (R), near-infrared (NIR), and two shortwave infrared bands (SWIR1 and SWIR2). Input data also included the ALOS Global Digital Surface Model (DSM) and MCD43A1 BRDF / Albedo data from October 30, 2017.

[0105] Figure 2 This is a schematic diagram of the overall technology for remote sensing image topographic radiometric correction according to an embodiment of the present invention.

[0106] Figure 3 This diagram illustrates the Landsat 8 OLI remote sensing image to be corrected, the corresponding ALOS DSM, and the shadow recognition results in an embodiment of the present invention. By comparing sub-images a and c, the correspondence between the shadow area and the dark area in the image can be clearly seen, demonstrating the accuracy of shadow discrimination.

[0107] Figure 4 This is a schematic diagram of the MCD43A1 and the BRDF weight parameter distribution after interpolation to the OLI scale (taking the blue band as an example) according to an embodiment of the present invention. By comparing the two lines, it can be intuitively seen that the interpolation process not only improves the spatial resolution of the data, but also maintains the original spatial distribution pattern and relative magnitude of the parameters, successfully transferring the prior BRDF knowledge to the high-resolution image to be corrected.

[0108] Figure 5a , Figure 5b This is a schematic diagram showing the band-wise display of the original reflectance, reconstructed average incident radiance, and resulting reflectance images according to an embodiment of the present invention. Through a horizontal comparison of the three columns in the same row, it can be clearly seen that: 1) the original image exhibits severe topographic effects; 2) the reconstructed average incident radiance image successfully depicts the differences in topographic illumination; 3) after correction, the resulting image effectively reduces the brightness and darkness of the topographic features, making the reflectance closer to the properties of the ground objects themselves. The band-wise display demonstrates the universality of this method across the multispectral range.

[0109] Figure 6This diagram illustrates a comparison of various color composite effects between the original and resulting reflectance images in an embodiment of the present invention. By directly comparing the up and down rows, it is possible to assess whether the corrected image effectively maintains the original color relationships and spectral characteristics of the ground features while improving brightness uniformity (i.e., avoiding color distortion). Successful correction should result in the second row image eliminating terrain shadows while presenting ground feature colors more realistically and naturally, and maintaining or enhancing the color distinction between different land cover types.

[0110] The following is in conjunction with the above appendix Figure 2-6 This describes the process of topographic radiometric correction for OLI remote sensing images.

[0111] Step 1: Perform radiometric calibration and FLAASH atmospheric correction on the OLI remote sensing image to obtain the original reflectance image; calculate the reflected radiance image based on the surface incident irradiance related parameters generated by 6S simulation; use MCD43A1 as the BRDF weight parameters for OLI remote sensing image pixel interpolation; construct a discrete direction sample set for the upper hemisphere, and calculate the Ross Thick volume scattering kernel function matrix and the Li-Sparse Reciprocal geometric optics kernel function matrix.

[0112] Step 2: Set the search radius to 30, calculate the maximum occlusion angle of the remote sensing image pixels in each azimuth direction and the pixel position that produces the maximum occlusion angle; and divide the remote sensing image into two regions: shaded and unshaded. Step 3: Set the number of belts Calculate the 1×25-dimensional real spherical harmonic function vector.

[0113] Step 4: Calculate the directional weight function value of each pixel in the remote sensing image in each azimuth direction.

[0114] Step 5: Based on the solar incidence angle, select training sample sets for the shaded and unshaded areas respectively.

[0115] Step 6: Calculate the empirical reflectance of the training samples based on the BRDF weight parameters and the volume scattering and geometric optics scattering kernel function matrices.

[0116] Step 7: Calculate the incident radiance of the training samples based on the reflected radiance (Step 1) and empirical reflectance of the training samples.

[0117] Step 8: Calculate the directional spherical harmonic basis functions of the training sample under all incident directions, and perform a weighted average by solid angle to obtain the terrain-modulated spherical harmonic basis function value of the sample.

[0118] Step 9: Based on the training sample sets of the shaded and unshaded areas respectively, use the least squares method to solve for the spherical harmonic illumination coefficient vectors of the shaded and unshaded areas.

[0119] Step 10: Calculate the topographic modulation spherical harmonic basis function corresponding to each pixel in the remote sensing image, and reconstruct the average incident radiance image by combining the spherical harmonic illumination coefficient vectors of the shadow area and the non-shadow area.

[0120] Step 11: Based on the reflected radiance image (Step 1) and the reconstructed average incident radiance image, calculate the topographically corrected true surface reflectance result image.

[0121] Compared with existing terrain correction methods, the technical method proposed in this invention comprehensively considers the accessibility of direct sunlight from different incident directions, the visibility of sky scattering, and the enhancement effect of reflection from adjacent terrain. It also reconstructs an intuitive image of the average incident radiance distribution, and the retrieved true surface reflectance image effectively reduces reflectance fluctuations caused by terrain undulations. The overall brightness consistency and uniformity are significantly improved. Its specific effects are mainly reflected in three aspects: (1) The minimum values ​​of each band in the corrected image are generally higher than those in the original image, indicating that the brightness of the shadow area has been effectively improved. The maximum values ​​of each band in the corrected image are lower than those in the original image, indicating that the bright areas have been moderately suppressed, avoiding brightness saturation caused by excessive reflection. In addition, the mean values ​​of each band in the corrected image are significantly higher than those in the original image, reflecting that the overall brightness level has been appropriately compensated. The standard deviation is smaller than that of the original image, indicating that the brightness difference and spatial fluctuation of the image have been effectively weakened, and the overall image is more uniform and stable. Furthermore, the information entropy of each band in the result image is higher than that in the original image, indicating that the spectral information of the shadow area has been restored to a certain extent.

[0122] (2) The corrected image maintains a high degree of consistency with the original image in overall tone, with no obvious color shift, indicating that the technique effectively maintains the spectral consistency of ground features while improving brightness. In addition, the overall brightness of the corrected image is improved, and the difference in brightness caused by terrain undulations is significantly reduced, resulting in a flatter and more balanced brightness distribution.

[0123] (3) For pixels on the sunny slope, the spectral curves before and after correction are quite similar, indicating that the correction of the sunny slope area basically maintains the original radiation characteristics. For pixels in the shaded area, the reflectance is significantly improved after correction, and its spectral curve is closer to that of the neighboring non-shaded pixels, indicating that the correction process effectively compensates for the brightness loss caused by terrain shading.

[0124] The remote sensing image topographic radiometric correction device provided by the present invention is described below. The remote sensing image topographic radiometric correction device described below and the remote sensing image topographic radiometric correction method described above can be referred to in correspondence with each other.

[0125] like Figure 7The image shown is a remote sensing image topographic radiometric correction device provided by the present invention, comprising: The terrain parameter calculation module 710 is used to acquire the remote sensing image to be corrected and its corresponding digital elevation model, and to calculate the terrain parameters of each pixel in the remote sensing image to be corrected based on the digital elevation model. The shadow region segmentation module 720 is used to divide the remote sensing image to be corrected into shadow regions and non-shadow regions according to the terrain parameters. The weight function determination module 730 is used to determine the directional weight function of each pixel in each discrete incident direction based on the terrain parameters. The terrain modulation spherical harmonic basis function determination module 740 is used to calculate the corresponding directional spherical harmonic basis function based on the directional weight function of each pixel itself; and to calculate the corresponding terrain modulation spherical harmonic basis function based on the directional spherical harmonic basis function of each pixel in all discrete incident directions. The training sample selection module 750 is used to select representative pixels as training samples in the shaded area and the non-shaded area respectively. The illumination coefficient solving module 760 is used to construct a linear model based on the reflection characteristics of the training samples and the terrain modulation spherical harmonic basis function, and solve for the spherical harmonic illumination coefficient that characterizes the overall three-dimensional illumination distribution of the shadow area and the non-shadow area. The average incident radiance reconstruction module 770 is used to reconstruct the average incident radiance of each pixel based on the topographic modulation spherical harmonic basis function of each pixel and the spherical harmonic illumination coefficient corresponding to its region. The true reflectance inversion module 780 is used to invert the true reflectance of the land surface after eliminating the influence of topography based on the reflected radiance observed by the sensor and the reconstructed average incident radiance.

[0126] Specifically, the functions of each module in the user account management system provided in this embodiment of the invention correspond one-to-one with the operation flow of each step in the above method-like embodiments, and the achieved effects are also the same. For details, please refer to the above embodiments, and this will not be repeated in this embodiment of the invention.

[0127] Figure 8 An example is a schematic diagram of the physical structure of an electronic device, such as... Figure 8As shown, the electronic device may include: a processor 810, a communication interface 820, a memory 830, and a communication bus 840, wherein the processor 810, the communication interface 820, and the memory 830 communicate with each other through the communication bus 840. The processor 810 can call logical instructions in the memory 830 to execute a remote sensing image topographic radiometric correction method. This method includes: acquiring the remote sensing image to be corrected and its corresponding digital elevation model, and calculating the topographic parameters of each pixel in the remote sensing image to be corrected based on the digital elevation model; dividing the remote sensing image to be corrected into shaded and unshaded areas according to the topographic parameters; determining the direction weight function of each pixel in each discrete incident direction according to the topographic parameters; calculating the corresponding directional spherical harmonic basis function based on the direction weight function of each pixel itself; and calculating the directional spherical harmonic basis function of each pixel in all discrete incident directions. The topographically modulated spherical harmonic basis function (THB) is calculated. Representative pixels are selected as training samples in both shaded and unshaded areas. Based on the reflection characteristics of the training samples and the THB, a linear model is constructed, and the spherical harmonic illumination coefficients characterizing the overall three-dimensional illumination distribution in both shaded and unshaded areas are solved. Based on the THB of each pixel and the corresponding spherical harmonic illumination coefficients of its region, the average incident radiance of each pixel is reconstructed. Based on the reflected radiance observed by the sensor and the reconstructed average incident radiance, the true surface reflectance after eliminating the influence of topography is obtained.

[0128] Furthermore, the logical instructions in the aforementioned memory 830 can be implemented as software functional units and, when sold or used as independent products, can be stored in a computer-readable storage medium. Based on this understanding, the technical solution of the present invention, in essence, or the part that contributes to the prior art, or a part of the technical solution, can be embodied in the form of a software product. This computer software product is stored in a storage medium and includes several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute all or part of the steps of the methods of the various embodiments of the present invention. The aforementioned storage medium includes various media capable of storing program code, such as USB flash drives, portable hard drives, read-only memory (ROM), random access memory (RAM), magnetic disks, or optical disks.

[0129] On the other hand, the present invention also provides a computer program product, which includes a computer program that can be stored on a non-transitory computer-readable storage medium. When the computer program is executed by a processor, the computer can execute the remote sensing image topographic radiometric correction method provided by the above methods. The method includes: acquiring the remote sensing image to be corrected and its corresponding digital elevation model, and calculating the topographic parameters of each pixel in the remote sensing image to be corrected based on the digital elevation model; dividing the remote sensing image to be corrected into shadowed and unshadowed regions according to the topographic parameters; determining the direction weight function of each pixel in each discrete incident direction according to the topographic parameters; and calculating the corresponding directionality based on the direction weight function of each pixel itself. Spherical harmonic basis functions (SHBF) are used to calculate the corresponding topographic modulated SHBF for each pixel in all discrete incident directions. Representative pixels are selected as training samples in both shaded and unshaded areas. A linear model is constructed based on the reflection characteristics of the training samples and the SHBF, and the spherical harmonic illumination coefficients characterizing the overall three-dimensional illumination distribution in both shaded and unshaded areas are solved. The average incident radiance of each pixel is reconstructed based on the SHBF for each pixel and the corresponding spherical harmonic illumination coefficients of its region. The true surface reflectance after eliminating the influence of topography is obtained by inverting the reflected radiance observed by the sensor and the reconstructed average incident radiance.

[0130] In another aspect, the present invention also provides a non-transitory computer-readable storage medium storing a computer program thereon, which, when executed by a processor, implements the remote sensing image topographic radiometric correction method provided by the methods described above. This method includes: acquiring a remote sensing image to be corrected and its corresponding digital elevation model, and calculating topographic parameters for each pixel in the remote sensing image to be corrected based on the digital elevation model; dividing the remote sensing image to be corrected into shaded and unshaded regions according to the topographic parameters; determining the direction weight function for each pixel in each discrete incident direction according to the topographic parameters; calculating the corresponding directional spherical harmonic basis function based on the direction weight function of each pixel; and calculating the direction weight function of each pixel in each discrete incident direction. The directional spherical harmonic basis functions of each pixel in all discrete incident directions are used to calculate the corresponding topographic modulated spherical harmonic basis functions. Representative pixels are selected as training samples in both shaded and unshaded areas. Based on the reflection characteristics of the training samples and the topographic modulated spherical harmonic basis functions, a linear model is constructed, and the spherical harmonic illumination coefficients characterizing the overall three-dimensional illumination distribution in both shaded and unshaded areas are solved. Based on the topographic modulated spherical harmonic basis functions of each pixel and the spherical harmonic illumination coefficients corresponding to its region, the average incident radiance of each pixel is reconstructed. Based on the reflected radiance observed by the sensor and the reconstructed average incident radiance, the true surface reflectance after eliminating the influence of topography is obtained by inversion.

[0131] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units; that is, they may be located in one place or distributed across multiple network units. Some or all of the modules can be selected to achieve the purpose of this embodiment according to actual needs. Those skilled in the art can understand and implement this without any creative effort.

[0132] Through the above description of the embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus necessary general-purpose hardware platforms, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solutions, in essence or the part that contributes to the prior art, can be embodied in the form of a software product. This computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, magnetic disk, optical disk, etc., including several instructions to cause a computer device (which may be a personal computer, server, or network device, etc.) to execute the methods of various embodiments or some parts of embodiments.

[0133] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.

Claims

1. A method for topographic radiometric correction of remote sensing images, characterized in that, include: Acquire the remote sensing image to be corrected and its corresponding digital elevation model, and calculate the terrain parameters of each pixel in the remote sensing image to be corrected based on the digital elevation model; The remote sensing image to be corrected is divided into shadowed and unshadowed areas based on the terrain parameters. The orientation weight function of each pixel in each discrete incident direction is determined based on the terrain parameters; Based on the directional weight function of each pixel itself, calculate its corresponding directional spherical harmonic basis function; Based on the directional spherical harmonic basis function of each pixel in all discrete incident directions, calculate its corresponding terrain-modulated spherical harmonic basis function. Representative pixels were selected as training samples from both the shaded and unshaded areas. Based on the reflection characteristics of the training samples and the terrain-modulated spherical harmonic basis function, a linear model is constructed, and the spherical harmonic illumination coefficients characterizing the overall three-dimensional illumination distribution of the shadowed and unshadowed areas are solved. Based on the topographically modulated spherical harmonic basis function of each pixel, and combined with the spherical harmonic illumination coefficient corresponding to its region, the average incident radiance of each pixel is reconstructed. Based on the reflected radiance observed by the sensor and the reconstructed average incident radiance, the true surface reflectance after eliminating the influence of topography is obtained by inversion.

2. The remote sensing image topographic radiometric correction method according to claim 1, characterized in that, The calculation of terrain parameters for each pixel in the remote sensing image to be corrected based on the digital elevation model includes: The surface normal vector of each pixel in the remote sensing image to be corrected is determined based on the digital elevation model. The terrain parameters for each pixel are calculated based on the surface normal vector. The terrain parameters include at least the slope, aspect, solar incidence angle, and maximum shading angle in multiple discrete azimuth directions.

3. The method for topographic radiometric correction of remote sensing images according to claim 1, characterized in that, The directional weight function for each pixel in each discrete incident direction is determined based on the terrain parameters, including: Based on whether each discrete incident direction is the direction of direct sunlight and whether the current pixel is in the shadow area, the corresponding weight calculation mode is selected; the weight calculation mode includes at least: direct sunlight and non-shadow mode, direct sunlight and shadow mode, non-direct sunlight and non-shadow mode, and non-direct sunlight and shadow mode. In non-shaded mode, the weighting function is jointly determined by the proportion of direct sunlight and the sky scattering visibility factor; In shadow mode, the weighting function is determined by the proportion of ambient reflectance, the neighbor reflectance visibility enhancement factor, and the neighbor reflectance ratio.

4. The method for topographic radiometric correction of remote sensing images according to claim 1, characterized in that, Representative pixels were selected as training samples from both the shaded and unshaded areas, including: Within the shaded area, a cumulative histogram is generated based on the solar incidence angle of all pixels. After removing extreme values, the remaining range of solar incidence angles is divided into several intervals. Within each interval, select the pixel whose brightness is closest to the average brightness of the pixels in the current interval as the training sample for the current interval. Training samples are selected from the non-shaded areas using the same principle as those selected from the shaded areas.

5. The method for topographic radiometric correction of remote sensing images according to claim 1, characterized in that, The linear model is constructed based on the reflection characteristics of the training samples and the terrain-modulated spherical harmonic basis function, and the spherical harmonic illumination coefficients characterizing the overall three-dimensional illumination distribution of the shadowed and unshadowed areas are solved, including: Radiometric calibration and atmospheric correction are performed on the remote sensing images to be corrected to obtain preliminary surface reflectance. Based on the radiative transfer model, the reflected radiance in the sensor observation direction is calculated according to the preliminary surface reflectance; Construct and compute a set of real spherical harmonic basis function vectors; For each training sample, the incident radiance of each training sample is calculated based on the reflected radiance, and the terrain-modulated spherical harmonic basis function corresponding to each training sample is calculated based on the real spherical harmonic basis function vector. Based on the training sample sets of the shaded and unshaded regions respectively, an overdetermined linear equation system is established with incident radiance as the dependent variable and the terrain-modulated spherical harmonic basis function as the independent variable. The overdetermined linear equations are solved using the least squares method to obtain the spherical harmonic illumination coefficient vectors of the shaded region and the unshaded region, respectively.

6. The method for topographic radiometric correction of remote sensing images according to claim 5, characterized in that, The step of calculating the incident radiance of each training sample based on the reflected radiance, and calculating the terrain-modulated spherical harmonic basis function corresponding to each training sample based on the real spherical harmonic basis function vector, includes: The real spherical harmonic basis function vector and the direction weight function of each training sample in each discrete incident direction are used; The directional spherical harmonic basis functions of each training sample in each discrete direction are calculated based on the directional weight function. For each training sample, the directional spherical harmonic basis functions of all incident directions are weighted and averaged to obtain the terrain-modulated spherical harmonic basis function corresponding to the training sample.

7. The method for topographic radiometric correction of remote sensing images according to claim 6, characterized in that, The empirical reflectance of each training sample in the sensor observation direction is calculated based on the bidirectional reflectance distribution function parameters, including: Based on the weight parameters of the kernel-driven model, and the pre-calculated volume scattering kernel function matrix and geometric optics kernel function matrix, the biaxial reflectivity of the training sample is calculated in the sensor observation direction and in any discrete incident direction in the upper hemisphere. The empirical reflectance of the training samples is obtained by weighting the dihedral reflectances of all discrete incident directions according to the solid angle.

8. The method for topographic radiometric correction of remote sensing images according to claim 1, characterized in that, The average incident radiance of each pixel is reconstructed by combining the terrain-modulated spherical harmonic basis function of each pixel with the spherical harmonic illumination coefficient of the region to which it belongs, including: For the pixels in the shadow region, the real spherical harmonic basis function vector is weighted using the directional weight function in each discrete incident direction to obtain the directional spherical harmonic basis function in each incident direction, and the terrain modulation spherical harmonic basis function of the pixel in the shadow region is obtained by weighted averaging according to the solid angle. For pixels in the non-shaded area, the real spherical harmonic basis function vector is weighted using the directional weight function in each discrete incident direction to obtain the directional spherical harmonic basis function in each incident direction, and the terrain modulation spherical harmonic basis function of the non-shaded area pixel is obtained by weighted averaging according to solid angle. The average incident radiance of each pixel in the shadow region is obtained by performing an inner product operation between the terrain modulation spherical harmonic basis function of the shadow pixel and the spherical harmonic illumination coefficient vector of the shadow region. The average incident radiance of each pixel in the non-shaded region is obtained by performing an inner product operation between the topographic modulation spherical harmonic basis function of the non-shaded pixel and the spherical harmonic illumination coefficient vector of the non-shaded region.

9. The method for topographic radiometric correction of remote sensing images according to claim 1, characterized in that, The process of retrieving the true surface reflectance after eliminating topographic effects based on the reflected radiance observed by the sensor and the reconstructed average incident radiance includes: The sensor-observed reflectance radiance corresponding to each band of the remote sensing image to be corrected is used as the numerator; Use the reconstructed average incident radiance of the corresponding band as the denominator; The ratio calculation is performed pixel by pixel and band by band to obtain the topographically corrected surface reflectance of each band; The surface reflectance after topographic correction for each band is combined to form the true surface reflectance after eliminating the influence of topographic undulations.

10. A remote sensing image topographic radiometric correction device, characterized in that, include: The terrain parameter calculation module is used to acquire the remote sensing image to be corrected and its corresponding digital elevation model, and to calculate the terrain parameters of each pixel in the remote sensing image to be corrected based on the digital elevation model. The shadow region segmentation module is used to divide the remote sensing image to be corrected into shadow regions and non-shadow regions according to the terrain parameters. The weighting function determination module is used to determine the directional weighting function of each pixel in each discrete incident direction based on the terrain parameters. The terrain-modulated spherical harmonic basis function determination module is used to calculate the corresponding directional spherical harmonic basis function based on the directional weight function of each pixel itself; and to calculate the corresponding terrain-modulated spherical harmonic basis function based on the directional spherical harmonic basis function of each pixel in all discrete incident directions. The training sample selection module is used to select representative pixels as training samples in the shaded area and the non-shaded area respectively. The illumination coefficient solution module is used to construct a linear model based on the reflection characteristics of the training samples and the terrain modulation spherical harmonic basis function, and solve for the spherical harmonic illumination coefficients that characterize the overall three-dimensional illumination distribution of the shadow area and the non-shadow area. The average incident radiance reconstruction module is used to reconstruct the average incident radiance of each pixel based on the terrain-modulated spherical harmonic basis function of each pixel and the spherical harmonic illumination coefficient corresponding to its region. The true reflectance inversion module is used to invert the true surface reflectance after eliminating the influence of topography based on the reflected radiance observed by the sensor and the reconstructed average incident radiance.

Images