CN109213964B - Satellite AOD product correction method fusing multi-source characteristic geographic parameters - Google Patents

Satellite AOD product correction method fusing multi-source characteristic geographic parameters Download PDF

Info

Publication number
CN109213964B
CN109213964B CN201810776814.2A CN201810776814A CN109213964B CN 109213964 B CN109213964 B CN 109213964B CN 201810776814 A CN201810776814 A CN 201810776814A CN 109213964 B CN109213964 B CN 109213964B
Authority
CN
China
Prior art keywords
aod
aeronet
data
satellite
variables
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN201810776814.2A
Other languages
Chinese (zh)
Other versions
CN109213964A (en
Inventor
邹滨
游介文
江佳
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Central South University
Original Assignee
Central South University
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Central South University filed Critical Central South University
Priority to CN201810776814.2A priority Critical patent/CN109213964B/en
Publication of CN109213964A publication Critical patent/CN109213964A/en
Application granted granted Critical
Publication of CN109213964B publication Critical patent/CN109213964B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
    • G06F17/10—Complex mathematical operations
    • G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Data Mining & Analysis (AREA)
  • General Physics & Mathematics (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • Mathematical Physics (AREA)
  • Computational Mathematics (AREA)
  • Mathematical Analysis (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Operations Research (AREA)
  • Probability & Statistics with Applications (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Algebra (AREA)
  • Evolutionary Biology (AREA)
  • Databases & Information Systems (AREA)
  • Software Systems (AREA)
  • General Engineering & Computer Science (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Image Processing (AREA)

Abstract

The invention discloses a satellite AOD product correction method fusing multi-source characteristic geographic parameters, which comprises the steps of firstly, constructing a satellite AOD product correction regression model based on a random forest algorithm, obtaining an AERONET AOD planar simulation space distribution result as satellite AOD data after primary correction; then comparing the satellite AOD with ground AERONET AOD real data one by one to obtain the deviation of the satellite AOD and the ground AERONET AOD real data, carrying out spatial interpolation estimation, and correcting the satellite AOD value subjected to preliminary correction based on the estimated deviation value; and finally, setting a correction threshold value, and comparing the corrected satellite AOD data with the ground AERONET AOD data site by site. The invention relates to a satellite AOD product correction method which is high in comprehensiveness, strong in referenceability, effective and accurate under the condition that AOD data of current satellite remote sensing inversion and a real AOD value have a certain difference, so that atmospheric PM is accurately researched by utilizing satellite AOD2.5And the effective concentration space-time distribution characteristics and the change rule and the evaluation of the atmospheric pollution condition of the continuous large-range geographic area provide references.

Description

Satellite AOD product correction method fusing multi-source characteristic geographic parameters
Technical Field
The invention relates to the field of satellite remote sensing inversion products, in particular to an aerosol optical thickness (AOD) correction method.
Background
The atmospheric aerosol is a multiphase system composed of various solid and liquid particles with the diameter of 0.001-100 mu m suspended in the atmosphere. The main sources of atmospheric aerosols include natural sources such as sand dust, salt particles, emissions of volcanic eruptions and smoke from forest combustion, and man-made sources such as smoke from the combustion of fossil and non-fossil fuels, transportation and various industrial emissions. In recent decades, the rapid development of global urbanization and industrialization has increased the man-made emission, and serious atmospheric pollution has seriously affected air quality and regional climate change, which becomes a global environmental problem. Most studied and focused on the Optical properties of aerosols at home and abroad are Aerosol Optical Depths (AODs). AOD is an important parameter for describing the reduction effect of aerosol on light, and PM is constructed by utilizing satellite AOD products2.5AOD relational model to estimate surface PM2.5The simulation of the concentration space-time distribution aims at solving the problem that the PM is difficult to be accurately revealed by the traditional ground observation means2.5The spatial and temporal distribution of pollution over a wide area and evolution characteristics provide a new solution.
The ground-based Aerosol remote sensing monitoring can obtain AOD accurate information on a point scale, and the ground-based AOD product widely accepted and used internationally at present is an Aerosol automatic observation Network (Aerosol Robotic Network, AeroNET) initiated and supported by the American NASA. The AERONET observation network covers 1497 sites worldwide, and a ground-based CIMEL solar photometer is utilized to obtain AODs with regional representatives in a global range. The observation error of the AERONET AOD is 0.01-0.02, and the accuracy of an AOD data product obtained by satellite remote sensing can be verified and evaluated by taking the observation value of the AERONET as a true value. However, compared with ground monitoring, the satellite remote sensing observation range is wider, the space-time limitation is smaller, and the AOD space-time distribution characteristics of a large-range geographic area can be obtained. However, aerosol products inverted by satellite remote sensing are affected by various factors such as ground surface coverage, atmospheric radiation, meteorological factors, instrument errors, inversion algorithms and products, and satellite AOD data and real AOD values have certain differences, so that result deviation of researching atmospheric pollution concentration space-time distribution characteristics by using direct satellite AOD products is caused, and a method for effectively performing space-time correction on satellite AOD data products is urgently needed.
The change of regional geographic characteristic elements has strong influence on the height of the AOD, land utilization, traffic roads, population and the like are closely related to artificial aerosol emission intensity, and meteorological factors such as wind speed, precipitation, relative humidity and the like can influence the diffusion and sedimentation of atmospheric aerosol. Therefore, the space-time distribution characteristic of the AOD can be simulated by analyzing the space distribution rule of the multi-source geographic characteristic parameters; meanwhile, the difference of regional geographical element distribution can further explain the difference between the AOD product of satellite remote sensing inversion and the AOD product of ground observation, and feasibility is provided for correcting satellite AOD data products. Therefore, the satellite AOD product is corrected by utilizing the deviation of the satellite AOD data and the ground AERONET real AOD data and fusing multi-source characteristic geographic parameters on the basis of a machine learning algorithm and a ground statistical theory, and the optimal, unbiased, accurate and reliable satellite AOD correction product can be obtained.
Disclosure of Invention
In order to effectively solve the problems of low precision and large error of an AOD product obtained by satellite remote sensing inversion, the invention provides a satellite AOD product correction method fusing multi-source characteristic geographic parameters.
In order to achieve the technical purpose, the technical scheme of the invention is that,
a satellite AOD product correction method fusing multi-source characteristic geographic parameters comprises the following steps:
step 1, establishing a sample training set by utilizing a modeling factor data set of a corresponding position of a matched ground AERONET monitoring station, then, taking the ground AERONET AOD data of the training set as dependent variables, taking the multi-source characteristic geographic parameters and the satellite AOD data as independent variables, randomly selecting corresponding dependent variables and independent variables from each sample as branch nodes of the decision tree to construct a decision tree corresponding to each sample, so that the decision trees form a random forest regression model, and the importance of independent variables in the model as characteristic variables is judged, then optimizing the quantity of the characteristic variables and the decision trees to complete the establishment of a random forest regression model, inputting satellite AOD data to be corrected and known multisource characteristic geographic parameters of corresponding places into the random forest regression model, obtaining an AERONET AOD surface-shaped simulation space distribution result as satellite AOD data after primary correction;
step 2, comparing the preliminarily corrected satellite AOD data obtained in the step 1 with ground station-by-station AERONET AOD real data to obtain the deviation of the two, performing common Krigin interpolation on the deviation to obtain the planar spatial distribution of the deviation between the ground AERONET AOD and the preliminarily corrected satellite AOD, and obtaining the secondarily corrected satellite AOD data according to the preliminarily corrected satellite AOD data and the deviation;
and 3, setting a correction threshold, comparing the satellite AOD data subjected to secondary correction with ground station-by-station AERONET AOD real data, then comparing the satellite AOD data with the correction threshold, returning to the step 2 if the satellite AOD data is larger than the correction threshold, and outputting a finally obtained satellite AOD correction product and a precision evaluation result until the currently obtained result is completely the same as the previous result or is smaller than the correction threshold.
In the method, in the first step, the modeling factor data set of the corresponding position of the matched ground AERONET monitoring station comprises ground AERONET AOD station data and multi-source characteristic geographic parameters, the multi-source characteristic geographic parameters comprise land utilization/coverage, road traffic, terrain elevation, wind speed, precipitation, relative humidity, temperature and air pressure, and the ground AERONET AOD station data serving as a dependent variable and the multi-source characteristic geographic parameters serving as an independent variable are matched one by one according to the monitoring stations to form the modeling factor data set.
In the method, in the first step, the establishing of the training set by the samples includes the following steps:
sampling partial observed values of an original sample D containing p characteristic variables by adopting a Bootstrap method, and randomly generating K training sets theta1、θ2、…θk。
In the method, in the first step, the forming of the random forest regression model comprises the following steps:
for samples in each training set, randomly selecting a fixed number n of variables from p variables as branch nodes of the decision tree to construct the decision tree, wherein n is<p to generate a corresponding decision tree H (X, θ) based on each training set1)}、{H(X,θ2)}、…{H(X,θk) And (v), wherein X is an independent variable, the K decision trees form a random forest regression model, a Gini coefficient is used for selecting a characteristic variable of the random forest regression model, and a Gini (v) calculation formula of the model at a node v is as follows:
Figure BDA0001730571490000041
in the formula (I), the compound is shown in the specification,
Figure BDA0001730571490000042
is the observed value, X, of the jth variable at node viKini at split node vInformation Gain (X)iV) is the difference in purity between the children of nodes v and v.
The method, the Gain (X) of the Kini informationiV) the calculation formula is as follows:
Gain(Xi,v)＝Gini(Xi,v)-wLGini(Xi,vL)-wRGini(Xi,vR)
in the formula, vLAnd vRAre the left and right child nodes, w, respectively, at node vLAnd wRRespectively, the proportion of characteristic variables belonging to left and right child nodes, at each node, randomly extracting mtry variables from p variables, and finally obtaining the characteristic variable with the maximum information gain for splitting the node v, wherein
Figure BDA0001730571490000043
In the method, in the first step, the importance determination of the feature variables in the model includes the following steps:
evaluating predictive feature variable X using predictive variable importance scoresiSelecting the characteristic variables with large contribution degree from the contribution to the random forest regression model, wherein the calculation formula is as follows:
Figure BDA0001730571490000044
in the formula (I), the compound is shown in the specification,
Figure BDA0001730571490000045
is in the random forest of ntree trees, is XiA split set of nodes.
In the method, in the first step, optimizing the number of feature variables and decision trees includes the following steps:
and traversing the mtry parameters from 1 to n, performing n times of tests, recording the error rate of each test, selecting the mtry value with the lowest error rate as the optimal mtry value, fixing the obtained optimal mtry value to the number ntree of the decision trees, traversing the ntree from 1 to n, performing n times of tests, recording the error rate of each test, and selecting the ntree value with the lowest error rate as the optimal ntree value.
In the second step, the deviation is interpolated by ordinary kriging, and the formula is as follows:
Figure BDA0001730571490000051
in the formula (I), the compound is shown in the specification,
Figure BDA0001730571490000052
is a point (x)0,y0) Estimated value, i.e. Z0＝Z(x0,y0)；λiAre weight coefficients and satisfy point (x) at the same time0,y0) Estimate of (c)
Figure BDA0001730571490000053
With the true value Z0A set of optimal coefficients with the smallest difference;
firstly, unbiased:
Figure BDA0001730571490000054
optimality:
Figure BDA0001730571490000055
in the method, in the third step, setting the correction threshold includes the following steps:
the correction threshold is expressed in terms of the Expected Error (EE), and is formulated as follows:
EE＝a+b×σAERONET
wherein σAERONETRepresenting the ground AERONET AOD, coefficients a and b being parameters for adjusting the precision of the corrected product, 0<a<1,0<b<1。
In the method, in the third step, the precision evaluation is realized by the following steps:
the accuracy evaluation index is expressed by a correlation coefficient R and a relative error delta, and the formula is as follows:
Figure BDA0001730571490000056
Figure BDA0001730571490000061
in the formula, x and y respectively represent an AERONET AOD value and a corrected satellite AOD value,
Figure BDA0001730571490000062
mean values of annual AERONET AOD and corrected satellite AOD are shown, respectively.
The invention has the technical effects that: (1) the method is a novel satellite AOD product correction method for coupling a machine learning algorithm and geostatistics theories and technologies under the conditions of large satellite AOD error and low precision; (2) meanwhile, the method is a relatively accurate satellite AOD data correction method, multi-source characteristic geographic parameters are fused, iterative correction is carried out on a satellite AOD data product based on the deviation between ground AERONET AOD real data and satellite AOD data until an optimal unbiased estimated satellite AOD product is obtained, and the repair result is accurate and reliable; (3) the method is also a relatively high-efficiency satellite AOD product correction method, the required technology and data can be obtained in real time, the algorithm is high-efficiency, and the satellite AOD product can be quickly corrected.
Drawings
FIG. 1 is a flowchart of a satellite AOD product correction method fusing multi-source characteristic geographic parameters, wherein (a) shows construction and simulation of an AOD random forest regression model fusing multi-source characteristic geographic parameters, (b) shows geostatistical modeling of a deviation value of a satellite AOD and a ground AERONET AOD, and (c) shows MODIS AOD correction products of iterative modeling solution optimal unbiased estimation.
FIG. 2 is a comparison graph of precision of a satellite AOD correction product fused with multi-source characteristic geographic parameters and an uncorrected satellite AOD product, wherein (a) is a precision graph of an original AOD product, and (b) is a precision graph of the AOD product corrected by the method.
Detailed Description
The following is a detailed description of a preferred embodiment of the invention, taken in conjunction with the accompanying drawings.
1. AOD random forest regression model construction and simulation fusing multi-source characteristic geographic parameters
The invention discloses a method for constructing and simulating an AOD random forest regression model fusing multi-source characteristic geographic parameters, which comprises the following steps:
step 1): and collecting and extracting multi-source characteristic geographic elements. In this example, the kyojin ji area is taken as an example, and geographical parameter sets related to regional atmospheric aerosol emission and diffusion are collected, mainly including weather daily value data (air temperature, air pressure, humidity, precipitation, wind speed), land cover/utilization type (construction land, grassland, woodland, water area), road traffic, population grids and terrain elevation data. For planar geographic elements, such as grid data of population grids, terrain elevations and the like, directly extracting grid attribute values corresponding to a ground AERONET monitoring station as corresponding characteristic variables; for point-like geographic elements, such as point-like meteorological elements, spatial distribution grid data of the elements are obtained by using a spatial interpolation method and then extracted; as for the road and land elements, because the space scale effect exists, the total length of the road in a certain range of the ground monitoring station and the area occupation ratio of various land are used as characteristic variables, and the area occupation ratio of various land and the total length of all roads in the range of 500-5000m around the station are respectively extracted and calculated by establishing a buffer zone.
Step 2): and (5) satellite AOD product data processing. In the example, taking a DT DB Combined AOD data product inverted by an MODIS sensor mounted on a tera satellite as an example, the product is downloaded in an NASA official website, and is subjected to projection, splicing and resampling, so that raster data with a spatial resolution of 10 × 10km is generated, and AOD data corresponding to the kyazine wing area is cut out. And finally, extracting satellite AOD data of corresponding points according to the spatial position of the ground AERONET monitoring station.
Step 3): and (3) processing data of the ground AERONET AOD product.
Wave band matching
The AeroNET can provide AOD data at wave bands of 1020, 870, 670, 500, 440nm and the like, while the satellite remote sensing inversion product MODIS only provides AOD data at 3 wave bands of 660, 550, 470 nm. The two AOD data products do not have corresponding identical wavebands. Therefore, the AERONET AOD data needs to be subjected to wave band interpolation to obtain AOD data matched with the MODIS product wave band.
On a wave band without the influence of water vapor, the spectral distribution of aerosol particles meets Junge distribution, and AOD exists between the AOD and the wavelength
Figure BDA0001730571490000085
Relation formula
σ(λ)＝βλ-α
Wherein σ (λ) represents an AOD value at a wavelength of λ; beta represents
Figure BDA0001730571490000086
The value of the atmospheric turbidity coefficient is influenced by the total number of aerosol particles, the particle spectrum distribution and the refractive index; alpha is
Figure BDA0001730571490000087
Wavelength index, the value of which is related to the aerosol particle radius. Assuming a band λ1，λ2The part is not influenced by water vapor and is provided with the following formula:
σ(λ1)＝βλ1 -α
σ(λ2)＝βλ2 -α
the formula above can be combined to obtain:
Figure BDA0001730571490000081
in the formula (I), the compound is shown in the specification,
Figure BDA0001730571490000082
is a wavelength lambda1-λ2And (3) wavelength index of (d). At a known lambda1，λ2And wavelength index
Figure BDA0001730571490000083
In the case of (2), wavelength λ1-λ2The value of the optical thickness of the aerosol at any wavelength λ in between can be interpolated by the following formula:
Figure BDA0001730571490000084
based on the above formula, the present example uses AERONETAOD data at 440 and 870nm and a wavelength index α between 440 and 870nm440-870And obtaining the AERONET AOD at 550nm by interpolation.
② space-time matching
AERONET AOD data is continuous observation data based on fixed time (typically 15min) intervals of monitored site locations, while MODIS Combined AOD values represent instantaneous monitoring values of 10km × 10km area, while satellite transit time is 30 minutes at 10 am of local time. Since the MODIS AOD is not the same as the AERONET AOD on both the temporal and spatial scales. If the AERONET observation value of the satellite transit time is directly taken to be compared with the MODIS AOD value of a single pixel, the comparison reliability is low, and the singular value phenomenon is easy to occur. In order to improve the comparability of two products, the AERONET site is taken as a center in the example, an abnormal value is removed by adopting MODISAOD of 5 × 5 pixels around the site, the average value is taken as a satellite AOD average value at the AERONET site, the time average value of ground AERONET observation data of 1.5h before and after the satellite passes through is taken as the AERONET observation average value, a certain screening condition is added, at least 5 or more effective pixels in the 5 × 5 pixels around the AERONET site need to be met, and at least 2 or more effective observation values of the AERONET surface observation data can be matched and compared within 1.5h before and after the satellite passes through.
Step 4): construction of random forest regression model
Constructing a Random Forest (Random Forest) regression model, taking ground AERONET AOD site data as a dependent variable, and fusing multi-source characteristic geographic parameters (meteorological elements and land utilization)Coverage, road traffic, terrain elevation, etc.) and MODIS DT DB Combined AOD products are independent variables. The main steps of random forest regression are as follows: sample data is randomly sampled: sampling partial observed values of an original sample D containing p characteristic variables by adopting a Bootstrap method (with replaced sampling), and randomly generating K training sets theta1、θ2、…θkAnd the data which is not extracted forms an out-of-bag data set (OBB) which can be used as a test sample set. Randomly selecting characteristic variables: for each sample in the training set, randomly selecting a fixed number n (n) from p variables<p) as branch nodes of the decision tree, thereby generating a corresponding decision tree { H (X, theta) for each training set1)}、{H(X,θ2)}、…{H(X,θk) And (6) forming a random forest regression model by the K decision trees, and finally obtaining a predicted value by the model in an average value mode.
In the example, a random forest regression model characteristic variable is selected by using a Gini coefficient, and Gini (v) calculation formula of the model at a node v is as follows:
Figure BDA0001730571490000091
in the formula (I), the compound is shown in the specification,
Figure BDA0001730571490000101
is the observed value of the jth variable at node v. XiGain (X) of kini information at split node viV) is the difference between the purities of the children of nodes v and v. The calculation formula is as follows:
Gain(Xi,v)＝Gini（Xi,v)-wLGini(Xi,vL)-wRGini(Xi,vR)
in the formula, vLAnd vRAre the left and right child nodes, w, respectively, at node vLAnd wRRespectively, the ratio of the characteristic variable belonging to the left and right child nodes. At each node, randomly decimating mtry among p variables
Figure BDA0001730571490000102
And (4) using the variables, namely the characteristic variable which finally obtains the maximum information gain, for splitting the node v.
This example uses predictor variable importance scores to evaluate predictor feature variable XiSelecting the characteristic variables with large contribution degree from the contribution to the random forest regression model, wherein the calculation formula is as follows:
Figure BDA0001730571490000103
in the formula (I), the compound is shown in the specification,
Figure BDA0001730571490000104
is in the random forest of ntree trees, is XiA split set of nodes.
In the process of constructing the random forest model, mtry and ntree are two important parameters which respectively represent the number of variables sampled randomly and the number of decision trees when branches are constructed. The mtry is n variables randomly selected from the P variables when the decision tree is constructed, namely the number of branch nodes of the decision tree, and the prediction error rate of the random forest model can be reduced by selecting a proper mtry value. The invention sets the mtry parameters from 1 to n through traversal, performs n times of tests, records the error rate of each test, and selects the mtry value with the lowest error rate as the optimal mtry value. The ntree represents the number of the decision trees, and generally speaking, when the ntree value is too small, the error rate of the model is high, and when the ntree value is too large, the complexity of the model is improved, and the efficiency of the model is reduced. In the example, the optimal mtry value obtained in the above way is fixed and is not changed, the ntree is traversed from 1 to n, n times of tests are carried out, the error rate of each test is recorded, and the ntree value with the lowest error rate is selected as the optimal ntree value.
Step 5): simulating AERONET AOD spatial distribution
And after the random forest regression model is completed, forming an input variable by the multi-source characteristic geographic parameters and the MODIS AOD data, inputting the input variable into the constructed random forest regression model, and obtaining an AERONET AOD planar simulation spatial distribution result as initially corrected MODIS AOD data.
2. Geostatistical modeling of deviation values of MODIS AOD and ground AERONET AOD
The geostatistical modeling of the deviation value of the MODIS AOD and the ground AERONET AOD adopted by the invention comprises the following steps:
step 1): solving deviation value of MODIS AOD and ground AERONET AOD
And (3) comparing the preliminarily corrected MODIS AOD data obtained in the step (1) with ground station-by-station AERONET AOD real data to obtain the deviation delta (t) of the two data.
Step 2): geostatistical modeling of deviations using kriging interpolation
And (3) performing Krigin interpolation on the deviation delta (t) between the ground station AERONET AOD and the MODIS AOD obtained in the step to obtain the planar spatial distribution of the deviation delta (t) between the ground AERONET AOD and the primarily corrected MODIS AOD. The general Kriging interpolation (OK) is a commonly used Kriging interpolation method, and the basic principle is shown in the following formula:
Figure BDA0001730571490000111
in the formula (I), the compound is shown in the specification,
Figure BDA0001730571490000112
is a point (x)0,y0) Estimated value, i.e. Z0＝Z(x0,y0)；λiAre weight coefficients and satisfy point (x) at the same time0,y0) Estimate of (c)
Figure BDA0001730571490000113
With the true value Z0The set of optimal coefficients with the smallest difference.
Firstly, unbiased:
Figure BDA0001730571490000114
optimality:
Figure BDA0001730571490000115
step 3): modified MODIS AOD product
And (3) obtaining a secondarily corrected MODIS AOD data product by using the primarily corrected MODIS AOD data obtained in the step (1) and adding the deviation delta (t) between the ground AERONET AOD simulated in the step (2) and the primarily corrected MODIS AOD.
3. MODIS AOD correction product for solving optimal unbiased estimation through iterative modeling
The invention adopts an MODIS AOD correction product for solving the optimal unbiased estimation by iterative modeling, which comprises the following steps:
step 1): and setting a correction threshold according to the actual product precision requirement. The correction threshold is expressed in terms of the Expected Error (EE), and is formulated as follows:
EE＝a+b×σAERONET
wherein σAERONETRepresenting a terrestrial AERONET AOD. The coefficients a and b are predefined constants (0)<a<1,0<b<1)。
In this example, by setting the coefficients a to 0.05 and b to 0.2, the correction threshold EE is 0.05+0.2
σAERONET。
Step 2): and (3) comparing the twice-corrected MODIS AOD data obtained in the step (2) with ground AERONETAOD data site by site, judging whether the difference value delta (t +1) of the two data is larger than a set threshold value, if so, repeating the step (2) until the current obtained result is completely the same as the previous result or smaller than the correction threshold value, and outputting the finally obtained MODIS AOD correction product and the precision evaluation result. Wherein, the accuracy evaluation index is expressed by a correlation coefficient R and a relative error (delta), and the formula is as follows:
Figure BDA0001730571490000121
Figure BDA0001730571490000122
in the formula, x and y respectively represent AERONET AOD value and corrected MODIS AOD value,
Figure BDA0001730571490000123
Figure BDA0001730571490000124
mean values of annual AERONET AOD and corrected MODIS AOD are shown, respectively.
In the example, when the cycle is stopped, t is 3, the correlation coefficient R of the finally obtained satellite AOD correction product and the ground AERONET AOD data is 0.98, the relative error delta is 10.12%, the correlation coefficient R is increased by 0.22 compared with the original AOD data, and the relative error delta is reduced by 34.01%, which shows that the precision of the satellite AOD data product can be effectively improved.
While the invention has been described in further detail with reference to specific preferred embodiments thereof, it will be understood by those skilled in the art that various changes in form and details may be made therein without departing from the spirit and scope of the invention as defined by the appended claims.

Claims (10)

1. A satellite AOD product correction method fusing multi-source characteristic geographic parameters is characterized by comprising the following steps:
step 1, establishing a sample training set by utilizing a modeling factor data set of a corresponding position of a matched ground AERONET monitoring station, then, taking the ground AERONET AOD data of the training set as dependent variables, taking the multi-source characteristic geographic parameters and the satellite AOD data as independent variables, randomly selecting corresponding dependent variables and independent variables from each sample as branch nodes of the decision tree to construct a decision tree corresponding to each sample, so that the decision trees form a random forest regression model, and the importance of independent variables in the model as characteristic variables is judged, then optimizing the quantity of the characteristic variables and the decision trees to complete the establishment of a random forest regression model, inputting satellite AOD data to be corrected and known multisource characteristic geographic parameters of corresponding places into the random forest regression model, obtaining an AERONET AOD surface-shaped simulation space distribution result as satellite AOD data after primary correction;
step 2, comparing the preliminarily corrected satellite AOD data obtained in the step 1 with ground station-by-station AERONET AOD real data to obtain the deviation of the two, performing common Krigin interpolation on the deviation to obtain the planar spatial distribution of the deviation between the ground AERONET AOD real data and the preliminarily corrected satellite AOD data, and obtaining the secondarily corrected satellite AOD data according to the preliminarily corrected satellite AOD data and the deviation;
and 3, setting a correction threshold, comparing the satellite AOD data subjected to secondary correction with ground station-by-station AERONET AOD real data, then comparing the satellite AOD data with the correction threshold, returning to the step 2 if the satellite AOD data is larger than the correction threshold, and outputting a finally obtained satellite AOD correction product and a precision evaluation result until the currently obtained result is completely the same as the previous result or is smaller than the correction threshold.
2. The method according to claim 1, wherein in step 1, the modeling factor data set of the corresponding position of the matched ground AERONET monitoring station comprises ground AERONET AOD station data and multi-source characteristic geographic parameters, wherein the multi-source characteristic geographic parameters comprise land utilization/coverage, road traffic, terrain elevation, wind speed, precipitation, relative humidity, temperature and air pressure, and the ground AERONET AOD station data serving as a dependent variable and the multi-source characteristic geographic parameters serving as independent variables are matched one by one according to the monitoring stations to form a modeling factor data set.
3. The method of claim 1, wherein the step 1 of establishing the training set by the samples comprises the steps of:
for an original sample D containing p characteristic variables, a Bootstrap method is adopted to carry out partial observation value of DSampling and randomly generating K training sets theta1、θ2、…θK。
4. The method as claimed in claim 3, wherein in step 1, the step of composing the random forest regression model comprises the steps of:
for samples in each training set, randomly selecting a fixed number n of variables from p characteristic variables as branch nodes of the decision tree to construct the decision tree, wherein n is<p to generate a corresponding decision tree H (X, θ) based on each training set1)}、{H(X,θ2)}、…{H(X,θK) And (v), wherein X is an independent variable, the K decision trees form a random forest regression model, a Gini coefficient is used for selecting a characteristic variable of the random forest regression model, and a Gini (v) calculation formula of the model at a node v is as follows:
Figure FDA0003135733880000021
in the formula (I), the compound is shown in the specification,
Figure FDA0003135733880000022
is the observed value, X, of the jth variable at node viGain (X) of kini information at split node viV) is the difference in purity between the children of nodes v and v.
5. The method of claim 4, wherein the kuni information Gain (X)iV) the calculation formula is as follows:
Gain(Xi,v)＝Gini(Xi,v)-wLGini(Xi,vL)-wRGini(Xi,vR)
in the formula, vLAnd vRAre the left and right child nodes, w, respectively, at node vLAnd wRRespectively, the characteristic variables belong to the proportion of the left and right subnodes, mtry variables are randomly extracted from p characteristic variables at each node, and finally the characteristic of the maximum information gain is obtainedThe variable is used for splitting node v, where mtry ═ round
Figure FDA0003135733880000031
6. The method according to claim 4, wherein the step 1 of determining the importance of the independent variable as the characteristic variable in the model comprises the following steps:
evaluating predictive feature variable X using predictive variable importance scoresiAnd (3) contributing to the random forest regression model, and selecting characteristic variables with large contribution degree from the random forest regression model, wherein the calculation formula is as follows:
Figure FDA0003135733880000032
in the formula (I), the compound is shown in the specification,
Figure FDA0003135733880000033
is in random forest of K trees, is XiA split set of nodes.
7. The method of claim 5, wherein the step 1 of optimizing the number of feature variables and decision trees comprises the steps of:
and traversing the number K of the decision trees, fixing the obtained optimal mtry value to be unchanged, traversing the number K from 1 to n, performing n times of tests, recording the error rate of each test, and selecting the K value with the lowest error rate as the optimal K value.
8. The method according to claim 4, wherein in step 2, the deviation is interpolated by ordinary kriging, and the formula is as follows:
Figure FDA0003135733880000034
in the formula (I), the compound is shown in the specification,
Figure FDA0003135733880000035
is a point (x)0,y0) Estimated value, i.e. Z0＝Z(x0,y0)；λiAre weight coefficients and satisfy point (x) at the same time0,y0) Estimate of (c)
Figure FDA0003135733880000036
With the true value Z0A set of optimal coefficients with the smallest difference;
firstly, unbiased:
Figure FDA0003135733880000037
optimality:
Figure FDA0003135733880000041
9. the method of claim 4, wherein the step 3 of setting the correction threshold comprises the steps of:
the correction threshold is expressed in terms of the expected error EE, and is given by the following equation:
EE＝a+b×σAERONET
wherein σAERONETRepresenting the ground AERONET AOD, coefficients a and b being parameters for adjusting the precision of the corrected product, 0<a<1,0<b<1。
10. The method according to claim 4, wherein in the step 3, the precision evaluation is realized by the following steps:
the accuracy evaluation index is expressed by a correlation coefficient R and a relative error delta, and the formula is as follows:
Figure FDA0003135733880000042
Figure FDA0003135733880000043
in the formula, x and y respectively represent an AERONET AOD value and a corrected satellite AOD value,
Figure FDA0003135733880000044
mean values of annual AERONET AOD and corrected satellite AOD are shown, respectively.
CN201810776814.2A 2018-07-13 2018-07-13 Satellite AOD product correction method fusing multi-source characteristic geographic parameters Active CN109213964B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN201810776814.2A CN109213964B (en) 2018-07-13 2018-07-13 Satellite AOD product correction method fusing multi-source characteristic geographic parameters

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN201810776814.2A CN109213964B (en) 2018-07-13 2018-07-13 Satellite AOD product correction method fusing multi-source characteristic geographic parameters

Publications (2)

Publication Number Publication Date
CN109213964A CN109213964A (en) 2019-01-15
CN109213964B true CN109213964B (en) 2021-08-17

Family

ID=64990105

Family Applications (1)

Application Number Title Priority Date Filing Date
CN201810776814.2A Active CN109213964B (en) 2018-07-13 2018-07-13 Satellite AOD product correction method fusing multi-source characteristic geographic parameters

Country Status (1)

Country Link
CN (1) CN109213964B (en)

Families Citing this family (21)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN109785569A (en) * 2019-01-28 2019-05-21 中科光启空间信息技术有限公司 A kind of forest fire monitoring method based on 3S technology
CN110398446A (en) * 2019-08-01 2019-11-01 成都信息工程大学 A kind of PM of combination AOD and meteorologic parameter2.5Concentration approximating method
CN110610258B (en) * 2019-08-20 2022-05-10 中国地质大学（武汉） Urban air quality refined estimation method and device fusing multi-source space-time data
CN111078678B (en) * 2019-12-18 2021-03-23 中国气象局乌鲁木齐沙漠气象研究所 Satellite precipitation data correction method based on multi-source information fusion and scale reduction
CN111104639B (en) * 2019-12-24 2022-06-10 福州大学 Point-surface fused time sequence PM2.5 spatial distribution estimation method
CN111323352B (en) * 2020-04-09 2021-08-24 中南大学 Regional PM2.5 remote sensing inversion model fusing fine particulate matter concentration data
CN111737850B (en) * 2020-05-15 2024-04-19 中国科学院空天信息创新研究院 Multi-source satellite AOD fusion method based on uncertainty on pixel scale
CN112016696B (en) * 2020-08-14 2022-10-04 武汉大学 PM integrating satellite observation and ground observation 1 Concentration inversion method and system
CN112131789B (en) * 2020-09-18 2023-11-21 中国人民解放军国防科技大学 Multispectral precipitation detection system and method based on random forest algorithm
CN112395808A (en) * 2020-11-18 2021-02-23 南京林业大学 Biomass remote sensing mapping method combining random forest and collaborative kriging
CN112926625B (en) * 2020-11-25 2023-12-22 北方工业大学 Deviation influence factor analysis method for satellite radiation data
CN112884342B (en) * 2021-03-10 2024-03-12 陕西九州遥感信息技术有限公司 Quality evaluation and cross calibration method for water-color satellite atmospheric roof radiation product
CN113065098A (en) * 2021-03-23 2021-07-02 北京航天创智科技有限公司 PM2.5 concentration detection method and device based on satellite remote sensing and ground gas system data
CN113077083B (en) * 2021-03-26 2023-01-20 云南电网有限责任公司电力科学研究院 Method for monitoring pollution condition based on satellite remote sensing technology
CN113297527B (en) * 2021-06-09 2022-07-26 四川大学 PM based on multisource city big data 2.5 Overall domain space-time calculation inference method
CN113534203B (en) * 2021-07-12 2022-08-23 长光卫星技术股份有限公司 On-orbit cross-radiation calibration method based on AERONET aerosol data
CN113609940A (en) * 2021-07-26 2021-11-05 南通智能感知研究院 Wide-coverage soil heavy metal remote sensing mixed mapping method
CN114357885B (en) * 2022-01-06 2024-03-01 河南大学 Photosynthetic effective radiation scattering proportion prediction method fusing multisource data
CN115753687A (en) * 2022-12-23 2023-03-07 西安工程大学 Atmospheric pollution remote sensing monitoring method, system and computer readable medium
CN117150894B (en) * 2023-08-25 2024-08-13 中国气象科学研究院 Correction method, system and device for satellite observation data
CN118673336B (en) * 2024-08-23 2024-10-29 水利部交通运输部国家能源局南京水利科学研究院 Secondary screening-based medium-and-long-term precipitation prediction method and system integrating multi-head self-attention mechanism

Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103034910A (en) * 2012-12-03 2013-04-10 北京农业信息技术研究中心 Regional scale plant disease and insect pest prediction method based on multi-source information
CN106933776A (en) * 2017-03-02 2017-07-07 宁波大学 A kind of method that MODIS AOD products missing data is repaired

Patent Citations (2)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN103034910A (en) * 2012-12-03 2013-04-10 北京农业信息技术研究中心 Regional scale plant disease and insect pest prediction method based on multi-source information
CN106933776A (en) * 2017-03-02 2017-07-07 宁波大学 A kind of method that MODIS AOD products missing data is repaired

Non-Patent Citations (2)

* Cited by examiner, † Cited by third party
Title
Fast Image Interpolation via Random Forests;Jun-Jie Huang et al.;《IEEE Transactions on Image Processing》;20151031;第24卷(第10期);第3232-3245页 *
基于卫星遥感和气象再分析资料的北京市PM2.5浓度反演研究;邵琦 等;《地理与地理信息科学》;20180531;第34卷(第3期);第32-38页 *

Also Published As

Publication number Publication date
CN109213964A (en) 2019-01-15

Similar Documents

Publication Publication Date Title
CN109213964B (en) Satellite AOD product correction method fusing multi-source characteristic geographic parameters
CN105181898B (en) Atmospheric pollution monitoring and management method as well as system based on high-density deployment of sensors
CN112905560B (en) Air pollution prediction method based on multi-source time-space big data deep fusion
CN109344865B (en) Data fusion method for multiple data sources
CN110673229A (en) Atmospheric pollutant diffusion track tracking method based on hotspot grid technology
Son et al. The effect of particulate matter on solar photovoltaic power generation over the Republic of Korea
Hong et al. Inferring vertical variability and diurnal evolution of O3 formation sensitivity based on the vertical distribution of summertime HCHO and NO2 in Guangzhou, China
CN110595968B (en) PM2.5 concentration estimation method based on geostationary orbit satellite
CN112016696B (en) PM integrating satellite observation and ground observation 1 Concentration inversion method and system
KR102706021B1 (en) Forest fire susceptibility mapping method and apparatus using artificial intelligence
CN110595960B (en) PM2.5 concentration remote sensing estimation method based on machine learning
CN115081557A (en) Night aerosol optical thickness estimation method and system based on ground monitoring data
CN115236005A (en) Soil heavy metal content inversion method integrating spectral and spatial characteristics
CN113156395B (en) Aerosol laser radar data fusion method and system
CN116223395A (en) Near-surface trace gas concentration inversion model and inversion method
CN114219345B (en) Secondary air quality prediction optimization method based on data mining
CN114707396A (en) All-time PM2.5Near real-time production method of concentration seamless lattice point data
CN117219183A (en) High coverage near ground NO in cloudy rain areas 2 Concentration estimation method and system
KR102002593B1 (en) Method and apparatus for analyzing harmful gas diffusion in a specific space
MBAYE ARMA model for short-term forecasting of solar potential: application to a horizontal surface of Dakar site.
CN115937714A (en) Net Primary Productivity Estimation Method Based on Neural Network and CASA Model
CN114357885B (en) Photosynthetic effective radiation scattering proportion prediction method fusing multisource data
CN113311509A (en) Sensitivity test method of MWHTS (metal wrap-through magnetic field resonance) to sea surface air pressure based on neural network
Talib et al. Relationship between Particulate Matter (PM2. 5) Concentration and Aerosol Optical Depth (AOD) throughout Peninsular Malaysia
CN114581043B (en) Cross-network modular multi-source multi-factor wave-condition analysis system

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant