GB2605726A

GB2605726A - No title - not published

Info

Publication number: GB2605726A
Application number: GB2209053.4A
Authority: GB
Inventors: Zhao Yuxin; Yang Shuo; Deng Xiong; Zhao Ting; Hao Rixu; Liu Yanlong; Zhao Hengde
Original assignee: Harbin Engineering University
Current assignee: Harbin Engineering University
Priority date: 2021-03-17
Filing date: 2022-03-15
Publication date: 2022-10-12
Also published as: GB2605726A8; GB202209053D0

Abstract

The disclosure provides a method for assimilating data based on a localized equivalent weights particle filter with statistical observations, including: obtaining a model integral initial background field; determining whether a statistical observation start time is reached, and accumulating observations to calculate an average value of statistical observations; calculating proposal densities according to the average value of the statistical observations to adjust ensemble particles; at an assimilation time, calculating weights of the particles by the equivalent weight method to adjust states of the particles; adjusting the particles by a resampling method to keep the number of the ensemble particles stable, and updating the states of the particles at observation positions; determining weights of the particles around assimilation observation corresponding positions by a localization function; updating weights of the particles according to localization weights, and updating states of surrounding particles; and calculating posterior estimation values of the states in the localized equivalent weights particle filter with statistical observations. The disclosure can effectively improve the data assimilation quality of a non-Gaussian grid model, and can be better applied to real-time data assimilation in a grid complex model, thereby improving the assimilation quality. The method may be applied in the field of atmospheric and oceanographic data assimilation.

Description

Method for Assimilating Data Based on Localized Equivalent Weights Particle Filter with Statistical Observations

Technical Field

[0001] The disclosure relates to a method for assimilating data based on a localized equivalent weights particle filter with statistical observations, and belongs to the field of atmospheric and oceanographic data assimilation.

Background

[0002] There are two ways to study ocean dynamics: one is to use numerical models, and the other is to directly observe the atmosphere and ocean. Numerical model simulation is mainly used for reflecting the characteristics of sea areas, and direct observations by satellites, etc. truly reflect the characteristics of ocean observations. Due to the massive acquisition of satellite remote sensing data, the ocean environment has particularity, and the ocean observations have defects in spatial distribution. Data assimilation is a research method that can organically combine two basic oceanographic research methods, namely numerical models and observations. Data assimilation refers to a method of continuously incorporating new observation during the dynamic operation of a dynamic model on the basis of considering the spatiotemporal distribution of data and the errors of an observation field and a background field. By continuously incorporating new observations into the model, the trajectory predicted by the model simulation can be gradually corrected to make the trajectory closer to a real trajectory, so that the prediction accuracy of the model simulation can be improved. The main purpose of data assimilation is to combine observations with theoretical model results to absorb the advantages of both, in order to obtain results closer to actual results. At present, data assimilation methods have been widely used in the fields of atmosphere, ocean and land surface to provide more accurate initial fields and optimize ocean model parameters for the prediction of ocean model states, thereby improving the climate prediction ability of ocean models.

[0003] Data assimilation algorithms, as a core of data assimilation, mainly rely on accurate observation data and reasonable numerical models. According to the correlations between assimilation algorithms and models, the data assimilation algorithms are divided into two categories: continuous data assimilation algorithms and sequential data assimilation algorithms. For example, the patent application with a patent application number of 201910038258.3 provides a coupled data assimilation and parameter optimization method based on an optimal observation time window. The patent uses a coupled data assimilation and parameter optimization method based on an optimal observation time window, and belongs to the technical fields of data assimilation, parameter optimization and numerical prediction of coupled climate model systems. Effective observation information is extracted to the maximum extent to fit the characteristic variability of a coupled model state, the time-varying characteristics of internal parameters of the model are ignored, and time average coefficients within the time window are introduced, thereby realizing more accurate estimation and optimization of model parameters, and strengthening the numerical prediction ability of the coupled model for the atmosphere and ocean. For example, the patent application with a patent application number of 202010013183.6 provides an adaptive localization method of satellite data assimilation in a vertical direction and an ensemble Kalman filter weather assimilation prediction method. The patent uses an adaptive localization method of satellite data assimilation in a vertical direction and an ensemble Kalman filter weather assimilation prediction method. In the adaptive localization method, according to any observation data and model variables given in an ensemble Kalman filter assimilation system, correlation coefficients between the observation data and the model variables are calculated. Then, the grouped correlation coefficients are used for estimating an original localization function of the observation data and the model variables. According to the profile of the correlation coefficients, the positions of satellite observations are estimated. The obtained adaptive localization parameters are used for predicting typhoons in a region model. Compared with a prediction result without using the disclosure, the error of the prediction result relative to observations is significantly reduced. At the same time, the use of the disclosure also significantly improves the prediction of the rapid intensification stage of typhoons. For example, the patent application with a patent application number of 201910430413.6 provides a water quality model particle filter assimilation method based on multi-source observation data. The patent includes the steps of constructing a two-dimensional water quality model; initializing the state variables and parameters of particles; generating boundary conditions of the particles; performing resampling to obtain new particles; calculating optimal estimation values of the simulated state variables and parameters of the two-dimensional water quality model; recursing the parameters of the particles from a time t to a time t+1; and updating the time, and continuing to generate boundary conditions of the particles until the operations at all times are completed, thereby realizing the particle filter assimilation of the two-dimensional water quality model. A particle filter algorithm is used for reasonably integrating the water quality multi-source observation data into the two-dimensional water quality model to dynamically update the parameters of the two-dimensional water quality model, thereby improving the simulation accuracy and prediction ability of the two-dimensional water quality model. [0004] The particle filter algorithm is an ensemble data assimilation method. Since the particle filter algorithm is not constrained by the assumption of model state quantity and error Gaussian distribution, the particle filter algorithm is suitable for any non-linear non-Gaussian dynamic system. A Monte Carlo sampling method is also used for approximating the posterior probability density distribution of the state quantity, which can better represent the change information of a non-linear system. The particle filter algorithm is simple and easy to implement. Furthermore, 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 that directly update the state values of the particles, the particle filter algorithm only updates the weights of the particles when updating the particles, and the state values actually represented by the particles remain unchanged, thereby avoiding the situation that the state values of the particles exceed its physical value range during the updating process. [0005] For the improvement of the particle filter, the equivalent weights particle filter method proposed by British professor VanLeeuwen uses the proposal density idea to effectively ameliorate the particle degradation and depletion problems in the traditional particle filter, and uses fewer ensemble particles to achieve the assimilation effect of traditional methods using more particles. Furthermore, the equivalent weights particle filter method based on statistical observations effectively improves the dependence on the future observation information, is better applied to the field of real-time data assimilation, and effectively improves the assimilation quality. However, the equivalent weights particle filter method lacks a corresponding localization scheme, so that the method is difficult to adapt to complex grid high-latitude models.

Summary

[0006] In order to provide an appropriate localization scheme for an equivalent weights particle filter, solve the limitations of applications of the method in complex grid modes, and make the method have better potential and value in practical applications, the disclosure provides a method for assimilating data based on a localized equivalent weights particle filter with statistical observations.

[0007] The objectives of the disclosure are realized by the method including the following steps: [0008] step 1: obtaining a model integral initial background field, characterized by: [0009] first introducing a model initial field into a model integral equation, integrating the model equation first to enable the model equation to reach a chaotic state to avoid the occurrence of a fluctuation problem of the model equation, and performing subsequent model state integration by taking a model variable reaching the chaotic state as an initial background field and taking the background field as an integration starting point; [0010] step 2: determining whether a statistical observation start time is reached, and accumulating observations to calculate an average value of statistical observations, characterized in that: [0011] a start time of statistical observation calculations is determined according to a given r value, and selecting an appropriate 7-can effectively improve the reliability of statistical observations, and can better guide ensemble particles to be close to historical observations at an assimilation time; when the statistical observations start, the historical observations at assimilation observation corresponding positions are accumulated to calculate an average value, calculating the average value can effectively avoid sudden jumps in the historical information of the corresponding observations, and the future observation information in traditional methods is replaced with the statistical average value; [0012] step 3: calculating proposal densities according to the average value of the statistical observations to adjust ensemble particles, and at a given assimilation time, calculating weights of the particles by an equivalent weights method to adjust states of the particles, characterized by [0013] calculating proposal densities of the ensemble particles at observation corresponding positions according to the average value of the historical observations at observation positions at a statistical assimilation time, selecting an optimal proposal density according to the proposal densities of the particles in the ensemble, determining a posterior probability density of each of the particles in the ensemble for statistical observations, adjusting the ensemble particles to be close to the average value of statistical historical observations, and furthermore, updating states of the particles according to the proposal density; [0014] step 4: at the given assimilation time, calculating weights of the particles by an equivalent weights method, and adjusting states of the particles at observation corresponding positions, characterized by: [0015] calculating proposal densities by statistical historical observations at observation corresponding positions, substituting the ensemble particles adjusted according to the proposal densities into a weight formula in the equivalent weights particle filter, and redetermining the weights of the particles in the ensemble according to the formula to ensure that the particles in the ensemble can obtain relatively close and optimal weights according to observations, so as to further adjust the states of the ensemble particles; [0016] step 5: adjusting the ensemble particles by a resampling method to keep the number of the particles stable, and updating the states of the particles at observation corresponding positions, characterized in that: [0017] the ensemble particles at observation corresponding positions are adjusted by a resampling method to mainly ensure that the ensemble particles with poor weight performance in the equivalent weights method are adjusted, and the particles with smaller weights are proposed to maintain the stability of the number of the ensemble particles; [0018] step 6: determining weights of the particles around assimilation observation positions by a localization function, characterized in that: [0019] after the states of the ensemble particles at observation corresponding positions at an assimilation time are updated, the localization function is used to continue to adjust the states of the ensemble particles within an observation influence radius; the selection of the localization function refers to localization schemes in localized particle filter methods; in the calculation process, the localization function is used for describing the position relationship between the weights of the ensemble particles in a region and the given observation, and the weights of the adjusted ensemble particles are determined with reference to localization parameters according to the existing assimilation observations; [0020] step 7: updating weights of the particles according to localization weights, and updating states of surrounding particles, characterized in that: [0021] after the weights of the influenced ensemble particles are determined according to assimilation observations, the states of the corresponding ensemble particles need to be adjusted according to the weights; before state adjustment, the weights of the ensemble particles need to be normalized; according to the normalized weight and the linear combination of sampled particles and prior particles, the states of the ensemble particles after localization are updated, so that the states of the ensemble particles in the vicinity of observations can be adjusted according to non-corresponding point observations; and [0022] step 8: calculating posterior estimation values of states in the localized equivalent weights particle filter with statistical observations, characterized in that: [0023] after the states of the ensemble particles around observations are adjusted by a localization scheme, all ensemble average values are reintroduced into the model equation, and the model integral is updated during the assimilation process to obtain final assimilation analysis results.

[0024] Step 2 specifically includes the following steps: determining whether a statistical observation start time is reached, and accumulating observations to calculate an average value of statistical observations, characterized by: [0025] determining a start time of statistical observation calculations according to a given value, where a calculation formula of I-is: [0026] [0027] where t ° represents a previous observation time that needs assimilation, E" represents an observation interval between two observations, and t represents a current time; in order to ensure the reliability of statistical observation results, when the observation interval that needs assimilation is larger, r is generally selected to be 0.8-0.9, and when the assimilation interval is smaller, T is generally selected to be 0.5-0.8; a statistical method is used for counting the observations before the start of an assimilation time to calculate an average value, (to is used for representing a time sequence of observation information in the observation interval at an observation position that needs assimilation, and assuming that t j is determined to be a statistical observation start time according to, the average value of statistical observations can be expressed as: E yi Ytn (t" [0028] [0029] Step 3 specifically includes the following steps: calculating proposal densities according to the average value of the statistical observations to adjust ensemble particles, and at a given assimilation time, calculating weights of the particles by an equivalent weights method to adjust states of the particles, characterized by: [0030] step 3.1: calculating probability densities of the particles at observation corresponding positions according to the average value of the statistical observations, characterized in that: [0031] in a Bayesian theory, the solution of a posterior probability density needs to rely on a prior probability density distribution and a likelihood probability solution, a particle filter method uses a conditional posterior probability density distribution based on the Bayesian n theory, and the posterior probability density distribution for the statistical observation Y can be expressed as: t 1 Pcli Iriz) ii(xi) p (x) -z_s i=1 (y" ) q.\-7 1 vn) [0032] [0033] where q( v) represents a proposal density; for the calculation of proposal densities, probability densities of the ensemble particles need to be calculated first and a model equation is used for integrating from a particle state xi at a time n-1; an expression for the probability densities of the ensemble particles is: [0034] P(x7lx:1-1) exp 2 y cri( [0035] where 9( V') represents a probability density from a particle state at a time n-1 to a time n, represents a state of the i th particle at the time n,' ( ) represents an integral result of the model equation at the time n-1, and represents a model error covariance expressed as: Af XQ=COV{X,,7}=-I( [0036] Al [0037] where Al represents a total number of the ensemble particles, and 7 represents an average value of the ensemble particles; and the probability densities of the particles are calculated according to the formula; [0038] step 3.2: calculating proposal densities of the particles according to the probability densities of the ensemble particles, characterized in that: [0039] based on the prior density calculated in the previous step in the model equation, proposal densities in equivalent weights particles can be further obtained, and an expression for the proposal densities calculated according to the prior density is: [0040] (1(.71 ocexp[--1 (x" ( -(x7 ')-le(y-V '(xn [0041] where q (x7 kr represents a proposal density taking an observation Y as a target, K" = OH T (11011 + I?) 1, and an Kil represents a relaxation stress matrix expressed as expression of K", a model error covariance (-2 and an observation error covariance R is: Al = COV {X, , =-1(x -)2 [0042] M [0043] where represents a state of the i th particle at the time n, H represents an observation operator which is generally 1, and Y represents observation information; in an equivalent weights particle filter based on statistical observations, the calculation of the proposal densities is not limited to the observation information at the assimilation time, and the proposal densities are based on the historical average value of statistical observations before the assimilation time, in order to adjust the positions of the ensemble particles to be close to the observation information before the assimilation time; and [0044] step 3.3: selecting an optimal proposal density, and calculating weights of the particles, characterized in that: [0045] to find the optimal proposal density, assuming that q) ), the weight of each of the particles in the ensemble based on statistical observations is expressed as: = p(y kJ) ) if 1,, q x, ,y [0046] [0047] where cp, represents a weight of the th particle in the ensemble, represents a state of the i th particle at a time, , -v represents statistical observation information at the time p(xl q i-1,y -1 represents a proposal density, Ix( -1 x: t represents a probability density, '1/4-p(yR lx;} represents a probability density of the particle state at the time for statistical observations, and finally, the weight of each of the particles through the proposal density is obtained; for adjusting the ensemble particles to be close to statistical observations through the proposal densities, the states of the ensemble particles can be expressed as: [0048] (x) +oc), (xt h(xt [0049] where represents an observation operator, a prior relaxation coefficient is 130, and Ber)(Y ensures that the ensemble particles are close to statistical observations, so the relaxation coefficient can be expressed as: [0050] B(r)= In-QHT R1 [0051] where r represents a statistical observation start threshold, and h represents a scale factor that controls a degree of relaxation to observations.

[0052] Step 4 specifically includes the following steps: at the given assimilation time, calculating weights of the particles by an equivalent weights method, and adjusting states of the particles at observation corresponding positions, including the following steps: [0053] step 4.1: calculating weights of the particles based on the ensemble particles adjusted according to proposal densities, characterized in that: [0054] the proposal densities enable each of the particles in the ensemble to be close to the observation information, so that each of the particles obtains nearly the same weight; an expression of the weight of each of the particles is: (y" x; )P(x:' (x" kir% [0055] [0056] where represents a cumulative multiplication result of proposal density weights in an assimilation time interval, and a minimum weight c' of the particles can be determined and is expressed as: [0057] -log o, = -log co, +1(y" -HI ()CM"T (HQHT R kY -(x7'))

I

[0058] where represents a prior weight of the th particle, and the prior weight is the weight of each of the particles in a proposal density calculation process; Y represents an observation vector; H represents an observation projection operator, and the projection operator is H =1 in a simple model; -r represents a state vector; a superscript T represents a matrix transposition; (2 represents a model error covariance matrix; R represents an observation error covariance matrix; assuming that a target weight C of the particles in the ensemble is obtained, in order to ensure that 80% of the particles in the particle ensemble can reach the calculated weight to avoid the occurrence of particle degradation and depletion problems: [0059] ' -log (971 = -log cor"' + (yn - ))1 (11(2117 R)1 (t - x:")) [0060] therefore, it is found that most of the particles in the ensemble in a state at the time n can maintain the same weight and some particles that do not reach the equivalent weights can be adjusted by a resampling method; and [0061] step 4.2: adjusting states of the particles according to weights, characterized in that: [0062] after weights of the ensemble particles are obtained, states of the particles are adjusted according to the weights, and then, the state at the time n can be expressed as: = f (x1+ a,KV -#f (x71) [0063] [0064] where Y represents an observation vector; H represents an observation projection K = OHT (HQHT +R)1, Q represents a model operator (H =1); x represents a state vector; error covariance, and R represents an observation error covariance; when the weights of the ensemble particles are approximately equal in the equivalent weights particle filter, the vector can be expressed as: [0065] a 1 blar [0066] where in two expressions a1=0.5x.TR1HKx and b = 0 5A-TR -log tio71 -log coTt x = y" -Iff represents a selected target weight, and represents a relative weight of the particles at the current time; in order to ensure a random effect of the states of the particles in the ensemble, random terms are added, and a final analytical equation is obtained: + a,K (y" -HI (x,"-1)) [0067] [0068] Step 6 specifically includes the following steps: determining weights of the particles around assimilation observation corresponding positions by a localization function, characterized in that: [0069] after the states of the ensemble particles at the assimilation observation corresponding positions are updated by the equivalent weights particle filter method and the resampling method, the localization function is used to continue to adjust the states of the ensemble particles within an observation influence radius r; the selection of the localization function refers to localization schemes in localized particle filter methods; in the calculation process, the localization function is used for describing the relationship between the weights of the ensemble particles in a region and a given observation position, and the weight in the particle filter represents a likelihood probability of observations; in the calculation, a localization operator IlY'x°141 is used, the operator is used for describing a relative position coefficient of the observation Y and the state of the ensemble particles, and a corresponding particle weight is expressed as: (Di4+ )1EY' [0070] '-= y, [0071] where the localization operatori[y,x,,t']is mainly used in local analysis to determine relative position information of the ensemble particles and observation information; when the observation -11 coincides with the X1of the ensemble particles, the maximum value 1 of the function is taken, and when the distance between the two exceeds the given influence radius T, the function value is 0, which means that the observation has no adjusting effect on the state of the ensemble particles; finally, vector weights of the ensemble particles are expressed as: --I/1117 =nexp f=1 rliy,-hi(x)r ',TY-111(0r la co, = exp [0072] [0073] where represents an observation error covariance, and h' represents a measurement operator; and therefore, it can be seen that the weights of the ensemble particles are not only related to the observation information, but also related to relative positions of observations and particles.

[0074] Step 7 specifically includes the following steps: updating weights of the particles according to localization weights, and updating states of surrounding particles, characterized in that: [0075] after the weights of the ensemble particles within the influence radius of the observation information at the assimilation time are obtained, the states of the corresponding ensemble particles need to be adjusted according to the weights; before state adjustment, the weights of the ensemble particles need to be normalized to ensure that the sum of the weights in the ensemble is 1; a normalized weight formula can be expressed as: [0076] V = 0 1 1,1 1=1 [0077] the states of the particles in the ensemble are readjusted according to the normalized weight, and the states of posterior particles can be expressed as: [0078] 0 -1) -.7x(Yo) [0079] where represents a posterior average value, kn represents the nth sampled particle, vectors rI and r2 can enable new particles to form sampled particles to be linearly combined with prior particles to finally realize posterior update of the states of localized particles, and calculation formulas of the vectors and 1.2 can be expressed as: o- [0080] L[lr ü) Ne -1, x, 0.) -L v, +c (x" -x, [0081] r, =c r 1\7,(I -1[x3,yi,r]) [0082] 1[x1, y,, ) ; in this process, the [0083] where o-I represents error variances of all observations up to posterior correlations between the states of the ensemble particles are ignored, but corresponding correlations are provided through corresponding sampling steps; for the lirn y y 1 r I *l[x] -> c -> 1, localization operator, when, and since when a posterior variance is approximately equal to a sampled particle variance, and in a similar way, urn,; , =0 can be obtained, which enables the posterior particles to obtain the states of the /[x -> 0 lim =1 1[x, y cc sampled particles; when J'y i, and, ; since when ii °, a =0 posterior variance is equal to a prior variance, and in a similar way, i" can be obtained; this sampling method provides a way to adjust the particles to fit a general Bayesian posterior solution in the vicinity of observations; because each of the sampled particles is combined with a prior particle, a posterior ensemble includes model states unique to ensemble particles, which avoids the collapse of an ensemble variance during observation assimilation. ;[0084] Step 8 specifically includes the following steps: calculating posterior estimation values of states in the localized equivalent weights particle filter with statistical observations, characterized by [0085] updating the states of all particles at the observation corresponding positions and within the observation influence radius according to the observation information, calculating a I Al = x ;M ;posterior estimation ensemble average value of the states: , substituting the updated posterior estimation ensemble average value as an initial value of an analytical model into a model integral equation again for the next prediction and assimilation, and repeating the above steps within the assimilation time when there are available observations to obtain a final analytical field which can be used as a data field for reflecting a current environment state. [0086] Compared with the prior art, the disclosure has the following beneficial effects: [0087] (1) By introducing a localization scheme into the equivalent weights particle filter of statistical observations, the assimilation quality can be effectively improved during sparse observations. The root mean square error of the assimilation result of the localized equivalent weights particle filter with statistical observations is better than the localization improvement of the localized particle filter method and traditional equivalent weights particle filter. [0088] (2) After the localization scheme is introduced, the equivalent weights particle filter of statistical observations can be better used in complex grids and high-dimensional modes, and the particle filter method has better potential in practical applications. ;Brief Description of Figures ;[0089] FIG. 1 shows a process of a localized equivalent weights particle filter with statistical observations. ;[0090] FIG. 2 shows comparison charts of root mean square errors of the localized equivalent weights particle filter with statistical observations. ;[0091] FIG. 3 shows a flow chart of data assimilation of a traditional equivalent weights particle filter. ;[0092] FIG. 4 shows a flow chart of data assimilation based on an equivalent weights particle filter of statistical observations. ;[0093] FIG. 5 shows a flow chart of data assimilation based on a localized equivalent weights particle filter with statistical observations. ;Detailed Description ;[0094] The disclosure will be further described in detail below with reference to the accompanying drawings and specific implementations. ;[0095] In order to more simply and clearly describe the specific implementation steps of an assimilation method based on a localized weight particle filter of statistical observations, a simple Lorenz-96 model is taken as an example for a brief description. This model can better reflect the influence of a localization method on assimilation results. ;[0096] Step 1: A model integral initial background field is obtained. ;[0097] The Lorenz-96 model is selected because it has 40 state variables, and then, a localization scheme can be better applied. In the Lorenz-96 model, first, an initial starting point of the model is input into a model equation for integrating 1 million steps to spin-up, so that the model variable reaches a chaotic state to avoid model fluctuations and deviations in an assimilation integration process of the model. The model variable reaching the chaotic state is taken as an initial background field of the model, and based on this background field of the model, initial values of the ensemble particles in particle filter are obtained through the model equation, thereby effectively improving the reliability of the states of the ensemble particles obtained by solving the model equation. ;[0098] Step 2: Whether a statistical observation start time is reached is determined, and observations are accumulated to calculate an average value of statistical observations. [0099] In data assimilation of equivalent weights particle filter of statistical historical observations, the future observation information of a traditional method is replaced with statistical observations, ensemble particles are guided to be close to historical statistical results of observations at the assimilation time to increase the number of effective particles in a sampling process, and an appropriate threshold r is selected to determine when to start statistical observations. The selection of V is related to a time interval of the observation information that needs assimilation. When r is too small, more useless observations will be introduced into the statistics, and statistical errors will be increased. When I-is too large, the statistical observation information will be too little, and proposal densities can not be correctly calculated to guide particles to be close to observations. A calculation formula of I is: [00100] [00101] where t ° represents a previous observation time that needs assimilation, t n represents an observation interval between two observations, and t represents a current time. In order to ensure the reliability of statistical observation results, when the observation interval that needs assimilation is larger, I is generally selected to be 0.8-0.9, and when the assimilation interval is smaller, I is generally selected to be 0.5-0.8. A statistical method is used for counting the observations before the start of an assimilation time to calculate an average value, y fyto t is used for representing a time sequence of observation information in the observation interval at an observation position that needs assimilation, and assuming that t is determined to be a statistical observation start time according to V, the average value of statistical observations can be expressed as: [00102] [00103] Step 3: Proposal densities are calculated according to the average value of the statistical observations to adjust ensemble particles, and at a given assimilation time, weights of the particles are calculated by an equivalent weights method to adjust states of the particles. ;[00104] Step 3.1: Probability densities of the particles at observation corresponding positions are calculated according to the average value of the statistical observations. ;[00105] In a Bayesian theory, the solution of a posterior probability density needs to rely on a prior probability density distribution and a likelihood probability solution, a particle filter method uses a conditional posterior probability density distribution based on the Bayesian n n theory, and the posterior probability density distribution for the statistical observation Y can be expressed as: p(x,,Lyi [00106] [00107] where CI (n ' Xin represents a proposal density. For the calculation of proposal densities, probability densities of the ensemble particles need to be calculated first X1 I and a model equation is used for integrating from a particle state at a time n-1. An expression for the probability densities of the ensemble particles is: P (Yu lx7) P lx7 2'1)n) p(j) )= [00108] P x exp -1(y f(r 9)T1' 2) [00109] where represents a probability density from a particle state at a time n-1 to a time n, represents a state of the th particle at the time n, (xi" 1) represents an integral result of the model equation at the time n-1, and Q represents a model error covariance expressed as: ;- ;[00110] M 1 [00111] where Al represents a total number of the ensemble particles, and r7 represents an average value of the ensemble particles. The probability densities of the particles are calculated according to the formula. ;[00112] Step 3.3: Proposal densities of the particles are calculated according to the probability densities of the ensemble particles. ;[00113] Based on the prior density calculated in the previous step in the model equation, proposal densities in equivalent weights particles can be further obtained, and an expression for the proposal densities calculated according to the prior density is: [00114] 1 -f(x;'-')-Kn (I; (xn-I)) [00115] where try(x) represents a proposal density taking an observation Y as a = OH T (HOH target, Kn represents a relaxation stress matrix expressed as, and an expression of K", a model error covariance (--) and an observation error covariance R is: ;AI ;=covfx y =-Dx, [00116] [00117] where x, represents a state of the i th particle at the time n, H represents an observation operator which is generally 1, and represents observation information. In an equivalent weights particle filter based on statistical observations, the calculation of the proposal densities is not limited to the observation information at the assimilation time, and the proposal densities are based on the historical average value of statistical observations before the assimilation time, in order to adjust the positions of the ensemble particles to be close to the observation information before the assimilation time. ;[00118] Step 3.4: An optimal proposal density is selected, and weights of the particles are calculated. ;[00119] The selection of proposal densities is considered to be an important criterion for controlling the positions of the particles and observation information and the calculation of weights of the particles. Selecting an appropriate and optimal proposal density can ensure the number of effective particles in sampling and simultaneously ensure the reliability of weights of the particles, and is also the most important part in the equivalent weights particle filter. To el(xi' IC' = -11 find the optimal proposal density, assuming that, the weight of each of the particles in the ensemble based on statistical observations is expressed as: [00120] (o.), = ti -1 (ye) x,e1) q x, x,,y 1r [00121] where (9, represents a weight of the 1th particle in the ensemble, represents a state of the th particle at a time U,Y represents statistical observation information at the p(x,' ) ,y time, represents a probability density, -1 represents a proposal p(y" ket,' ) density, represents a probability density of the particle state at the time for statistical observations, and finally, the weight of each of the particles through the proposal density is obtained. For adjusting the ensemble particles to be close to statistical observations through the proposal densities, the states of the ensemble particles can be expressed as: [00122] (xi) =f + (r)(y -h(x,1-1)) [00123] where represents an observation operator, a prior relaxation coefficient is BO, and B(r)(yl -(x.f1) ensures that the ensemble particles are close to statistical observations, so the relaxation coefficient can be expressed as: [00124] BO= In-QHT 1?-1 [00125] where T represents a statistical observation start threshold, and h represents a scale factor that controls a degree of relaxation to observations. ;[00126] Step 4: At the given assimilation time, weights of the particles are calculated by an equivalent weights method, and states of the particles at observation corresponding positions are adjusted. ;[00127] Step 4.1: Weights of the particles are calculated based on the ensemble particles adjusted according to the proposal densities. ;[00128] The proposal densities are calculated according to the observation information in assimilation to adjust the positions of the particles in the ensemble, and at the assimilation time, the particles are closer to the observation information to ensure that most of the particles obtain equal weights in a posterior probability density function at the assimilation time. The proposal densities enable each of the particles in the ensemble to be close to the observation information, so that each of the particles obtains nearly the same weight. An expression of the weight of each of the particles is: [00129] ro - lx7) V Ix 71 q(xn 0 'est [00130] where represents a cumulative multiplication result of proposal density weights in an assimilation time interval, and a minimum weight CD of the particles can be determined and is expressed as: [00131] -log ao, = -log ore' (y" ))] (11Q117 + 1?)1 (y"- [00132] ' 2 f (i)), where (9rn represents a prior weight of the i th particle, and the prior weight is the weight of each of the particles in a proposal density calculation process; Y represents an observation vector; H represents an observation projection operator, and the projection operator is H = I in a simple model; -I. represents a state vector; a superscript T represents a matrix transposition; (2 represents a model error covariance matrix; and R represents an observation error covariance matrix. Assuming that a target weight C, of the particles in the ensemble is obtained, in order to ensure that 80% of the particles in the particle ensemble can reach the calculated weight to avoid the occurrence of particle degradation and depletion problems: [00133] -logo staTger = + -1 (/' -Ht (x," 1)) (HQ1-11 + R) -Hf 1)) , [00134] therefore, it is found that most of the particles in the ensemble in a state at the time n can maintain the same weight, and some particles that do not reach the equivalent weights can be adjusted by a resampling method. ;[00135] Step 4.2: States of the particles are adjusted according to weights. ;[00136] After weights of the ensemble particles are obtained, states of the particles are adjusted according to the weights, and then, the state at the time n can be expressed as: [00137] a,K(yn -fif (x)), [00138] where Y represents an observation vector; H represents an observation R)1 K = OH' (HQH projection operator ( H = I);X represents a state vector; and, Q represents a model error covariance, and R represents an observation error covariance. When the weights of the ensemble particles are approximately equal in the equivalent weights particle filter, the vector a' can be expressed as: [00139] a = 1-V1-11, a, , [00140] where in two expressions a O =.5x,1 12-11-1Kxand b, = 0.5x: x, -log co:"'n -log col' x = Y" (x,' 1) weight, and LU represents a relative weight of the particles at the current time. In order to ensure a random effect of the states of the particles in the ensemble, random terms are added, and a final analytical equation is obtained: n-1+ a,K (y" -ny (x;')) [00141] represents a selected target [00142] Step 5: The ensemble particles are adjusted by a resampling method to keep the number of the particles stable, and the states of the particles at observation corresponding positions are updated. ;[00143] The resampling method is used for removing particles with smaller weights in the ensemble particles at the observation corresponding positions in the equivalent weights particle filter, simultaneously ensuring the stability of the total number of the particles in the ensemble, and replicating particles with larger weights. The equivalent weights particle filter method can ensure that most of the particles are preserved. However, in order to prevent a very small number of the particles from meeting the requirement of an equivalent weights in an equivalent weights process, the resampling method is used for performing the final state adjustment to keep the number of the particles stable. ;[00144] Step 6: Weights of the particles around assimilation observation corresponding positions are determined by a localization function. ;[00145] After the states of the ensemble particles at the assimilation observation corresponding positions are updated by the equivalent weights particle filter method and the resampling method, the localization function is used to continue to adjust the states of the ensemble particles within an observation influence radius r. The selection of the localization function refers to localization schemes in localized particle filter methods. In the calculation process, the localization function is used for describing the relationship between the weights of the ensemble particles in a region and a given observation position, and the weight in the particle filter represents a likelihood probability of observations. In the calculation, a localization operator r] is used, the operator is used for describing a relative position coefficient of the observation Y and the state xi of the ensemble particles, and a corresponding particle weight is expressed as: [00146] co; = Pcv x l[y where the localization operator is mai [00147] nly used in local analysis to determine relative position information of the ensemble particles and observation information. When the observation Y coincides with the xi of the ensemble particles, the maximum value 1 of the function is taken, and when the distance between the two exceeds the given influence radius 1, the function value is 0, which means that the observation has no adjusting effect on the state of the ensemble particles. Finally, vector weights of the ensemble particles are expressed as: [ dty,-hi ()}2/2 of = exp --Li rily, (x)12 [00148] I = 121 exp [-, r ;- ;[00149] where represents an observation error covariance, and h-1 represents a measurement operator. Therefore, it can be seen that the weights of the ensemble particles are not only related to the observation information, but also related to relative positions of observations and particles. ;[00150] Step 7: Weights of the particles are updated according to localization weights, and states of surrounding particles are updated. ;[00151] After the weights of the ensemble particles within the influence radius of the observation information at the assimilation time are obtained, the states of the corresponding ensemble particles need to be adjusted according to the weights. Before state adjustment, the weights of the ensemble particles need to be normalized to ensure that the sum of the weights in the ensemble is 1. A normalized weight formula can be expressed as: ;WI = ;[00152] [00153] The states of the particles in the ensemble are readjusted according to the normalized weight, and the states of posterior particles can be expressed as: [00154] , vci) -011 (I, -1) -OH =x ± ri 0 lXA: -X)± r, 0 (X" -X) [00155] where X represents a posterior average value, k, represents the nth sampled particle, and vectors and r2 can enable new particles to form sampled particles to be linearly combined with prior particles to finally realize posterior update of the states of localized particles, wherein calculation formulas of the vectors rl and 112 can be expressed as: [00156] = (1 la 1 A) -0, 2 -xi)j N, - C N,(1-1[x y cr]) ;C-- ;[00157] [00158] [00159] where represents error variances of all observations up to. In this process, the posterior correlations between the states of the ensemble particles are ignored, but corresponding correlations are provided through corresponding sampling steps. For the I v y,* 1 c 0 lim Since when /[x 1 localization operator, when, and cO 1.j a posterior variance is approximately equal to a sampled particle variance, and in a similar way, limr,7=0 can be obtained, which enables the posterior particles to obtain the states of the lim r =I l[x, , r] -> 0 c -> l[x" y r] = 0, a sampled particles. When, and. Since when urn =0 posterior variance is equal to a prior variance, and in a similar way, can be obtained.

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

[00160] Step 8: Posterior estimation values of states in the localized equivalent weights particle filter with statistical observations are calculated.

[00161] The states of all particles at the observation corresponding positions and within the observation influence radius are updated according to the observation information, a posterior estimation ensemble average value M,=, of the states is calculated, the updated posterior estimation ensemble average value as an initial value of an analytical model is substituted into the model integral equation again for the next prediction and assimilation, and the above steps are repeated within the assimilation time when there are available observations to obtain a final analytical field which can be used as a data field for reflecting a current environment state.

[00162] The disclosure provides a data assimilation technology based on a localized equivalent weights particle filter with statistical observations. Compared with the traditional data assimilation technology based on an equivalent weights particle filter of statistical observations, the disclosure has the following significant characteristics: A localization scheme adapted to applications of the equivalent weights particle filter method of statistical observations is proposed to effectively solve the limitations of this method for applications of complex high-latitude grid models. At the assimilation time, a statistical observation method is used for updating the states of the ensemble particles at observation corresponding positions; then, according to the position information and localization parameters of assimilation observations, the observations are used for adjusting the weights of available ensemble particles that are sensitive to the observations around the observation positions; and then, according to the weights, the states of the ensemble particles in the vicinity of the observations are adjusted, thereby improving the utilization rate of the observations and improving the assimilation quality. The method proposed in this patent can effectively improve the assimilation ability of the traditional equivalent weights particle filter method of statistical observations under complex grid models and sparse observation conditions, can effectively improve the assimilation quality, and also can improve the application prospects of the statistical equivalent weights particle filter method.

Claims

What is claimed is: 1. A method for assimilating data based on a localized equivalent weights particle filter with statistical observations, characterized by comprising the following steps: step 1: obtaining a model integral initial background field, characterized by: first introducing a model initial field into a model integral equation, integrating the model equation first to enable the model equation to reach a chaotic state, and performing subsequent model state integration by taking a model variable reaching the chaotic state as an initial background field and taking the background field as an integration starting point; step 2: determining whether a statistical observation start time is reached, and accumulating observations to calculate an average value of statistical observations, characterized by determining a start time of statistical observation calculations according to a given / value, wherein a calculation formula of T is: T (tj t0) -t0) wherein t ° represents a previous observation time that needs assimilation, t n represents an observation interval between two observations, and t represents a current time; in order to ensure the reliability of statistical observation results, when the observation interval that needs assimilation is larger, / is generally selected to be 0.8-0.9, and when the assimilation interval is smaller, T is generally selected to be 0.5-0.8; a statistical method is used for counting the observations before the start of an assimilation time to calculate an average value, Yf = frt°, , yr') is used for representing a time sequence of observation information in the observation interval at an observation position that needs assimilation, and assuming that t j is determined to be a statistical observation start time according to T, the average value of statistical observations can be expressed as: yr, (t" -t' + 1 step 3: calculating proposal densities according to the average value of the statistical observations to adjust ensemble particles, and at a given assimilation time, calculating weights of the particles by an equivalent weights method to adjust states of the particles, characterized by: step 3.1: calculating probability densities of the particles at observation corresponding positions according to the average value of the statistical observations, characterized in that: in a Bayesian theory, the solution of a posterior probability density needs to rely on a prior probability density distribution and a likelihood probability solution, a particle filter method uses a conditional posterior probability density distribution based on the Bayesian n theory, and the posterior probability density distribution for the statistical observation Y can be expressed as: n n -1) (rn kbz) t P V kit') P ) AT p (yn) (xn " tin lx7t-vil) wherein represents a proposal density; for the calculation of proposal densities, probability densities of the ensemble particles need to be calculated first, and a model equation is used for integrating from a particle state x, at a time n-1; an expression for the probability densities of the ensemble particles is: p(xr lx) x exp [--21 -( I)) P(x wherein ' represents a probability density from a particle state at a time n-1 to a time n, represents a state of the th particle at the time n, -f ( ' ) represents an integral result of the model equation at the time n-1, and Q represents a model error covariance expressed as: m wherein Al represents a total number of the ensemble particles, and.7 represents an average value of the ensemble particles; and the probability densities of the particles are calculated according to the formula; step 3.2: calculating proposal densities of the particles according to the probability densities of the ensemble particles, characterized in that: based on the prior density calculated in the previous step in the model equation, proposal densities in equivalent weights particles can be further obtained, and an expression for the proposal densities calculated according to the prior density is: X.7,) 61(x 7f: e X p ( 1)) K -H (xn 1)) Q -(x7 '- -H (x17 I))I -qwherein 11: )1 represents a proposal density taking an observation Y as a K" = OH T (IVITT I, and an target, K" represents a relaxation stress matrix expressed as expression of K" , a model error covariance and an observation error covariance R is: t4 R = COV 37} = -E(X, y)2 M I-1 wherein represents a state of the th particle at the time n, H represents an observation operator which is generally 1, and Y represents observation information; in an equivalent weights particle filter based on statistical observations, the calculation of the proposal densities is not limited to the observation information at the assimilation time, and the proposal densities are based on the historical average value of statistical observations before the assimilation time, in order to adjust the positions of the ensemble particles to be close to the observation information before the assimilation time; and step 3.3: selecting an optimal proposal density, and calculating weights of the particles, characterized in that: to find the optimal proposal density, assuming that q (x11)= P(x7kil,Y"), the weight of each of the particles in the ensemble based on statistical observations is expressed as: (a), )1') = p (y1) f t-" ) X, , wherein (0, represents a weight of the i th particle in the ensemble, represents a state of the th particle at a time, Y represents statistical observation information at the pfrtime I' , r' r r x, , y represents a probability density, - represents a proposal p(y" density, represents a probability density of the particle state at the time for statistical observations, and finally, the weight of each of the particles through the proposal density is obtained; for adjusting the ensemble particles to be close to statistical observations through the proposal densities, the states of the ensemble particles can be expressed as: (x,)'I =,f(x, )+B (r)(y -wherein h Xit represents an observation operator, a prior relaxation coefficient is B(r), and ensures that the ensemble particles are close to statistical observations, so the relaxation coefficient can be expressed as: B (r) = r OHT wherein T represents a statistical observation start threshold, and b represents a scale factor that controls a degree of relaxation to observations; step 4: at the given assimilation time, calculating weights of the particles by an equivalent weights method, and adjusting states of the particles at observation corresponding positions, characterized by: step 4.1: calculating weights of the particles based on the ensemble particles adjusted according to proposal densities, characterized in that: the proposal densities enable each of the particles in the ensemble to be close to the observation information, so that each of the particles obtains nearly the same weight; an expression of the weight of each of the particles is: sl P(Ym Ix) P 11-t-') (x" Y") rest wherein represents a cumulative multiplication result of proposal density weights in an assimilation time interval, and a minimum weight c°, of the particles can be determined and is expressed as: = -log Orr 1 (y" Hf (xn I)) (HQI -1' 1?) (312 Hf 7 1)) wherein represents a prior weight of the th particle, and the prior weight is the weight of each of the particles in a proposal density calculation process; -represents an observation vector; H represents an observation projection operator, and the projection operator is H = I in a simple model; X represents a state vector; a superscript T represents a matrix transposition; -Q represents a model error covariance matrix; R represents an observation error covariance matrix; assuming that a target weight C1 of the particles in the ensemble is obtained, in order to ensure that 80% of the particles in the particle ensemble can reach the calculated weight to avoid the occurrence of particle degradation and depletion problems: = -log a)7 +4/ -10x; (HQI + R) (y" -( )) therefore, it is found that most of the particles in the ensemble in a state at the time n can maintain the same weight, and some particles that do not reach the equivalent weights can be adjusted by a resampling method; and step 4.2: adjusting states of the particles according to weights, characterized in that: after weights of the ensemble particles are obtained, states of the particles are adjusted according to the weights, and then, the state at the time n can be expressed as: f ( 1)* ceiK (y" -Iff wherein -v represents an observation vector; H represents an observation projection K = OH T (110HT ±R)1 Q represents a model operator (H =1); represents a state vector error covariance, and R represents an observation error covariance; when the weights of the ensemble particles are approximately equal in the equivalent weights particle filter, the vector a, can be expressed as: a =1-'11-I a wherein in two expressions a^ = 0.5x R 11Kx and 17, = 0.5x{R'1x log re' -log orn x = y' - (x;')cat nrget represents a selected target weight, and represents a relative weight of the particles at the current time; in order to ensure a random effect of the states of the particles in the ensemble, random terms are added, and a final analytical equation is obtained: 1-1 a,K (yr' - (x)) step 5: adjusting the ensemble particles by a resampling method to keep the number of the particles stable, and updating the states of the particles at observation corresponding positions; step 6: determining weights of the particles around assimilation observation corresponding positions by a localization function, characterized in that: after the states of the ensemble particles at the assimilation observation corresponding positions are updated by the equivalent weights particle filter method and the resampling method, the localization function is used to continue to adjust the states of the ensemble particles within an observation influence radius r; the selection of the localization function refers to localization schemes in localized particle filter methods; in the calculation process, the localization function is used for describing the relationship between the weights of the ensemble particles in a region and a given observation position, and the weight in the particle filter represents a likelihood probability of observations; in the calculation, a localization operator I[y'x is used, the operator is used for describing a relative position coefficient of the observation and the state of the ensemble particles, and a corresponding particle weight is expressed as: CO x,,r1+ 1, wherein the localization operator i[y,x,,uiis mainly used in local analysis to determine relative position information of the ensemble particles and observation information; when the observation -V coincides with the of the ensemble particles, the maximum value 1 of the function is taken, and when the distance between the two exceeds the given influence radius T, the function value is 0, which means that the observation has no adjusting effect on the state of the ensemble particles; finally, vector weights of the ensemble particles are expressed as: = exp [y1,, r]_/PI 1[Y1,x1,d{Y,-17,(0}2 I /2 = exp 2 1=1 wherein cri represents an observation error covariance, and h represents a measurement operator; and therefore, it can be seen that the weights of the ensemble particles are not only related to the observation information, but also related to relative positions of observations and particles; step 7: updating weights of the particles according to localization weights, and updating states of surrounding particles, characterized in that: after the weights of the ensemble particles within the influence radius of the observation information at the assimilation time are obtained, the states of the corresponding ensemble particles need to be adjusted according to the weights; before state adjustment, the weights of the ensemble particles need to be normalized to ensure that the sum of the weights in the ensemble is 1; a normalized weight formula can be expressed as: the states of the particles in the ensemble are readjusted according to the normalized weight, and the states of posterior particles can be expressed as: 0) ()),-1) -X + -x 1+r -I) -X(' ) ) wherein "c represents a posterior average value, kn represents the nth sampled particle, vectors r1 and 12 can enable new particles to form sampled particles to be linearly combined with prior particles to finally realize posterior update of the states of localized particles, and calculation formulas of the vectors ri and 12 can be expressed as: (v 1) -CYO 2 - C./ (X,); - )] r =c r 2.1 1,1 * N,(1-1[x2,yEr]) = l[x y JP- o-wherein represents error variances of all observations up to Y,; in this process, the posterior correlations between the states of the ensemble particles are ignored, but corresponding correlations are provided through corresponding sampling steps; for the limr l[x y r] c -> 0 >0 since when nxi'yi 'r]= 1 localization operator, when, and, a posterior variance is approximately equal to a sampled particle variance, and in a similar way, lim fr; =0 can be obtained, which enables the posterior particles to obtain the states of the lim r.j =1 I[x 1] -> 0 c -> oc l[t y r] -° sampled particles; when y, and ") ; since when, a urnr, j = 0 posterior variance is equal to a prior variance, and in a similar way, can be obtained; -r this sampling method provides a way to adjust the particles to fit a general Bayesian posterior solution in the vicinity of observations; because each of the sampled particles is combined with a prior particle, a posterior ensemble comprises model states unique to the ensemble particles, which avoids the collapse of an ensemble variance during observation assimilation; and step 8: calculating posterior estimation values of states in the localized equivalent weights particle filter with statistical observations, characterized by: updating the states of all particles at the observation corresponding positions and within the observation influence radius according to the observation information, calculating a 1 m posterior estimation ensemble average value M 7=1 of the states, substituting the updated posterior estimation ensemble average value as an initial value of an analytical model into the model integral equation again for the next prediction and assimilation, and repeating the above steps within the assimilation time when there are available observations to obtain a final analytical field which can be used as a data field for reflecting a current environment state.