A method for joint correction of remote sensing images of atmosphere, terrain and geometry based on radiation transfer model

Abstract

The present invention provides a remote sensing image atmospheric terrain geometry joint correction method based on a radiation transfer model, which belongs to the field of remote sensing image processing technology. By establishing a ground-air coupled radiation transfer model, and breaking through the high-precision digital surface model and high-precision automatic registration technology of hyperspectral images, along the sun-undulating surface target and background-sensor geometric path, high-precision iterative calculation of components such as direct sunlight, sky scattering, ground-air coupled radiation, and surrounding background reflection is performed, and then the reflection characteristics of regular grid ground pixels are reconstructed, and finally a high-precision hyperspectral surface reflectivity with integrated radiation and geometry correction is obtained. The present invention effectively overcomes the problems of over-correction and under-correction.

Description

Technical Field

[0001] The present invention belongs to the technical field of remote sensing image processing, and specifically relates to a remote sensing image atmospheric terrain geometry joint correction method based on a radiation transfer model. The present invention can be used in an airborne or satellite-borne data preprocessing process, thereby serving the application of data acquired by a remote sensor, and is a supporting technology in the field of earth observation. Background Art

[0002] Correction of atmospheric and terrain effects is a key step in the remote sensing data preprocessing chain. With the improvement of satellite sensor performance, high-precision preprocessing methods can meet the new applications of remote sensing in many fields. Hyperspectral imaging is affected by Rayleigh scattering and Mie scattering of atmospheric molecules and aerosols. Terrain features also affect the propagation, reflection, absorption and scattering of radiation. In this context, atmospheric correction algorithms based on radiation transfer models have gradually replaced empirical methods of atmospheric correction. In mountainous environments, images must also be corrected for terrain effects. How to eliminate interference information in hyperspectral data and accurately invert the true reflectivity of objects is a key issue in the quantification of hyperspectral remote sensing applications.

[0003] At present, related research work at home and abroad is mainly based on two ideas: one is to first use terrain correction, and then perform atmospheric correction on the radiance data that eliminates the terrain effect; the other is the idea of ​​joint atmospheric-terrain correction. In previous studies, terrain correction and atmospheric correction were mostly performed separately.

[0004] Although terrain correction and atmospheric correction are performed successively, it is the current routine operation process. Considering them as two independent processes cannot truly reflect the radiation transmission mechanism of the coupling between the atmosphere and the surface, and it is easy to over-correct or under-correct the reflectivity in some areas. Therefore, the reflectivity inversion of ground objects is developing towards an integrated model trend of atmosphere-topography joint correction. According to the idea of ​​atmosphere-topography joint correction, although some methods have been used to establish an atmosphere-topography coupling model, the model construction is not detailed enough. Because the precise registration of remote sensing images and terrain models is not considered, the calculation of some radiation is not accurate enough, and the estimation of the initial reflectivity value is too rough, resulting in under-correction and over-correction of some undulating areas, and it is impossible to obtain high-precision reflectivity results. Summary of the invention

[0005] In order to solve the problem of high-precision inversion of surface reflectance caused by complex lighting conditions in undulating areas, the present invention proposes a remote sensing image atmospheric terrain geometry joint correction method based on a radiation transfer model. By establishing a ground-atmosphere coupled radiation transfer model and breaking through the high-precision digital surface model and high-precision automatic registration technology of hyperspectral images, along the sun-undulating surface target and background-sensor geometric path, the atmospheric terrain radiation transfer equation is combined with the medium spectral resolution atmospheric transmittance calculation model (MODTRAN), and components such as direct sunlight, sky scattering, ground-atmosphere coupled radiation and surrounding background reflection are iteratively calculated with high precision, and then the reflection characteristics of regular grid ground pixels are reconstructed, and finally the hyperspectral surface reflectance with integrated radiation and geometry correction is obtained.

[0006] In order to achieve the above object, the present invention adopts the following technical scheme:

[0007] A method for joint correction of remote sensing image atmosphere and terrain geometry based on a radiation transfer model comprises the following steps:

[0008] Step 1: Obtain a hyperspectral remote sensing radiance image and a digital elevation model / digital surface model, and obtain DEM / DSM data as input data; the digital elevation model is represented by DEM, and the digital surface model is represented by DSM;

[0009] Step 2: Ignore the terrain undulations, consider the contribution of the target and background to the radiance entering the sensor, and solve for the initial surface reflectance ;

[0010] Step 3: Use the rational polynomial coefficient model as the basic method of geometric correction to generate the remote sensing image input geometry file, realize the geographic positioning of the remote sensing image, and realize the automatic registration of DEM / DSM data and remote sensing images to obtain the latitude, longitude and elevation values ​​of each pixel;

[0011] Step 4: According to the DEM / DSM data obtained in step 1, the slope, aspect, sky visibility factor, terrain visibility factor and shadow coefficient are calculated, and they are sampled to each pixel position of the radiance image to accurately extract the terrain parameters of each pixel;

[0012] Step 5: Set the imaging time, pixel location and elevation, and atmospheric parameters, and use the radiation transfer model MODTRAN to calculate the direct solar radiation. , atmospheric scattered radiation , Program Radiation And transmittance ,in is the row number where the pixel is located, The column number where the pixel is located, is the ground elevation at the pixel location, is the wavelength;

[0013] Step 6: Based on the initial reflectivity result of step 2 , calculate the reflection of surrounding objects , Earth-atmosphere coupled radiation and atmospheric proximity effect , and calculate the intermediate reflectivity value ;

[0014] Step 7: Calculate the reflectivity of step 6 repeatedly for a set number of times until the final reflectivity result converges, thereby obtaining the final high-precision surface reflectivity result under the undulating surface.

[0015] Furthermore, in step 2, the initial surface reflectivity The specific calculation method is:

[0016]

[0017] in, L is the entrance pupil radiance, is the target pixel reflectance, is the average reflectivity of the background objects, is the downward hemispheric albedo of the atmosphere, For the process radiation, and The entrance pupil radiance is composed of three components, namely the part of the target pixel reflection that directly enters the sensor, the part of the radiation reflected by the background objects that enters the sensor, and the path radiation. The entrance pupil radiance is composed of three components, namely the part of the target pixel reflection that directly enters the sensor, the part of the radiation reflected by the background objects that enters the sensor, and the path radiation.

[0018] Furthermore, the specific calculation method of the terrain parameters in step 4 is as follows:

[0019]

[0020]

[0021] in, are the slope and the aspect, are the rates of elevation change in the north-south and east-west directions, respectively;

[0022]

[0023] in, They are the sky visibility factor and the terrain visibility factor, For this location The maximum zenith angle in the direction, The value of is between 0° and 360°; the projected shadow is calculated by the depth buffer algorithm in computer graphics, and a decimal between 0 and 1 is used to represent the proportion of the shadow area. When the coefficient is 0, the pixel is in the shadow; when the coefficient is 1, the pixel is in the sunlight area.

[0024] Furthermore, in step 6, the intermediate reflectivity value is calculated The method is as follows:

[0025]

[0026] Furthermore, in step 7, the reflectivity intermediate value As the new initial value, the reflectivity is iteratively calculated by setting the number of times until the final reflectivity result converges, that is, the final high-precision surface reflectivity result under the undulating surface is obtained.

[0027] Beneficial effects achieved by the present invention:

[0028] (1) The influence of terrain undulation and complex lighting conditions on reflectivity inversion was considered, and a ground-atmosphere coupled radiation transfer model was established. Stereo positioning was performed through the RPC model. A breakthrough was made in the high-precision automatic registration technology of DSM / DEM data and hyperspectral images, and a joint correction method for atmospheric terrain geometry was proposed.

[0029] (2) The initial reflectivity value is solved with the help of a flat surface model, and then each radiation component is iteratively calculated to finally obtain a high-precision surface reflectivity product, effectively overcoming the problems of over-correction and under-correction. BRIEF DESCRIPTION OF THE DRAWINGS

[0030] Figure 1 It is a flow chart of a method for joint correction of remote sensing image atmosphere and terrain geometry based on a radiation transfer model of the present invention;

[0031] Figure 2 Schematic diagram of the radiation amount of the sun-target-sensor path;

[0032] Figure 3a , Figure 3b This is a comparison chart before and after the joint correction of atmosphere and terrain geometry; Figure 3a is the GF-5 radiance image, Figure 3b This is the result after joint correction of atmosphere and terrain geometry. DETAILED DESCRIPTION

[0033] like Figure 1 As shown, the method for joint correction of remote sensing image atmosphere terrain geometry based on radiation transfer model of the present invention comprises the following steps:

[0034] Step 1: Obtain hyperspectral remote sensing radiance images and digital elevation models (DEM) / digital surface models (DSM) as input data;

[0035] Step 2: Ignore the terrain undulations, consider the contribution of the target and background to the radiance entering the sensor, and solve for the initial surface reflectance , using the following formula:

[0036]

[0037] in, L is the entrance pupil radiance, is the target pixel reflectance, is the average reflectivity of the background objects, is the downward hemispheric albedo of the atmosphere, For the process radiation, and As coefficient, it is considered that the pupil radiance consists of three components, namely, the part reflected by the target pixel directly entering the sensor, the part reflected by the background object entering the sensor, and the path radiation. By using MODTRAN to calculate atmospheric parameters such as water vapor and aerosol content, and using MODTRAN to simulate the radiance data under known atmospheric parameters, and substituting it into the above formula, the radiance under different atmospheric parameters can be obtained. and , considering the proximity effect, we can find , the initial value of the surface reflectivity can be obtained.

[0038] Step 3: Use the rational polynomial coefficients (RPC) model as the basic method of geometric correction to produce the remote sensing image input geometry file (IGM), experiment with the geographic positioning of remote sensing images, and realize the automatic registration of DEM / DSM data and remote sensing images to obtain the latitude, longitude and elevation values ​​of each pixel.

[0039] Step 4: According to the DEM / DSM data obtained in step 1, the slope, aspect, sky visibility factor, terrain visibility factor and shadow coefficient are solved, and the solved parameters are sampled to each pixel of the radiance image to accurately extract the terrain parameters of each pixel, including:

[0040]

[0041]

[0042] in, are the slope and the aspect, are the rates of elevation change in the north-south and east-west directions, respectively;

[0043]

[0044] in, They are the sky visibility factor and the terrain visibility factor, For this location The maximum zenith angle in the direction, The value of is between 0° and 360°; the projected shadow is calculated by the depth buffer algorithm in computer graphics, and a decimal between 0 and 1 is used to represent the proportion of the shadow area. When the coefficient is 0, the pixel is in the shadow; when the coefficient is 1, the pixel is in the sunlight area.

[0045] Step 5: Follow Figure 2 In the sun-target-sensor path shown in the figure, the radiation in the sun-target-sensor path is solved. The radiation in the figure includes: direct solar radiation a, atmospheric scattered radiation b, ground-atmosphere coupled radiation c, surrounding ground object reflection d, target reflected radiation E, path radiation F and atmospheric proximity effect G.

[0046] Set the imaging time, pixel location and elevation, and atmospheric parameters, and use MODTRAN to calculate the direct solar radiation , atmospheric scattered radiation , Program Radiation And transmittance .in is the row number where the pixel is located, The column number where the pixel is located, is the ground elevation at the pixel location, is the wavelength;

[0047] Direct radiation refers to the component of solar radiation that reaches the top of the atmosphere and then passes through the atmosphere directly to the target object; atmospheric scattered radiation refers to the component of solar radiation that reaches the target object after interacting with the atmosphere inside the atmosphere when passing through the atmosphere.

[0048] Step 6: Based on the initial surface reflectivity of step 2 , calculate the reflection of surrounding objects , Earth-atmosphere coupled radiation and atmospheric proximity effect .

[0049] The ground-atmosphere coupled radiation means that when a part of the energy passes through the atmosphere, the radiation interacts multiple times between the atmosphere and the undulating surface, and reaches the target object after scattering and reflection; the reflection of the surrounding objects means that the target pixel receives part of the radiation reflected by the surrounding undulating pixels to the target object. The path radiation is the radiation component that reaches the sensor entrance pupil after the solar radiation interacts with the atmosphere in the atmosphere, and some radiation components are scattered once or multiple times; the atmospheric proximity effect is due to the fact that a small part of the radiation components reflected by the surrounding objects will be reflected to the target object, and another small part will interact with the atmosphere to reach the sensor entrance pupil.

[0050] Use the reflectivity calculation formula to calculate the reflectivity median value :

[0051]

[0052] Step 7: Set the reflectivity middle value As a new initial value, the reflectivity is calculated in step 6 repeatedly by setting the number of times until the final reflectivity result converges, and the final high-precision surface reflectivity result under the undulating surface can be obtained.

[0053] In order to verify the effect of the designed atmosphere-topography geometry joint correction method, the Guanting Reservoir in Hebei Province was selected as the test area to carry out correction effect analysis. Figure 3a , Figure 3b The comparison of the true color composite image before and after the atmospheric terrain geometry joint correction of the Gaofen-5 AHSI image data is given. Combined with the image close-up, it can be seen that the scene becomes flat after correction, effectively eliminating the dependence on slope. This effect is particularly obvious in the box-marked area, because the influence of the surrounding terrain is accurately corrected, making the Pearl Lake in the middle area of ​​the image clearly distinguishable.

Claims

1. A method for joint correction of remote sensing imagery atmosphere and terrain geometry based on a radiation transfer model, characterized in that: The following steps are involved: Step 1: Obtain a hyperspectral remote sensing radiance image and a digital elevation model / digital surface model to obtain DEM / DSM data as input data; Digital elevation model is represented by DEM, and digital surface model is represented by DSM; Step 2: Ignore the terrain undulation, consider the contribution of the target and background to the radiance entering the sensor, and solve for the initial surface reflectivity ρ0; Step 3: Use the rational polynomial coefficient model as the basic method of geometric correction to perform geometric correction on the hyperspectral remote sensing radiance image, realize the geographic positioning of the remote sensing image, and realize the automatic registration of DEM / DSM data and remote sensing image to obtain the latitude, longitude and elevation values ​​of each pixel; Step 4: According to the DEM / DSM data obtained in step 1, the slope, aspect, sky visibility factor, terrain visibility factor and shadow coefficient are calculated, and they are sampled to each pixel position of the radiance image to accurately extract the terrain parameters of each pixel; Step 5: Set the imaging time, pixel location and elevation, and atmospheric parameters, and use the radiation transfer model to calculate the direct solar radiation. Atmospheric diffuse radiation Cheng Radiation L atm (x, y, z, λ) and transmittance T(z, λ), where x is the row number of the pixel, y is the column number of the pixel, z is the ground elevation at the pixel location, and λ is the wavelength; Step 6: Calculate the surrounding ground object reflection according to the initial reflectivity result ρ0 Earth-atmosphere coupled radiation and atmospheric proximity effect L a_adj (x, y, z, λ), and calculate the intermediate reflectivity value ρ n , put ρ n The value of is assigned to ρ0; Step 7: Calculate the intermediate reflectivity value ρ by looping step 6 for a set number of times n , until the reflectivity result converges, the final high-precision surface reflectivity result under the undulating surface can be obtained.

2. The method for joint correction of remote sensing image atmosphere and terrain geometry based on radiation transfer model as claimed in claim 1, characterized in that: In step 2, the specific calculation method of the initial surface reflectivity is: Where L is the entrance pupil radiance, ρ is the target pixel reflectivity, and ρ e is the average reflectivity of the background objects, S is the downward hemispherical albedo of the atmosphere, and L a is the path radiation, A and B are coefficients, and the entrance pupil radiance is composed of three components, namely the part reflected by the target pixel directly entering the sensor, the part reflected by the background object entering the sensor, and the path radiation.

3. The method for joint correction of remote sensing image atmosphere and terrain geometry based on radiation transfer model as claimed in claim 1, characterized in that: The specific calculation method of the terrain parameters in step 4 is as follows: Among them, Slo and Asp are slope and slope direction respectively, f x 、f y are the rates of elevation change in the north-south and east-west directions, respectively; Among them, V sky 、V terrain They are the sky visibility factor and the terrain visibility factor, where H φ It is the maximum zenith angle of the pixel in the φ direction, and the value of φ is between 0 and 2π. The projected shadow is calculated by the depth buffer algorithm in computer graphics, and a decimal between 0 and 1 is used to represent the proportion of the shadow area. When the coefficient is 0, the pixel is in the shadow; when the coefficient is 1, the pixel is in the sunlight area.

4. The method for joint correction of remote sensing image atmosphere and terrain geometry based on radiation transfer model as claimed in claim 1, characterized in that: In step 6, calculate the reflectivity median value ρ n The method is as follows: Among them, L(x,y,z,λ) is the entrance pupil radiance and T(z,λ) is the atmospheric transmittance.

5. The method for joint correction of remote sensing image atmosphere and terrain geometry based on radiation transfer model as claimed in claim 1, characterized in that: In step 7, the reflectivity intermediate value ρ n As the new initial value, the reflectivity is iteratively calculated by setting the number of times until the final reflectivity result converges, that is, the final high-precision surface reflectivity result under the undulating surface is obtained.

Images