WO2022194117A1

WO2022194117A1 - Statistical observation localized equivalent-weights particle filter-based data assimilation method

Info

Publication number: WO2022194117A1
Application number: PCT/CN2022/080805
Authority: WO
Inventors: 赵玉新; 杨硕; 邓雄; 赵廷; 郝日栩; 刘延龙; 赵恒德
Original assignee: 哈尔滨工程大学
Priority date: 2021-03-17
Filing date: 2022-03-15
Publication date: 2022-09-22
Also published as: CN113051529A; CN113051529B

Abstract

The present invention provides a statistical observation localized equivalent-weights particle filter-based data assimilation method, comprising: acquiring a mode integration initial background field; determining whether a statistical observation start moment is reached, and accumulating observations to calculate a statistical observation mean value; calculating a proposal density according to the statistical observation mean value to adjust ensemble particles; at a given assimilation moment, calculating weights of the particles by using an equivalent-weights method, and adjusting states of the particles; using a resampling method to adjust the ensemble particles to keep the number of the particles stable, and updating the states of the particles at positions corresponding to the observations; using a localization function to determine the weights of the particles surrounding a position corresponding to an assimilation observation; according to a localization weight, updating the weights of the particles, and updating the states of the surrounding particles; and calculating a state posteriori estimation value of the statistical observation localized equivalent-weights particle filter. The present invention can effectively improve the data assimilation quality of a non-Gaussian gridding mode, and can be better applied to real-time data assimilation in a complex gridding mode, thereby improving the quality of assimilation.

Description

A Data Assimilation Method Based on Statistical Observation Localized Average Weight Particle Filter

technical field

The invention relates to a data assimilation method based on statistical observation and localized equalization weight particle filtering, belonging to the field of atmospheric and ocean data assimilation.

Background technique

There are two ways to study ocean dynamics, one is to use numerical models to study, and the other is to directly observe the atmosphere and ocean. Numerical model simulation is mainly used to reflect the characteristics of sea areas, and direct observations such as satellites truly reflect the characteristics of ocean observations. Due to the massive acquisition of satellite remote sensing data, the particularity of the marine environment and the spatial distribution of ocean observations are insufficient. Data assimilation is a research method that can organically combine numerical models and observations, two basic means of oceanographic research. Data assimilation refers to a method of continuously incorporating new observational data in the dynamic operation of the dynamic model on the basis of considering the spatiotemporal distribution of the data and the errors of the observational field and background field. By continuously integrating new observation data into the model, the trajectory predicted by the model simulation can be gradually corrected to make it closer to the real trajectory, and the accuracy of the model simulation prediction can be improved. The main purpose of data assimilation is to combine observational data with theoretical model results, absorbing the advantages of both, in order to obtain results that are closer to reality. Forecasting of model states provides more accurate initial fields and optimizes ocean model parameters to improve the climate forecasting capabilities of ocean models.

As the core of data assimilation, data assimilation algorithm mainly relies on accurate observation data and reasonable numerical model. According to the correlation between the assimilation algorithm and the model, the data assimilation algorithm is divided into two categories: continuous data assimilation algorithm and sequential data assimilation algorithm. For example, the patent application of Patent Application No. 201910038258.3 uses the coupled data assimilation and parameter optimization method based on the optimal observation time window. This patent uses the coupled data assimilation and parameter optimization method based on the optimal observation time window. Technical field of parameter optimization and numerical prediction. Extract the effective observation information to the greatest extent to fit the characteristic variability of the coupled model state, ignore the time-varying characteristics of the internal parameters of the model, and introduce the time-averaged coefficients within the time window to achieve more accurate estimation and optimization of model parameters. Atmospheric and ocean numerical prediction capabilities of coupled models. For example, the patent application No. 202010013183.6 of the patent application satellite data assimilation in the vertical direction of the adaptive localization method and ensemble Kalman filter weather assimilation forecast method, the patent uses the satellite data assimilation in the vertical direction of the adaptive localization method and ensemble card Mann filter weather assimilation forecasting method. The adaptive localization method calculates the correlation coefficient between the observation data and the model variable according to any observation data and model variables given in the ensemble Kalman filter assimilation system; and then uses the grouped correlation coefficient to estimate the relationship between the observation data and the model variable. The original localization function; according to the profile of the correlation coefficient, the position of the satellite observation is estimated, and the obtained adaptive localization parameter is used to predict the typhoon in the regional model. The prediction result is compared with the prediction result without using the present invention, The error relative to the observation is obviously reduced, and at the same time, the use of the present invention also significantly improves the forecast of the typhoon in the stage of rapid intensification. For example, the patent application No. 201910430413.6 is a water quality model particle filter assimilation method based on multi-source observation data. The patent constructs a two-dimensional water quality model; initializes the state variables and parameters of the particles; generates the boundary conditions of the particles; resamples to obtain a new set of particles; The parameters are recursively pushed from time t to time t+1; at update time, the boundary conditions of particles continue to be generated, until the operation is completed at all times, and the particle filter assimilation of the two-dimensional water quality model is realized. The multi-source observation data of water quality is reasonably integrated into the two-dimensional water quality model by using the particle filter algorithm, and the parameters of the two-dimensional water quality model are dynamically updated, which improves the simulation accuracy and prediction ability of the two-dimensional water quality model.

Particle filter algorithm is an ensemble data assimilation method. Since particle filter algorithm is not constrained by the assumption of model state quantity and error Gaussian distribution, it is suitable for any nonlinear non-Gaussian dynamic system. It also uses the Monte Carlo sampling method to approximate the posterior probability density distribution of the state quantity, which can better represent the change information of the nonlinear system. The particle filter algorithm is simple and easy to implement. At the same time, compared with the current mainstream Kalman filter series algorithms, the particle filter algorithm does not have complex operations such as matrix transposition and inversion, so the calculation efficiency is higher. Compared with the Kalman filter series algorithms, which directly update the state value of particles, the particle filter algorithm only updates the weight of the particle when updating the particle, and the actual state value of the particle remains unchanged, which can avoid the update process. There is a situation where the particle state value is outside its physical value range.

For the improvement of particle filtering, the equal-weight particle filtering method proposed by British professor Van Leeuwen uses the proposed density idea to effectively improve the particle degradation and particle depletion problems in traditional particle filtering. It can achieve more traditional methods by using fewer aggregated particles. The assimilation effect of particles. At the same time, the weighted particle filter method based on statistical observation can effectively improve its dependence on future observation information, making it better applied to the field of real-time data assimilation and effectively improving the assimilation quality. However, the equal-weight particle filter method lacks a corresponding localization scheme, which makes it difficult to adapt to the complex gridded high-latitude patterns.

SUMMARY OF THE INVENTION

The purpose of the present invention is to provide a suitable localization scheme for the equal-weight particle filter, solve the limitation of the application of the method in the complex grid mode, and make the method have better potential and value of practical application. Provided is a method for data assimilation based on localized weighted particle filtering based on statistical observations.

The object of the present invention is achieved in this way:

Step 1: Obtain the initial background field of the mode integral

First, the initial field of the model is introduced into the model integral equation, and the model equation is integrated first to make the model equation reach a chaotic state. This method can avoid the occurrence of the fluctuation problem of the model equation. The model variable that reaches the chaotic state is used as the initial background field, and the background field As the starting point of the integration, the subsequent mode state integration is carried out

Step 2: Determine whether the start time of statistical observation is reached, and accumulate the observations to obtain the mean value of the statistical observation

According to the given τ value, the starting moment of the statistical observation calculation is determined, and selecting the appropriate τ can effectively improve the reliability of the statistical observation, and can better guide the ensemble particles to approach the historical observation at the time of assimilation. When the statistical observation starts, the historical observation at the corresponding position of the assimilation observation is accumulated to obtain the mean value. The calculation of the mean value can effectively avoid the sudden jump in the historical information of the corresponding observation, and the statistical mean value is used to replace the future observation information in the traditional method.

Step 3: Calculate the proposed density according to the mean value of the statistical observations. At a given assimilation time, the particle weight is calculated using the average weight method and the particle state is adjusted.

Calculate the proposed density of the set of particles corresponding to the observation position according to the historical mean of the observation position at the moment of statistical assimilation, select the optimal proposed density according to the proposed density of the particles in the set, and determine the posterior probability density of each particle in the set for the statistical observation, At the same time, the ensemble particles are adjusted to be close to the statistical historical observation mean, and the particle state is updated according to the proposed density.

Step 4: At a given assimilation time, use the average weight method to calculate the particle weight, and adjust the observed particle state at the corresponding position

Calculate the proposed density using the statistical historical observation of the corresponding position of the observation. According to the set particles adjusted by the proposed density, the adjusted particles are brought into the weight formula in the equal-weight particle filter, and the weight of the particles in the set is re-determined according to the formula to ensure that the set is in the set. According to the observation, the particles of , can obtain a relatively close and optimal weight, and further adjust the state of the aggregated particles.

Step 5: Use the resampling method to adjust the aggregated particles to keep the number of particles stable, and update the state of the particles at the corresponding position of the observation

The resampling method is used to adjust the aggregated particles at the corresponding position of the observation, mainly to ensure that the aggregated particles with poor weight performance in the equal-weight method are adjusted, and the particles with smaller weights are proposed to maintain the stability of the aggregated particle number.

Step 6: Use the localization function to determine the particle weights around the assimilated observation location

After updating the state of the set of particles at the corresponding position observed at the time of assimilation, the localization function is used to continue to adjust the state of the set of particles within the influence radius of the observation. For the selection of the localization function, refer to the localization scheme in the localized particle filter method, and use the localization function in the calculation process to describe the positional relationship between the weight of the aggregated particles in the area and the given observation. The assimilation observations of , refer to the localization parameters to determine the weights of the adjusted ensemble particles.

Step 7: Update the particle weight according to the localization weight, and update the surrounding particle state;

After the weight of the affected aggregate particles is determined according to the assimilation observation, the state of the corresponding aggregate particle needs to be adjusted by this weight. Before the state adjustment, the weight of the aggregate particle needs to be normalized. The linear combination of the prior particles updates the state of the aggregated particles after localization, so that the observed nearby aggregated particles can be adjusted according to the non-corresponding point observations.

Step 8: Calculate the state posterior estimate of the localized weighted particle filter for statistical observations

After using the localization scheme to adjust the observed state of the surrounding ensemble particles, the overall ensemble mean is reintroduced into the model equation, and the model integral is updated during the assimilation process to obtain the final assimilation analysis result.

The second step is specifically: judging whether the start time of the statistical observation is reached, and accumulating the observations to obtain the mean value of the statistical observation;

According to the given value of τ, the starting moment of the statistical observation calculation is determined. The calculation formula of τ is:

where t ⁰ is the previous observation time to be assimilated, t ⁿ is the observation interval between two observations, and t ^j represents the current time; in order to ensure the reliability of the statistical observation results, when the observation interval to be assimilated is large, τ is generally The selection is 0.8-0.9. When the assimilation interval is small, τ is generally selected as 0.5-0.8; the statistical method is used to calculate the average value of the observations before the assimilation time.

represents the time series of the observation information in the observation interval at the observation position to be assimilated. Assuming that ^tj is the start time of the statistical observation according to τ, the mean value of the statistical observation can be expressed as:

The third step is specifically: calculating the proposed density adjustment set particle according to the mean value of statistical observation, at a given assimilation moment, using the average weight method to calculate the particle weight, and adjust the particle state;

Step 3.1: Calculate the particle probability density at the corresponding position of the observation according to the mean value of the statistical observation;

In Bayesian theory, the solution of the posterior probability density needs to rely on the prior probability density distribution and the likelihood probability solution. The particle filter method uses the conditional posterior probability density distribution based on Bayesian theory.

The posterior probability density distribution of , can be expressed as

in

In order to propose the density, for the calculation of the proposed density, it is necessary to first calculate the probability density of the aggregated particles, using the mode equation, from the state of the particle at time n-1

Integrating, the probability density expression of the aggregated particles is:

represents the probability density from the particle state at time n-1 to time n,

represents the state of the ith particle at time n,

Represents the integral result of the mode equation at time n-1, and Q represents the mode error covariance. The expression is:

M represents the total number of aggregated particles,

Represents the mean value of aggregate particles; the probability density of particles is calculated by this formula;

Step 3.3: Obtain the particle proposal density from the probability density of the aggregated particles

Based on the prior density calculated in the previous step in the model equation, the proposed density in the uniformly weighted particle can be further obtained. The proposed density obtained from the prior density is expressed as:

formula

Represents the proposed density targeting observation y, K ⁿ represents the relaxation stress matrix whose expression is K ⁿ =QH ^T (HQH ^T +R) ^-1 , K ⁿ and the model error covariance Q and the observation error covariance R expression for:

Indicates the state of the ith particle at time n, H indicates that the observation operator generally takes a value of 1, and y indicates the observation information; in the particle filter based on the statistical observation average weight, the calculation of the proposed density is not limited to the observation information at the time of assimilation, The proposed density is based on the historical mean of statistical observations before the assimilation time, in order to adjust the ensemble particle position close to the observation information before the assimilation time;

Step 3.4: Select the optimal proposal density and calculate particle weights

Finding the optimal proposal density hypothesis

The weight of each particle in the set based on statistical observations is expressed as:

where ω _i represents the weight of the ith particle in the set,

represents the state of the ith particle at time t ^j ,

represents the statistical observation information at time t ^j ,

represents the probability density,

represents the proposal density,

Represents the probability density of the particle state at time t ^j for the statistical observation, and finally obtains the weight of each particle passing through the proposed density; for the proposed density to adjust the aggregated particles to approach the statistical observation, the aggregated particle state can be expressed as:

in

is the observation operator, and the prior relaxation coefficient is B(τ),

Make sure that the ensemble particles are close to the statistical observations, so the relaxation coefficient can be expressed as:

B(τ)=bτQH ^T R ^-1

where τ is the statistical observation start threshold, and b is the scaling factor that controls the degree of relaxation to observations;

Step 4 is specifically as follows: at a given assimilation moment, using the average weight method to calculate the particle weight, and adjust the observed particle state at the corresponding position;

Step 4.1: Calculate the weight of the particle on the basis of the proposed density-adjusted set of particles;

The proposed density brings each particle in the ensemble close to the observed information, so that it gets approximately the same weight; the weight expression for each particle is:

in

Represents the cumulative multiplication result of the proposed density weight in the assimilation time interval, and the minimum weight ω _i of the particle can be determined as:

in

Indicates the prior weight of the i-th particle, which is the weight of each particle in the proposed density calculation process; y denotes the observation vector; H denotes the observation projection operator, which is H=1 in the simple mode; x represents the state vector; the superscript T represents the matrix transposition; Q represents the mode error covariance matrix; R represents the observation error covariance matrix; assuming that the particle target weight C _i in the set is obtained, in order to ensure that 80% of the particles in the particle set can reach the Calculated weights to avoid particle degradation and depletion problems:

Thus, it is found that for the state of n time

Most of the particles in the set can maintain the same weight, and some particles that do not reach the average weight can be adjusted by resampling;

Step 4.2: Adjust particle state according to weight

After obtaining the aggregate particle weight, adjust the particle state according to the weight, then the state at time n can be expressed as:

In the formula, y represents the observation vector; H represents the observation projection operator (H=1); x represents the state vector; K=QH ^T (HQH ^T +R) ^-1 , Q is the mode error covariance, and R is the observation error covariance , when the weights of the set particles are approximately equal in the equal-weight particle filter, the vector α _i can be expressed as:

middle

and

in these two expressions

represents the selected target weight,

Represents the relative weight of the particle at the current moment; in order to ensure the random effect of the particle state in the set, a random term is added, and the final analysis equation is obtained:

Step 6 is specifically: using the localization function to determine the particle weight around the corresponding position of the assimilation observation;

After using the equal-weight particle filtering method and the resampling method to update the aggregate particle state at the corresponding position of the assimilation observation, use the localization function to continue to adjust the aggregate particle state within the influence radius r of the observation; for the selection of the localization function, refer to the localization function. The localization scheme in the particle filter method uses the localization function to describe the relationship between the weight of the aggregated particles in the area and the given observation position during the calculation process. The weight in the particle filter represents the likelihood probability of the observation. In the calculation Using the localization operator l[y, _xi ,r], the operator is used to describe the relative position coefficient of the observation y and the state of the set particle _xi , and the corresponding particle weight is expressed as:

ω _i,j =p(y|x _i,j )l[y,x _j ,r]+[1-l[y,x _j ,r]]

Among them, the localization operator l[y, _xi ,r] is mainly used in local analysis to judge the relative position information of the set particle and the observation information; when the observation y and the set particle x _i coincide, the maximum value of this function is 1 , when the distance between the two exceeds the given influence radius r, the function value is 0, indicating that the observation has no adjustment effect on the state of this set of particles; the vector weight of the final set of particles is expressed as:

in

is the observation error covariance, h _j is the measurement operator; it can be seen that the weight of the aggregated particles is not only related to the observation information, but also related to the relative position of the observation and the particle;

Step 7 is specifically: update the particle weight according to the localized weight, and update the surrounding particle state;

After obtaining the weight of the aggregated particles within the influence radius of the observation information at the time of assimilation, the state of the corresponding aggregated particle needs to be adjusted by this weight. Before the state adjustment, the weight of the aggregated particle needs to be normalized to ensure that the sum of the weights in the aggregate is 1 , the normalized weight formula can be expressed as:

According to the normalized weight, the particle state in the set is re-adjusted, and the posterior particle state can be expressed as:

in

Represents the posterior mean, k _n is the nth sampled particle, where vectors r ₁ and r ₂ can linearly combine new particles to form sampled particles and prior particles, and finally realize the posterior update of the localized particle state, where the vector The calculation formulas of r ₁ and r ₂ can be expressed as:

r _2,j =c _j r _1,j

formula

is the error variance of all observations up to y _i , in which the posterior correlation between ensemble particle states is ignored, but the corresponding correlation is provided through the corresponding sampling step, where for the localization operator when l When [x _j , y _j , r] → 1, c _j → 0, and we have

Because the posterior variance is approximately equal to the sampling particle variance when l[x _j , y _j , r]=1, the same can be obtained

This enables the posterior particle to obtain the state of the sampled particle, when l[x _j ,y _j ,r]→0, c _j →∞,

Because the posterior variance is equal to the prior variance when l[x _j , y _j , r]=0, the same can be obtained

This sampling method provides a way to tune the particles to fit a general Bayesian posterior solution in the vicinity of the observation. Because each sampled particle is combined with a prior particle, the posterior ensemble contains the model state unique to the ensemble particle, This avoids the collapse of ensemble variance during observation assimilation;

The eighth step is specifically: calculating the state a posteriori estimated value of the localized weighted particle filter for statistical observation;

According to the observation information and the corresponding position of the new observation and the state of all particles within the influence radius of the observation, calculate the posterior estimated set mean of the state:

Take the updated posterior estimated set mean as the initial value of the analysis model, and bring it back into the model integral equation for the next step of prediction and assimilation. Repeat the above steps within the assimilation time when there are available observations to obtain the final analysis field. A field can act as a data field that reflects the current state of the environment.

Compared with the prior art, the beneficial effects of the present invention are:

(1) The localization scheme is introduced into the statistical observation weighted particle filter, which can effectively improve the assimilation quality when the observation is sparse. The root mean square error of the assimilation result of the statistical observation localization weighted particle filter is better than the localized particle filter method. Localization improvement of traditional equal-weight particle filter;

(2) After the introduction of the localization scheme, the statistical observation weighted particle filter can be better applied to the complex gridded medium and high-dimensional patterns, so that the particle filter method has better practical application potential.

Description of drawings

Figure 1 shows the process of localized weighted particle filtering for statistical observation;

Figure 2 is a comparison chart of the root mean square error of the localized average weight particle filter for statistical observation;

Fig. 3 is the data assimilation flow chart of traditional equal weight particle filter;

Fig. 4 is a flow chart of data assimilation based on statistical observation average weight particle filter;

Figure 5 is a flow chart of data assimilation based on localized average weight particle filter based on statistical observations.

Detailed ways

The present invention will be described in further detail below with reference to the accompanying drawings and specific embodiments.

In order to more simply and clearly describe the specific implementation steps of the localized weighted particle filter assimilation method based on statistical observations, a simple Lorenz-96 model is taken as an example for a brief description, which can better reflect the effect of the localization method on the assimilation results. influences.

Step 1: Obtain the initial background field of the mode integral

The Lorenz-96 mode is selected because it has 40 state variables, which can better apply the localization scheme. In the Lorenz-96 mode, input the initial starting point of the mode into the model equation and integrate 1 million steps for spin-up to make the mode The variable reaches the chaotic state, to avoid the mode fluctuation deviation introduced in the assimilation integration process, and the mode variable that reaches the chaotic state is used as the initial background field of the mode. Based on this mode background field, the initial value of the particle filter set particles is obtained through the mode equation. , which can effectively improve the reliability of the ensemble particle state obtained by solving the mode equation.

In the data assimilation of statistical historical observation weighted particle filter, the statistical observation is used to replace the future observation information of the traditional method, and the aggregated particles are guided to be close to the historical statistical results observed at the time of assimilation, and the number of effective particles in the sampling process is increased, and an appropriate threshold τ is selected to judge the appropriate Start statistical observations. The choice of τ is related to the time interval of the observation information that needs to be assimilated. When τ is too small, more useless observations will be introduced into the statistics and the statistical error will be increased. When τ is too large, the statistical observation information will be too small. It is not possible to correctly calculate the proposed density to guide the particle closer to the observation. The calculation formula of τ is:

where t ⁰ is the previous observation time that needs to be assimilated, t ⁿ is the observation interval between two observations, and t ^j represents the current time. In order to ensure the reliability of the statistical observation results, when the observation interval to be assimilated is large, τ is generally selected as 0.8-0.9, and when the assimilation interval is small, τ is generally selected as 0.5-0.8. Use statistical methods to calculate the mean of the observations before the assimilation time, using

Step 3.1: Calculate the probability density of particles at the corresponding position of the observation according to the mean value of the statistical observation

The posterior probability density distribution of , can be expressed as

in

Integrating, the probability density expression of the aggregated particles is:

represents the state of the ith particle at time n,

M represents the total number of aggregated particles,

represents the aggregate particle mean. The probability density of particles is calculated by this formula.

formula

Indicates the state of the ith particle at time n, H indicates that the observation operator generally takes a value of 1, and y indicates the observation information. In the weighted particle filter based on statistical observation, the calculation of the proposed density is not limited to the observation information at the time of assimilation. The proposed density can use the historical mean value of the statistical observation before the time of assimilation. The purpose is to adjust the position of the set particles to be close to the observation information before the time of assimilation. .

Step 3.4: Select the optimal proposal density and calculate particle weights

The selection of the proposed density is considered to be an important criterion for controlling the position of particles and observation information and the calculation of particle weights. Selecting the appropriate and optimal proposed density can ensure the number of valid particles in sampling and the reliability of particle weights, which is also an equal-weight particle filter. The most important part in finding the optimal proposal density hypothesis

where ω _i represents the weight of the ith particle in the set,

represents the state of the ith particle at time t ^j ,

represents the statistical observation information at time t ^j ,

represents the probability density,

represents the proposal density,

Represents the probability density of the particle state at time ^tj for the statistical observation, and finally obtains the weight of each particle passing through the proposed density. For the proposed density-adjusted ensemble particle approaching statistical observations, the ensemble particle state can be expressed as:

in

is the observation operator, and the prior relaxation coefficient is B(τ),

B(τ)=bτQH ^T R ^-1

where τ is the statistical observation start threshold, and b is the scaling factor that controls the degree of relaxation to observations.

Step 4.1: Calculate the weights of particles based on the proposed density-adjusted aggregate particles

According to the observation information in the assimilation, the position of the particles in the proposed density adjustment set is calculated, and the particles are closer to the observation information at the assimilation time, so as to ensure that most of the particles get equal weights in the posterior probability density function at the assimilation time. The proposed density brings each particle in the ensemble close to the observational information, so that it gets approximately the same weight. The weight expression for each particle is:

in

Indicates the prior weight of the i-th particle, which is the weight of each particle in the proposed density calculation process; y denotes the observation vector; H denotes the observation projection operator, which is H=1 in the simple mode; x represents the state vector; the superscript T represents the matrix transpose; Q represents the mode error covariance matrix; R represents the observation error covariance matrix. Assuming that the particle target weight C _i in the set is obtained, in order to ensure that 80% of the particles in the particle set can reach the calculated weight and avoid the occurrence of particle degradation and depletion problems:

Thus, it can be found that for the state of n time

Most of the particles in the set can maintain the same weight, and some particles that do not reach the average weight can be adjusted by resampling.

Step 4.2: Adjust particle state according to weight

middle

and

in these two expressions

represents the selected target weight,

Represents the relative weight of the particle at the current moment. In order to ensure the random effect of the particle state in the set, a random term is added, and the final analytical equation is obtained:

The resampling method is used to remove the particles with smaller weights in the set of particles at the corresponding position observed in the equal-weight particle filter, while ensuring the stability of the total number of particles in the set, and copy the particles with larger weights. The equal-weight particle filtering method itself ensures that most of the particles are preserved, but in order to prevent a very small number of particles that do not meet the requirements of equal weight in the process of equal-weighting, the resampling method is used to perform final state adjustment to keep the number of particles stable.

Step 6: Use the localization function to determine the particle weight around the corresponding position of the assimilation observation

After using the equal-weight particle filtering method and the resampling method to update the collective particle state at the corresponding position of the assimilation observation, the localization function is used to continue to adjust the collective particle state within the influence radius r of the observation. For the selection of the localization function, refer to the localization scheme in the localized particle filter method. In the calculation process, the localization function is used to describe the relationship between the weight of the aggregated particles in the area and the given observation position, and the weight in the particle filter is used. Indicates the likelihood probability of the observation. The localization operator l[y,x _i ,r] is used in the calculation. This operator is used to describe the relative position coefficient of the observation y and the set particle state x _i , and the corresponding particle weight. Expressed as:

ω _i,j =p(y|x _i,j )l[y,x _j ,r]+[1-l[y,x _j ,r]]

Among them, the localization operator l[y, x _i , r] is mainly used in local analysis to judge the relative position information of aggregated particles and observation information. When the observation y coincides with the set particle _xi , the maximum value of this function is 1. When the distance between the two exceeds the given influence radius r, the function value is 0, which means that the observation has no adjustment effect on the state of the set particle. . The vector weight of the final set of particles is expressed as:

in

is the observation error covariance, and h _j is the measurement operator. It can be seen that the weight of the aggregated particles is not only related to the observation information, but also related to the relative position of the observation and the particle.

in

Represents the posterior mean, k _n is the nth sampled particle, where vectors r ₁ and r ₂ can linearly combine new particles to form sampled particles and prior particles, and finally achieve a posteriori update of the localized particle state. The calculation formula of the vectors r ₁ and r ₂ can be expressed as:

r _2,j =c _j r _1,j

formula

This allows the posterior particle to obtain the state of the sampled particle. When l[x _j ,y _j ,r]→0, c _j →∞,

This sampling method provides a way to tune the particles to fit a general Bayesian posterior solution in the vicinity of the observation. Because each sampled particle is combined with a prior particle, the posterior ensemble contains the model state unique to the ensemble particle, This avoids the collapse of ensemble variance during observation assimilation.

The present invention proposes a data assimilation technology based on localized average weight particle filtering based on statistical observation. Compared with the traditional statistical observation weighted particle filter data assimilation technology, the significant feature of the present invention is that a localization scheme adapted to its application is proposed for the statistical observation weighted particle filter method, which effectively solves the problem of the method for complex high-latitude grids. The limitation of the application of the assimilation mode is that at the time of assimilation, the statistical observation method is used to update the state of the aggregated particles at the corresponding position of the observation, and then according to the position information and localization parameters of the assimilation observation, the observation is used to adjust the surrounding of the observation position, and the aggregated particles are available for observation sensitive. According to the weight of the observation, the state of the nearby set of particles is adjusted to improve the utilization rate of the observation and improve the quality of assimilation. The method proposed in this patent can effectively improve the assimilation ability of the traditional statistical observation weighted particle filter method under complex grid patterns and sparse observation conditions, can effectively improve the assimilation quality, and at the same time improve the application prospect of the statistical average weighted particle filter method.

Claims

A method for data assimilation based on localized weighted particle filtering based on statistical observation, characterized in that it comprises the following steps:

Step 1: Obtain the initial background field of the mode integral;

First, the initial field of the mode is introduced into the mode integral equation, and the mode equation is integrated first, so that the mode equation reaches the chaotic state, and the mode variable that reaches the chaotic state is used as the initial background field, and the background field is used as the starting point of integration to carry out the subsequent mode state integration;

Step 2: Judging whether the start time of statistical observation is reached, and accumulating the observation to obtain the mean value of statistical observation;

According to the given value of τ, the starting moment of the statistical observation calculation is determined. The calculation formula of τ is:

where t 0 is the previous observation time to be assimilated, t n is the observation interval between two observations, and t j represents the current time; in order to ensure the reliability of the statistical observation results, when the observation interval to be assimilated is large, τ is generally The selection is 0.8-0.9. When the assimilation interval is small, τ is generally selected as 0.5-0.8; the statistical method is used to calculate the average value of the observations before the assimilation time.
represents the time series of the observation information in the observation interval at the observation position to be assimilated. Assuming that tj is the start time of the statistical observation according to τ, the mean value of the statistical observation can be expressed as:

Step 3: Calculate the proposed density adjustment set particle according to the mean value of statistical observation. At a given assimilation time, use the average weight method to calculate the particle weight and adjust the particle state;

Step 3.1: Calculate the particle probability density at the corresponding position of the observation according to the mean value of the statistical observation;

In Bayesian theory, the solution of the posterior probability density needs to rely on the prior probability density distribution and the likelihood probability solution. The particle filter method uses the conditional posterior probability density distribution based on Bayesian theory.
The posterior probability density distribution of , can be expressed as

in
In order to propose the density, for the calculation of the proposed density, it is necessary to first calculate the probability density of the aggregated particles, using the mode equation, from the state of the particle at time n-1
Integrating, the probability density expression of the aggregated particles is:

represents the probability density from the particle state at time n-1 to time n,
represents the state of the ith particle at time n,
Represents the integral result of the mode equation at time n-1, and Q represents the mode error covariance. The expression is:

M represents the total number of aggregated particles,
Represents the mean value of aggregate particles; the probability density of particles is calculated by this formula;

Step 3.3: Obtain the particle proposal density from the probability density of the aggregated particles

Based on the prior density calculated in the previous step in the model equation, the proposed density in the uniformly weighted particle can be further obtained. The proposed density obtained from the prior density is expressed as:

formula
Represents the proposed density targeting observation y, K n represents the relaxation stress matrix whose expression is K n =QH T (HQH T +R) -1 , K n and the model error covariance Q and the observation error covariance R expression for:

Indicates the state of the ith particle at time n, H indicates that the observation operator generally takes a value of 1, and y indicates the observation information; in the particle filter based on the statistical observation average weight, the calculation of the proposed density is not limited to the observation information at the time of assimilation, The proposed density is based on the historical mean of statistical observations before the assimilation time, in order to adjust the ensemble particle position close to the observation information before the assimilation time;

Step 3.4: Select the optimal proposal density and calculate particle weights

Finding the optimal proposal density hypothesis
The weight of each particle in the set based on statistical observations is expressed as:

where ω i represents the weight of the ith particle in the set,
represents the state of the ith particle at time t j ,
represents the statistical observation information at time t j ,
represents the probability density,
represents the proposal density,
Represents the probability density of the particle state at time t j for the statistical observation, and finally obtains the weight of each particle passing through the proposed density; for the proposed density to adjust the aggregated particles to approach the statistical observation, the aggregated particle state can be expressed as:

in
is the observation operator, and the prior relaxation coefficient is B(τ),
Make sure that the ensemble particles are close to the statistical observations, so the relaxation coefficient can be expressed as:

B(τ)=bτQH T R -1

where τ is the statistical observation start threshold, and b is the scaling factor that controls the degree of relaxation to observations;

Step 4: At a given assimilation moment, use the average weight method to calculate the particle weight, and adjust the observed particle state at the corresponding position;

Step 4.1: Calculate the weight of the particle on the basis of the proposed density-adjusted set of particles;

The proposed density brings each particle in the ensemble close to the observed information, so that it gets approximately the same weight; the weight expression for each particle is:

in
Represents the cumulative multiplication result of the proposed density weight in the assimilation time interval, and the minimum weight ω i of the particle can be determined as:

in
Indicates the prior weight of the i-th particle, which is the weight of each particle in the proposed density calculation process; y denotes the observation vector; H denotes the observation projection operator, which is H=1 in the simple mode; x represents the state vector; the superscript T represents the matrix transposition; Q represents the mode error covariance matrix; R represents the observation error covariance matrix; assuming that the particle target weight C i in the set is obtained, in order to ensure that 80% of the particles in the particle set can reach the Calculated weights to avoid particle degradation and depletion problems:

Thus, it is found that for the state of n time
Most of the particles in the set can maintain the same weight, and some particles that do not reach the average weight can be adjusted by resampling;

Step 4.2: Adjust particle state according to weight

After obtaining the aggregate particle weight, adjust the particle state according to the weight, then the state at time n can be expressed as:

In the formula, y represents the observation vector; H represents the observation projection operator (H=1); x represents the state vector; K=QH T (HQH T +R) -1 , Q is the mode error covariance, and R is the observation error covariance , when the weights of the set particles are approximately equal in the equal-weight particle filter, the vector α i can be expressed as:

middle
and
in these two expressions
represents the selected target weight,
Represents the relative weight of the particle at the current moment; in order to ensure the random effect of the particle state in the set, a random term is added, and the final analysis equation is obtained:

Step 5: Use the resampling method to adjust the aggregated particles to keep the number of particles stable, and update the state of the particles corresponding to the observed position;

Step 6: Use the localization function to determine the particle weight around the corresponding position of the assimilation observation;

After using the equal-weight particle filtering method and the resampling method to update the aggregate particle state at the corresponding position of the assimilation observation, use the localization function to continue to adjust the aggregate particle state within the influence radius r of the observation; for the selection of the localization function, refer to the localization function. The localization scheme in the particle filter method uses the localization function to describe the relationship between the weight of the aggregated particles in the area and the given observation position during the calculation process. The weight in the particle filter represents the likelihood probability of the observation. In the calculation Using the localization operator l[y, xi ,r], the operator is used to describe the relative position coefficient of the observation y and the state of the set particle xi , and the corresponding particle weight is expressed as:

ω i,j =p(y|x i,j )l[y,x j ,r]+[1-l[y,x j ,r]]

Among them, the localization operator l[y, xi ,r] is mainly used in local analysis to judge the relative position information of the set particle and the observation information; when the observation y and the set particle x i coincide, the maximum value of this function is 1 , when the distance between the two exceeds the given influence radius r, the function value is 0, indicating that the observation has no adjustment effect on the state of this set of particles; the vector weight of the final set of particles is expressed as:

in
is the observation error covariance, h j is the measurement operator; it can be seen that the weight of the aggregated particles is not only related to the observation information, but also related to the relative position of the observation and the particle;

Step 7: Update the particle weight according to the localization weight, and update the surrounding particle state;

After obtaining the weight of the aggregated particles within the influence radius of the observation information at the time of assimilation, the state of the corresponding aggregated particle needs to be adjusted by this weight. Before the state adjustment, the weight of the aggregated particle needs to be normalized to ensure that the sum of the weights in the aggregate is 1 , the normalized weight formula can be expressed as:

According to the normalized weight, the particle state in the set is re-adjusted, and the posterior particle state can be expressed as:

in
Represents the posterior mean, k n is the nth sampled particle, where vectors r 1 and r 2 can linearly combine new particles to form sampled particles and prior particles, and finally realize the posterior update of the localized particle state, where the vector The calculation formulas of r 1 and r 2 can be expressed as:

r 2,j =c j r 1,j

formula
is the error variance of all observations up to y i , in which the posterior correlation between ensemble particle states is ignored, but the corresponding correlation is provided through the corresponding sampling step, where for the localization operator when l When [x j , y j , r] → 1, c j → 0, and we have
Because the posterior variance is approximately equal to the sampling particle variance when l[x j , y j , r]=1, the same can be obtained
This enables the posterior particle to obtain the state of the sampled particle, when l[x j ,y j ,r]→0, c j →∞,
Because the posterior variance is equal to the prior variance when l[x j , y j , r]=0, the same can be obtained
This sampling method provides a way to tune the particles to fit a general Bayesian posterior solution in the vicinity of the observation. Because each sampled particle is combined with a prior particle, the posterior ensemble contains the model state unique to the ensemble particle, This avoids the collapse of ensemble variance during observation assimilation;

Step 8: Calculate the state a posteriori estimate of the localized weighted particle filter for statistical observation;

According to the observation information and the corresponding position of the new observation and the state of all particles within the influence radius of the observation, calculate the posterior estimated set mean of the state:
The updated posterior estimated set mean is used as the initial value of the analysis model, and is re-introduced into the model integral equation for the next step of prediction and assimilation. Repeat the above steps within the assimilation time when there are available observations to obtain the final analysis field. A field can act as a data field that reflects the current state of the environment.