CN-121998877-A - Remote sensing image terrain radiation correction method, device and equipment

CN121998877ACN 121998877 ACN121998877 ACN 121998877ACN-121998877-A

Abstract

The invention relates to the field of remote sensing image processing and provides a remote sensing image topographic radiation correction method, device and equipment, wherein the method comprises the steps of obtaining a remote sensing image to be corrected, calculating topographic parameters of each pixel, dividing the remote sensing image to be corrected into a shadow area and a non-shadow area, and determining a direction weight function of each pixel; the method comprises the steps of calculating a terrain modulation spherical harmonic basis function corresponding to each pixel, selecting training samples in a shadow area and a non-shadow area respectively, constructing a linear model, solving spherical harmonic illumination coefficients representing overall three-dimensional illumination distribution of the shadow area and the non-shadow area, reconstructing average incident radiance of each pixel by combining the spherical harmonic illumination coefficients corresponding to the areas of each pixel based on the terrain modulation spherical harmonic basis function, and inverting to obtain the ground surface true reflectivity. The invention solves the problem of inaccurate correction of remote sensing image radiation distortion under complex terrain in the prior art, and realizes accurate inversion of true reflectivity of the earth surface.

Inventors

WANG QIAN
ZHANG YANJUN
WANG CHENG
HAN QI

Assignees

天津师范大学

Dates

Publication Date: 20260508
Application Date: 20251230

Claims (10)

1. A method for correcting terrain radiation of a remote sensing image, comprising: acquiring a remote sensing image to be corrected and a digital elevation model corresponding to the remote sensing image to be corrected, and calculating the topography parameter 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 a shadow area and a non-shadow area according to the topographic parameters; Determining a direction weight function of each pixel in each discrete incidence direction according to the terrain parameters; Calculating a corresponding directional spherical harmonic basis function based on the direction weight function of each pixel; calculating a corresponding terrain modulation spherical harmonic basis function based on the directional spherical harmonic basis functions of each pixel in all discrete incidence directions; respectively selecting representative pixels in the shadow area and the non-shadow area as training samples; Constructing a linear model based on the reflection characteristic of the training sample and the terrain modulation spherical harmonic basis function, and solving a spherical harmonic illumination coefficient representing the overall three-dimensional illumination distribution of the shadow region and the non-shadow region; Reconstructing average incident radiance of each pixel by combining the spherical harmonic illumination coefficients corresponding to the areas of the spherical harmonic fundamental functions based on the terrain modulation spherical harmonic fundamental functions of each pixel; And inverting to obtain the ground surface true reflectivity after eliminating the terrain influence according to the reflected radiance observed by the sensor and the reconstructed average incident radiance.
2. The method according to claim 1, wherein calculating the topographic parameters of each pixel in the remote sensing image to be corrected based on the digital elevation model comprises: determining the earth 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 according to the surface normal vector, wherein the topographic parameters at least comprise gradient, slope direction, sun incidence angle and maximum shielding angle in a plurality of discrete azimuth directions.
3. The method of claim 1, wherein determining a direction weight function for each discrete direction of incidence for each pixel based on the topographical parameters comprises: Selecting a corresponding weight calculation mode according to whether each discrete incidence direction is a direct solar radiation direction and whether a current pixel is positioned in a shadow area, wherein the weight calculation mode at least comprises a direct solar radiation and shadow mode, a non-direct solar radiation and shadow mode and a non-direct solar radiation and shadow mode; In the non-shadow mode, the weight function is determined by the direct solar component proportion and the sky scattering visibility factor; in shadow mode, the weighting function is determined by a combination of the proportion of ambient reflected irradiance, the adjacent reflected visibility enhancement factor, and the adjacent reflectance ratio.
4. The method of claim 1, wherein selecting representative pixels as training samples in the shadow and non-shadow areas, respectively, comprises: in the shadow area, generating an accumulated histogram based on the sun incidence angles of all pixels, and dividing the residual sun incidence angle range into a plurality of intervals after eliminating extreme values; In each interval, selecting a pixel closest to the brightness average value of the pixels in the current interval as a training sample of the current interval; And selecting the training samples in the non-shadow area by adopting the same principle as that of selecting the training samples in the shadow area.
5. The method of claim 1, wherein constructing a linear model based on the reflection characteristics of the training samples and the terrain modulation spherical harmonic basis functions and solving for spherical harmonic illumination coefficients characterizing the overall three-dimensional illumination distribution of the shadow region and the non-shadow region comprises: Performing radiometric calibration and atmospheric correction on the remote sensing image to be corrected to obtain a preliminary earth surface reflectivity; Calculating to obtain the reflected radiance in the observation direction of the sensor according to the preliminary surface reflectivity based on a radiation transmission model; constructing and calculating a group of real spherical harmonic basis function vectors; For each training sample, calculating the incident radiance of each training sample according to the reflected radiance, and calculating a terrain modulation spherical harmonic basis function corresponding to each training sample according to the real spherical harmonic basis function vector; respectively based on training sample sets of the shadow area and the non-shadow area, establishing an overdetermined linear equation set taking incident radiance as a dependent variable and taking the terrain modulation spherical harmonic basis function as an independent variable; And solving the overdetermined linear equation set by adopting a least square method to respectively obtain a spherical harmonic illumination coefficient vector of the shadow region and a spherical harmonic illumination coefficient vector of the non-shadow region.
6. The method of claim 5, wherein calculating the incident radiance of each training sample from the reflected radiance and calculating the terrain-modulated spherical harmonic basis function corresponding to each training sample from the real spherical harmonic basis function vector comprises: Utilizing the real spherical harmonic basis function vector and the direction weight function of each training sample in each discrete incident direction; Calculating a directional spherical harmonic basis function of each training sample in each discrete direction according to the direction weight function; And carrying out weighted average on directional spherical harmonic basis functions of all incidence directions of each training sample to obtain a terrain modulation spherical harmonic basis function corresponding to the training sample.
7. The method of claim 6, wherein calculating the empirical reflectivity of each training sample in the direction of sensor observation based on the bi-directional reflectance distribution function parameters comprises: Based on the weight parameters of the nuclear drive model, a pre-calculated volume scattering kernel function matrix and a geometrical optics kernel function matrix, calculating the two-way reflectivity of the training sample in any discrete incidence direction of the sensor observation direction and the upper hemisphere; and (3) carrying out weighted average on the two-way reflectivity in all the discrete incidence directions according to the solid angle to obtain the empirical reflectivity of the training sample.
8. The method according to claim 1, wherein reconstructing the average incident radiance of each pixel based on the terrain modulation spherical harmonic basis function of each pixel in combination with the spherical harmonic illumination coefficient corresponding to the region to which the spherical harmonic basis function belongs comprises: For the pixels of the shadow area, weighting the real spherical harmonic basis function vector by using the direction weight function of the pixels in each discrete incidence direction to obtain directional spherical harmonic basis functions of each incidence direction, and obtaining the terrain modulation spherical harmonic basis functions of the pixels of the shadow area according to solid angle weighted average; For the pixels of the non-shadow area, weighting the real spherical harmonic basis function vector by using the direction weight function of the pixels in each discrete incidence direction to obtain directional spherical harmonic basis functions of each incidence direction, and obtaining the terrain modulation spherical harmonic basis functions of the pixels of the non-shadow area according to solid angle weighted average; Performing inner product operation on the spherical harmonic basis function of the shadow pixel topography modulation and the spherical harmonic illumination coefficient vector of the shadow region to obtain average incident radiance of each pixel of the shadow region; And carrying out inner product operation on the spherical harmonic basis function of the non-shadow pixel terrain modulation and the spherical harmonic illumination coefficient vector of the non-shadow region to obtain the average incident radiance of each pixel of the non-shadow region.
9. The method according to claim 1, wherein the inversion of the reflected radiance and the reconstructed average incident radiance according to the sensor observation to obtain the true reflectivity of the earth surface after eliminating the influence of the terrain comprises: Taking the sensor observation reflection radiance corresponding to each wave band of the remote sensing image to be corrected as a molecule; taking the average incident radiance reconstructed in the corresponding wave band as a denominator; performing ratio operation pixel by pixel and wave band by wave band to obtain the surface reflectivity after terrain correction of each wave band; And combining the earth surface reflectivities after the topography correction of each wave band to form the earth surface true reflectivities after the relief influence of the topography is eliminated.
10. A remote sensing image terrain radiation correction device, comprising: The terrain parameter calculation module is used for acquiring a remote sensing image to be corrected and a digital elevation model corresponding to the remote sensing image to be corrected, and calculating the terrain parameter of each pixel in the remote sensing image to be corrected based on the digital elevation model; The shadow region dividing module is used for dividing the remote sensing image to be corrected into a shadow region and a non-shadow region according to the terrain parameter; the weight function determining module is used for determining a direction weight function of each pixel in each discrete incidence direction according to the terrain parameter; The terrain modulation spherical harmonic basis function determining module is used for calculating a corresponding directional spherical harmonic basis function based on the direction weight function of each pixel, and calculating a corresponding terrain modulation spherical harmonic basis function based on the directional spherical harmonic basis function of each pixel in all discrete incidence directions; The training sample selection module is used for selecting representative pixels in the shadow area and the non-shadow area respectively as training samples; the illumination coefficient solving module is used for constructing a linear model based on the reflection characteristic of the training sample and the terrain modulation spherical harmonic basis function and solving a spherical harmonic illumination coefficient representing the integral three-dimensional illumination distribution of the shadow area and the non-shadow area; the average incident radiance reconstruction module is used for reconstructing the average incident radiance of each pixel by combining the spherical harmonic illumination coefficients corresponding to the area to which the spherical harmonic fundamental function belongs based on the terrain modulation spherical harmonic fundamental function of each pixel; and the real reflectivity inversion module is used for inverting to obtain the ground surface real reflectivity after eliminating the topographic influence according to the reflected radiance observed by the sensor and the reconstructed average incident radiance.

Description

Remote sensing image terrain radiation correction method, device and equipment Technical Field The invention relates to the field of remote sensing image processing, in particular to a remote sensing image topographic radiation correction method, device and equipment. Background The remote sensing technology becomes a core support of a modern earth observation system, and provides an irreplaceable technical means for acquiring earth surface information in a large range and multiple phases. Particularly in mountainous areas and hilly areas with complex terrains, the remote sensing image can provide key data such as surface coverage, vegetation states, environmental changes and the like, and supports quantitative application in various fields such as resource exploration, ecological monitoring, disaster assessment and the like. However, due to the influence of the topography fluctuation, the surface radiation signal received by the remote sensing sensor can generate obvious geometric and radiation distortion, and the distortion directly restricts the quantitative application precision of the remote sensing data in the complex topography area. To eliminate the terrain effects, various terrain radiation correction methods have been developed. Early solutions, such as band-ratio methods and empirical statistical models, often relied on specific assumptions (such as proportional to the change in reflectance of different bands) or required flat area references, were less ubiquitous in complex terrain and were vulnerable to spectral information. Classical models based on lambertian assumptions are widely used, but often result in undercorrected or overcorrected shadow areas due to neglecting the non-lambertian reflection properties of the earth's surface and the incident-observation geometry. The subsequent Minnaert model introduces an empirical index to simulate the non-lambertian effect, but the physical meaning of the parameter is not clear, and the dichroism reflection characteristic under the complex illumination condition is difficult to truly describe. Some of the more advanced models began to consider radiation component decomposition or adjacent terrain reflection, but the contribution to sky scattering or adjacent reflection still had a simplified or uniform processing approach that failed to accurately quantify the difference in radiation in different directions of incidence. Therefore, how to construct a terrain radiation correction method which can strictly follow a physical radiation transmission mechanism and efficiently process multi-directional illumination effects under complex terrains so as to accurately recover the true spectral characteristics of the earth surface becomes a key technical bottleneck for improving the quantitative application level of remote sensing data in mountain areas, and is an important subject to be solved in the current remote sensing pretreatment field. Disclosure of Invention The invention provides a remote sensing image terrain radiation correction method, device and equipment, which solve the problem of inaccurate remote sensing image radiation distortion correction under complex terrain in the prior art and realize accurate inversion of true reflectivity of the earth surface. The invention provides a remote sensing image topographic radiation correction method, which comprises the following steps: acquiring a remote sensing image to be corrected and a digital elevation model corresponding to the remote sensing image to be corrected, and calculating the topography parameter 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 a shadow area and a non-shadow area according to the topographic parameters; Determining a direction weight function of each pixel in each discrete incidence direction according to the terrain parameters; Calculating a corresponding directional spherical harmonic basis function based on the direction weight function of each pixel; calculating a corresponding terrain modulation spherical harmonic basis function based on the directional spherical harmonic basis functions of each pixel in all discrete incidence directions; respectively selecting representative pixels in the shadow area and the non-shadow area as training samples; Constructing a linear model based on the reflection characteristic of the training sample and the terrain modulation spherical harmonic basis function, and solving a spherical harmonic illumination coefficient representing the overall three-dimensional illumination distribution of the shadow region and the non-shadow region; Reconstructing average incident radiance of each pixel by combining the spherical harmonic illumination coefficients corresponding to the areas of the spherical harmonic fundamental functions based on the terrain modulation spherical harmonic fundamental functions of each pixel; And inverting to obtain the