Aviation hyperspectral image surface reflectance inversion method

By combining the radiation transmission model and the SCS-C terrain correction model, the organic combination of atmospheric and terrain correction is achieved, solving the problems of low inversion accuracy and difficult parameter acquisition in the prior art, and improving the inversion accuracy of surface reflectivity of aerial hyperspectral images.

CN119782709APending Publication Date: 2025-04-08BEIJING RES INST OF URANIUM GEOLOGY
View PDF 0 Cites 0 Cited by

Patent Information

Application Number
CN202411958757.1
Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
Filing Date
2024-12-30
Publication Date
2025-04-08

AI Technical Summary

Technical Problem

The prior art is difficult to organically combine atmospheric correction and topographic correction, which leads to difficulty in inversion of high-precision surface reflectivity of aeronautical hyperspectral images and difficult to obtain correction parameters.

Method used

The radiation transmission model is used to perform atmospheric correction, combined with the SCS-C terrain correction model, the atmospheric water vapor content and aerosol optical thickness are inverted through the APDA ratio method, and the organic combination of atmospheric and terrain correction is achieved.

Benefits of technology

The inversion accuracy of surface reflectivity of aeronautical hyperspectral remote sensing images is improved, over-correction problems are reduced, and the parameter acquisition process is simplified.

✦ Generated by Eureka AI based on patent content.

Smart Images

  • Figure CN119782709A_ABST
    Figure CN119782709A_ABST
Patent Text Reader

Abstract

The invention belongs to the field of remote sensing, and particularly relates to an aviation hyperspectral image surface reflectance inversion method, which comprises the following steps: acquiring aviation hyperspectral remote sensing image data and corresponding digital elevation model (DEM) data; the method comprises the following steps: preprocessing aviation hyperspectral remote sensing image data, and converting a DN value of a sensor into radiance; based on the preprocessed aviation hyperspectral remote sensing image data, adopting an APDA ratio method to invert the atmospheric water vapor content; calculating the optical thickness of the aerosol; solving the atmospheric radiation transfer equation of the visible light-short wave infrared spectrum to obtain an intermediate reflectivity; preprocessing the obtained digital elevation model (DEM) data; calculating to obtain a solar effective incident angle; calculating an empirical parameter C value in the SCS-C terrain correction model; and based on the intermediate reflectivity, performing terrain correction calculation in combination with an SCS-C model to obtain the final real surface reflectivity. According to the invention, high-precision inversion of the surface reflectance of the aviation hyperspectral remote sensing image is realized.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] The invention belongs to the field of remote sensing, and in particular relates to a method for inverting surface reflectance of aerial hyperspectral images. Background Art

[0002] The wide application of hyperspectral remote sensing technology in the fields of mineral resource exploration, ecological environment monitoring, military target detection and identification, etc., fully demonstrates the advantage of imaging spectroscopy technology in obtaining geometric, radiation and spectral multi-dimensional information "integration of images and spectra". In mountainous areas with undulating terrain, aerial hyperspectral images are affected by both atmospheric and topographic factors. High-precision inversion of surface reflectivity is the basis for the application of hyperspectral data. The main purpose of atmospheric correction is to reduce the impact of atmospheric absorption, scattering and reflection on the received signal of the remote sensor, and terrain correction aims to reduce the spectral differences of similar objects on sunny and shady slopes caused by terrain effects.

[0003] At present, there is a lack of suitable and efficient methods to organically combine atmospheric correction and terrain correction for the inversion of surface reflectance of aerial hyperspectral images. Existing methods have obvious problems of over-correction and difficulty in obtaining correction parameters. The present invention carries out atmospheric correction based on the radiation transfer model, and its parameters are clear and easy to obtain. At the same time, the SCS-C model is introduced for terrain correction, which better overcomes the over-correction problem of traditional terrain correction models. Summary of the invention

[0004] The purpose of the present invention is to provide a method for inverting the surface reflectance of aerial hyperspectral images. The method introduces the SCS-C terrain correction model on the basis of atmospheric correction based on the radiation transfer model to achieve high-precision inversion of the surface reflectance of aerial hyperspectral remote sensing images.

[0005] The technical solution to achieve the purpose of the present invention is:

[0006] A method for inverting surface reflectance of aerial hyperspectral images, the method comprising:

[0007] Step 1: Obtain aerial hyperspectral remote sensing image data and corresponding digital elevation model DEM data;

[0008] Step 2: Preprocess the aerial hyperspectral remote sensing image data obtained in step 1 and convert the DN value of the sensor into radiance;

[0009] Step 3: Based on the aerial hyperspectral remote sensing image data preprocessed in step 2, the APDA ratio method is used to invert the atmospheric water vapor content;

[0010] Step 4: Calculate aerosol optical depth;

[0011] Step 5: Based on the atmospheric water vapor content and aerosol optical depth obtained in steps 3 and 4, combined with the MODTRAN model and the metadata of the image, the atmospheric radiation transfer equation in the visible-shortwave infrared spectrum is solved to obtain other parameters such as atmospheric path radiation La, atmospheric transmittance Gb, atmospheric geometry coefficient Gt and atmospheric hemispherical albedo S, and then the intermediate reflectivity is obtained;

[0012] Step 6: Preprocess the digital elevation model (DEM) data obtained in step 1

[0013] Step 7: Based on the digital elevation model DEM data processed in step 6, terrain related factors are calculated to obtain slope and aspect data, and the effective solar incidence angle is further calculated from the slope and aspect;

[0014] Step 8: Calculate the empirical parameter C value in the SCS-C terrain correction model;

[0015] Step 9: Based on the intermediate reflectivity obtained in step 5, the terrain correction calculation is performed in combination with the SCS-C model to obtain the final true surface reflectivity.

[0016] The preprocessing of the aerial hyperspectral remote sensing image data obtained in step 1 in step 2 includes geometric correction and radiation calibration. The calculation formula for converting the DN value of the sensor into radiance is:

[0017] L=gain×DN+offset (1)

[0018] Where L is the radiance, DN is the pixel value of band λ, gain and offset are the gain and offset of band λ respectively.

[0019] The formula for calculating the water vapor content PW using the APDA algorithm in step 3 is:

[0020]

[0021] Among them, L atm,m L is the atmospheric radiation of the absorption channel. atm,r is the atmospheric radiation of the reference channel, λ r Indicates the central wavelength of the reference channel, L m and L r are the radiances of the absorption and reference channels before correction, respectively; LIR([x],[y])| a represents the linear regression line passing through the point (x, y) at the point x = a; α, β, γ are fitting parameters, which are related to the water vapor content and aerosol type; R APDA is the APDA channel ratio, which can be obtained according to the algorithm definition formula; PW is the water vapor content.

[0022] The step 3 comprises:

[0023] Step 3.1: Generate water vapor content using the MODTRAN radiation transfer model and look up the parameters in table LUT1. The parameters in the lookup table include: water vapor content PW, path radiation of water vapor absorption channel and reference channel, APDA channel ratio R APDA ;

[0024] Step 3.2: Use formula (3) to fit the ratio R of all APDA channels in LUT1 APDA and the corresponding water vapor content PW, and store the fitting parameters α, β, γ;

[0025] Step 3.3: Set the initial water vapor content PW1, and obtain the process radiation of the water vapor absorption channel and the reference channel from LUT1;

[0026] Step 3.4: Calculate the APDA channel ratio R of the pixel using formula (2) APDA , and calculate the water vapor content value PW2 of the pixel according to formula (4);

[0027] Step 3.5: According to PW2 combined with LUT1, the new water vapor absorption channel and the path radiation value of the reference channel are obtained to replace the previously obtained path radiation;

[0028] Step 3.6: Use formula (2) again to calculate the APDA channel ratio R APDA , according to formula (4), the APDA ratio is converted into the final water vapor content PW3.

[0029] The step 4 comprises:

[0030] Step 4.1: Select the location of the dark pixel based on the NDVI value and the apparent reflectance of the near-infrared band;

[0031] Step 4.2: Based on the average water vapor value obtained by water vapor content inversion, call MODTRAN to establish the visibility lookup table LUT2;

[0032] Step 4.3: Calculate the surface reflectance of the red and blue bands, as well as the parameters corresponding to different visibility VIS in the visibility lookup table LUT2, including: atmospheric radiation L a , atmospheric hemispheric albedo S, radiation correction parameter G b , G t ;

[0033] Step 4.4: Substitute the parameters obtained in step 4.3 into the atmospheric radiation transfer equation, combine the lookup table LUT2 to obtain the visibility VIS value that satisfies the equation, and then interpolate to obtain the visibility VIS corresponding to the red and blue bands. Calculate the average value of the two to be the visibility VIS at the dark pixel;

[0034] Step 4.5: Based on the relationship between visibility at 550 nm and aerosol optical depth, calculate the aerosol optical depth.

[0035] The atmospheric radiation transfer equation in step 5 is:

[0036]

[0037] Among them, ρ is the surface reflectivity of the pixel, ρ e is the average surface reflectivity of the pixel and its neighboring pixels, S is the downward hemispherical albedo of the atmosphere, and L a is the atmospheric radiation, G b and G t is a coefficient that depends on the atmospheric transmittance and geometry;

[0038] The calculation formula of the average surface reflectance of the pixel and the nearby pixels is:

[0039]

[0040] Among them, L e It is the average value of the pixel and the spatial radiation image around the pixel, and can be calculated using the original image.

[0041] The step 5 comprises:

[0042] Step 5.1: Based on the average visibility, call MODTRAN to create the surface reflectance lookup table LUT3.

[0043] Step 5.2: Based on the water vapor content PW of the pixel, the atmospheric correction parameters of the pixel are obtained by spatial interpolation of the lookup table LUT3: atmospheric path radiation L a , atmospheric hemispheric albedo S, radiation correction parameter G b , G t ;

[0044] Step 5.3: According to formula (6) and the atmospheric correction parameters obtained in step 5.2, the average surface reflectance ρ of the pixel and its neighboring pixels is calculated. e ;

[0045] Step 5.4: For the parameters obtained in steps 5.1-5.3, solve the surface reflectance ρ of the pixel according to formula (5), and this surface reflectance is called the intermediate reflectance.

[0046] The calculation formulas for the slope, slope direction, and solar effective incident angle in step 7 are:

[0047]

[0048]

[0049] Among them, θ t and Asp represent the slope and aspect respectively, f x and f y are the elevation changes in the north-south and east-west directions, i t is the effective incident angle of the sun, θ S is the solar zenith angle, is the solar azimuth.

[0050] The calculation formula of the empirical parameter C value in the SCS-C terrain correction model in step 8 is as follows:

[0051] DN = acosi t +b (10)

[0052] C = b / a

[0053] Among them, b is the intercept of the regression equation, a is the slope of the equation, and C is the empirical parameter.

[0054] The calculation formula of the true reflectivity of the surface in step 9 is:

[0055]

[0056] Among them, ρ T is the slope surface reflectivity, that is, the intermediate reflectivity obtained in step 5, ρ H is the final true reflectivity of the surface.

[0057] The beneficial technical effects of the present invention are:

[0058] 1. The present invention provides a method for inverting the surface reflectance of aerial hyperspectral images. Based on the atmospheric correction using the MODTRAN model, the SCS-C terrain correction model is introduced to organically combine the atmospheric correction and the terrain correction, thereby improving the inversion accuracy of the surface reflectance of aerial hyperspectral remote sensing images.

[0059] 2. In the surface reflectance inversion method of aerial hyperspectral image provided by the present invention, the atmospheric correction process adopts the APDA atmospheric preprocessing differential absorption method to invert the water vapor content, which reduces the interference of other atmospheric components on the water vapor absorption signal. At the same time, the dark pixel method is used to invert the aerosol optical thickness, which improves the estimation accuracy of parameters such as atmospheric transmittance, atmospheric path radiation and atmospheric hemispherical albedo in terrain undulating areas.

[0060] 3. In the surface reflectance inversion method of aerial hyperspectral image provided by the present invention, the terrain correction process adopts the SCS-C correction model, which takes into account the changes of the canopy to the terrain, can improve the over-correction problem, and the model parameters are easy to obtain and the calculation efficiency is high. BRIEF DESCRIPTION OF THE DRAWINGS

[0061] Figure 1 This is a flow chart of a method for inverting surface reflectance of aerial hyperspectral images provided by the present invention. DETAILED DESCRIPTION

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

[0063] The invention provides an aerial hyperspectral image surface reflectance inversion method, which starts from the theoretical basis of atmosphere correction and terrain correction, inverts atmospheric water vapor content and aerosol optical thickness based on aerial hyperspectral data, adopts the MODTRAN model to calculate parameters such as atmospheric transmittance, atmospheric path radiation and atmospheric hemispherical albedo, and calculates terrain correction parameters based on the atmospheric radiation transfer equation of the visible light-shortwave infrared spectrum and the atmospheric parameters simulated by MODTRAN, and at the same time calculates terrain correction parameters based on the SCS-C model, so as to realize the aerial hyperspectral surface reflectance inversion coupled with atmosphere correction and terrain correction.

[0064] like Figure 1 As shown, the specific steps include:

[0065] Step 1: Obtain aerial hyperspectral remote sensing image data and corresponding digital elevation model DEM data;

[0066] Step 2: Preprocess the aerial hyperspectral remote sensing image data obtained in step 1 and convert the DN value of the sensor into radiance;

[0067] The preprocessing of the aerial hyperspectral remote sensing image data obtained in step 1 includes geometric correction and radiation calibration. The calculation formula for converting the DN value of the sensor into radiance is:

[0068] L=gain×DN+offset (1)

[0069] Where L is the radiance, DN is the pixel value of band λ, gain and offset are the gain and offset of band λ respectively.

[0070] Step 3: Based on the aerial hyperspectral remote sensing image data preprocessed in step 2, the atmospheric water vapor content is inverted using the APDA ratio method (atmospheric preprocessing differential absorption method);

[0071] The calculation of water vapor content PW using the APDA algorithm is completed according to the following two formulas:

[0072]

[0073] Among them, R APDA is the APDA channel ratio, L atm,m L is the atmospheric radiation of the absorption channel. atm,r is the atmospheric radiation of the reference channel, λr Indicates the central wavelength of the reference channel, L m and L r are the radiances of the absorption and reference channels before correction, respectively. LIR([x],[y])| a represents the linear regression line passing through the point (x, y) at the point x = a. α, β, γ are fitting parameters, which are related to the water vapor content and aerosol type; PW is the water vapor content.

[0074] The specific implementation includes the following steps:

[0075] Step 3.1: Generate water vapor content using the MODTRAN radiation transfer model and look up the parameters in table LUT1. The parameters in the lookup table include: water vapor content PW, path radiation of water vapor absorption channel and reference channel, APDA channel ratio R APDA ;

[0076] Step 3.2: Use formula (3) to fit the ratio R of all APDA channels in LUT1 APDA and the corresponding water vapor content PW, and store the fitting parameters α, β, γ;

[0077] Step 3.3: Set the initial water vapor content PW1, and obtain the process radiation of the water vapor absorption channel and the reference channel from LUT1;

[0078] Step 3.4: Calculate the APDA channel ratio R of the pixel using formula (2) APDA , and calculate the water vapor content value PW2 of the pixel according to formula (4);

[0079] Step 3.5: According to PW2 combined with LUT1, the new water vapor absorption channel and the path radiation value of the reference channel are obtained to replace the previously obtained path radiation;

[0080] Step 3.6: Use formula (2) again to calculate the APDA channel ratio R APDA , according to formula (4), the APDA ratio is converted into the final water vapor content PW3.

[0081] Step 4: Calculate aerosol optical depth;

[0082] The specific implementation includes the following steps:

[0083] Step 4.1: Select the location of the dark pixel based on the NDVI value and the apparent reflectance of the near-infrared band;

[0084] Step 4.2: Based on the average water vapor value obtained by water vapor content inversion, call MODTRAN to establish the visibility lookup table LUT2;

[0085] Step 4.3: Calculate the surface reflectance of the red and blue bands, as well as the parameters corresponding to different visibility VIS in the visibility lookup table LUT2, including: atmospheric radiation L a , atmospheric hemispheric albedo S, radiation correction parameter G b , G t ;

[0086] Step 4.4: Substitute the parameters obtained in step 4.3 into the atmospheric radiation transfer equation, combine the lookup table LUT2 to obtain the visibility VIS value that satisfies the equation, and then interpolate to obtain the visibility VIS corresponding to the red and blue bands. Calculate the average value of the two to be the visibility VIS at the dark pixel;

[0087] Step 4.5: Based on the relationship between visibility at 550 nm and aerosol optical depth, calculate the aerosol optical depth.

[0088] Step 5: Based on the atmospheric water vapor content and aerosol optical depth obtained in steps 3 and 4, combined with the MODTRAN model and image metadata, the atmospheric radiation transfer equation in the visible-shortwave infrared spectrum is solved to obtain the atmospheric path radiation L a , atmospheric transmittance G b , the coefficient of atmospheric geometry G t and other parameters such as atmospheric hemispherical albedo S, and then the intermediate reflectivity is obtained;

[0089] The atmospheric radiation transfer equation is:

[0090]

[0091] Among them, ρ is the surface reflectivity of the pixel, ρ e is the average surface reflectivity of the pixel and its neighboring pixels, S is the downward hemispherical albedo of the atmosphere, and L a is the atmospheric radiation, G b and G t is a coefficient that depends on the transmittance of the atmosphere and the geometric conditions.

[0092] It indicates the radiation brightness of the solar radiation that enters the surface of the earth after passing through the atmosphere and then directly enters the sensor. It indicates the radiance of the surface pixels entering the sensor due to the reflection of the atmosphere.

[0093] In the formula, ρ and ρ e The main reason for the difference is the “neighboring pixel effect” caused by atmospheric scattering, ρ e It can be estimated by the following equation:

[0094]

[0095] Among them, Le It is the average value of the pixel and the spatial radiation image around the pixel, and can be calculated using the original image.

[0096] The specific implementation includes the following steps:

[0097] Step 5.1: Based on the average visibility, call MODTRAN to create the surface reflectance lookup table LUT3.

[0098] Step 5.2: Based on the water vapor content PW of the pixel, the atmospheric correction parameters of the pixel are obtained by spatial interpolation of the lookup table LUT3: atmospheric path radiation L a , atmospheric hemispheric albedo S, radiation correction parameter G b , G t ;

[0099] Step 5.3: According to formula (6) and the atmospheric correction parameters obtained in step 5.2, the average surface reflectance ρ of the pixel and its neighboring pixels is calculated. e .

[0100] Step 5.4: For the parameters obtained in steps 5.1-5.3, solve the surface reflectance ρ of the pixel according to formula (5), and this surface reflectance is called the intermediate reflectance;

[0101] Step 6: Preprocess the digital elevation model (DEM) data obtained in step 1

[0102] The preprocessing of the digital elevation model DEM data obtained in step 1 includes operations such as clipping and smoothing;

[0103] Step 7: Based on the digital elevation model DEM data processed in step 6, terrain related factors are calculated to obtain slope and aspect data, and the effective solar incidence angle is further calculated from the slope and aspect;

[0104] The specific calculation method is as follows:

[0105]

[0106]

[0107] Among them, θ t and Asp represent the slope and aspect respectively, f x and f y are the elevation changes in the north-south and east-west directions, i t is the effective incident angle of the sun, θ s is the solar zenith angle, is the solar azimuth.

[0108] Step 8: Calculate the empirical parameter C value in the SCS-C terrain correction model;

[0109] The specific calculation method is as follows:

[0110] DN = acosi t +b (10)

[0111] A regression equation is established between the pixel DN value and the cosine value of the effective solar incidence angle, where b is the intercept of the regression equation and a is the slope of the equation, and the empirical parameter C=b / a is obtained.

[0112] Step 9: Based on the intermediate reflectivity obtained in step 5, the terrain correction calculation is performed in combination with the SCS-C model to obtain the final true surface reflectivity.

[0113] The specific calculation method is as follows:

[0114]

[0115] Among them, ρ T is the slope surface reflectivity, that is, the intermediate reflectivity obtained in step 5, ρ H is the final true reflectivity of the surface.

[0116] The present invention is described in detail above with reference to the accompanying drawings and embodiments, but the present invention is not limited to the above embodiments, and various changes can be made within the knowledge of ordinary technicians in the field without departing from the purpose of the present invention. The contents not described in detail in the present invention can adopt the existing technology.

Claims

1. A method for inverting surface reflectance from aerial hyperspectral images, characterized in that: The method comprises: Step 1: Obtain aerial hyperspectral remote sensing image data and corresponding digital elevation model DEM data; Step 2: Preprocess the aerial hyperspectral remote sensing image data obtained in step 1 and convert the DN value of the sensor into radiance; Step 3: Based on the aerial hyperspectral remote sensing image data preprocessed in step 2, the APDA ratio method is used to invert the atmospheric water vapor content; Step 4: Calculate aerosol optical depth; Step 5: Based on the atmospheric water vapor content and aerosol optical depth obtained in steps 3 and 4, combined with the MODTRAN model and the metadata of the image, the atmospheric radiation transfer equation in the visible-shortwave infrared spectrum is solved to obtain the atmospheric path radiation La, atmospheric transmittance Gb, atmospheric geometry coefficient Gt and atmospheric hemispherical albedo S, and then the intermediate reflectivity is obtained; Step 6: Preprocess the digital elevation model (DEM) data obtained in step 1 Step 7: Based on the digital elevation model DEM data processed in step 6, terrain related factors are calculated to obtain slope and aspect data, and the effective solar incidence angle is further calculated from the slope and aspect; Step 8: Calculate the empirical parameter C value in the SCS-C terrain correction model; Step 9: Based on the intermediate reflectivity obtained in step 5, the terrain correction calculation is performed in combination with the SCS-C model to obtain the final true surface reflectivity.

2. The method for inverting surface reflectance of aerial hyperspectral images according to claim 1, characterized in that: The preprocessing of the aerial hyperspectral remote sensing image data obtained in step 1 in step 2 includes geometric correction and radiation calibration. The calculation formula for converting the DN value of the sensor into radiance is: L=gain×DN+offset (1) Where L is the radiance, DN is the pixel value of band λ, gain and offset are the gain and offset of band λ respectively.

3. The method for inverting surface reflectance of aerial hyperspectral images according to claim 1, characterized in that: The formula for calculating the water vapor content PW using the APDA algorithm in step 3 is: Among them, L atm,m L is the atmospheric radiation of the absorption channel. atm,r is the atmospheric radiation of the reference channel, λ r Indicates the central wavelength of the reference channel, L m and L r are the radiances of the absorption and reference channels before correction, respectively; LIR([x],[y])| a represents the linear regression line passing through the point (x, y) at the point x = a; α, β, γ are fitting parameters, which are related to the water vapor content and aerosol type; R APDA is the APDA channel ratio, which can be obtained according to the algorithm definition formula; PW is the water vapor content.

4. The method for inverting surface reflectance of aerial hyperspectral images according to claim 3, characterized in that: The step 3 comprises: Step 3.1: Generate water vapor content using the MODTRAN radiation transfer model and look up the parameters in table LUT1. The parameters in the lookup table include: water vapor content PW, path radiation of water vapor absorption channel and reference channel, APDA channel ratio R APDA ; Step 3.2: Use formula (3) to fit the ratio R of all APDA channels in LUT1 APDA and the corresponding water vapor content PW, and store the fitting parameters α, β, γ; Step 3.3: Set the initial water vapor content PW1, and obtain the process radiation of the water vapor absorption channel and the reference channel from LUT1; Step 3.4: Calculate the APDA channel ratio R of the pixel using formula (2) APDA , and calculate the water vapor content value PW2 of the pixel according to formula (4); Step 3.5: According to PW2 combined with LUT1, the new water vapor absorption channel and the path radiation value of the reference channel are obtained to replace the previously obtained path radiation; Step 3.6: Use formula (2) again to calculate the APDA channel ratio R APDA , according to formula (4), the APDA ratio is converted into the final water vapor content PW3.

5. The method for inverting surface reflectance of aerial hyperspectral images according to claim 1, characterized in that: The step 4 comprises: Step 4.1: Select the location of the dark pixel based on the NDVI value and the apparent reflectance of the near-infrared band; Step 4.2: Based on the average water vapor value obtained by water vapor content inversion, call MODTRAN to establish the visibility lookup table LUT2; Step 4.3: Calculate the surface reflectance of the red and blue bands, as well as the parameters corresponding to different visibility VIS in the visibility lookup table LUT2, including: atmospheric radiation L a , atmospheric hemispheric albedo S, radiation correction parameter G b , G t ; Step 4.4: Substitute the parameters obtained in step 4.3 into the atmospheric radiation transfer equation, combine the lookup table LUT2 to obtain the visibility VIS value that satisfies the equation, and then interpolate to obtain the visibility VIS corresponding to the red and blue bands. Calculate the average value of the two to be the visibility VIS at the dark pixel; Step 4.5: Based on the relationship between visibility at 550 nm and aerosol optical depth, calculate the aerosol optical depth.

6. The method for inverting surface reflectance of aerial hyperspectral images according to claim 1, characterized in that: The atmospheric radiation transfer equation in step 5 is: Among them, ρ is the surface reflectivity of the pixel, ρ e is the average surface reflectivity of the pixel and its neighboring pixels, S is the downward hemispherical albedo of the atmosphere, and L a is the atmospheric radiation, G b and G t is a coefficient that depends on the atmospheric transmittance and geometry; The calculation formula of the average surface reflectance of the pixel and the nearby pixels is: Among them, L e It is the average value of the pixel and the spatial radiation image around the pixel, which can be calculated using the original image.

7. The method for inverting surface reflectance of aerial hyperspectral images according to claim 6, characterized in that: The step 5 comprises: Step 5.1: Based on the average visibility, call MODTRAN to create the surface reflectance lookup table LUT3. Step 5.2: Based on the water vapor content PW of the pixel, the atmospheric correction parameters of the pixel are obtained by spatial interpolation of the lookup table LUT3: atmospheric path radiation L a , atmospheric hemispheric albedo S, radiation correction parameter G b , G t ; Step 5.3: According to formula (6) and the atmospheric correction parameters obtained in step 5.2, the average surface reflectance ρ of the pixel and its neighboring pixels is calculated. e ; Step 5.4: For the parameters obtained in steps 5.1-5.3, solve the surface reflectance ρ of the pixel according to formula (5), and this surface reflectance is called the intermediate reflectance.

8. The method for inverting surface reflectance of aerial hyperspectral images according to claim 1, characterized in that: The calculation formulas for the slope, slope direction, and solar effective incident angle in step 7 are: Among them, θ t and Asp represent the slope and aspect respectively, f x and f y are the elevation changes in the north-south and east-west directions, i t is the effective incident angle of the sun, θ S is the solar zenith angle, is the solar azimuth.

9. The method for inverting surface reflectance of aerial hyperspectral images according to claim 1, characterized in that: The calculation formula of the empirical parameter C value in the SCS-C terrain correction model in step 8 is as follows: DN=so t +b (10) C = b / a Among them, b is the intercept of the regression equation, a is the slope of the equation, and C is the empirical parameter.

10. The method for inverting surface reflectance of aerial hyperspectral images according to claim 1, characterized in that: The calculation formula of the true reflectivity of the surface in step 9 is: Among them, ρ T is the slope surface reflectivity, that is, the intermediate reflectivity obtained in step 5, ρ H is the final true reflectivity of the surface.