阅读说明：本技术 一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法 (Aerosol optical thickness remote sensing inversion method considering surface non-Lambert characteristics ) 是由 田信鹏 高志强 于 2019-09-05 设计创作，主要内容包括：本发明涉及一种考虑地表非朗伯特性的气溶胶光学厚度遥感反演方法,包括获取遥感数据进行预处理；利用遥感图像地表反射率数据基于先验知识约束完成半经验核驱动模型的每日核系数估算；通过离散余弦变换的惩罚最小二乘估计时间序列平滑滤波算法构建单位时间内核系数地表先验知识数据集；建立气溶胶反演通用查找表；利用半经验核驱动模型结合核系数地表先验知识数据集,确定非朗伯体辐射传输前向模型的地表反射率参数；利用上述结果反演获取气溶胶光学厚度。本发明可获取高精度的空间分辨率为500m的气溶胶光学厚度产品,能提供空间相对更为连续的气溶胶分布,显示了在支持区域污染精细管控、重点城市污染传输通道监测及污染物溯源等领域的应用潜力。(The invention relates to an aerosol optical thickness remote sensing inversion method considering surface non-Lambert characteristics, which comprises the steps of obtaining remote sensing data for preprocessing; utilizing the remote sensing image earth surface reflectivity data to complete the estimation of the daily nuclear coefficient of the semi-empirical nuclear driving model based on the prior knowledge constraint; constructing a unit-time kernel coefficient earth surface prior knowledge data set by a punished least square estimation time sequence smoothing filtering algorithm of discrete cosine transform; establishing an aerosol inversion general lookup table; determining the earth surface reflectivity parameter of a non-Lambertian radiation transmission forward model by utilizing a semi-empirical kernel driving model in combination with a kernel coefficient earth surface prior knowledge data set; and (5) obtaining the optical thickness of the aerosol by utilizing the result inversion. The invention can obtain aerosol optical thickness products with high-precision spatial resolution of 500m, can provide aerosol distribution with relatively more continuous space, and shows application potential in the fields of supporting regional pollution fine control, monitoring of key urban pollution transmission channels, tracing of pollutants and the like.)

1. An aerosol optical thickness remote sensing inversion method considering surface non-Lambert characteristics is characterized by comprising the following steps:

Acquiring remote sensing data and preprocessing the remote sensing data;

constructing a ground surface prior knowledge data set, namely utilizing remote sensing image ground surface reflectivity data of a long-time sequence to complete daily kernel coefficient estimation of a semi-empirical kernel driving model based on prior knowledge condition constraint;

And (3) constructing a lookup table: establishing a general atmospheric parameter data set, and establishing an MODIS sensor aerosol inversion lookup table file through integral operation on data in the atmospheric parameter data set according to an MODIS spectral response function;

determining the surface reflectivity parameters, namely determining four surface reflectivity parameters in a non-Lambertian radiation transmission forward model by utilizing an RT L S semi-empirical nuclear driving model and combining a constructed nuclear coefficient surface priori knowledge data set SANIFD;

Aerosol optical thickness inversion: and (3) according to the imaging angle information, the earth surface reflectivity parameters and the atmospheric parameter information in the lookup table, combining the apparent reflectivity data obtained by preprocessing, and obtaining the optical thickness of the aerosol by inversion based on a non-Lambert body radiation transmission forenomial model.

2. The aerosol optical thickness remote sensing inversion method considering surface non-Lambert characteristics as claimed in claim 1, wherein the estimation of the semi-empirical kernel-driven model kernel coefficient is completed based on prior knowledge condition constraint, and the cost function is as follows: j (f) ═ k (Kf- ρ) ^TΣ^T(Kf-ρ)+(f-f₀)^TM^T(f-f₀)；

Where K is the three kernel matrices, f is the corresponding three kernel coefficient vectors, ρ is the directional reflectivity vector, Σ is the covariance matrix of the reflectivity error, f is the three kernel coefficient vector ₀And M is a mean vector and a covariance matrix of the kernel coefficient prior knowledge;

The above cost letter The least squares solution expression of numbers can be written as: f ═ K ^TΣ^-1K+M^-1)^-1(K^TΣ^-1ρ+M^-1f₀) The estimation of the nuclear factor can be completed by this formula.

3. The aerosol optical thickness remote sensing inversion method considering surface non-Lambertian characteristics as claimed in claim 1, wherein the kernel coefficients of the semi-empirical kernel-driven model comprise: isotropic scattering f _isoAnd volume scattering f _volAnd geometrical optical scattering f _geoThree kernel coefficients.

4. The method for inverting the aerosol optical thickness remote sensing considering the non-Lambertian characteristics of the earth's surface as recited in claim 1, wherein the unit time is every several days.

5. the aerosol optical thickness remote sensing inversion method considering surface non-Lambert characteristics, as claimed in claim 1, is characterized in that a nuclear coefficient surface prior knowledge data set is constructed by a DCT-PLS time series smoothing filter algorithm, and the specific implementation formula is as follows:

In the formula, | | | represents the euclidean norm; Representing a time series variable to be smoothed; w represents the optimal smoothing filter weight; s is a parameter for controlling the smoothness of the filtering, and as the value of s increases, the smoother will be, for the calculation of one-dimensional data, the penalty term P can typically be calculated using the L aplace differential.

6. The remote sensing inversion method for aerosol optical thickness considering surface non-Lambert characteristics as claimed in claim 1, wherein said general atmospheric parameter data set records the following data:

Different observation geometry data: 18 sun zenith angles, 15 observation zenith angles and 19 opposite azimuth angles;

Aerosol optical thickness AOD condition (in the range of 0.05-3.0) 16 aerosol optical thicknesses;

Atmospheric parameter information in the spectral range of 0.3-2.5 μm, including t _dd(θ_s)、t_dd(θ_v)、t_dh(θ_s)、t_hd(θ_v)、And rho ₀。

7. The aerosol optical thickness remote sensing inversion method considering surface non-Lambertian characteristics as claimed in claim 1, wherein the four surface reflectance parameters in the non-Lambertian radiation transmission forward model comprise: r is _dd、r_dh、r_hd、r_hhThe calculation formula is as follows:

Wherein, K _volAnd K _geoRespectively volume scattering nuclei and geometric optical nuclei, f _iso、f_volAnd f _geoRespectively representing the nuclear coefficient values of isotropic scattering, volume scattering and geometric optical scattering, and relating to the wave band lambda; r is the surface reflectance value, R when calculated separately _dd、r_dh、r_hd、r_hhThe reflection characteristics of the earth surface are shown as earth surface direction-direction reflectivity, direction-hemisphere reflectivity, hemisphere-direction reflectivity and hemisphere-hemisphere reflectivity.

8. The aerosol optical thickness remote sensing inversion method considering surface non-Lambertian characteristics as claimed in claim 1, wherein the non-Lambertian radiation transmission forenomial model is:

Wherein the content of the first and second substances,

Matrix T (theta) _s)、T(θ_v) And Are respectively defined as:

T(θ_s)＝[t_dd(θ_s) t_dh(θ_s)]

In the formula, theta _s、θ_v、Respectively a solar zenith angle, a sensor zenith angle, a solar azimuth angle, a sensor azimuth angle and a relative azimuth angle; rho _TOAIs the apparent reflectance; rho ₀Scattering reflectivity for an atmospheric path, representing an atmospheric path radiation term; The spherical albedo of the atmosphere; t and r represent transmittance and reflectance, respectively, and subscripts h and d represent hemisphere (diffuse) and direction (direct); thus, t _hd(θ_v) Representing extinction of the atmosphere to earth surface diffusion radiation in a satellite observation direction, namely the uplink diffusion transmittance of the atmosphere; t is t _dd(θ_v) Representing extinction of the atmosphere on the earth surface reflected direct radiation in the satellite observation direction, namely the upward direct transmittance of the atmosphere; likewise, t _dh(θ_s) Characterizing the extinction of atmospheric air to the diffuse sky radiation in the incident direction, i.e. the down-atmospheric diffuse transmittance; t is t _dd(θ_s) The extinction of the atmosphere on the incident direction to the direct solar radiation is achieved, namely the downward direct transmittance of the atmosphere; r is _dd、r_dh、r_hd、r_hhThe reflection characteristics of the earth surface are shown as earth surface direction-direction reflectivity, direction-hemisphere reflectivity, hemisphere-direction reflectivity and hemisphere-hemisphere reflectivity.

9. The aerosol optical thickness remote sensing inversion method considering surface non-Lambert characteristics, as claimed in claim 1, characterized in that the preprocessing of the remote sensing data includes cloud detection, ice/snow and water body identification, and an improved dynamic threshold method is adopted, specifically expressed as:

1) Firstly, using a land surface priori knowledge data set to calculate and obtain the land surface reflectance value of each wave band: Then, the cloud pixel identification is realized by using different wave band cloud detection dynamic threshold models, and the formula is expressed as follows:

2) Because the geometric optical kernel and the volume scattering kernel in the adopted semi-empirical kernel driving model are not suitable for the reflectivity calculation of ice/snow and water surface types, the ice/snow and the water need to be subjected to mask removal in the aerosol inversion, which is specifically as follows:

2 nd and 5 th wave band apparent reflectivity value after remote sensing image preprocessing According to the formula Calculating and identifying ice/snow; surface reflectance values of 1 st and 2 nd wave bands preprocessed by using remote sensing images According to the formula And calculating and identifying the water body.

Technical Field

The invention relates to the field of atmospheric remote sensing, in particular to an aerosol optical thickness remote sensing inversion method considering surface non-Lambert characteristics.

Background

Aerosol, an important component of the earth's climate system, is one of three major active components in the atmosphere, whose quantity and properties vary greatly on both time and space scales, making it difficult to describe in a fixed fashion. Aerosol Optical Depth (AOD) is one of important Optical parameters for representing extinction effect of Aerosol particles, and has important indication effect on aspects of terrestrial gas system radiation balance, global climate change, atmospheric environment, human health and the like (Zhanying et al, 2013; Wangzhong et al, 2016). Because atmospheric aerosol has large space-time variability, satellite remote sensing becomes an effective means for providing reliable aerosol distribution information. The atmospheric layer top radiation information received by the satellite sensor is the final result of the co-action and mutual coupling of the earth surface and the atmosphere, and the atmospheric contribution is very limited compared with the earth surface, especially in high-reflectivity areas of cities. The remote sensing quantification process is to distinguish two unknowns from the single information, and the essence of the remote sensing quantification process is a ill-posed inversion problem. In actual aerosol remote sensing, the aerosol type and the earth surface reflectivity parameter need to be determined through priori knowledge, and the uncertainty of the aerosol type and the earth surface reflectivity parameter is also a main source of aerosol inversion errors. The development and improvement of the aerosol inversion algorithm are also the two aspects of the assumption around the aerosol type and the estimation of the earth surface reflectivity.

in recent years, many scholars at home and abroad make much effort on the aspect of estimation of surface direction reflection in Aerosol inversion, Thomas et al propose an ORAC (Oxford-RA L Aerosol and Cloud) algorithm to realize Aerosol parameter inversion of ATSR-2 and AATSR data, the algorithm uses BRDF/Albedo products to construct a reflectivity prior knowledge base to realize ground gas separation, Xue et al propose an Aerosol collaborative inversion algorithm considering BRDF effect based on AATSR satellite down-point and forward two-angle observation data, TSGuing et al propose a two-star collaborative inversion algorithm of surface reflectivity by using MODIS data, further improve the accuracy of the air Aerosol on the surface, on the basis of the assumption that BRDF is stable, the TSGuing et al still realize the ground surface reflection ratio and the non-proportional inversion algorithm of the ground surface reflectivity, and basically solve the problem of the non-proportional inversion algorithm of the non-atmospheric reflection index.

Disclosure of Invention

In order to solve the uncertainty of the estimation of the earth surface reflectivity caused by the assumption of the earth surface lambertian reflection in the traditional aerosol algorithm, the invention provides an aerosol inversion algorithm coupled with the earth surface two-way reflection characteristic, and finally, an aerosol inversion algorithm with stronger universality and capable of being operated in a business mode is constructed.

In order to achieve the purpose, the technical scheme adopted by the invention is as follows: an aerosol optical thickness remote sensing inversion method considering surface non-Lambert characteristics comprises the following steps:

Acquiring remote sensing data and preprocessing the remote sensing data;

constructing a ground surface prior knowledge data set, namely utilizing remote sensing image ground surface reflectivity data of a long-time sequence to complete daily nuclear coefficient estimation of a semi-empirical nuclear driving model based on prior knowledge condition constraint;

And (3) constructing a lookup table: establishing a general atmospheric parameter data set, and establishing an MODIS sensor aerosol inversion lookup table file through integral operation on data in the atmospheric parameter data set according to an MODIS spectral response function;

determining the surface reflectance parameters, namely determining four surface reflectance parameters in a non-Lambert body radiation transmission forward model by utilizing an RT L S semi-empirical nuclear driving model and combining a constructed nuclear coefficient surface priori knowledge data set SANIFD;

Aerosol optical thickness inversion: and (3) according to the imaging angle information, the earth surface reflectivity parameters and the atmospheric parameter information in the lookup table, combining the apparent reflectivity data obtained by preprocessing, and obtaining the optical thickness of the aerosol by inversion based on a non-Lambert body radiation transmission foreterm model.

The estimation of the kernel coefficient of the semi-empirical kernel-driven model is completed based on the prior knowledge condition constraint, and the adopted cost function Comprises the following steps: j (f) ═ k (Kf- ρ) ^TΣ^T(Kf-ρ)+(f-f₀)^TM^T(f-f₀)；

Where K is the three kernel matrices, f is the corresponding three kernel coefficient vectors, ρ is the directional reflectivity vector, Σ is the covariance matrix of the reflectivity error, f is the three kernel coefficient vector ₀And M is a mean vector and a covariance matrix of the kernel coefficient prior knowledge;

The least squares solution expression of the above cost function can be written as:

f＝(K^TΣ^-1K+M^-1)^-1(K^TΣ^-1ρ+M^-1f₀) The estimation of the nuclear factor can be completed by this formula.

The kernel coefficients of the semi-empirical kernel driving model include: isotropic scattering f _isoAnd volume scattering f _volAnd geometrical optical scattering f _geoThree kernel coefficients.

The unit time is every several days.

the method comprises the following steps of constructing a kernel coefficient earth surface priori knowledge data set through a DCT-PLS time sequence smoothing filtering algorithm, wherein a specific implementation formula is as follows:

In the formula, | | | represents the euclidean norm; Representing a time series variable to be smoothed; w represents the optimal smoothing filter weight; s is a parameter for controlling the smoothness of the filtering, and as the value of s increases, the smoother will be, for the calculation of one-dimensional data, the penalty term P can typically be calculated using the L aplace differential.

The general atmospheric parameter dataset records the following data:

Different observation geometry data: 18 sun zenith angles, 15 observation zenith angles and 19 opposite azimuth angles;

Aerosol optical thickness AOD condition (in the range of 0.05-3.0) 16 aerosol optical thicknesses;

Atmospheric parameter information in the spectral range of 0.3-2.5 μm, including t _dd(θ_s)、t_dd(θ_v)、t_dh(θ_s)、t_hd(θ_v)、And rho ₀。

The four surface reflectivity parameters in the non-lambertian radiation transmission forward model include: r is _dd、r_dh、r_hd、 r_hhThe calculation formula is as follows:

Wherein, K _volAnd K _geoRespectively volume scattering nuclei and geometric optical nuclei, f _iso、f_volAnd f _geoRespectively representing isotropic scattering, volume scattering and geometric optical scattering nuclear coefficient values which are related to a wave band lambda; r is the surface reflectance value, R when calculated separately _dd、r_dh、r_hd、r_hhThe reflection characteristics of the earth surface are shown as earth surface direction-direction reflectivity, direction-hemisphere reflectivity, hemisphere-direction reflectivity and hemisphere-hemisphere reflectivity.

The non-lambertian radiation transmission antecedent model is:

Wherein the content of the first and second substances,

Matrix T (theta) _s)、T(θ_v) And Are respectively defined as:

T(θ_s)＝[t_dd(θ_s)t_dh(θ_s)]

In the formula, theta _s、θ_v、Respectively a solar zenith angle, a sensor zenith angle, a solar azimuth angle, a sensor azimuth angle and a relative azimuth angle; rho _TOAIs the apparent reflectance; rho ₀Scattering reflectivity for an atmospheric path, representing an atmospheric path radiation term; The spherical albedo of the atmosphere; t and r represent transmittance and reflectance, respectively, and subscripts h and d represent hemisphere (diffuse) and direction (direct); thus, t _hd(θ_v) Representing extinction of the atmosphere to earth surface diffusion radiation in a satellite observation direction, namely the uplink diffusion transmittance of the atmosphere; t is t _dd(θ_v) Representing extinction of the atmosphere on the earth surface reflected direct radiation in the satellite observation direction, namely the upward direct transmittance of the atmosphere; likewise, t _dh(θ_s) Characterizing the extinction of atmospheric air to the diffuse sky radiation in the incident direction, i.e. the down-atmospheric diffuse transmittance; t is t _dd(θ_s) The extinction of the atmosphere on the incident direction to the direct solar radiation is achieved, namely the downward direct transmittance of the atmosphere; r is _dd、r_dh、 r_hd、r_hhThe reflection characteristics of the earth surface are shown as earth surface direction-direction reflectivity and direction-hemisphere reflectivity Hemispherical-directional reflectivity and hemispherical-hemispherical reflectivity.

The method comprises the following steps of preprocessing remote sensing data, including cloud detection, ice/snow and water body identification, and adopting an improved dynamic threshold method, wherein the method is specifically represented as follows:

1) Firstly, using a land surface priori knowledge data set to calculate and obtain the land surface reflectance value of each wave band: Then, the cloud pixel identification is realized by using different wave band cloud detection dynamic threshold models, and the formula is expressed as follows:

2) Because the geometric optical kernel and the volume scattering kernel in the adopted semi-empirical kernel driving model are not suitable for the reflectivity calculation of ice/snow and water surface types, the ice/snow and the water need to be subjected to mask removal in the aerosol inversion, which is specifically as follows:

2 nd and 5 th wave band apparent reflectivity value after remote sensing image preprocessing According to the formula Calculating and identifying ice/snow; using remote sensing images Pretreated surface reflectance values of 1 st and 2 nd wave bands According to the formula And calculating and identifying the water body.

The invention has the beneficial effects that:

The invention can obtain aerosol optical thickness products with high-precision spatial resolution of 500m, can provide aerosol distribution information with relatively more continuous space, and shows application potential in the fields of supporting regional pollution fine control, monitoring of key urban pollution transmission channels, tracing of pollutants and the like. To the accomplishment of the foregoing and related ends, one or more aspects of the invention comprise the features hereinafter fully described and particularly pointed out in the claims.

Drawings

FIG. 1 is a schematic flow chart of an aerosol optical thickness inversion method of the present invention.

FIG. 2 is a schematic diagram of a process for constructing a prior knowledge data set of the earth's surface according to the present invention.

Figure 3 is an exemplary diagram of the data layers of the aerosol product of the present invention.

Figure 4 is a graph comparing the optical thickness of the aerosol obtained by the present invention with the spatial distribution of a MODIS aerosol product.

FIG. 5 is a cross-comparison verification chart of AOD inversion results and MODIS aerosol products (a: DT algorithm, b: DB algorithm, c: DTDB algorithm, d: MAIAC algorithm) on the AERONET site of the Xianghe vegetation.

FIG. 6 is a cross-comparison verification chart of AOD inversion results and MODIS aerosol products (a: DT algorithm, b: DB algorithm, c: DTDB algorithm, d: MAIAC algorithm) at AERONET sites in Beijing city.

Detailed Description

The present invention now will be described more fully hereinafter with reference to the accompanying drawings, in which some of the details that are pertinent to the practice of the invention will be described. This invention may, however, be embodied in many different forms and these embodiments are provided so that this invention will be thorough and complete, and will fully convey the scope of the invention to those skilled in the art.

The method adopts a rapid and accurate radiation transmission model considering the anisotropic reflection characteristic of the earth surface, the model divides the radiation process into an uplink component and a downlink component, and each component comprises a direct radiation part and a diffusion part. The division can clearly distinguish and determine each reflectivity component of the incident and surface reflection radiation, and the apparent reflectivity calculation formula is as follows:

Wherein the content of the first and second substances,

Matrix T (theta) _s)、T(θ_v) And Are respectively defined as:

T(θ_s)＝[t_dd(θ_s) t_dh(θ_s)]

In the formula, theta _s、θ_v、Respectively the sun zenith angle and the sensor A zenith angle, a solar azimuth angle, a sensor azimuth angle and a relative azimuth angle; rho _TOAIs the apparent reflectance; rho ₀Scattering reflectivity for an atmospheric path, representing an atmospheric path radiation term; The spherical albedo of the atmosphere; t and r represent transmittance and reflectance, respectively, and subscripts h and d represent hemisphere (diffuse) and direction (direct); thus, t _hd(θ_v) Representing extinction of the atmosphere to earth surface diffusion radiation in a satellite observation direction, namely the uplink diffusion transmittance of the atmosphere; t is t _dd(θ_v) Representing extinction of the atmosphere on the earth surface reflected direct radiation in the satellite observation direction, namely the upward direct transmittance of the atmosphere; likewise, t _dh(θ_s) Characterizing the extinction of atmospheric air to the diffuse sky radiation in the incident direction, i.e. the down-atmospheric diffuse transmittance; t is t _dd(θ_s) The extinction of the atmosphere on the incident direction to the direct solar radiation is achieved, namely the downward direct transmittance of the atmosphere; r is _dd、r_dh、 r_hd、r_hhThe reflection characteristics of the earth surface are shown as earth surface direction-direction reflectivity, direction-hemisphere reflectivity, hemisphere-direction reflectivity and hemisphere-hemisphere reflectivity.

since the estimation of the surface parameters and the atmospheric parameters are independent, the estimation of each reflectivity parameter of the surface can be realized by a semi-empirical kernel driving model, the atmospheric parameters can be acquired from the L UT constructed in advance by radiation transmission codes according to the observation geometrical conditions, and the acquisition method of the two parameters is explained in detail in the following inversion process.

Fig. 1 is a schematic flow chart of an aerosol optical thickness inversion method according to the present invention. The process is as follows:

(1) obtaining remote sensing data, wherein the remote sensing data required in the algorithm construction and product production process mainly comprise 1) apparent reflectivity data (MOD02 HKM) with MODIS L evel-1 spatial resolution of 500m, 2) MODIS L evel-1 geometric position data (MOD03), 3) surface reflectivity product data (MOD09GA) with MODIS L evel-2 spatial resolution of 500m and BRDF/albedo product data (MCD43A1), 4) vegetation index product data (MOD13A1) with MODIS L evel-3 spatial resolution of 500m, and 5) MODIS L evel-2 aerosol product data (MOD04, MCD19A2), wherein all the data can be freely downloaded from L AADS DAAC (web of Lahttps: Laeseoss.

(2) Data processing: preprocessing MOD02HKM data such as geometric correction and cloud/ice and snow detection is carried out, and apparent reflectivity data of an aerosol inversion region can be obtained.

(3) Constructing a data set of the prior knowledge of the earth surface: in the algorithm, a ground surface priori knowledge data set is three nuclear system values of isotropic scattering, volume scattering and geometric optical scattering in a semi-empirical nuclear driving model, and the construction process comprises two steps: 1) utilizing MOD09GA data of a long-time sequence and based on a priori knowledge conditional constraint algorithm to realize f in a semi-empirical core driving model (formula (8)) _iso、f_volAnd f _geothree kernel coefficient estimation (daily products), 2) utilizing the day-by-day kernel coefficient data estimated in the first step, and based on the punishment least Square estimation of Discrete Cosine Transform (DCT-PLS) time series smoothing filter algorithm to synthesize and construct f of finally every 8 days _iso、f_volAnd f _geoThe construction process of the nuclear coefficient Surface prior knowledge data set (Surface ANIsotropic FactorDatabase, SANIFD) is shown in FIG. 2.

(4) And (3) constructing a lookup table: a general atmospheric parameter data set is constructed based on an atmospheric radiation transmission program (MODTRAN 5), and an aerosol inversion lookup table file is constructed by combining a spectral response function of an MODIS sensor.

(5) Determination of surface reflectivity: and (4) determining each surface reflectivity parameter in the non-Lambert body radiation transmission forward model by combining a semi-empirical kernel driving model based on the surface prior knowledge data set constructed in the step (3).

(6) Aerosol optical thickness inversion: and (3) according to the imaging angle information, the earth surface reflectivity and the atmospheric parameter information in the lookup table, combining the apparent reflectivity data obtained in the step (2), and obtaining the optical thickness of the aerosol in an inversion mode.

The specific inversion process is as follows:

1) Geometric correction and apparent reflectivity calculation. Firstly, a longitude and latitude data set in MOD03 data is used as a control point, and actual longitude and latitude corresponding to MOD02HKM are calculated through pixel-by-pixel interpolation; then, geometric correction is performed for each scan band data; and finally, splicing all the scanning band correction results to finish geometric correction, and simultaneously realizing radiometric calibration of data to obtain a real apparent reflectivity data image.

2) And (5) cloud detection. The invention adopts an improved dynamic threshold cloud detection algorithm to carry out cloud detection on MODIS data, and the basic process is as follows: firstly, using a land surface priori knowledge data set to calculate and obtain the land surface reflectance values of all wave bands in formulas (2) to (4) Then, simulating the apparent emissivity corresponding to the cloud pixel by using a different-waveband cloud detection dynamic threshold model Finally, real apparent reflectivity data is utilized And comparing the cloud pixel with the simulated apparent reflectivity to realize the identification of the cloud pixel, wherein the formula is expressed as follows:

3) Ice/snow and water body identification. As the geometric optical kernel and the volume scattering kernel in the adopted semi-empirical kernel driving model are not suitable for the reflectivity calculation of ice/snow and water body surface types, the surface type mask needs to be removed in the aerosol inversion. The algorithm used utilizes the apparent reflectivity of MODIS bands 2 and 5 to identify ice/snow (equation (6)). Similarly, the apparent reflectivities of the 1 st and 2 nd wave bands (formula (7)) are used for identifying the water body image element:

4) And (5) constructing a lookup table. In order to realize the aerosol inversion of different sensors, a General Atmospheric parameter data set (GAtmParas) is constructed by using a MODTRAN 5 Atmospheric radiation transmission program, and the data set records the Atmospheric parameter information of different observation geometries (18 solar zenith angles, 15 observation zenith angles and 19 opposite azimuth angles) and a spectral range of 0.3-2.5 μm under the AOD condition (16 aerosol optical thicknesses in a range of 0.05-3.0), including t _dd(θ_s)、t_dd(θ_v)、 t_dh(θ_s)、t_hd(θ_v)、And rho ₀Six atmospheric parameters. According to the MODIS spectral response function, the atmospheric parameter data set is utilized, and the construction of an aerosol inversion lookup table of the MODIS sensor can be completed through integral operation.

5) And (4) determining the surface reflectivity. In actual calculation, r in formula (1) can be calculated using a semi-empirical core-driven model (formula (8)) _dd、r_dh、r_hd、r_hhFour reflectivities.

Wherein R is the surface reflectance value, K _volAnd K _geoRespectively volume scattering nuclei and geometric optical nuclei, f _iso、 f_volAnd f _geoThe nuclear coefficient values of isotropic scattering, volume scattering and geometric optical scattering are respectively expressed and are related to the wave band lambda.

furthermore, experiments have shown that equation (9) can be well fitted when calculating the albedo of black and white space using the RossThick volume scattering kernel and the L iSearser geometric optical kernel integration, and therefore, the present invention also uses this equation to determine the reflectivity of each earth surface, with polynomial coefficients as shown in Table 1:

TABLE 1 polynomial h of formula (9) _k(θ)＝g_0k+g_1kθ²+g_2kθ³Coefficient value

6) Production of aerosol products. Firstly, according to the imaging geometric condition of remote sensing data, searching and calculating corresponding atmospheric parameters through the lookup table constructed in the step (4); then, by utilizing the earth surface reflectivity information determined in the step (5) and combining the formula (1), calculating and obtaining an apparent reflectivity vector under each aerosol optical thickness condition; and (3) finally, calculating the corresponding optical thickness of the aerosol by using the observed value of the real apparent reflectivity calculated in the step (1). An example of an aerosol optical thickness product produced by the present invention is shown in figure 3.

In order to evaluate the performance of the inversion method, the optical aerosol thickness obtained by the method of the present invention and the spatial distribution of each aerosol product of MODIS in beijing area are analyzed by cross-contrast, as shown in fig. 4. It can be seen that the AODs obtained by the algorithms have good consistency in spatial distribution, and are all represented by high AOD load in southwest urban areas and low AOD value in vegetation areas. The new algorithm provided by the invention can obtain AOD results in regions with low vegetation reflectivity, high urban reflectivity and the like, has good spatial continuity and spatial resolution up to 500m, and can provide more detailed information than that of MODIS aerosol products.

In addition, in order to evaluate the accuracy of the inversion algorithm more precisely, the aerosol results obtained by the new algorithm of the invention are respectively compared and analyzed with MODIS C6.1 DT, DB, DTDB algorithm products and C6 MAIAC algorithm products by using ground observation data of two AERONET sites of Xiang river (XHS) and Beijing (BJS) in 2013-2016, as shown in FIGS. 5 and 6. The aerosol and MAIAC obtained by the new algorithm have higher precision in cities and vegetation, and the specific gravity (PWE) in Expected errors at a BJS/XHS site is 74.9/82.1% and 73.1/81.2% respectively, so that the high estimation of the MODIS DT, DB and DTDB algorithm on the AOD can be effectively improved, the AOD reflection precision is improved, particularly, the vegetation area performance is more obvious, and the specific gravity (PAE) in the Expected errors is reduced by about 15%. Compared with the MAIAC algorithm, the new algorithm has a higher correlation coefficient and a smaller RMSE value with the measured data. In addition, the AOD spatial resolution obtained by the new algorithm is 500m, so that more detailed information than that of MODIS products can be provided, and the method is more beneficial to environmental monitoring and the like of urban areas.

In conclusion, the aerosol optical thickness information with high precision and high spatial resolution in the city can be obtained by the method, and the application potential in the fields of supporting the fine control of regional pollution, monitoring of key urban pollution transmission channels, tracing of pollutants and the like is displayed.

While the foregoing is directed to the preferred embodiment of the present invention, it will be understood by those skilled in the art that various changes and modifications may be made without departing from the spirit and scope of the invention as defined in the appended claims.