Aerosol and gas profile inversion method and system based on MAX-DOAS
By constructing a Bayesian inversion framework using the MCMC Bayesian method, the problems of low vertical resolution and local optima in the MAX-DOAS inversion algorithm are solved, achieving high-resolution aerosol and gas profile inversion and improving inversion accuracy and reliability.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Applications(China)
- Current Assignee / Owner
- HEFEI INSTITUTE OF PHYSICAL SCIENCE CHINESE ACADEMY OF SCIENCES
- Filing Date
- 2025-12-30
- Publication Date
- 2026-05-12
AI Technical Summary
Existing MAX-DOAS-based inversion algorithms suffer from low vertical resolution, local optima, large linearization errors, and incomplete uncertain estimations, resulting in large inversion errors and failing to meet the needs of urban-scale pollution monitoring and refined research.
A Bayesian inversion framework is constructed using the MCMC Bayesian method. The degree of agreement between observations and the forward model is quantified by the likelihood function. The posterior probability distribution of aerosol and gas profiles is inverted using the MCMC sampling algorithm. The measurement error, smoothing error and total error are mathematically strictly separated by the Bayesian error decomposition algorithm. A two-layer MCMC inversion engine is used for global search.
This method achieves an improvement in the resolution of 100m aerosol profile inversion, enhances inversion accuracy and reliability, provides complete uncertainty quantification information, solves the local optimum problem of gradient descent method, and ensures the global optimum and reliability of inversion results.
Smart Images
Figure CN122016676A_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of environmental monitoring technology, and in particular to a method and system for aerosol and gas profile inversion based on MAX-DOAS. Background Technology
[0002] In today's rapidly developing economy, human activities such as industrialization and urbanization have had a detrimental impact on the environment. Aerosols and trace gases are important components of atmospheric pollutants, playing a crucial role in atmospheric physical processes. MAX-DOAS (Multiaxial Differential Absorption Spectroscopy) is a ground-based remote sensing technology that acquires oblique column concentration information of aerosols and trace gases in the atmosphere by measuring solar scattered light at different angles. This information is then matched with a vertical distribution inversion algorithm to achieve gas profile inversion, enabling the detection of atmospheric pollutants.
[0003] Currently, the inversion algorithms based on MAX-DOAS are mainly divided into table lookup methods and optimal estimation methods, both of which still have some obvious shortcomings: 1. Existing MAX-DOAS inversion algorithms generally have low vertical resolution, around 200m, resulting in insufficient near-surface atmospheric stratification and failing to meet the needs of urban-scale pollution monitoring and refined research, particularly in distinguishing the complex structures within the boundary layer. Furthermore, due to the low vertical resolution, smoothing errors account for a large proportion of the inversion error, reaching over 70%, significantly limiting the accuracy of the inversion algorithm.
[0004] 2. Existing MAX-DOAS aerosol and gas profile inversion algorithms based on the optimal estimation method basically all use the Gauss-Newton method to iteratively solve the value function, but they have the following serious drawbacks: 1) Local Optimal Solution: When the optimal estimation method is used to solve for the optimal solution, it will use the gradient descent method. When there are many atmospheric layers, singular values will be generated and it will get stuck in the local optimal solution, making it impossible to continue to calculate the global optimal solution. Moreover, the convergence result of the gradient method is strongly dependent on the initial guess value. Different initial values may lead to completely different inversion results, and there is a lack of objective initial value selection criteria.
[0005] 2) Linearization error: The optimal estimation method uses Taylor's formula for expansion, assuming that the forward model can be linearized near the prior point. However, the atmospheric radiation transfer process has strong nonlinear characteristics. When the linearization error is large, the Gauss-Newton method iterative process is prone to divergence or false convergence, causing the traditional convergence criterion (gradient magnitude less than the threshold) to fail in nonlinear systems. The algorithm may stop at the saddle point or oscillate between multiple local extrema.
[0006] 3) Incomplete uncertainty estimation: Existing gas inversion error analysis assumes that the posterior distribution is Gaussian, but the actual posterior distribution often has non-Gaussian characteristics such as multimodal, skewed, and thick tails. This leads to serious distortion of the confidence interval, and the actual coverage is far lower than the calculated result. It is impossible to accurately identify the main source of error, and the error prediction results under different observation conditions are unreliable.
[0007] Patent document CN113834792 proposes a 50-meter resolution trace gas profile inversion method based on MAX-DOAS. It obtains the aerosol profile through a first inversion calculation using a Monte Carlo sampling algorithm. The aerosol profile, combined with the gas column concentration and prior gas profile, is then used in a second inversion calculation using the Monte Carlo sampling algorithm to obtain the trace gas profile. This invention uses the Monte Carlo method to replace the Gauss-Newton iterative solution. However, it employs a Monte Carlo random sampling method, randomly generating a weighting factor α in the zero-to-one interval, calculating the value function for each α value, and sorting them to find the optimal solution. The search space is a one-dimensional scalar parameter, and based on a preset number of samplings, there is no clear convergence diagnostic criterion. Therefore, a certain degree of inversion error still exists. Summary of the Invention
[0008] The technical problem to be solved by this invention is how to reduce the inversion error of aerosol and trace gas profile inversion based on MAX-DOAS inversion algorithm.
[0009] To solve the above-mentioned technical problems, the present invention adopts the following technical solution: a method for inverting aerosol and gas profiles based on MAX-DOAS, comprising the following steps: S1. Input the original MAX-DOAS spectral data, obtain the O4 oblique column concentration and trace gas oblique column concentration through spectral fitting, construct a Bayesian inversion framework, quantify the degree of agreement between the observation and the forward model through the likelihood function, use the MCMC sampling algorithm to invert the posterior probability distribution of the aerosol extinction profile, and extract the posterior mean from the posterior probability distribution of the aerosol extinction profile as the optimal aerosol extinction profile. S2. Fix the optimal aerosol extinction profile obtained in step S1, combine it with the trace gas oblique column concentration, and invert the posterior probability distribution of the gas concentration profile through the conditional likelihood function, the prior probability distribution of gas concentration and the MCMC sampling algorithm. Extract the posterior mean from the posterior probability distribution of the trace gas concentration profile as the optimal trace gas concentration profile. S3. Using the Bayesian error decomposition algorithm, combined with the posterior probability distribution of the aerosol extinction profile obtained in step S1 and the posterior probability distribution of the trace gas concentration profile obtained in step S2, the measurement error, smoothing error, and total error are strictly separated mathematically, and the inversion progress is monitored in real time and anomaly diagnosis is performed.
[0010] As a further optimized technical solution, the main flow of the MCMC sampling algorithm used in steps S1 and S2 is as follows: Input the oblique column concentration and prior profile into the Bayesian inversion framework to construct the posterior probability distribution; Using the posterior probability as the target distribution, the MCMC sampling algorithm is used to randomly walk in the parameter space, and the radiative transfer model is called to calculate the forward simulation value. The new state is accepted or rejected based on the degree of agreement between the simulation and the observation. After iterative convergence, the posterior sample set is obtained. The inversion profile is obtained by extracting the posterior mean, standard deviation and confidence interval from the posterior sample set through statistical analysis.
[0011] As a further optimized technical solution, the posterior probability is used as the target distribution. The MCMC sampling algorithm is used to randomly walk in the parameter space, and the radiative transfer model is called to calculate the forward simulation value. The acceptance or rejection of the new state is determined based on the degree of agreement between the simulation and observation. After iterative convergence, the posterior sample set is obtained. The specific process of extracting the posterior mean, standard deviation, and confidence interval through statistical analysis of the posterior sample set to obtain the inversion profile is as follows: Initially, several walkers are set up, and the initial state is generated by randomly perturbing them near the prior mean. The number of sampling steps and the number of burn-in steps are set. Perform a stretching movement for each walker, randomly select another walker and propose a new state, and call the radiative transfer model to calculate the likelihood function; The posterior probability and acceptance probability are calculated using the likelihood function and prior probability distribution. Random numbers between 0 and 1 are generated, and the relationship between the random numbers and the acceptance probability is determined. The calculation is repeated according to the set number of sampling steps until convergence. The convergent posterior sample is retained, the posterior statistic is calculated as the inversion profile, the standard deviation is used as the uncertainty, and the quantiles are extracted as the confidence interval.
[0012] As a further optimized technical solution, the MCMC sampling algorithms in steps S1 and S2 construct a two-layer MCMC inversion engine. The posterior probability distribution of the aerosol extinction profile is obtained through the first-layer MCMC sampling algorithm, and the posterior probability distribution of the trace gas concentration profile is obtained through the second-layer MCMC sampling algorithm.
[0013] As a further optimized technical solution, the specific execution steps of the two-layer MCMC inversion engine include: The O4 oblique column concentration and aerosol extinction prior profile are input into the first layer MCMC sampling algorithm of the two-layer MCMC inversion engine to construct a Bayesian inversion framework. The posterior probability distribution of the aerosol extinction profile is obtained by quantifying the degree of agreement between the observation and the forward model through the likelihood function. The posterior mean is extracted from the posterior probability distribution of the aerosol extinction profile as the optimal aerosol extinction profile. Next, the optimal aerosol extinction profile obtained in the first layer is used as a fixed parameter. Combined with the trace gas oblique column concentration data, it is input into the second layer MCMC sampling algorithm. Based on Bayes' theorem, combined with the conditional likelihood function and the prior probability distribution of gas concentration, the posterior probability distribution of the gas concentration profile is obtained through the MCMC sampling algorithm. The MCMC sampling chain is then subjected to burn-in processing and convergence diagnosis. The posterior mean of the posterior probability distribution of the gas concentration profile is calculated as the optimal trace gas concentration profile. The vertical profiles of aerosol and gas are obtained through vertical integration.
[0014] As a further optimized technical solution, the first-layer MCMC sampling algorithm specifically includes: The O4 oblique column concentration and aerosol extinction prior profile are input into the forward modeling module, and the weighting function of air quality factor and aerosol extinction profile of each layer is calculated by the SCIATRAN radiative transfer model. Based on the calculated weight function and the prior probability distribution of aerosols, the posterior probability distribution function of the aerosol extinction profile is constructed according to Bayes' theorem. The Metropolis-Hastings algorithm is used for random sampling. Samples are drawn from the posterior probability distribution of the aerosol extinction profile through the acceptance-rejection criterion. After burn-in processing and Gelman-Rubin convergence diagnosis, the posterior sample set of the aerosol extinction profile is obtained. The posterior mean is calculated from the posterior sample set of aerosol extinction profiles as the estimate of the optimal aerosol extinction profile. The posterior covariance matrix is calculated to quantify the uncertainty. The optimal aerosol extinction profile is substituted into the SCIATRAN radiative transfer model to verify the consistency between the simulated O4 oblique column concentration and the observation. The optimal aerosol extinction profile and its complete uncertainty quantification results are output.
[0015] As a further optimized technical solution, the second-layer MCMC sampling algorithm specifically includes: The optimal aerosol extinction profile obtained in the first layer is used as a fixed parameter. The trace gas oblique column concentration, the fixed aerosol extinction profile and the trace gas prior profile are input into the forward modeling module. The sensitivity weighting function of the trace gas oblique column concentration to the gas concentration of each layer is calculated by the SCIATRAN radiative transfer model. Based on the sensitivity weighting function and the prior profile of trace gas concentration, a conditional posterior probability distribution is constructed. The Stretch Move ensemble sampling algorithm combined with the Metropolis-Hastings criterion is used for sampling. Through the collaborative optimization of multiple sampling chains, the posterior sample set of the trace gas concentration profile is obtained. Extract the posterior mean from the posterior sample of the trace gas concentration profile as the optimal trace gas concentration profile, calculate the posterior covariance to quantify the uncertainty, perform vertical integration between the optimal trace gas concentration profile and the atmospheric density profile, calculate the vertical column concentration of trace gas, and output the trace gas concentration profile and its uncertainty quantification results.
[0016] As a further optimized technical solution, the posterior probability distribution of the aerosol extinction profile and the optimal aerosol extinction profile inversion process include: Constructing a Bayesian inversion framework:
[0017] Where y is the observation vector, which is the O4 column concentration obtained from MAX-DOAS observation; x is the state vector, which is the inversion profile to be obtained in the end. Let be the posterior probability distribution of the aerosol extinction profile, and let represent the probability distribution of x obtained after observing the observation vector y. Let x be the likelihood function, representing the probability of observing the observation vector y given the true state vector x. The prior probability distribution is the state vector x that is inferred from existing experience; This is the normalization constant; According to the above formula, the likelihood function Adopting a Gaussian distribution:
[0018] in, Map the state vector x to the observation space, and use the SCIATRAN radiative transfer model in the above formula; For the covariance of the observed data error, is the inverse matrix of the error covariance, and is the weight matrix of the residuals; the larger the error, the smaller the weight. Based on the above formula, construct the likelihood function. Logarithmic form:
[0019] in, Let i be the aerosol observation value at the i-th elevation angle. These are the model observations corresponding to the elevation angle. This is the observation error; Based on the Bayesian inversion framework described above, the constructed prior probability distribution adopts a multivariate Gaussian distribution:
[0020] Among them, the prior mean This indicates that a typical aerosol extinction profile already exists. Let be the prior covariance matrix, representing the uncertainty and vertical correlation of the profile; Based on the above algorithm, calculate the logarithmic form of the posterior probability distribution of the aerosol extinction profile:
[0021] in This is the weighted average of the residuals from aerosol observations.
[0022] Based on the above algorithm, calculate the posterior mean of the posterior probability distribution of the aerosol extinction profile:
[0023] in, is the posterior mean of the aerosol extinction profile, represents the optimal aerosol extinction profile, and N is the total number of sampling steps in the first layer of MCMC.
[0024] As a further optimized technical solution, the posterior probability distribution of the gas concentration profile and the inversion process of the optimal trace gas concentration profile include: Calculate the posterior probability distribution of the gas concentration profile:
[0025] in, The final gas profile distribution needs to be obtained through inversion. The posterior mean of the aerosol extinction profile obtained from the inversion step is given. To obtain the gas profile that needs to be inverted, To observe the concentration of trace gases in the oblique column, Let be the likelihood probability function for the gas. Let be the prior probability distribution of the gas; Based on the above algorithm, calculate the posterior mean of the posterior probability distribution of the trace gas concentration profile:
[0026] in, Let be the posterior mean of the trace gas concentration profile, represent the optimal trace gas concentration profile, and M be the total number of sampling steps in the second layer MCMC.
[0027] This invention also provides an aerosol and gas profile inversion system based on MAX-DOAS, comprising: The posterior probability distribution inversion module for aerosol extinction profiles is used to input the original MAX-DOAS spectral data, obtain the O4 oblique column concentration and trace gas oblique column concentration through spectral fitting, construct a Bayesian inversion framework, quantify the degree of agreement between observations and the forward model through the likelihood function, invert the posterior probability distribution of aerosol extinction profiles using the MCMC sampling algorithm, and extract the posterior mean from the posterior probability distribution of aerosol extinction profiles as the optimal aerosol extinction profile. The posterior probability distribution inversion module of the gas concentration profile is used to fix the optimal aerosol extinction profile obtained by the above module, and combine the trace gas oblique column concentration to invert the posterior probability distribution of the gas concentration profile through the conditional likelihood function, the prior probability distribution of gas concentration and the MCMC sampling algorithm. The posterior mean is extracted from the posterior probability distribution of the trace gas concentration profile as the optimal trace gas concentration profile. The monitoring module uses the Bayesian error decomposition algorithm, combined with the posterior probability distribution of the obtained aerosol extinction profile and the obtained trace gas concentration profile, to achieve strict mathematical separation of measurement error, smoothing error and total error, monitor the inversion progress in real time and perform anomaly diagnosis.
[0028] Beneficial Effects: This invention achieves 100m aerosol profile inversion for the first time, improving the resolution from 200m to 100m and significantly enhancing inversion accuracy. It employs the MCMC Bayesian method instead of traditional optimal estimation, fundamentally solving the problem of gradient descent algorithms getting trapped in local optima. By using an ensemble sampler to perform a global search in the high-dimensional parameter space, it can explore all reasonable solution space regions, avoiding strong dependence on initial values and ensuring convergence to the global optimum or a near-global optimum solution set. Based on the Bayesian inversion framework, this invention naturally obtains the posterior probability distribution of the inversion results, providing complete uncertainty quantification information. Furthermore, this invention provides confidence intervals, correlation matrices, and MCMC convergence diagnostic indicators for each height level, offering a more scientific and comprehensive basis for result reliability assessment and improving the accuracy and reliability of the inversion results.
[0029] The following differences exist between this document and the prior art document with publication number CN113834792: (1) Different sampling algorithms: The comparison file uses the Monte Carlo random sampling method, which randomly generates a weight factor α in the interval between zero and one, calculates the value function for each α value and sorts them to find the optimal solution, and the search space is a one-dimensional scalar parameter; the present invention uses the Stretch Move ensemble sampling algorithm, which runs sixty-four to one hundred and twenty-eight independent sampling chains for collaborative sampling, and the suggested steps are based on the positional relationship between the chains (x* k = x j + z(x k - x jBy using the Metropolis-Hastings criterion for probabilistic acceptance judgment and searching in a high-dimensional profile parameter space, the accuracy of the inversion results is further improved.
[0030] (2) Different mathematical frameworks: The comparative document is based on the value function minimization framework, which uses a weight factor α to linearly interpolate between the prior profile and the model weight function, and the solution is in the form of x. a = VCD × [α×K / ||K|| + (1-α)×x p / ||x p ||];This invention is based on a complete Bayesian inference framework, constructs a posterior probability density function P(x|y) ∝ P(y|x) ×P(x), obtains the posterior probability distribution through MCMC sampling, and the solution is in the form of the statistics of the posterior sample set.
[0031] (3) Different output information: The comparison file outputs the best line type and the corresponding profile point estimate, and calculates the uncertainty through the error propagation formula; the output of this invention includes the posterior mean (optimal estimate), the complete posterior sample set, and three independent components based on Bayesian error decomposition: measurement error, smoothing error and model error.
[0032] (4) Different convergence diagnosis: The comparison file is based on a preset number of samplings and has no clear convergence diagnosis criteria; this invention uses Gelman-Rubin multi-chain convergence diagnosis. It employs multiple diagnostic mechanisms, including statistical measures, effective sample size assessment, and autocorrelation analysis, to provide quantitative convergence judgment and automatic early termination functions. Attached Figure Description Figure 1 This is a flowchart of the aerosol and gas profile inversion method based on MAX-DOAS according to an embodiment of the present invention; Figure 2 This is a flowchart of the main MCMC algorithm described in an embodiment of the present invention; Figure 3 This is a flowchart illustrating the computation process of the MCMC algorithm described in an embodiment of the present invention. Detailed Implementation To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.
[0033] See Figure 1 This embodiment presents a method for inverting aerosol and gas profiles based on MAX-DOAS using the MCMC (Markov Chain Monte Carlo) algorithm, comprising the following steps: S1. Input the original MAX-DOAS spectral data, obtain the O4 oblique column concentration and trace gas oblique column concentration through spectral fitting, construct a Bayesian inversion framework, quantify the degree of agreement between the observation and the forward model through the likelihood function, use the MCMC sampling algorithm to invert the posterior probability distribution of the aerosol extinction profile, and extract the posterior mean from the posterior probability distribution of the aerosol extinction profile as the optimal aerosol extinction profile. S2. Fix the optimal aerosol extinction profile obtained in step S1, combine it with the trace gas oblique column concentration, and invert the posterior probability distribution of the gas concentration profile through the conditional likelihood function, the prior probability distribution of gas concentration and the MCMC sampling algorithm. Extract the posterior mean from the posterior probability distribution of the trace gas concentration profile as the optimal trace gas concentration profile. S3. Using the Bayesian error decomposition algorithm, combined with the posterior probability distribution of the aerosol extinction profile obtained in step S1 and the posterior probability distribution of the trace gas concentration profile obtained in step S2, the measurement error, smoothing error, and total error are strictly separated mathematically, and the inversion progress is monitored in real time and anomaly diagnosis is performed.
[0034] In the above steps, the MCMC sampling algorithms in steps S1 and S2 construct a two-layer MCMC inversion engine. The posterior probability distribution of the aerosol extinction profile is obtained through the first-layer MCMC sampling algorithm, and the posterior probability distribution of the trace gas concentration profile is obtained through the second-layer MCMC sampling algorithm.
[0035] The specific execution steps of the two-layer MCMC inversion engine include: The O4 oblique column concentration and aerosol extinction prior profile are input into the first layer MCMC sampling algorithm of the two-layer MCMC inversion engine to construct a Bayesian inversion framework. The posterior probability distribution of the aerosol extinction profile is obtained by quantifying the degree of agreement between the observation and the forward model through the likelihood function. The posterior mean is extracted from the posterior probability distribution of the aerosol extinction profile as the optimal aerosol extinction profile. Next, the optimal aerosol extinction profile obtained in the first layer is used as a fixed parameter. Combined with the trace gas oblique column concentration data, it is input into the second layer MCMC sampling algorithm. Based on Bayes' theorem, combined with the conditional likelihood function and the prior probability distribution of gas concentration, the posterior probability distribution of the gas concentration profile is obtained through the MCMC sampling algorithm. The MCMC sampling chain is then subjected to burn-in processing and convergence diagnosis. The posterior mean of the posterior probability distribution of the gas concentration profile is calculated as the optimal trace gas concentration profile. The vertical profiles of aerosol and gas are obtained through vertical integration.
[0036] The first-layer MCMC sampling algorithm specifically includes: The O4 oblique column concentration and aerosol extinction prior profile are input into the forward modeling module, and the weighting function of air quality factor and aerosol extinction profile of each layer is calculated by the SCIATRAN radiative transfer model. Based on the calculated weight function and the prior probability distribution of aerosols, the posterior probability distribution function of the aerosol extinction profile is constructed according to Bayes' theorem. The Metropolis-Hastings algorithm is used for random sampling. Samples are drawn from the posterior probability distribution of the aerosol extinction profile through the acceptance-rejection criterion. After burn-in processing and Gelman-Rubin convergence diagnosis, the posterior sample set of the aerosol extinction profile is obtained. The posterior mean is calculated from the posterior sample set of aerosol extinction profiles as the estimate of the optimal aerosol extinction profile. The posterior covariance matrix is calculated to quantify the uncertainty. The optimal aerosol extinction profile is substituted into the SCIATRAN radiative transfer model to verify the consistency between the simulated O4 oblique column concentration and the observation. The optimal aerosol extinction profile and its complete uncertainty quantification results are output.
[0037] The second-layer MCMC sampling algorithm specifically includes: The optimal aerosol extinction profile obtained in the first layer is used as a fixed parameter. The trace gas oblique column concentration, the fixed aerosol extinction profile and the trace gas prior profile are input into the forward modeling module. The sensitivity weighting function of the trace gas oblique column concentration to the gas concentration of each layer is calculated by the SCIATRAN radiative transfer model. Based on the sensitivity weighting function and the prior profile of trace gas concentration, a conditional posterior probability distribution is constructed. The Stretch Move ensemble sampling algorithm combined with the Metropolis-Hastings criterion is used for sampling. Through the collaborative optimization of multiple sampling chains, the posterior sample set of the trace gas concentration profile is obtained. Extract the posterior mean from the posterior sample of the trace gas concentration profile as the optimal trace gas concentration profile, calculate the posterior covariance to quantify the uncertainty, perform vertical integration between the optimal trace gas concentration profile and the atmospheric density profile, calculate the vertical column concentration of trace gas, and output the trace gas concentration profile and its uncertainty quantification results.
[0038] The following is a detailed explanation: See Figure 2 In a specific example, the main flow of the MCMC sampling algorithm used in steps S1 and S2 is as follows: Input the oblique column concentration (O4 oblique column concentration or trace gas oblique column concentration) and the prior profile (aerosol extinction prior profile or trace gas prior profile) into the Bayesian inversion framework to construct the posterior probability distribution; Using the posterior probability as the target distribution, the MCMC sampling algorithm is used to randomly walk in the parameter space, and the radiative transfer model is called to calculate the forward simulation value. The new state is accepted or rejected based on the degree of agreement between the simulation and the observation. After iterative convergence, the posterior sample set is obtained. The inversion profile is obtained by extracting the posterior mean, standard deviation and confidence interval from the posterior sample set through statistical analysis.
[0039] It should be noted that various traditional inversion algorithms use gradient descent to find extrema, which is prone to getting trapped in local optima in high-dimensional parameter spaces and is highly sensitive to the choice of initial values. This embodiment innovatively uses the MCMC method based on the Bayesian inversion framework to solve the MAX-DOAS inversion problem. By randomly walking in the parameter space, it can explore the global parameter space across local extrema, fundamentally solving the local extrema trap problem.
[0040] In addition, traditional methods rely on the linearization approximation of the forward model near the prior point. When the aerosol optical thickness is high or the gas concentration is large, the strong nonlinear characteristics of radiative transfer lead to the accumulation of linearization errors, which may cause iterative divergence or convergence to an incorrect solution. In contrast, the MCMC method directly calls the complete SCIATRAN radiative transfer model to calculate the forward simulation value without any linearization assumptions, and can accurately handle nonlinear effects such as multiple scattering and absorption saturation in radiative transfer.
[0041] Furthermore, traditional methods assume that the posterior distribution is Gaussian, and can only provide single-point estimates and approximate uncertainties based on linearization, failing to characterize the true form of the posterior distribution. In contrast, the MCMC method directly constructs a complete sample set of the posterior distribution through sampling, which can accurately reflect the non-Gaussian characteristics of the distribution, such as multimodality, skewness, and heavy tails, providing a rigorous statistical basis for error decomposition and uncertainty quantification.
[0042] See Figure 3 In a specific example, the process involves using the posterior probability as the target distribution, randomly walking through the parameter space using the MCMC sampling algorithm, calculating the forward simulated value using the radiative transfer model, and deciding whether to accept or reject the new state based on the degree of agreement between the simulation and observation. After iterative convergence, a posterior sample set is obtained. The specific procedure for extracting the posterior mean, standard deviation, and confidence interval through statistical analysis of the posterior sample set to obtain the inversion profile is as follows: Initially, several walkers are set up, and the initial state is generated by randomly perturbing them near the prior mean. The number of sampling steps and the number of burn-in steps are set. Perform a stretching movement for each walker, randomly select another walker and propose a new state, and call the radiative transfer model to calculate the likelihood function; The posterior probability and acceptance probability are calculated using the likelihood function and prior probability distribution. Random numbers between 0 and 1 are generated, and the relationship between the random numbers and the acceptance probability is determined. The calculation is repeated according to the set number of sampling steps until convergence. The convergent posterior sample is retained, the posterior statistic is calculated as the inversion profile, the standard deviation is used as the uncertainty, and the quantiles are extracted as the confidence interval. It should be noted that this embodiment employs an intelligent caching optimization mechanism, which significantly reduces the number of model calls by reusing forward simulation results of similar profiles, thus compressing the single-time inversion time to a shorter duration. Furthermore, the updates of each walker during MCMC sampling are independent, making it suitable for parallel computing. This fully utilizes multi-core CPU resources, further improving computational efficiency. This embodiment, through caching optimization and parallelization, combined with its advantages of global sampling and complete uncertainty quantification, significantly enhances the algorithm's robustness and result reliability.
[0043] In a specific example, the detailed process of the MAX-DOAS aerosol and gas profile inversion method based on the MCMC algorithm is as follows: Constructing a Bayesian inversion framework:
[0044] Where y is the observation vector, which is the O4 column concentration obtained from MAX-DOAS observation; x is the state vector, which is the inversion profile to be obtained in the end. Let be the posterior probability distribution of the aerosol extinction profile, and let represent the probability distribution of x obtained after observing the observation vector y. Let x be the likelihood function, representing the probability of observing the observation vector y given the true state vector x. The prior probability distribution is the state vector x that is inferred from existing experience; This is the normalization constant; According to the above formula, the likelihood function Adopting a Gaussian distribution:
[0045] in, Map the state vector x to the observation space, and use the SCIATRAN radiative transfer model in the above formula; For the covariance of the observed data error, is the inverse matrix of the error covariance, and is the weight matrix of the residuals; the larger the error, the smaller the weight. Based on the above formula, construct the likelihood function. Logarithmic form:
[0046] in, Let i be the aerosol observation value at the i-th elevation angle. These are the model observations corresponding to the elevation angle. This is the observation error; Based on the Bayesian inversion framework described above, the constructed prior probability distribution adopts a multivariate Gaussian distribution:
[0047] Among them, the prior mean This indicates that a typical aerosol extinction profile already exists. Let be the prior covariance matrix, representing the uncertainty and vertical correlation of the profile; Based on the above algorithm, calculate the logarithmic form of the posterior probability distribution of the aerosol extinction profile:
[0048] in This is the weighted average of the residuals from aerosol observations.
[0049] Based on the above algorithm, calculate the posterior mean of the posterior probability distribution of the aerosol extinction profile:
[0050] in, is the posterior mean of the aerosol extinction profile, represents the optimal aerosol extinction profile, and N is the total number of sampling steps in the first layer of MCMC.
[0051] Based on the above algorithm, calculate the posterior probability distribution of the gas concentration profile:
[0052] in, The final gas profile distribution needs to be obtained through inversion. The posterior mean of the aerosol extinction profile obtained from the above inversion steps is given. To obtain the gas profile that needs to be inverted, To observe the concentration of trace gases in the oblique column, Let be the likelihood probability function for the gas. Let be the prior probability distribution of the gas.
[0053] Based on the above algorithm, calculate the posterior mean of the posterior probability distribution of the trace gas concentration profile:
[0054] in, Let be the posterior mean of the trace gas concentration profile, represent the optimal trace gas concentration profile, and M be the total number of sampling steps in the second layer MCMC.
[0055] According to the above sampling algorithm, the MCMC sampling adopts an integrated affine invariant sampler (Emcee), which includes initializing K random walkers, K≥2n+2, where n is the parameter dimension; using the stretch move method to calculate the posterior probability of the new state; deciding whether to accept or reject the new state according to the Metropolis-Hastings criterion; iterating for 50 to 200 steps until convergence to obtain the posterior distribution sample; To determine whether to accept the new state, a formula for calculating the acceptance rate is constructed:
[0056] in, For the proposed new state, Let z be the current state, z be the stretching factor, and n be the parameter dimension.
[0057] The real-time inversion detection indicators include: The required acceptance rate for aerosol inversion is 0.2-0.7, and for gas inversion it is 0.15-0.65. Based on the Bayesian error decomposition algorithm, the total error of this inversion algorithm satisfies:
[0058]
[0059]
[0060] in, To account for measurement error, To smooth the error, and with each error component orthogonal to the others in the vertical direction; gain matrix , The posterior covariance matrix is... This is the transpose of the Jacobian matrix. The observation error covariance matrix originates from instrument noise, spectral fitting errors, etc.; the average kernel matrix... K, is the identity matrix, and K is the Jacobian matrix, used to describe the sensitivity of different height layers under different observation angles.
[0061] In summary, the embodiments of this invention achieve high-resolution aerosol and gas profile inversion, improving the aerosol resolution from 200m to 100m and the gas resolution from 200m to 50m, significantly enhancing the accuracy of the inversion. This invention employs a novel, more accurate solution algorithm, using the MCMC method based on a Bayesian inversion framework, avoiding the local optimum problem caused by traditional optimal estimation methods, while simultaneously improving the accuracy and reliability of the computational solution.
[0062] The present invention also provides a MAX-DOAS-based aerosol and gas profile inversion system corresponding to the above-mentioned MAX-DOAS-based aerosol and gas profile inversion method, comprising: The posterior probability distribution inversion module for aerosol extinction profiles is used to input the original MAX-DOAS spectral data, obtain the O4 oblique column concentration and trace gas oblique column concentration through spectral fitting, construct a Bayesian inversion framework, quantify the degree of agreement between observations and the forward model through the likelihood function, invert the posterior probability distribution of aerosol extinction profiles using the MCMC sampling algorithm, and extract the posterior mean from the posterior probability distribution of aerosol extinction profiles as the optimal aerosol extinction profile. The posterior probability distribution inversion module of the gas concentration profile is used to fix the optimal aerosol extinction profile obtained by the above module, and combine the trace gas oblique column concentration to invert the posterior probability distribution of the gas concentration profile through the conditional likelihood function, the prior probability distribution of gas concentration and the MCMC sampling algorithm. The posterior mean is extracted from the posterior probability distribution of the trace gas concentration profile as the optimal trace gas concentration profile. The monitoring module uses the Bayesian error decomposition algorithm, combined with the posterior probability distribution of the obtained aerosol extinction profile and the obtained trace gas concentration profile, to achieve strict mathematical separation of measurement error, smoothing error and total error, monitor the inversion progress in real time and perform anomaly diagnosis.
[0063] The specific steps performed in each module are the same as those described above for the aerosol and gas profile inversion method based on MAX-DOAS.
[0064] The above embodiments are only used to illustrate the technical solutions of the present invention, and are not intended to limit it. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some of the technical features. Such modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the spirit and scope of the technical solutions of the embodiments of the present invention.
Claims
1. A method for inverting aerosol and gas profiles based on MAX-DOAS, characterized in that: Includes the following steps: S1. Input the original MAX-DOAS spectral data, obtain the O4 oblique column concentration and trace gas oblique column concentration through spectral fitting, construct a Bayesian inversion framework, quantify the degree of agreement between the observation and the forward model through the likelihood function, use the MCMC sampling algorithm to invert the posterior probability distribution of the aerosol extinction profile, and extract the posterior mean from the posterior probability distribution of the aerosol extinction profile as the optimal aerosol extinction profile. S2. Fix the optimal aerosol extinction profile obtained in step S1, combine it with the trace gas oblique column concentration, and invert the posterior probability distribution of the gas concentration profile through the conditional likelihood function, the prior probability distribution of gas concentration and the MCMC sampling algorithm. Extract the posterior mean from the posterior probability distribution of the trace gas concentration profile as the optimal trace gas concentration profile. S3. Using the Bayesian error decomposition algorithm, combined with the posterior probability distribution of the aerosol extinction profile obtained in step S1 and the posterior probability distribution of the trace gas concentration profile obtained in step S2, the measurement error, smoothing error, and total error are strictly separated mathematically, and the inversion progress is monitored in real time and anomaly diagnosis is performed.
2. The aerosol and gas profile inversion method based on MAX-DOAS as described in claim 1, characterized in that: The main flow of the MCMC sampling algorithm used in steps S1 and S2 is as follows: Input the oblique column concentration and prior profile into the Bayesian inversion framework to construct the posterior probability distribution; Using the posterior probability as the target distribution, the MCMC sampling algorithm is used to randomly walk in the parameter space, and the radiative transfer model is called to calculate the forward simulation value. The new state is accepted or rejected based on the degree of agreement between the simulation and the observation. After iterative convergence, the posterior sample set is obtained. The inversion profile is obtained by extracting the posterior mean, standard deviation and confidence interval from the posterior sample set through statistical analysis.
3. The aerosol and gas profile inversion method based on MAX-DOAS as described in claim 2, characterized in that: The process of using the posterior probability as the target distribution, randomly walking in the parameter space through the MCMC sampling algorithm, calculating the forward simulation value by calling the radiative transfer model, and deciding whether to accept or reject the new state based on the degree of agreement between the simulation and observation, and obtaining the posterior sample set after iterative convergence; and extracting the posterior mean, standard deviation and confidence interval through statistical analysis of the posterior sample set to obtain the inversion profile is as follows: Initially, several walkers are set up, and the initial state is generated by randomly perturbing them near the prior mean. The number of sampling steps and the number of burn-in steps are set. Perform a stretching movement for each walker, randomly select another walker and propose a new state, and call the radiative transfer model to calculate the likelihood function; The posterior probability and acceptance probability are calculated using the likelihood function and prior probability distribution. Random numbers between 0 and 1 are generated, and the relationship between the random numbers and the acceptance probability is determined. The calculation is repeated according to the set number of sampling steps until convergence. The convergent posterior sample is retained, the posterior statistic is calculated as the inversion profile, the standard deviation is used as the uncertainty, and the quantiles are extracted as the confidence interval.
4. The aerosol and gas profile inversion method based on MAX-DOAS as described in claim 1, characterized in that: The MCMC sampling algorithms in steps S1 and S2 construct a two-layer MCMC inversion engine. The posterior probability distribution of the aerosol extinction profile is obtained through the first-layer MCMC sampling algorithm, and the posterior probability distribution of the trace gas concentration profile is obtained through the second-layer MCMC sampling algorithm.
5. The aerosol and gas profile inversion method based on MAX-DOAS as described in claim 4, characterized in that: The specific execution steps of the two-layer MCMC inversion engine include: The O4 oblique column concentration and aerosol extinction prior profile are input into the first layer MCMC sampling algorithm of the two-layer MCMC inversion engine to construct a Bayesian inversion framework. The posterior probability distribution of the aerosol extinction profile is obtained by quantifying the degree of agreement between the observation and the forward model through the likelihood function. The posterior mean is extracted from the posterior probability distribution of the aerosol extinction profile as the optimal aerosol extinction profile. Next, the optimal aerosol extinction profile obtained in the first layer is used as a fixed parameter. Combined with the trace gas oblique column concentration data, it is input into the second layer MCMC sampling algorithm. Based on Bayes' theorem, combined with the conditional likelihood function and the prior probability distribution of gas concentration, the posterior probability distribution of the gas concentration profile is obtained through the MCMC sampling algorithm. The MCMC sampling chain is then subjected to burn-in processing and convergence diagnosis. The posterior mean of the posterior probability distribution of the gas concentration profile is calculated as the optimal trace gas concentration profile. The vertical profiles of aerosol and gas are obtained through vertical integration.
6. The aerosol and gas profile inversion method based on MAX-DOAS as described in claim 5, characterized in that: The first-layer MCMC sampling algorithm specifically includes: The O4 oblique column concentration and aerosol extinction prior profile are input into the forward modeling module, and the weighting function of air quality factor and aerosol extinction profile of each layer is calculated by the SCIATRAN radiative transfer model. Based on the calculated weight function and the prior probability distribution of aerosols, the posterior probability distribution function of the aerosol extinction profile is constructed according to Bayes' theorem. The Metropolis-Hastings algorithm is used for random sampling. Samples are drawn from the posterior probability distribution of the aerosol extinction profile through the acceptance-rejection criterion. After burn-in processing and Gelman-Rubin convergence diagnosis, the posterior sample set of the aerosol extinction profile is obtained. The posterior mean is calculated from the posterior sample set of aerosol extinction profiles as the estimate of the optimal aerosol extinction profile. The posterior covariance matrix is calculated to quantify the uncertainty. The optimal aerosol extinction profile is substituted into the SCIATRAN radiative transfer model to verify the consistency between the simulated O4 oblique column concentration and the observation. The optimal aerosol extinction profile and its complete uncertainty quantification results are output.
7. The aerosol and gas profile inversion method based on MAX-DOAS as described in claim 5, characterized in that: The second-layer MCMC sampling algorithm specifically includes: The optimal aerosol extinction profile obtained in the first layer is used as a fixed parameter. The trace gas oblique column concentration, the fixed aerosol extinction profile and the trace gas prior profile are input into the forward modeling module. The sensitivity weighting function of the trace gas oblique column concentration to the gas concentration of each layer is calculated by the SCIATRAN radiative transfer model. Based on the sensitivity weighting function and the prior profile of trace gas concentration, a conditional posterior probability distribution is constructed. The StretchMove ensemble sampling algorithm combined with the Metropolis-Hastings criterion is used for sampling. Through the collaborative optimization of multiple sampling chains, the posterior sample set of the trace gas concentration profile is obtained. Extract the posterior mean from the posterior sample of the trace gas concentration profile as the optimal trace gas concentration profile, calculate the posterior covariance to quantify the uncertainty, perform vertical integration between the optimal trace gas concentration profile and the atmospheric density profile, calculate the vertical column concentration of trace gas, and output the trace gas concentration profile and its uncertainty quantification results.
8. The aerosol and gas profile inversion method based on MAX-DOAS as described in claim 1, characterized in that: The posterior probability distribution of the aerosol extinction profile and the inversion process of the optimal aerosol extinction profile include: Constructing a Bayesian inversion framework: Where y is the observation vector, which is the O4 column concentration obtained from MAX-DOAS observation; x is the state vector, which is the inversion profile to be obtained in the end. Let be the posterior probability distribution of the aerosol extinction profile, and let represent the probability distribution of x obtained after observing the observation vector y. Let x be the likelihood function, representing the probability of observing the observation vector y given the true state vector x. The prior probability distribution is the state vector x that is inferred from existing experience; This is the normalization constant; According to the above formula, the likelihood function Adopting a Gaussian distribution: in, Map the state vector x to the observation space, and use the SCIATRAN radiative transfer model in the above formula; For the covariance of the observed data error, is the inverse matrix of the error covariance, and is the weight matrix of the residuals; the larger the error, the smaller the weight. Based on the above formula, construct the likelihood function. Logarithmic form: in, Let i be the aerosol observation value at the i-th elevation angle. These are the model observations corresponding to the elevation angle. This is the observation error; Based on the Bayesian inversion framework described above, the constructed prior probability distribution adopts a multivariate Gaussian distribution: Among them, the prior mean This indicates that a typical aerosol extinction profile already exists. Let be the prior covariance matrix, representing the uncertainty and vertical correlation of the profile; Based on the above algorithm, calculate the logarithmic form of the posterior probability distribution of the aerosol extinction profile: in This is the weighted average of the residuals from aerosol observations. Based on the above algorithm, calculate the posterior mean of the posterior probability distribution of the aerosol extinction profile: in, is the posterior mean of the aerosol extinction profile, represents the optimal aerosol extinction profile, and N is the total number of sampling steps in the first layer of MCMC.
9. The aerosol and gas profile inversion method based on MAX-DOAS as described in claim 8, characterized in that: The posterior probability distribution of the gas concentration profile and the inversion process of the optimal trace gas concentration profile include: Calculate the posterior probability distribution of the gas concentration profile: in, The final gas profile distribution needs to be obtained through inversion. The posterior mean of the aerosol extinction profile obtained from the inversion step is given. To obtain the gas profile that needs to be inverted, To observe the concentration of trace gases in the oblique column, Let be the likelihood probability function for the gas. Let be the prior probability distribution of the gas; Based on the above algorithm, calculate the posterior mean of the posterior probability distribution of the trace gas concentration profile: in, Let be the posterior mean of the trace gas concentration profile, represent the optimal trace gas concentration profile, and M be the total number of sampling steps in the second layer MCMC.
10. A MAX-DOAS-based aerosol and gas profile inversion system, characterized in that: include: The posterior probability distribution inversion module for aerosol extinction profiles is used to input the original MAX-DOAS spectral data, obtain the O4 oblique column concentration and trace gas oblique column concentration through spectral fitting, construct a Bayesian inversion framework, quantify the degree of agreement between observations and the forward model through the likelihood function, invert the posterior probability distribution of aerosol extinction profiles using the MCMC sampling algorithm, and extract the posterior mean from the posterior probability distribution of aerosol extinction profiles as the optimal aerosol extinction profile. The posterior probability distribution inversion module of the gas concentration profile is used to fix the optimal aerosol extinction profile obtained by the above module, and combine the trace gas oblique column concentration to invert the posterior probability distribution of the gas concentration profile through the conditional likelihood function, the prior probability distribution of gas concentration and the MCMC sampling algorithm. The posterior mean is extracted from the posterior probability distribution of the trace gas concentration profile as the optimal trace gas concentration profile. The monitoring module uses the Bayesian error decomposition algorithm, combined with the posterior probability distribution of the obtained aerosol extinction profile and the obtained trace gas concentration profile, to achieve strict mathematical separation of measurement error, smoothing error and total error, monitor the inversion progress in real time and perform anomaly diagnosis.