Aerosol optical depth inversion method

Abstract

Provided in the present invention is an aerosol optical depth (AOD) inversion method. The method comprises: S1, acquiring a multi-spectral image and a corresponding MODIS image, and performing pre-processing; S2, performing radiation calibration on the multi-spectral image, so as to obtain an apparent reflectance of each waveband; S3, removing cloud pixels and cloud shadow pixels, so as to obtain an image to be subjected to inversion; S4, calculating a normalized vegetation index of each pixel in said image, and performing classification according to a dark pixel area and a bright pixel area; S5, estimating a red-blue waveband surface reflectance ratio in the dark pixel area by using a dark target method, and constructing a monthly red-blue waveband surface reflectance ratio library in the bright pixel area according to the MODIS image; S6, establishing an AOD inversion lookup table; and S7, iteratively calculating red-blue waveband surface reflectance theoretical ratios under different AOD values by using the AOD lookup table, and calculating differences between the red-blue waveband surface reflectance theoretical ratios and values in the monthly red-blue waveband surface reflectance ratio library, and using, as an inversion result, the AOD value when the difference is the minimum. The present invention can simultaneously acquire high-resolution AOD inversion results of dark pixels and bright pixels, so as to accurately describe a change in the aerosol in an urban area.

Description

A method for inverting aerosol optical thickness Technical Field

[0001] The present invention relates to the field of atmospheric remote sensing technology, and in particular to an inversion method for aerosol optical thickness. Background Art

[0002] Atmospheric aerosols are dispersed particles suspended in the air. They include primary aerosols (directly emitted into the atmosphere as particles) and secondary aerosols (converted from major atmospheric pollutants), with diameters ranging from 0.001 to 100 μm. Changes in atmospheric aerosol concentrations directly affect human health and air quality, and influence the Earth-atmosphere budget and climate change through direct or indirect effects of radiative forcing. Therefore, monitoring the spatial distribution and changes of aerosols is crucial.

[0003] Aerosol Optical Depth (AOD) is a key parameter for describing aerosol properties and represents atmospheric turbidity. Traditional aerosol monitoring is mostly ground-based. However, due to limited station deployment, it is difficult to obtain large-scale spatial distribution information of aerosols. In recent years, satellite remote sensing technology, with its high timeliness and large-scale observation capabilities, has enabled the acquisition of spatially continuous aerosol data, showing great potential for monitoring macro-environmental and pollution distribution. The key to AOD inversion is determining surface reflectivity information to facilitate ground-atmosphere decoupling. Inversion algorithms include dark target algorithms, deep blue algorithms, structure function methods, polarization algorithms, and multi-angle algorithms. The structure function method requires a "clean day" image as a benchmark and has high requirements for geometric correction; the polarization algorithm can only be applied to invert fine particle aerosols; and the multi-angle algorithm requires specific sensor support. All three algorithms are difficult to implement in operational applications. Dark target algorithms have been continuously improved and applied in recent years, but they still have significant limitations. They require support from the shortwave infrared band and are not suitable for bright pixel areas (such as cities and deserts). The Deep Blue algorithm can invert the AOD in bright pixel areas, but its accuracy is lower than that of the Dark Target algorithm. Moreover, the Deep Blue algorithm's inversion results are limited by the spatial resolution of the external surface reflectance data.

[0004] Traditional aerosol inversion methods cannot maintain the same spatial resolution for both dark and bright regions. Therefore, how to simultaneously invert the aerosol optical depth in both dark and bright pixel regions without sacrificing the spatial resolution of the original data remains a challenge in current aerosol optical depth inversion research.

[0005] Summary of the Invention

[0006] In view of the above problems, the present invention provides an inversion method for aerosol optical depth. By constructing a monthly red-blue band surface reflectance ratio library and adopting the red-blue band ratio method, high spatial resolution AOD inversion results including dark and bright pixels are obtained simultaneously. The method has good spatial continuity and can more accurately describe the spatial changes of aerosols in urban areas.

[0007] The technical solutions of the present invention are as follows:

[0008] S1 obtains multispectral images and corresponding MODIS images and preprocesses the multispectral images;

[0009] S2 performs radiometric calibration on the pre-processed multispectral image to obtain apparent reflectance data of each band, including a blue band, a green band, a red band, and a near-infrared band;

[0010] S3 removes cloud pixels and cloud shadow pixels from the pre-processed multispectral image to obtain the image to be inverted;

[0011] S4 calculates the normalized vegetation index of each pixel in the image to be inverted, and classifies each pixel in the image to be inverted according to the normalized vegetation index according to a preset threshold value, into a dark pixel area and a bright pixel area;

[0012] S5 estimates the surface reflectance ratio of the first red and blue bands using the dark target method in dark pixel areas. In bright pixel areas, based on the surface reflectance data of each band in MODIS imagery, a monthly red and blue band surface reflectance ratio library is constructed using the sub-minimum synthesis method through the Google Earth Engine platform.

[0013] S6 simulates the apparent reflectivity detected by different satellite sensors under different conditions based on the atmospheric radiation transfer model and establishes an AOD lookup table;

[0014] S7 uses the AOD lookup table to iteratively calculate the theoretical ratio of the red and blue band surface reflectance at different AOD values, and subtracts the theoretical ratio from the first red and blue band surface reflectance ratio and the monthly red and blue band surface reflectance ratio library, respectively, and takes the AOD value in the AOD lookup table corresponding to the theoretical ratio of the red and blue band surface reflectance with the smallest difference as the inversion result of the multispectral image;

[0015] Specifically, in step S5, in the dark pixel area, the dark target algorithm is used to estimate the surface reflectance ratio of the first red and blue bands, including the following steps:

[0016] In the dark pixel area of ​​the image to be inverted, the surface reflectance of the red and blue bands is estimated using the shortwave infrared band apparent reflectance through a linear relationship. The calculation formula is:

[0017] Among them, ρ BLUE is the surface reflectivity in the blue light band, ρ RED is the surface reflectivity in the red light band, ρ SWIRW It is the TOA apparent reflectivity in the 2.1μm short infrared band;

[0018] Divide the surface reflectance of the blue band by the surface reflectance of the red band to obtain the surface reflectance ratio of the first red and blue bands The value is fixed at 2.

[0019] Specifically, in step S5, in the bright pixel area, a sub-minimum synthesis method is used to construct a monthly red and blue band surface reflectance ratio library, including the following steps:

[0020] Obtain MODIS images through the Google Earth Engine platform;

[0021] MODIS images are preprocessed, and multiple images of the same day in a long time series are synthesized into one image using the minimum value synthesis method, and the area that matches the time and space of the multispectral image is cropped;

[0022] Obtain the surface reflectance data for each month and use the sub-minimum composite method to construct the monthly red and blue band surface reflectance ratio library ρ(mon,wrs): ρ(mon,wrs)=2ndMin_composite(ρ1(mon,wrs),...ρ n (mon,wrs))

[0023] Among them, ρ(mon,wrs) represents the construction of the monthly red and blue band surface reflectance ratio image library; ρ1, ρ n An image representing the ratio of red and blue band surface reflectance from a month of MODIS imagery.

[0024] Specifically, step S4 includes:

[0025] S41 calculates the Normalized Difference Vegetation Index (NDVI) of each pixel based on the apparent reflectance of the red and near-infrared bands of each pixel in the multispectral image.

[0026] S42 marks the pixels with NDVI>0.55 on the image to be inverted as dark pixel areas;

[0027] S43 marks the pixels with NDVI≦0.55 and the corresponding pixels on the image to be inverted as bright pixel areas.

[0028] 9. Specifically, in step S7, iterative calculations are performed at different AOD values, and the variation of different AOD values ​​is within the range of 0 to 2.0.

[0029] Specifically, the preprocessing in step S1 includes orthophotography and geometric correction.

[0030] Specifically, step S3 includes:

[0031] S31 obtains the band value of each pixel of each multispectral image by calculation;

[0032] S32 uses long-time series of remote sensing samples to construct a clear sky background field containing a blue light band. Based on the clear sky background field containing a blue light band and the blue light band reflectance of each pixel in the multispectral image, the cloud detection threshold is calculated.

[0033] S33 performs cloud detection on the image to be inverted using a cloud detection threshold, removes cloud and cloud shadow pixels, and obtains the image to be inverted.

[0034] Specifically, step S8 includes: using AERONET and SONET site data to verify the accuracy of the inverted AOD value, and the evaluation indicators include root mean square error and expected error, and the calculation formula is: EE=±(0.05+0.2×τ0)

[0035] Where RMSE is the root mean square error, EE is the expected error, τ is the inverted AOD value, τ0 is the station AOD value, and n is the number of verification samples.

[0036] The beneficial effects of the present invention are:

[0037] (1) The present invention divides MODIS images into dark pixel areas and bright pixel areas, and adopts different calculation methods for different areas. It can simultaneously obtain high-spatial-resolution AOD inversion results for both dark and bright pixels, with good spatial continuity, and can more accurately describe the spatial variation of aerosols in urban areas.

[0038] (2) In bright pixel areas, the present invention uses the sub-minimum synthesis method to construct a ratio library. This method assumes that the surface reflectance of most ground objects remains unchanged within a month, eliminates negative and invalid values ​​in the image, and is less affected by clouds, mountains, buildings and their shadows, thereby improving the accuracy of the inversion results. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments. Obviously, the drawings described below are only some embodiments of the present invention. For ordinary technicians in this field, other drawings can be obtained based on these drawings without paying any creative work.

[0040] FIG1 is a schematic diagram of a method according to an embodiment of the present invention;

[0041] FIG2 is a technical flow chart of an embodiment of the present invention. DETAILED DESCRIPTION

[0042] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. All other embodiments obtained by ordinary technicians in this field based on the embodiments of the present invention are within the scope of protection of the present invention.

[0043] To achieve the above object, the present invention provides an aerosol optical thickness inversion method, referring to Figures 1 and 2 . The method comprises the following steps:

[0044] S1 obtains multispectral images and corresponding MODIS images and preprocesses the multispectral images;

[0045] In this embodiment of the present invention, the multispectral imagery is derived from the 2022 imagery of the Guangdong-Hong Kong-Macao Greater Bay Area captured by the Landsat-8 Operational Land Imager (OLI) sensor. The imagery is orthorectified and geometrically corrected using a rational function model (RFM) using the satellite imagery's inherent rational polynomial coefficients (RPCs). The MODIS imagery is MOD09A1 data, which is atmospherically corrected surface reflectance data containing information in the red and blue bands.

[0046] S2 performs radiometric calibration on the pre-processed multispectral image to obtain apparent reflectance data of each band, including a blue band, a green band, a red band, and a near-infrared band;

[0047] In the embodiment of the present invention, calibration information is first obtained, including the maximum and minimum radiance values, the maximum and minimum DN values ​​of the pixels, the state parameters of the sensor, and the solar geometry information, including the solar zenith angle and the solar azimuth angle. The radiometric calibration of the satellite sensor detection signal converts the DN value of the image into the apparent reflectance of each band. Specifically, the apparent reflectance is calculated according to the calibration formula and the solar geometry information. The calibration formula is: L λ =(LMAX λ -LMIN λ ) / (Q calmax -Q calmin )×(Q cal -Q calmin )+LMIN λ

[0048] in:

[0049] L λ: Spectral radiation brightness obtained by the sensor, unit is W / (m 2 sr.μm);

[0050] Q cal : DN value corresponding to the pixel;

[0051] Q calmax and Q calmin : are the maximum and minimum DN values ​​of pixels respectively;

[0052] LMAX λ and LMIN λ : are the spectral radiances of the maximum and minimum DN values ​​of the corresponding pixels, in W / (m 2 sr.μm).

[0053] The conversion formula from spectral radiance to apparent reflectance is: ρ λ =π.L λ .d 2 / (ESUN λ .cosθ s )

[0054] in:

[0055] ρ λ : apparent reflectance, unitless;

[0056] π: Pi, approximately equal to 3.14159, no unit;

[0057] d: distance between the Earth and the Sun, astronomical unit;

[0058] ESUN λ : Average zenithal solar irradiance, unit W / (m 2 .μm);

[0059] θ s : Solar zenith angle, unit: degree.

[0060] S3 removes cloud pixels and cloud shadow pixels from the pre-processed multispectral image to obtain the image to be inverted;

[0061] In an embodiment of the present invention, step S3 includes:

[0062] S31 obtains the band value of each pixel of each multispectral image by calculation;

[0063] S32 uses long-time series of remote sensing samples to construct a clear sky background field containing a blue light band. Based on the clear sky background field containing a blue light band and the blue light band reflectance of each pixel in the multispectral image, the cloud detection threshold is calculated.

[0064] S33 uses the cloud detection threshold to detect clouds on the multispectral image, removes cloud and cloud shadow pixels, and obtains the image to be inverted.

[0065] In an embodiment of the present invention, a group of historical remote sensing samples is first obtained, which is a long time series. Cloud pixels are identified in the images of the remote sensing samples, and cloud shadow pixels are identified based on the flooding algorithm. Cloud and cloud shadow pixels are eliminated, and clear sky pixels are screened out to obtain T clear sky images. For pixels at the same position in these T clear sky images, the minimum value of the blue light band reflectance is obtained, and all minimum values ​​are arranged in the corresponding pixel arrangement order to obtain a clear sky background field containing the blue light band.

[0066] The band values ​​of each pixel (x, y) of each remote sensing data are calculated based on metadata, including: blue light band reflectance Band1 (x,y) , Green band reflectivity Band2 (x,y) , red light band reflectivity Band3 (x,y) , near-infrared band reflectivity Band4 (x,y) and the brightness temperature Band15 (x,y) .

[0067] The surface information Band1clear-sky(x,y) of the blue light band of each pixel in the clear sky background field is obtained. The blue light band reflectance of each pixel (x,y) of the remote sensing data is subtracted from Band1clear-sky(x,y) to obtain the difference. The difference is subjected to multi-scale retinal enhancement with color rejuvenation to obtain the reflectance difference. The cloud detection threshold is extracted from the reflectance difference image using the Otsu algorithm.

[0068] The cloud detection threshold is used to determine whether each pixel in the remote sensing data is a cloud pixel, and the clouds are identified. At the same time, the flooding algorithm is used to identify cloud shadows, and the clouds and cloud shadows in the remote sensing data are identified and removed.

[0069] S4 calculates the normalized difference vegetation index (NDVI) of each pixel in the image to be inverted, and classifies each pixel in the image to be inverted into a dark pixel area and a bright pixel area according to the normalized difference vegetation index according to a preset threshold;

[0070] In the embodiment of the present invention, step S4 includes:

[0071] S41 calculates the normalized vegetation index of each pixel based on the apparent reflectance of the red band and near-infrared band of each pixel in the multispectral image:

[0072] Among them, ρ RED is the apparent reflectivity in the red light band, ρ NIRis the apparent reflectance in the near-infrared band, and NDVI is the normalized difference vegetation index;

[0073] S42 marks the pixels with NDVI>0.55 on the image to be inverted as dark pixel areas;

[0074] S43 marks the pixels with NDVI≦0.55 and the corresponding pixels on the image to be inverted as bright pixel areas.

[0075] The present invention divides the inverted image into different categories according to a preset NDVI threshold. The data in the first category have NDVI values ​​greater than 0.55, indicating that the data in this category contain a large portion of densely vegetated areas, meaning that the scenes in these data can be considered dark pixel areas. The data in the second category have NDVI values ​​less than or equal to 0.55, indicating that the data in this category contain a large portion of complex surfaces, meaning that the scenes in these data can be considered bright pixel areas.

[0076] S5 estimates the surface reflectance ratio of the first red and blue bands using the dark target method in dark pixel areas. In bright pixel areas, based on the surface reflectance data of each band in MODIS imagery, a monthly red and blue band surface reflectance ratio library is constructed using the sub-minimum synthesis method through the Google Earth Engine platform.

[0077] In an embodiment of the present invention, in step S5, estimating the surface reflectance ratio of the red and blue bands using a dark target algorithm in the dark pixel area includes the following steps:

[0078] In the dark pixel area of ​​the image to be inverted, the surface reflectance of the red and blue bands is estimated using the shortwave infrared band apparent reflectance through a linear relationship. The calculation formula is:

[0079] Among them, ρ BLUE is the surface reflectivity in the blue light band, ρ RED is the surface reflectivity in the red light band, ρ SWIRW It is the TOA apparent reflectivity in the 2.1μm short infrared band;

[0080] Divide the surface reflectance of the blue band by the surface reflectance of the red band to obtain the surface reflectance ratio of the first red and blue bands The value is fixed at 2.

[0081] In bright pixel areas such as cities and bare soil, the red, blue and shortwave infrared bands used by the dark target algorithm usually do not satisfy a linear relationship, so the dark target algorithm cannot be used to invert aerosol optical depth for non-dark pixel surfaces. In bright pixel areas, the change in surface reflectance over time does not have as obvious seasonal and interannual changes as in dark pixel areas. Therefore, based on the assumption that the surface reflectance in bright pixel areas changes less over a certain period of time, the minimum synthesis technique can be used to construct a monthly-scale surface reflectance ratio library for the red and blue bands. Due to the influence of tall buildings in built-up areas, the minimum value of the surface reflectance may be incorrect, so the sub-minimum synthesis (second minimum synthesis) method is used to construct the surface reflectance ratio library.

[0082] Google Earth Engine (GEE) is a cloud-based platform for large-scale geospatial analysis and natural environment monitoring. It leverages Google's massive computing power to process petabytes of remote sensing imagery data. GEE archives a wide range of remote sensing imagery data, including MODIS, Landsat, and Sentinel.

[0083] In an embodiment of the present invention, in step S5, in the bright pixel area, a sub-minimum synthesis method is used to construct a monthly red and blue band surface reflectance ratio library, which includes the following steps:

[0084] Obtain MODIS images through the Google Earth Engine platform;

[0085] MODIS images are preprocessed, and multiple images of the same day in a long time series are synthesized into one image using the minimum value synthesis method, and the area that matches the time and space of the multispectral image is cropped;

[0086] Obtain the surface reflectance data for each month and use the sub-minimum composite method to construct the monthly red and blue band surface reflectance ratio library ρ(mon,wrs): ρ(mon,wrs)=2ndMin_composite(ρ1(mon,wrs),...ρ n (mon,wrs))

[0087] Among them, ρ(mon,wrs) represents the construction of the monthly red and blue band surface reflectance ratio image library; ρ1, ρ n An image representing the ratio of red and blue band surface reflectance from a month of MODIS imagery.

[0088] This paper uses atmospherically corrected surface reflectance data from the MOD09A1 archive of the Global Earth Observation (GEE). First, a monthly surface reflectance dataset for the Guangdong-Hong Kong-Macao Greater Bay Area is obtained through the GEE ImageCollection. Cloud and cloud shadow pixels are removed from the dataset using the map function using the image quality assurance (QA) file. A monthly red-blue band surface reflectance ratio library is constructed using subminimum synthesis techniques. Finally, the library is exported to Google Drive for download and used for AOD inversion in bright pixel areas.

[0089] Since MOD09A1 has a long revisit period and is often affected by clouds and cloud shadows, the minimum composite technique will result in many missing values ​​in the surface reflectance library. Therefore, we use the sub-minimum composite technique based on GEE to construct the ratio of the red and blue band surface reflectance library and downscale the resolution from 500m to 4km. The implementation principle is as follows: ρ(mon,wrs)=2ndMin_composite(ρ1(mon,wrs),...ρ n (mon,wrs))

[0090] Among them, ρ(mon,wrs) represents the construction of the monthly red and blue band surface reflectance ratio image library; ρ1, ρ n Represents the surface reflectance ratio image of the red and blue bands of the MOD09A1 image of the Guangdong-Hong Kong-Macao Greater Bay Area every month.

[0091] S6 simulates the apparent reflectivity detected by different satellite sensors under different conditions based on the atmospheric radiation transfer model and establishes an AOD lookup table;

[0092] In this embodiment of the present invention, the 6S atmospheric radiation transfer model is used to establish an AOD inversion lookup table, and its parameters are set as follows:

[0093] 1 Satellite zenith angle: multispectral image metadata acquisition, unit is degree;

[0094] 1. Solar zenith angle: Multispectral image metadata acquisition, unit is degree;

[0095] 1 relative azimuth: multispectral image metadata acquisition, unit is degree;

[0096] 17 AODs: 0.001, 0.01, 0.05, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.15, 1.4, 1.65, 2.0;

[0097] 2 atmosphere types: mid-latitude summer and mid-latitude winter;

[0098] 1 aerosol type: continental aerosol;

[0099] Output parameters: atmospheric radiation reflectivity, total transmittance of incident light from the top of the atmosphere to the ground, total transmittance in the direction of entering the satellite sensor's field of view, and hemispherical reflectivity downward from the bottom of the atmosphere;

[0100] A lookup table is generated for each multispectral image, the angle of which is fixed and can therefore be obtained from the multispectral metadata.

[0101] S7 uses the AOD lookup table to iteratively calculate the theoretical ratio of the red and blue band surface reflectance at different AOD values, and subtracts the theoretical ratio from the first red and blue band surface reflectance ratio and the monthly red and blue band surface reflectance ratio library, respectively, and takes the AOD value in the AOD lookup table corresponding to the theoretical ratio of the red and blue band surface reflectance with the smallest difference as the inversion result of the multispectral image;

[0102] When the land surface is assumed to be a uniform Lambertian surface, the apparent reflectance of the top layer of the atmosphere measured by satellite sensors can be expressed as:

[0103] in:

[0104] θs: solar zenith angle;

[0105] θv: satellite zenith angle;

[0106] relative azimuth;

[0107] apparent reflectance;

[0108] Atmospheric radiation reflectivity;

[0109] T g : total gas molecule transmittance;

[0110] T(θ s ) is the total transmittance of incident light from the top of the atmosphere to the ground;

[0111] T(θ v ): total transmittance in the direction of upward into the field of view of the satellite sensor;

[0112] S: downward hemispherical reflectivity of the bottom atmosphere, which can be calculated by the 6S radiation transfer model;

[0113] Surface reflectivity.

[0114] In the embodiment of the present invention, in step S7, iterative calculation is performed at different AOD values, and the variation of the different AOD values ​​is within the range of 0 to 2.0, with a step length of 0.01.

[0115] Based on the apparent reflectivity calculated above and the three parameters output by the 6S atmospheric radiation transfer model, the calculation is performed with an initial AOD value of 0.001 and a step size of 0.01. The red and blue band surface reflectivity under different AOD conditions is iteratively calculated using the following formula. The solution formula for the red and blue band surface reflectivity can be obtained by transforming the above formula:

[0116] According to the theoretical surface reflectance of the red and blue bands obtained above, calculate their ratio, and make the difference between it and the red and blue band surface reflectance ratio and the values ​​in the monthly red and blue band surface reflectance ratio database. The AOD value of the current iteration with the minimum absolute value of the difference is taken as the inversion result of the multispectral image;

[0117] The verification data comes from AERONET and SONET. SONET is an observation network established by the Chinese Academy of Sciences in collaboration with multiple Chinese universities and research institutes in typical areas of China. SONET stations are more evenly distributed across China, covering multiple provinces across the country. There are currently 23 stations. In the embodiment of the present invention, Level 1.5 observation data from the Guangdong-Hong Kong-Macao Greater Bay Area is selected to verify the accuracy of the AOD inversion results.

[0118] S8 uses the ground observation data of the site to verify the accuracy of the inversion results.

[0119] In the embodiment of the present invention, the accuracy of the inverted AOD values ​​is verified using AERONET and SONET site data. The evaluation indicators include the root mean square error (RMSE) and the expected error (EE), and the calculation formula is as follows: EE=±(0.05+0.2×τ0)

[0120] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any changes or substitutions that can be easily conceived by a person skilled in the art within the technical scope disclosed in the present invention should be included in the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be based on the scope of protection of the claims.

Claims

1. An aerosol optical thickness inversion method, characterized in that: The following steps are involved: S1 obtains multispectral images and corresponding MODIS images, and preprocesses the multispectral images; S2 performs radiation calibration on the pre-processed multispectral image to obtain apparent reflectance data of each band, wherein the band includes a blue light band, a green light band, a red light band, and a near infrared band; S3 removes cloud pixels and cloud shadow pixels from the pre-processed multispectral image to obtain the image to be inverted; S4 calculates the normalized vegetation index of each pixel in the image to be inverted, and classifies each pixel in the image to be inverted according to the normalized vegetation index according to a preset threshold value, into a dark pixel area and a bright pixel area; S5 estimates the surface reflectance ratio of the first red and blue bands using the dark target method in dark pixel areas; in bright pixel areas, based on the surface reflectance data of each band in MODIS images, the sub-minimum synthesis method is used to construct a monthly red and blue band surface reflectance ratio library through the Google Earth Engine platform; S6 simulates the apparent reflectivity detected by different satellite sensors under different conditions based on the atmospheric radiation transfer model and establishes an AOD lookup table; S7 uses the AOD lookup table to iteratively calculate the theoretical ratio of red and blue band surface reflectance at different AOD values, and subtracts the theoretical ratio from the first red and blue band surface reflectance ratio and the monthly red and blue band surface reflectance ratio library, and takes the AOD value of the AOD lookup table corresponding to the theoretical ratio of red and blue band surface reflectance with the smallest difference as the inversion result of the multispectral image; S8 uses the site ground observation data to verify the accuracy of the inversion results.

2. The aerosol optical thickness inversion method according to claim 1, characterized in that: In step S5, in the dark pixel area, the surface reflectance ratio of the first red and blue bands is estimated using a dark target algorithm, including the following steps: In the dark pixel area of ​​the image to be inverted, the short-wave infrared band apparent Reflectivity estimates the surface reflectivity of the red and blue bands. The calculation formula is: Among them, ρ BLUE is the surface reflectivity in the blue light band, ρ RED is the surface reflectivity in the red light band, ρ SWIRW It is the TOA apparent reflectivity in the 2.1μm short infrared band; Divide the surface reflectance of the blue light band by the surface reflectance of the red light band to obtain the surface reflectance ratio of the first red and blue bands The value is fixed at 2.

3. The aerosol optical thickness inversion method according to claim 2, characterized in that: In step S5, in the bright pixel area, a monthly red and blue band surface reflectance ratio library is constructed using the sub-minimum synthesis method, including the following steps: Obtain MODIS images through the Google Earth Engine platform; MODIS images are preprocessed, multiple images of the same day in a long time series are synthesized into one image by the minimum value synthesis method, and the area that matches the time and space of the multispectral image is cropped; Obtain the surface reflectance data for each month, and use the sub-minimum synthesis method to construct the monthly red and blue band surface reflectance ratio library ρ(mon,wrs): ρ(mon,wrs)=2ndMin_composite(ρ1(mon,wrs),...ρ n (mon,wrs)) Among them, ρ(mon,wrs) represents the construction of a monthly red and blue band surface reflectance ratio image library; ρ1, ρ n An image representing the ratio of red and blue band surface reflectance from MODIS images for one month.

4. The aerosol optical thickness inversion method according to claim 1, characterized in that: Step S4 includes: S41 calculates the normalized difference vegetation index NDVI of each pixel based on the apparent reflectance of the red band and near infrared band of each pixel in the multispectral image; S42 marks the pixels with NDVI>0.55 and the corresponding pixels on the image to be inverted as dark pixel areas; S43 marks the pixels with NDVI ≤ 0.55 as bright pixels on the image to be inverted. Pixel area.

5. The aerosol optical thickness inversion method according to claim 1, characterized in that: The preprocessing in step S1 includes orthophoto and geometric precision correction.

6. The aerosol optical thickness inversion method according to claim 1, characterized in that: Step S3 includes: S31 obtains the band value of each pixel of each multispectral image; S32 uses long-time series of remote sensing samples to construct a clear sky background field containing a blue light band, and calculates the cloud detection threshold based on the clear sky background field containing a blue light band and the blue light band reflectance of each pixel of the multispectral image; S33 performs cloud detection on the image to be inverted using the cloud detection threshold, removes cloud and cloud shadow pixels, and obtains the image to be inverted.

7. The aerosol optical thickness inversion method according to claim 1, characterized in that: In step S7, iterative calculation is performed at different AOD values, and the variation of different AOD values ​​is within the range of 0 to 2.

0.

8. The aerosol optical thickness inversion method according to claim 1, characterized in that: Step S8 includes: using AERONET and SONET site data to verify the accuracy of the AOD value obtained by inversion, and the evaluation indicators include root mean square error and expected error. The calculation formula is: EE=±(0.05+0.2×τ0) Among them, RMSE is the root mean square error, EE is the expected error, τ is the inverted AOD value, τ0 is the site AOD value, and n is the number of verification samples.

Images