CN104569952B - Aerosol optical thickness temporal-spatial distribution inversion method without mid-infrared channel sensor - Google Patents

Abstract

The invention discloses an aerosol optical thickness temporal-spatial distribution inversion method without a mid-infrared channel sensor. The method comprises the following steps: firstly eliminating cloud pixels, water body pixels and ice and ice and snow pixels in channels, then estimating the surface reflectance of the two channels by using a clear sky zenith reflectance image of a time sequence; meanwhile, establishing a lookup table in an offline manner on the basis of an aerosol mode, the combination of different incidence and observation attitudes and the surface reflectance; finally inverting the aerosol optical thickness temporal-spatial distribution result on the basis of the zenith reflectance data of the time sequence and the determined surface reflectance according to the block adjustment value. According to the method disclosed by the invention, the influences of an observation zenith angle and the length of a sliding time window are simultaneously considered in a multi-day combined observation method, and the precision of the obtained surface reflectance result is particularly better in East Asian regions with serious air pollution.

Description

Inversion method of spatio-temporal distribution of aerosol optical depth without mid-infrared channel sensor

technical field

The invention relates to a method for retrieving the time-space distribution of aerosol optical thickness based on a sensor without a mid-infrared channel, in particular to a method for retrieving the time-space distribution of aerosol optical thickness in East Asia.

Background technique

The main difficulties in retrieval of aerosols over land come from calibration, cloud removal, accurate separation of ground-atmosphere contributions, and determination of aerosol models. Among them, accurate separation of aerosol contributions from sensor signals requires accurate estimation of surface reflectance. The estimation error of the surface reflectance is considered to be one of the main sources of error in the inversion of aerosol optical depth.

At present, the more popular algorithm for terrestrial aerosol optical depth inversion using reflectance spectral intensity is the DT (Dark Target) method. Kaufman et al. obtained a certain linear relationship between the surface reflectance of the dark surface with dense vegetation in the three channels of red, blue and mid-infrared through a large number of aircraft test data, and applied this relationship to the inversion of MODIS land aerosol optical depth (Kaufman et al. et al., 1997). Levy et al. (2007) updated the MODIS terrestrial aerosol retrieval algorithm, using the surface reflectance relationship of three channels of 0.47 μm, 0.66 μm and 2.12 μm (the relationship between the surface reflectance of different channels varies with the scattering angle and vegetation index) to retrieve the aerosol optical depth, and at the same time obtain the surface reflectance and aerosol thickness ratio.

For sensors without mid-infrared channel settings, the DT method is powerless, such as the CCD sensor carried by the environmental star, the AVHRR sensor carried by the NOAA series satellites, and the VISSR sensor carried by the FY-2 series geostationary satellites. An early solution is to use multiple images of the same area in different time phases, select the day when the image only has the influence of background aerosols as a sunny day, and assume that the surface reflectance of the sunny day and the polluted day are unchanged, and use two days of images difference in inversion of aerosol properties (Fraser et al. 1984; Kaufman et al. 1990). Or another approach directly assumes that the surface albedo of a dark target is a constant value (Soufflet et al. 1997). Subsequent research improved the above method, assuming that the surface albedo was constant during this period, and selected clear days based on the time series of AVHRR images. Since AVHRR is a wide-range sensor, the satellite observation attitude changes greatly during the 8-day revisit period. Therefore, the improved method also considers the influence of observation attitude (mainly the satellite zenith angle) on the surface emissivity estimation (Hauser et al. 2005; Riffler et al. 2010). However, applying the above method to the air-polluted East Asian region often leads to larger inversion errors because it is more difficult to find clear days.

Especially considering that the air pollution in East Asia is much more serious than that in Europe and North America, the probability of finding a sunny day is small, so it is necessary to increase the length of the sliding time window for selecting sunny days. However, the longer the time window is, the more difficult it is to satisfy the assumption that the surface albedo is constant, and the greater the variation range of the albedo may be.

Contents of the invention

In view of the problems existing in the above-mentioned prior art, the present invention introduces the theory of block adjustment in the field of surveying and mapping for the first time, and extends the spatial scale to the time scale. Based on the time series remote sensing data, the aerosol in East Asia is inverted by using the block adjustment theory of time-spatial scale Spatiotemporal distribution results of optical thickness. At the same time, the effects of the satellite zenith angle and the length of the sliding time window on determining the surface reflectivity are considered.

In order to achieve the above-mentioned purpose, when using the aerosol optical thickness space-time distribution inversion method of the present invention to determine the surface reflectance using a sensor without a mid-infrared channel, the length of the satellite sliding time window should not be too long to avoid changes in the surface reflectance, nor should it be too short to avoid No sunny day can be found, so the length of the sliding time window is set to 30 - 61 days, the current pixel is set to be the current processing pixel, and the current satellite zenith angle is the satellite zenith angle of the current pixel. The specific steps are as follows:

(1) Subtract 10º from the current satellite zenith angle as the current lower limit of the angle, and add 10º as the current upper limit of the angle;

(2) Compare the current upper and lower limits of the angle with ±90º respectively, and if it exceeds, set the exceeded current upper or lower limit to 90º or -90º;

(3) Initially set the length of the current sliding time window to a minimum of 30 days;

(4) Select the zenith reflectivity of channel 1 whose satellite zenith angle is between the upper and lower limits of the current angle in the current sliding time window, and the number of selected zenith reflectivities must be greater than 7;

(5) When the conditions are met:

1) The subdark albedo is not equal to the zenith albedo for the current day; or

2) The subdark reflectance is equal to the zenith reflectance of the current day, and the current time period is equal to 61 days,

Use the sub-dark reflectance of the selected zenith reflectance as the background reflectance of the current pixel;

(6) When the background reflectance of the current pixel is not obtained, the current time period is increased by 2 days, and the current day remains unchanged; repeat steps (4), (5) and (6) until the background reflectance of the current pixel is obtained or The current time period exceeds 61 days;

(7) If the background reflectance of the current pixel is still not obtained, skip to step (1), change the upper and lower limits of the angle, and cycle from step (1) to step (7) until the background reflectance of the current pixel is obtained, or The upper limit exceeds 90º, and the lower limit exceeds -90º;

(8) If the background reflectance of channel 1 of the current day is found, the background reflectance of channel 2 of the corresponding sunny day is the background reflectance of channel 2 of the current day;

(9) Since there are still aerosol effects on the background day, it is necessary to remove the background aerosol effect of the obtained image. Assume that the aerosol optical thickness value at 550nm on the background day is 0.05, and the aerosol model is continental. Correction is performed to obtain the surface reflectance of channel 1 and channel 2.

Further, the aerosol optical depth inversion method, the specific steps are as follows:

1) Eliminate the cloud pixels, water body pixels and ice and snow pixels in the channel 1 and channel 2, and only keep the sunny pixels, and then filter out the ocean water bodies for the sunny pixels, and keep the land and inland water body pixels; for East Asia In the spring and winter of the region, keep the pixels to further filter out the ice and snow pixels;

2) Establish a look-up table offline based on the aerosol model, different combinations of incident and observation attitudes, and surface reflectivity;

3) Based on the time series zenith albedo data and the determined surface albedo, according to the block adjustment value, invert the spatio-temporal distribution of aerosol optical thickness; linearly interpolate the lookup table to the actual inversion input, when interpolating, first Find the geometric attitude angle of the actual inversion input, the zenith albedo corresponding to the given input near the surface albedo, and then follow the order of the albedo, relative azimuth, satellite zenith angle, and solar zenith The interpolation is carried out at the first stage, and the interpolation result of this stage is returned to the upper stage to continue interpolation until the corresponding aerosol optical depth is obtained when the actual geometric observation and surface reflectance are obtained.

Further, the sensors used are AVHRR sensors, CCD sensors, and VISSR sensors.

Description of drawings

Figure 1 is a comparison map of the surface reflectance of Beijing Station;

Figure 2 is a comparison chart of surface albedo at Xinglong Station;

Fig. 3 is a comparison chart of the present invention and AERONET site ground-based observation in East Asia;

Fig. 4 is the comparative scatter density diagram of the results of the present invention and MYD04 aerosol optical depth products in East Asia;

The AVHRR aerosol optical thickness map inverted by the method of the present invention in Fig. 5a;

Figure 5b is the MODIS dark target aerosol optical depth map;

Figure 5c is a true-color composite image of 1km MYD02 data;

Figure 5d is a comparison of AODHRR AOD values, MODIS dark target AOD values and AERONET AOD observations matched with their respective transit times.

detailed description

Specific embodiments of the present invention will be described below in conjunction with the accompanying drawings.

This embodiment takes the AVHRR sensor as an example to describe the aerosol optical depth inversion process and the method for determining the surface reflectance in the process, as well as the comparative verification of the obtained surface reflectance and aerosol optical depth.

Surface reflectivity determination method

For the determination of the surface albedo, due to the short length of the sliding time window, no sunny day will be found, and the length of the time window is too long, which makes it difficult to satisfy the assumption of constant surface albedo. This algorithm assumes that the minimum length of the sliding time window is 30 days, and the longest is no more than 61 days. The current pixel is the current processing pixel, and the current satellite zenith angle is the satellite zenith angle of the current pixel. Specific steps are as follows:

(1) Subtract 10º from the current satellite zenith angle as the current lower limit of the angle, and add 10º as the current upper limit of the angle;

(2) Compare the current upper and lower limits of the angle with ±90º respectively, if it exceeds, set the exceeded current upper (lower) limit to 90º (-90º);

(3) Initially set the length of the current sliding time window to a minimum of 30 days;

(4) Select the zenith reflectivity of channel 1 whose satellite zenith angle is between the upper and lower limits of the current angle in the current sliding time window, and the number of selected zenith reflectivities must be greater than 7;

(5) When one of the following conditions is met, the sub-dark reflectance of the selected zenith reflectance is used as the background reflectance of the current pixel: 1) The sub-dark reflectance is not equal to the current day's zenith reflectance; 2 ) The subdark reflectance is equal to the zenith reflectance of the current day, and the current time period is equal to 61 days;

(6) If the background reflectance of the current pixel is not obtained, the current time period is increased by 2 days, and the current day remains unchanged, and steps (4), (5) and (6) are repeated until the background reflectance of the current pixel or the current the period exceeds 61 days;

(7) If the background reflectance of the current pixel is still not obtained, skip to step (1), change the upper and lower limits of the angle, and cycle from steps (1) to (7) until the background reflectance of the current pixel is obtained, or the upper limit Exceeds 90º, and the lower limit exceeds -90º;

(8) If the background reflectance of channel 1 of the current day is found, the background reflectance of channel 2 of the corresponding sunny day is the background reflectance of channel 2 of the current day;

(9) Remove the background aerosol effect of the obtained image. Assuming that the aerosol optical depth value at 550nm on the background day is 0.05, and the aerosol model is continental, the background aerosol effect is corrected to obtain the surface reflectance of channel 1 and channel 2.

Different from other algorithms, this algorithm takes into account the influence of the length of the sliding time window and the satellite zenith angle when determining the surface reflectance.

2. Aerosol optical depth retrieval method

The cloud pixels, water body pixels, and ice and snow pixels of the 630nm and 850nm channels are firstly eliminated, and then the surface reflectance of these two channels is estimated by using the time series of clear-sky zenith reflectance images. Simultaneously, a look-up table is built offline based on aerosol patterns, different combinations of incident and observed attitudes, and surface reflectivity. Finally, based on the time series zenith albedo data and the determined surface albedo, according to the block adjustment value, the spatio-temporal distribution of aerosol optical depth is inverted. Specific steps are as follows:

1) Eliminate cloud pixels, water body pixels and ice and snow pixels. According to the CLAVR cloud results that come with AVHRR data, cloud pixels, some cloud pixels, and some clear pixels are filtered out, and only clear pixels are kept; then the ocean water body is filtered out for sunny pixels, and land and inland water body images are kept Yuan; especially for the spring and winter in East Asia, the reserved pixels need to further filter out the ice and snow pixels.

Among them, the ice and snow mask comes from a product developed by NOAA (ftp://www.orbit.nesdis.noaa.gov/pub/smcd/emb/snow/global_mult_snow_ice/).

2) Determine the surface reflectance according to the above method for determining the surface reflectance.

3) Build a lookup table offline based on aerosol patterns, different combinations of incident and observed attitudes, and surface reflectivity.

Based on 6S (Second Simulation of a Satellite Signal in the Solar Spectrum-Vector) software, the required lookup table is established according to the specific settings of each parameter in Table 1. Although the calculation of the lookup table is time-consuming, it can be compensated by the final fast inversion of the aerosol optical depth.

Table 1. The settings of each parameter of the lookup table in aerosol retrieval

variable name Enter the number enter wavelength 2 0.63 μm, 0.85 μm Sun zenith angle (degrees) 8 0, 12, 24, 36, 48, 54, 60, 66 Satellite zenith angle (degrees) 11 0, 8, 14, 20, 24, 30, 36, 42, 48, 54, 60 Relative azimuth (degrees) 15 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120, 132, 144, 160, 180 Aerosol Optical Depth 5 0.001,0.2,0.5,1.0,2.0 Surface reflectivity 4 0.01, 0.12, 0.23, 0.34 Atmospheric Mode 4 Mid-latitude summer, Mid-latitude winter, Subarctic summer, Tropical

4) Aerosol optical depth inversion: the look-up table calculates the zenith reflectance under the given input aerosol optical depth, surface reflectance and geometric observation conditions, but the geometric observation conditions of the pixels in the actual inversion and the obtained Surface reflectance is often not a given input set, so a lookup table needs to be interpolated to the actual inversion input. This method uses linear interpolation. Since the unknown quantity to be obtained is the aerosol optical thickness, when interpolating, it is necessary to first find the geometric attitude angle of the actual inversion input and the zenith reflectance corresponding to the given input near the surface reflectance, and then according to the surface reflectance, relative azimuth angle, satellite zenith angle, and solar zenith angle, interpolate each level in turn, and return the interpolation result of this level to the previous level to continue interpolation until the corresponding aerosol optics when the actual geometric observation and surface reflectance are obtained thickness. The zenith reflectance values are known for both channels, and unlike other methods, the aerosol optical depth is determined from adjusted values rather than least squares.

3. Comparison of surface reflectance results

The surface reflectance results were compared with three other surface reflectance results under the same incidence and observation attitude conditions, and the four results used the same AVHRR dataset preprocessed and corrected for gas absorption. 1) Use AVHRR zenith reflectivity, ground-based AERONET aerosol optical depth products within AVHRR transit time ±30 minutes, ground-based AERONET spectral distribution and birefringence index products within AVHRR transit time ±3 hours to invert the surface reflectance of AERONET sites ; 2) MODIS BRDF/Albedo product (MCD43C2) reprocesses the incident and observation attitudes to AVHRR, and because it has time series results, it can provide time-consistent result references. 3) Apply the European Riffler et al. (2010) algorithm to the same AVHRR dataset, which only considers the effect of the satellite zenith angle.

The comparison of surface reflectance measured by Beijing Station as shown in Figure 1 and Xinglong Station as shown in Figure 2, the time series shown is from January 31, 2009 to May 31, 2010, where red is the inversion method of the present invention Results; blue is the European algorithm result; green is the surface reflectance after MCD43C2 processing; purple is the surface reflectance after AERONET site processing. Through the comparison of the two sites, it is found that the surface albedo inversion result of the present invention is closer to the surface albedo processed by AERONET and MCD43, showing the advantages of this algorithm compared with the European algorithm.

4. Aerosol optical depth inversion results and comparative verification

As shown in Figure 3, this algorithm is applied to East Asia (18°N-54°N, 100°E-136°E), and the inversion time period is from January 31, 2009 to May 31, 2010. And compared with the ground-based observations of AERONET station.

The aerosol optical depth inversion results of the method of the present invention at 550nm in the figure, the aerosol optical depth inversion results of the European algorithm are compared with the MODIS aerosol optical depth products. The thin solid line is the 1:1 line, and the thick solid line, dotted line, and dashed line are the fitting straight lines of MODIS, European algorithm, and this algorithm, respectively. The text in the figure shows the fitting formula, correlation coefficient ( R ) and RMSE error.

As shown in Figure 4, the scatter density map of the comparison between the results of this algorithm and the MYD04 aerosol optical depth product in East Asia on February 27, 2009. On February 27, 2009, the comparative scatter density diagram of the results of the method of the present invention and the product of MYD04 aerosol optical depth in East Asia. Data were sorted and binned into 0.03 x 0.03 bins. The thin and thick solid lines are the 1:1 line and the linear fit line, respectively. The text in the figure expresses the linear fitting formula, correlation coefficient ( R ), RMSE error and the number of points involved in the fitting (N).

As shown in Figures 5a-5d, the mapping and comparison of aerosol optical depth in East Asia on March 26, 2009: Figure 5a is the AVHRR aerosol optical depth map retrieved by the method of the present invention; Figure 5b is the MODIS dark target aerosol optical depth Fig. 5c is a true-color composite map of 1km MYD02 data; Fig. 5d is the comparison of AVHRR aerosol optical depth values, MODIS dark target aerosol optical depth values and AERONET aerosol optical depth different border crossing times in the same region).

One of the innovative points of the method proposed by the present invention is that the influence of the observed zenith angle and the length of the sliding time window is considered in the multi-day joint observation method, especially for the East Asian region with serious air pollution, the accuracy of the obtained surface reflectance results is higher it is good. The second innovation point is that the block adjustment value is used instead of the least squares method as the criterion for the aerosol optical depth inversion process. Compared with AERONET ground-based observations, it is found that the aerosol optical depth results obtained by inversion are more accurate than the algorithm that does not consider the above factors.

The present invention is a method for retrieving the time-space distribution of aerosol optical thickness for sensors without a mid-infrared channel. It has the technical characteristics of following natural laws. This method solves the technical problem of retrieving aerosol optical thickness and does not belong to the patent law 25 The rules and methods of intellectual activities mentioned in Article 2.

Claims (3)

1. the aerosol optical depth spatial and temporal distributions inversion method of no mid-infrared channel sensor is it is characterised in that using in no When infrared channel sensor determines Reflectivity for Growing Season, sliding time window length is 30-61 days, and current pixel is currently processed picture Unit, present satellites zenith angle is the satellite zenith angle of current pixel, comprises the following steps that：

(1) present satellites zenith angle is deducted 10 ° of current lower limits as angle, plus 10 ° of current upper bound as angle；

(2) the current bound of angle is compared with ± 90 ° respectively, if it was exceeded, the current upper bound exceeding or lower limit are set For 90 ° or -90 °；

(3) current sliding time window length is initially set the shortest 30 days；

(4) select current sliding time window in satellite zenith angle between the passage 1 of current angular bound zenith reflectivity, And the zenith reflectivity number selected must be more than 7；

(5) when meeting condition：

1) secondary dark reflectivity is not equal to the zenith reflectivity when the day before yesterday；Or

2) secondary dark reflectivity is equal to the zenith reflectivity when the day before yesterday, and current slot is equal to 61 days, will be anti-for the zenith selected Penetrate the background reflectivity as current pixel for time dark reflectivity of rate；

(6) without the background reflectivity obtaining current pixel, current slot increases by 2 days, constant when the day before yesterday；Circulation step (4), (5) and (6) are until obtaining the background reflectivity of current pixel or current slot more than 61 days；

(7) if not obtaining the background reflectivity of current pixel yet, skipping to step (1), changing the bound of angle, from step (1) to step (7) circulation until obtaining the background reflectivity of current pixel, or the upper limit is more than 90 °, and lower limit exceedes -90 °；

(8) if finding the background reflectivity of the passage 1 when the day before yesterday, the background reflectivity of the passage 2 of corresponding sunny day is Background reflectivity when the passage 2 of the day before yesterday；

(9) obtained by removing, the background aerosol effect of image is it is assumed that aerosol optical depth value at 550nm for the background day is 0.05, aerosol model is continent type, background aerosol effect is corrected, obtains passage 1 and the earth surface reflection of passage 2 Rate.

2. aerosol optical depth spatial and temporal distributions inversion method as claimed in claim 1 is it is characterised in that aerosol optical depth Inversion method comprises the following steps that：

1) reject described passage 1 and passage 2 medium cloud pixel, water body pixel and ice and snow pixel, only retain sunny pixel, then to fine Bright pixel filters ocean water body, retains land and Inland Water pixel；For winter in spring, retain pixel and further filter out ice and snow picture Unit；

2) look-up table is set up offline based on aerosol model, different incidence and observation attitude integration and Reflectivity for Growing Season；

3) Reflectivity for Growing Season based on time series zenith reflectivity data and determination, according to block adjustment value, inverting gas is molten Glue optical thickness spatial and temporal distributions result；Look-up table is carried out with linear interpolation input to practical inversion, during interpolation, first look for reality Given input corresponding zenith reflectivity near the geometry attitude angle of inverting input, Reflectivity for Growing Season, then according to earth's surface Reflectivity, relative bearing, satellite zenith angle, the order of solar zenith angle, successively to often grading row interpolation, and inserting this grade Value result returns upper level and continues interpolation, until it is thick to obtain corresponding aerosol optical when actual geometry observation and Reflectivity for Growing Season Degree.

3. aerosol optical depth spatial and temporal distributions inversion method as claimed in claim 1 or 2 is it is characterised in that red in described nothing Outer tunnel sensor is AVHRR sensor, ccd sensor, VISSR sensor.