CN113360854B - Data assimilation method based on adaptive covariance expansion - Google Patents
Data assimilation method based on adaptive covariance expansion Download PDFInfo
- Publication number
- CN113360854B CN113360854B CN202110911181.3A CN202110911181A CN113360854B CN 113360854 B CN113360854 B CN 113360854B CN 202110911181 A CN202110911181 A CN 202110911181A CN 113360854 B CN113360854 B CN 113360854B
- Authority
- CN
- China
- Prior art keywords
- expansion
- analysis
- observation
- members
- expansion factor
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 139
- 230000003044 adaptive effect Effects 0.000 title claims abstract description 12
- 238000004458 analytical method Methods 0.000 claims abstract description 63
- 239000013598 vector Substances 0.000 claims abstract description 37
- 239000011159 matrix material Substances 0.000 claims abstract description 31
- 230000010339 dilation Effects 0.000 claims abstract description 22
- 230000008569 process Effects 0.000 claims abstract description 19
- 238000001914 filtration Methods 0.000 claims abstract description 18
- 230000010354 integration Effects 0.000 claims abstract description 7
- 230000006870 function Effects 0.000 description 25
- 230000004807 localization Effects 0.000 description 14
- 238000012360 testing method Methods 0.000 description 13
- 238000002474 experimental method Methods 0.000 description 12
- 230000008859 change Effects 0.000 description 11
- 238000004364 calculation method Methods 0.000 description 7
- 238000009826 distribution Methods 0.000 description 7
- 238000004422 calculation algorithm Methods 0.000 description 6
- 230000000694 effects Effects 0.000 description 6
- 230000006872 improvement Effects 0.000 description 6
- 238000005070 sampling Methods 0.000 description 5
- 230000035945 sensitivity Effects 0.000 description 3
- 230000007423 decrease Effects 0.000 description 2
- 230000007246 mechanism Effects 0.000 description 2
- 238000006467 substitution reaction Methods 0.000 description 2
- 230000008961 swelling Effects 0.000 description 2
- 241000282414 Homo sapiens Species 0.000 description 1
- 239000000654 additive Substances 0.000 description 1
- 230000000996 additive effect Effects 0.000 description 1
- 238000012271 agricultural production Methods 0.000 description 1
- 230000004075 alteration Effects 0.000 description 1
- HTIQEAQVCYTUBX-UHFFFAOYSA-N amlodipine Chemical compound CCOC(=O)C1=C(COCCN)NC(C)=C(C(=O)OC)C1C1=CC=CC=C1Cl HTIQEAQVCYTUBX-UHFFFAOYSA-N 0.000 description 1
- 230000009286 beneficial effect Effects 0.000 description 1
- 230000000739 chaotic effect Effects 0.000 description 1
- 238000012512 characterization method Methods 0.000 description 1
- 230000019771 cognition Effects 0.000 description 1
- 238000012937 correction Methods 0.000 description 1
- 238000002790 cross-validation Methods 0.000 description 1
- 238000011161 development Methods 0.000 description 1
- 238000011156 evaluation Methods 0.000 description 1
- 239000012530 fluid Substances 0.000 description 1
- 238000004519 manufacturing process Methods 0.000 description 1
- 238000005259 measurement Methods 0.000 description 1
- 238000012986 modification Methods 0.000 description 1
- 230000004048 modification Effects 0.000 description 1
- 238000005457 optimization Methods 0.000 description 1
- 238000005293 physical law Methods 0.000 description 1
- 230000001737 promoting effect Effects 0.000 description 1
- 238000009877 rendering Methods 0.000 description 1
- 238000011160 research Methods 0.000 description 1
- 238000010187 selection method Methods 0.000 description 1
- 230000009466 transformation Effects 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/18—Complex mathematical operations for evaluating statistical data, e.g. average values, frequency distributions, probability functions, regression analysis
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01W—METEOROLOGY
- G01W1/00—Meteorology
- G01W1/10—Devices for predicting weather conditions
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/16—Matrix or vector computation, e.g. matrix-matrix or matrix-vector multiplication, matrix factorization
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Data Mining & Analysis (AREA)
- Computational Mathematics (AREA)
- Theoretical Computer Science (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Software Systems (AREA)
- Databases & Information Systems (AREA)
- Algebra (AREA)
- General Engineering & Computer Science (AREA)
- Operations Research (AREA)
- Life Sciences & Earth Sciences (AREA)
- Environmental & Geological Engineering (AREA)
- Bioinformatics & Cheminformatics (AREA)
- Bioinformatics & Computational Biology (AREA)
- Evolutionary Biology (AREA)
- Computing Systems (AREA)
- Probability & Statistics with Applications (AREA)
- Atmospheric Sciences (AREA)
- Biodiversity & Conservation Biology (AREA)
- Ecology (AREA)
- Environmental Sciences (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
Abstract
The invention discloses a data assimilation method based on adaptive covariance expansion, which comprises the following steps: obtaining atmospheric observation values based ont‑1, performing mode integration on the analysis value to obtain analysis timetThe forecast field of (1); according to the forecast field set, estimating to obtaintEnsemble prediction error covariance matrix P for a time instanttAnd a dilation factor vector; updating the set members of the set Kalman filtering in the initial assimilation process by the expansion factor vectors to increase the set variance of the new set Kalman filtering to form new set members, and performing iterative updating on the new set members by using a set Kalman filtering method to obtain final analysis set members; pattern prediction is performed using the members of the analysis set as initial fields. According to the invention, through adjustment of the expansion factor, the updated prediction error covariance matrix conforms to the statistical relationship of the prediction error covariance, the innovation amount and the observation error covariance, the calculated expansion factor is more reasonable, and the assimilation performance is obviously improved.
Description
Technical Field
The invention belongs to the technical field of weather forecasting, and particularly relates to a data assimilation method based on adaptive covariance expansion.
Background
From ancient times to present, no matter how society develops and changes towards generations, the natural exploration of people never stops. From prehistoric mankind to modern civilization, people are constantly trying to discover and master the weather development law so as to guide work and life of themselves. Accurate weather forecast has important influence on agricultural production, military activities and daily life of human beings. Abbe proposed at the beginning of the twentieth century that the laws of physics can be used for weather forecasting, who recognizes that predicting atmospheric conditions can be regarded as a mathematical physics problem at the beginning, the basic idea is to obtain the true state of the atmosphere at some future time by solving a controlled partial differential equation that describes the hydrodynamic and thermodynamic characteristics of the atmosphere, based on the currently observed weather conditions. This is also the theoretical basis for Numerical Weather forecasting (NWP). Under such a typical mathematical physics problem framework, three most basic conditions need to be satisfied to obtain accurate weather forecast numerical results: firstly, a state value which reflects the atmosphere at the initial moment of mode forecast as accurately as possible; secondly, the numerical mode of the atmospheric and marine physical motion states can be accurately described; three are the boundary conditions required for the numerical mode. It can be seen that the error between the atmospheric forecast status obtained by NWP and the actual atmospheric status is mainly derived from three aspects: the method comprises the following steps of firstly, predicting an error between a background field at an initial moment and an actual atmospheric state, secondly, predicting a mode error, and accurately expressing the error generated by the actual physical process in the weather evolution process by the current atmospheric mode; and thirdly, the error of the given boundary condition from the actual state. It is worth noting that the chaotic nature of the weather system determines that even a small initial error may grow rapidly with the increasing duration of the forecast, eventually leading to significant deviations in the forecast results. It is therefore clearly not practical to consider numerical prediction as an accurate initial boundary value problem. While improving the mode accuracy, it is very important how to fully and effectively utilize the current available observation data to correct the system state at the forecast initial time to obtain a more accurate forecast initial field and boundary field. Data assimilation is an effective way to effectively fuse observations and patterns to obtain better initial values and reasonable boundary values.
It is simply understood that data assimilation exists because the factor value mode and the observation system cannot independently accomplish accurate estimation of the system state. The numerical mode aims to reveal the intrinsic physical mechanism of the system and give the change rule of the system along with the time based on the mechanism and the initial state of the system. The establishment of such a model itself requires the estimation of the relevant parameters of the model in dependence on the observed data. For a real weather system, people cannot predict the real system state in advance, so that a 'real' numerical mode capable of representing the atmospheric evolution state cannot be obtained, and only approximation can be performed through the existing priori knowledge. The observation represents the 'true value' of the system at a certain moment in a certain sense, but the real observation is rare, the space-time distribution is extremely uneven, and various errors such as instrument errors and representation errors exist, so that the real state of the system cannot be completely replaced by the observation. This requires that the observed information be fused with the patterns in a reliable way, i.e. often referred to as data assimilation.
Data assimilation is an important means for effectively fusing a numerical mode and observation, and the method not only can provide an accurate and reasonable initial value for mode prediction, deepens cognition on the physical law of atmospheric state evolution, but also can perform cross validation on an observation system. Therefore, it has been a very important research area in numerical weather forecasting.
The error covariance is constructed by the prediction set in the ensemble Kalman filtering, so that the huge calculation amount caused by integrating the ensemble error covariance in the traditional Kalman filtering method is reduced, and meanwhile, the method can be applied to a nonlinear system. The method is simple in concept, easy to implement and beneficial to parallelism, and becomes an assimilation method with wide application. Because the set membership in the prediction set is often smaller than the mode dimension, a sampling error introduced, a system error existing in the mode, a characterization error existing in observation acquisition, an instrument error and other errors exist, the set variance constructed by using the prediction set is often underestimated. This directly leads to the situation that the filter divergence occurs when the ensemble kalman filter method is applied to some high-dimensional systems.
Covariance dilation is one of the most dominant methods in current ensemble filtering methods to prevent filter divergence due to ensemble error covariance being underestimated. The main idea is to improve the variance of the set error by various methods, so as to compensate the problem that the covariance of the prediction error is underestimated due to the mode error, the sampling error and other unknown source errors. The expansion method mainly comprises three expansion methods, namely additive expansion, multiplicative expansion and prior loose expansion. Of which multiplicative expansion is most widely used. The initial multiplicative dilation algorithm dilates the covariance by defining a dilation factor. In 2009, Anderson indicated that in the global weather forecast system, north america were observed more densely, and variance could be underestimated due to the existence of pattern errors and sampling errors, but in south america where observation is sparser, if the same dilation method as in north america is used, the variance of pattern variables in south america would also increase, leading to inconsistency with climatic conditions, and more seriously, the physical balance of the pattern field would be destroyed, thus rendering pattern results unusable. It follows that it is important to take into account spatial variations in the expansion process. For multiplicative expansion, a more general representation is to consider the expansion factor as a expansion matrix D. And applying the matrix D to the background error covariance to obtain an expanded background error covariance matrix, wherein the key point of multiplicative expansion is to determine an expansion operator matrix D.
Disclosure of Invention
In view of this, the present invention provides a adaptive covariance expansion method (hereinafter referred to as an SDACI method) based on innovation statistics and spatial dependence, which defines an expansion factor vector having the same dimension as a mode variable, and solves a target function by constructing the function, and the main idea of the function is that an updated prediction error covariance matrix can maximally conform to a statistical relationship about a prediction error covariance, an innovation amount, and an observation error covariance in the innovation statistics by adjusting the expansion factor.
The invention discloses a data assimilation method based on adaptive covariance expansion, which comprises the following steps:
obtaining atmospheric observation values based ont-1, performing mode integration on the analysis value to obtain analysis timetThe forecast field of (1);
according to the forecast set, estimating to obtaintEnsemble prediction error covariance matrix P for a time instanttAnd a dilation factor vector;
substituting the expansion factor vector into the analysis process to obtain each analysis set member;
Weather forecasting is performed using the analysis set members.
Further, the analysis time is calculated according to the following formulatThe forecast field of (1):
wherein,is composed oft-1 at the first momenti The number of members of an analysis set,is composed oftAt the first momentiThe members of each of the forecast sets,mas to the number of members of the set,M t-1is thatt-a mode operator at time 1.
Further, it is calculated according to the following formulaEnsemble prediction error covariance matrixP t :
Further, the inflation factor vector is calculated according to the following formula:
constructing an objective function according to the following formula, and solving the minimum value of the objective function to finally obtain an expansion factor vector:
,for a given statistical prediction error covariance matrix at time t,d is an innovation quantity which represents the deviation between the background field and the observation field in the assimilation process, and H is a tangential operator of a nonlinear observation operator, namely for any state vector x0, the innovation quantity comprises the following components:
wherein the swellingVector of expansion factorIs one dimension ofnEach element representing an expansion factor corresponding to a mode variable.
Further, the expansion factor vector is usedSubstituting the set analysis member updating process of the set Kalman filtering based on random disturbance to solve each analysis set member:
Where K is the Kalman gain matrix, yoRepresenting a time series of observation vectors,is a non-linear observation operator, subscripttRepresentstAt the moment of time, the time of day,to representtAt the first momentiAnd (4) disturbing the observation error.
The method converts the determination of the expansion factor into the optimization process of the objective function so as to determine the established parameters in the expansion method, fully considers the influence of the spatial distribution on the expansion factor, and can effectively improve the assimilation performance compared with the expansion factor method only considering the time distribution.
Drawings
FIG. 1 is a flow chart of a data assimilation method based on adaptive covariance expansion according to the present invention;
FIG. 2 analysis of mass comparison plots for three expansion methods with set membership of 10;
FIG. 3 analysis of mass comparison plots for three expansion methods with set membership of 20;
FIG. 4 analysis of mass comparison plots for three expansion methods with set membership of 30;
FIG. 5 is a graph of the fixed expansion method analyzing the trend of performance with set membership;
FIG. 6 shows the Wu method analyzing the trend of performance with the number of members in a set;
FIG. 7 is a graph of the trend of performance versus set membership for the present invention;
FIG. 8 is a graph comparing the performance of the present invention and the Wu process;
FIG. 9 is a graph showing the performance change of the assimilation results obtained by the three expansion methods when the standard deviation of the observation error is 0.5;
FIG. 10 is a graph showing the variation in performance of assimilation results obtained by three expansion methods with an observation error standard deviation of 2.0;
FIG. 11 analysis performance of the fixed dilation method with set membership of 20 as a function of standard deviation of observed error;
FIG. 12 is a graph of the analytical performance of the Wu method as a function of standard deviation of observation error for a set membership of 20;
FIG. 13 is a graph of the analytical performance of the present invention as a function of standard deviation of observation error for a set membership of 20;
FIG. 14 is a graph of the scale of the improvement in analytical performance of the present invention relative to the Wu method at standard deviations of 0.5 and 2.0 for observed error.
Detailed Description
The invention is further described with reference to the accompanying drawings, but the invention is not limited in any way, and any alterations or substitutions based on the teaching of the invention are within the scope of the invention.
The innovation quantity is the difference value between an observed field and a background field, and can measure the deviation between the observed field and the background field, which is a very important variable in assimilation analysis. The invention provides a self-adaptive expansion method based on innovation amount statistics.
In the framework of linear system estimation theory, regardless of the forecasting process and neglecting time factors, the analysis field x of a single analysis moment a Can be given by equation 1:
wherein,a background field is represented by a field of the background,the increment of the analysis is represented by,
d is an innovation quantity representing the observation field y for the gain matrix in the analysis processoAndthe difference between them. x is the number of b In the form of a prior background field,the nonlinear observation operator has the main function of projecting the mode state to an observation space, and is essentially the combination of a series of mode point-to-observation point interpolation operators and mode control variable-to-observation transformation operators. And H is a non-linear observation operatorTangential operators of, i.e. for, arbitrary state vectors x0The method comprises the following steps:
wherein
According to the above definition, the innovation quantity d can then be expressed as:
wherein x istIn order to be a true field of the pattern,andrespectively, indicating observation error and background error.
The expected estimate of the amount of innovation d is:
=0 (6)
similarly, the covariance estimate of the innovation is expressed as:
considering that the covariance statistics of the innovation d and the background error covariance and observation error covariance satisfy the above equation, it can be considered to estimate the inflation factor using this relationship.
Given at Statistical prediction error covariance matrix P of time instants f By using the idea of multiplicative expansion method, a certain expansion factor is givenλActing on P f Then the covariance of the background error after dilation isλP f . Assuming the amount of innovation is
The error covariance of the random vector d can be expressed by the following second order regression equation
further obtain
Will be provided withλExpanding the vector with one dimension same as the dimension of the pattern spaceλ= (λ 1 , λ 2 , . . . , λ n )TAnd an objective function is provided based on the formula (11) to measure the quality of the expansion adjustment, and the optimal expansion factor vector is provided by solving the objective function. The objective function is set as follows:
λis one dimension ofn Each element representing an expansion factor corresponding to a mode variable,andrespectively representi Expansion factorλ i Lower and upper bounds. As can be seen from the definition of the objective function, the essence of the objective function is to find a vector of expansion factorsλAfter correction ofThe statistical conclusion according with the innovation amount can be updated. In special cases whenλ i The same value is taken, that is, the spatial distribution of the expansion factor is not considered, which is the method proposed by Wu (Wu method for short).
To solve the problem (12), Ψ Inf: (12) is first calculatedλ) For each ofλ i To find
therefore, it is not only easy to use
From this it can be derived
According to the definition of Q and A, the method is easy to obtainI.e. only need to findIs ready to obtainAs can be seen from the definition of A:
while
This can be deduced:
thus, equation 19 can be simplified
Then
From the above analysis, the objective function Ψ Inf (is calculatedλ) Value of (D) and gradient function thereofSee algorithms 0.1 and 0.2 for detailed descriptions of specific steps.
Algorithm 0.1 Objective function calculation based on innovation statistics
Knowing the dilation factor vector λ, the set of prediction fieldsens_pObservation vector y o Observation error covariance matrix R, observation operator H, set membershipmNumber of mode statesn
1: initializing an objective function value costfun = 0;
3:for i=1 to m do
4:for j=1 to n do
6:end for
7:end for
9:Q ← ddT;
10:Q ← Q – R;
11:for i =1 to m do
14: end for
15:C ← HXf;
16:A ← CCT;
17:Q ← Q - A;
18:Q ← QQT;
19:costfun ← Tr(Q);
20:return costfun.
algorithm 0.2 Objective function gradient calculation based on innovation statistics
Knowing the dilation factor vector λ, the set of prediction fieldsens_pObservation vector y o Observation error covariance matrix R, observation operator H, set membershipmNumber of mode statesn
1: initializing objective function value gradient = 0;
3:for i = 1 to m do
4:for j = 1 to n do
6:end for
7: end for
9:Q ← ddT;
10:Q ← Q – R;
11:for i = 1 to m do
14:end for
15:B ← Xf(HXf)T;
16:C ← HXf;
17:A ← CCT;
18:Q ← A – Q;
19:Q ← QQT;
20:D ← BQH;
21:for i = 1 to n do
22:gradient[i] ← 2λi -1ei TDei;
23:end for
24:return gradient.
therefore, as shown in fig. 1, the data assimilation method based on adaptive covariance expansion of the present invention mainly comprises the following steps:
s10: based ont-1, performing mode integration on the analysis value to obtain analysis timetThe forecast field of (1);
wherein,is composed oft-1 at the first momenti The number of members of an analysis set,is composed oftAt the first momentiThe members of each of the forecast sets,mis the number of members of the set.
S20: according to the forecast set, estimating to obtaintEnsemble prediction error covariance matrix P for a time instanttAnd a dilation factor vector;
To findObtaining an expansion factor vector, wherein R is a diagonal matrix,the background error covariance matrix after dilation:,,andrespectively represent the ith expansion factorLower and upper bounds of (1);
s30: substituting the expansion factor vector into the set analysis member updating process of the set Kalman filtering based on random disturbance to solve the members of each analysis set。
Where K is the Kalman gain matrix, yoRepresenting a time series of observation vectors,is a non-linear observation operator, subscripttRepresentstAt the moment of time, the time of day,to representtAt the first momentiAnd (4) disturbing the observation error.
The above formula is an analysis set member calculation method obtained in the process of updating the set analysis members of the EnKF based on random disturbance, and if different set kalman filtering algorithms such as a classical kalman filtering method, a deterministic set kalman filtering method, and the like are adopted, the step S30 may be calculated according to a corresponding implementation manner, which is not limited in this embodiment.
Pattern prediction is performed using the members of the analysis set.
The numerical weather forecast refers to a method for predicting an atmospheric motion state and a weather phenomenon in a certain time period by solving a fluid mechanics and thermodynamics equation set describing a weather evolution process through a numerical calculation method under a certain initial value and edge value condition according to an atmospheric actual condition. The numerical model can be simply viewed as a set of nonlinear equations and their numerical solutions. In this numerical solution, the atmospheric state value set is often classified as a column vector in the mathematical problem, the column vector is a mode state vector, and is denoted as x, and the mode state vector is an abstract representation of meteorological elements in the numerical prediction mode.
The data assimilation aims to adopt a certain analysis method to fuse the prior knowledge obtained in a numerical system and observation data information which can be obtained currently at different times, different sources and different resolutions together by a scientific method, and generate a corrected system state at an analysis moment on the basis of meeting physical constraint conditions, so that an accurate initial field is provided for atmospheric (ocean) numerical prediction.
Each set member vector x in the set data assimilation method represents a mode state variable value, a plurality of set members are used for simulating possible distribution of the mode state variable, and a set average value is used for determining a final value of the mode state variable. Therefore, the dimension of one set member represents the number of corresponding mode variables in the numerical forecasting mode, and the variables of the set member represent the mode variables (namely, the common meteorological elements such as temperature, pressure, humidity and wind power) in the forecasting mode. From the perspective of the business forecasting system, the assimilated analysis field actually refers to a meteorological element field obtained by analyzing the mode variables by using different assimilation methods, and provides an initial field for mode forecasting. The mode prediction based on the initial field is common knowledge in the art, and the description of the embodiment is omitted.
To test the performance of the present invention, the present section performed numerical experiments and analyses using the Lorenz-96 model. This is a very widely used numerical model in the theoretical study of the assimilation method. In the experiment of the present invention, the mode initial value was set to
The integration scheme uses a fourth order Runge-Kutta method with time steps of 0.05 dimensionless units, which roughly corresponds to a real-time integration step of 6 hours, assuming a typical time scale for atmospheric dissipation of 5 days. The assimilation method uses the ensemble square root filter method EnSRF. And (3) taking an initial value of the true mode according to the formula (27), operating for 120000 steps to generate, wherein the simulated observation is full observation, namely each mode variable can be obtained through direct observation, and the simulated observation is generated by adding normal random disturbance with the mean value of 0 and the covariance matrix of R to the true state of each lattice point. R is a diagonal array, each observation error is independent, and the standard deviation settings of the observation errors are respectively 0.5,1.0 and 2.0.
The Lorenz-96 mode is a forced dissipative mode whose forced term (Force) is formed by a variableF In the experimental process, the invention selects two test schemes of a perfect mode and an imperfect mode, wherein the perfect mode represents no mode error, namely, the forced item of Lorenz-96 in the assimilation operation processF = 8, imperfect mode test process set differently by process of assimilationF To achieve this, settings are 4, 5, 6, 7, 9, 10, 11 and 12, respectively. The total number of mode integration steps is 100000 steps, the synchronization is carried out once every 4 steps averagely, and the spin-up time is 3000 steps. The invention adopts a fixed covariance localization scheme, namely the localization radius parameter does not change along with the change of time, the localization function is a GC function, and the localization radius is determined through experiments. Meanwhile, in order to test the sensitivity of the invention to the set members, 10, 20 and 30 sets of sample numbers are selected for testing respectively. The evaluation criteria for performance were evaluated by analyzing the root mean square error (RMSE, equation (28))
In order to ensure the reliability of the experiment, 10 times of experiments are carried out in each group, and the average value of 10 experiments is taken as the final result.
The present invention will employ a range of differencesF To simulate different degrees of pattern error. Truth value is composed ofF Production of = 8, assimilation testF Is gotThe values are 4, 5, 6, 7, 8, 9, 10, 11 and 12, respectively. Since the selection of the localization radius directly affects the assimilation effect, the localization parameters need to be determined first. Table 1, table 2, and table 3 show the parameter settings used when obtaining the optimal RMSE using a fixed expansion scheme under different standard deviations of observation errors, different set membership, and different mode error values, respectively, the optimal parameters are selected by manual adjustment, and the local radius value range is 2∼40, the value step length is 2, and the expansion factor is selected according toIs selected by a formula, whereinρThe expansion factor is forgetting factor, and has a value range of 0.1-1.0, a value step length of 0.05, and a value range of 1.0∼3.16. Thus for a fixed set membership the number of experimental groups manually adjusted was 20 x 19 = 380.
From the perspective of optimal parameter setting, as the number of sets increases, the sensitivity of the expansion factor to the number of set members is not high relative to the localization radius, that is, under the same mode error value, the value of the optimal value of the expansion factor does not change much, but the whole trend is descending. As can be seen from Table 1, when the standard deviation of the observed error is 0.5 and the forcing term is 4, the optimal expansion coefficient is not changed due to the increase of the set membership, and is 3.16; taking the set membership numbers 10 and 20 as examples, 5 groups in 9 groups of expansion factor data have the same value, and the value change of the other 4 groups is very small; when the set membership is 20 and 30, the optimal values of the 9 groups of data are completely the same. Similar situations also occur when the standard deviation of the observation error takes values of 1.0 and 2.0. This further shows that the covariance dilation scheme mainly solves the problem of insufficient error estimation, and the set membership mainly causes the variation of random sampling error, so that the set membership has less influence on the value of dilation factor compared with mode error. In addition, as can be seen from tables 1, 2 and 3, when the set membership is fixed, asFThe closer to the forcing term that generates the true value(s) ((F = 8), the smaller the expansion factor,otherwise, the larger the size. This shows that the larger the mode error is, the greater the force of covariance adjustment is, and the other side also reflects that covariance expansion is an effective method for solving the problem of variance underestimation caused by introduction of various errors. In addition, the value of the standard deviation of the observation error also influences the value of the expansion factor, and under the condition that other settings are the same, the larger the standard of the observation error is, the smaller the coefficient of the expansion factor is.
From the change of localization radius, when the pattern error is fixed, in most cases, the more the number of set members, the maximum of the optimal localization radius is. When in useF = 4,ρ=0.When the number of the set members is increased from 10 to 30 at 5, the optimal localization radius is increased from 10 to 20, and when the number of the set members is increased from 10 to 20F = 8, the optimal localization radius increases from 30 to 36. This means that as the number of members of the set increases, the sampling error becomes smaller and smaller, so that the pseudo-correlation between distant grid points decreases. Meanwhile, the larger the standard deviation of the observation error is, the more the localization radius value is increased. The algorithm of the present invention is aimed at determining the expansion factor and also paying more attention to its effect on RMSE, so in the following tests the settings of the localization radii were fixed using the settings of tables 1, 2 and 3.
TABLE 1 number of different sets of membership with standard deviation of observation error of 0.5 hours, assimilation experiment optimal parameter setting
TABLE 2 number of different sets of membership with standard deviation of observation error of 1.0 hour, assimilation experiment optimal parameter settings
TABLE 3 number of different sets of membership with standard deviation of observation error of 2.0 hour, assimilation experiment optimal parameter settings
The set membership refers to how large set samples are adopted by the ensemble Kalman filtering method to complete estimation of a prediction error covariance matrix, and is an important parameter influencing the performance of the ensemble Kalman filtering method. Generally speaking, the larger the set membership, the better the assimilation result, so different covariance expansion methods should be adopted to conform to the rule. Given the dimension size of the Lorenz-96 pattern, the present invention selected set membership of 10, 20, and 30 samples to test the sensitivity of the present invention to set membership.
FIGS. 2 to 7 show the results of analysis using different dilation factor determination schemes for the set membership ens of 10, 20, and 30, respectively, and for the standard deviation of observation error of 1.0. The former is analyzed in terms of set membership, and the latter is considered in terms of the method of inflation factor. The fixed expansion factor means that the expansion factor is manually adjusted by an exhaustive method, the adjustment range is 1.0-3.16, under the condition that the fixed localization radius and the initial set are the same, 19 groups of assimilation experiments are carried out, and the expansion factor corresponding to the minimum RMSE is determined to be the optimal expansion factor scheme value (green dot transverse line) from the assimilation experiments. This result is used to illustrate the optimal performance that can be achieved by manual empirical tuning under given circumstances, which corresponds to a global optimal solution that is sought over a range of expansion factor values. Because the real value of the expansion factor cannot be obtained, the expansion factor is regarded as a test reference value, the performance of the self-adaptive expansion factor method is expected to be close to the standard, and the closer the performance of the self-adaptive expansion factor method is to the standard, the better the effect of the self-adaptive expansion factor method is considered. Table 4 gives details of the values of RMSE for the fixed dilation method, the Wu method, and the present invention for different set membership and forcing settings with an observed error standard deviation of 1.0.
TABLE 4
From the test results, it can be seen that the scheme (fig. 2,3, 4) with a fixed expansion factor is optimal for the same values of the members of the set, and neither the Wu method nor the present invention can achieve the same performance as the optimal expansion scheme. However, the present invention is closer to the result of the optimal fixed expansion method than the Wu method regardless of the change in the forcing term, which reflects to some extent the effectiveness of the adaptive expansion scheme of the present invention.
In addition, in order to more intuitively understand the degree of performance improvement of the invention relative to the Wu method, the following indexes are used for measurement:
in the formula,RMSE wu for the RMSE value calculated for the Wu method,RMSE sdaci the RMSE values for the present invention. From fig. 8, the analysis performance of the present invention is significantly improved compared to the Wu method. The smaller the set member, the larger the pattern error, and the more obvious the effect. When the number of members in the set is 10, the analysis performance of the method is improved by about 17.08 percent on average compared with the Wu method, and the minimum value appears inF 1.30% for = 8, the maximum occurring whenF When = 12, 24.94% is reached; when the set membership is 20, the analysis performance of the method is averagely improved by about 13.93 percent compared with that of a Wu method, and the value of a forcing term when the minimum value and the maximum value occur is consistent with that when the set membership is 10, the values respectively reach 0.69 percent and 21.83 percent; as the number of members in the set is increased to 30, the improvement degree of the analysis performance of the method is further reduced compared with that of the Wu method, the average value reaches 12.86%, and the minimum value and the maximum value are 0.98% and 20.48% respectively. This shows that the present invention can more effectively adjust variance compared to Wu method, thereby reducing the problem of variance underestimation caused by mode error, and has a very good promoting effect on improving assimilation quality.
From experimental results, the analysis results vary with the mode error no matter how the set membership and the dilation method varyThe trends were approximately the same. As the set membership increases, the analytical performance increases, i.e., the RMSE decreases (FIGS. 2,3, 4). From the perspective of the forced term setting, all schemes areFThe best analysis results were obtained when = 8, and their performance gap was at a minimum at this time. Regardless of the expansion scheme, the analysis performance is gradually improved along with the increase of the number of members of the set, but the performance improvement effect is less obvious when the number of the members of the set is larger. When the number of members in the set is 10, the Wu method and the invention areF In the case of a value of = 8, the analytical RMSE is significantly higher than the optimal solution, but the difference is almost negligible when the number of set members is increased to 20, 30. As the mode error increases, the assimilation performance obtained by different expansion factor methods has obvious difference, namely the RMSE is gradually increased when the difference between the forcing term and 8 is larger. As can be seen from Table 4, when the number of members of the setm = 10, RMSE for Wu method and invention is 0.7327 and 0.6158, respectively, when forced entry is set to 7, RMSE reaches 1.3858 and 1.0402, respectively, when forced entry is increased to 12. In all tests, whenF = 4 orF RMSE reached maximum when = 12. This shows that the dilation factor scheme is very necessary for systems with large mode errors, and the optimal dilation method can bring about improvement of analysis performance under the condition that other parameter settings are the same.
The standard deviation of the observed error is a main parameter for describing the error distribution of the observation entering the assimilation system and is also an important factor influencing the assimilation performance. In order to test and compare the analysis performance of the SDACI method under the condition of different standard deviations of observation errors, the invention tests the analysis performance of the method adopting different expansion factors under the conditions that the set membership is 20 and the standard deviations of the observation errors are 0.5 and 2.0, and the selection method of the expansion factor method is consistent with the previous method.
Fig. 9 and 10 and fig. 11, 12 and 13 show the analysis performance of the three fixed expansion methods under different values of standard deviation of observation error. The overall analytical performance of the fixed expansion method is still the best of the three fixed expansion methods, and secondly, the analytical performance of the Wu method is the worst. Such asAs shown in fig. 9 and 10, the difference between the Wu method and the experimental results of the present invention and the fixed swelling method is significantly reduced as the standard deviation of the observation error increases. Moreover, when the value of the forcing term is 8, i.e., in the case of perfect mode, the Wu method and the analysis value of the present invention are comparable to the fixed expansion factor method, when the standard deviation of the observation error is 2.0,F where = 7, the analytical field RMSE of the SDACI method of the invention is already substantially identical to that of the fixed expansion method, whenF = 9, the analytical results according to the invention exceed the fixed swelling method. The increase of standard deviation of observation errors means the increase of observation errors entering the system, and the analysis performance of the invention in the situation can be so close to that of a fixed expansion method, and the effectiveness of the invention in the aspect of covariance expansion is proved to a certain extent.
The change of standard deviation of observation errors does not affect the performance comparison relationship of the Wu method and the analysis condition of the invention. The analytical performance of the present invention is superior to the Wu method regardless of the change in the forcing term, i.e., the change in the pattern error, when observing the standard deviation of the errorρ= 0.5, the analysis performance of the method is improved by about 12.05 percent compared with the Wu method on averageFThe maximum lifting ratio reached 21.31% when the ratio was 10 (FIG. 14), and the RMSE value obtained by the Wu method was 0.5047 (Table 5) when the RMSE value was 0.6414 according to the analysis of the Wu methodF When the lifting ratio is not less than 5, the lifting ratio is minimum and reaches 3.33%, and the RMSE values of the Wu method and the analysis of the invention respectively reach 0.6628 and 0.6407; standard deviation of observation errorρ= 2.0, the analytical performance improvement of the present invention averaged 6.73% over the Wu method, when the forcing term is appliedF When = 12, the lift ratio reached the maximum of 14.33%. At this point, the Wu method and analytical RMSE values for the present invention are 1.7398 and 1.4904, respectively; the minimum value occurring atF In the case of = 8, 0.88%, and the RMSE value analyzed by the Wu method was 0.8972, the RMSE value obtained by the present invention was 0.8893.
TABLE 5 fixed swell method, Wu method and RMSE results of the invention
From the experimental results of fig. 11, 12 and 13, it can be seen that, according to the three expansion factor methods, as the covariance of the observation errors increases, that is, the observation errors entering the assimilation system increase, no matter what expansion factor scheme is adopted, the RMSE variation trend of the analysis values is the same, that is, the larger the observation errors are, the lower the analysis performance is, which accords with the fundamental characteristics of the assimilation system, and good observation plays a decisive role in the assimilation quality.
Based on the Lorenz-96 mode, a large number of set assimilation experiments under the conditions of different set members, different observation error covariances and different forcing terms show that compared with the similar Wu method, the method provided by the invention considers the spatial distribution of the expansion factors, so that the expansion factors obtained by calculation are more reasonable, and the assimilation performance is obviously improved.
The above embodiment is an embodiment of the present invention, but the embodiment of the present invention is not limited by the above embodiment, and any other changes, modifications, substitutions, combinations, and simplifications which do not depart from the spirit and principle of the present invention should be regarded as equivalent replacements within the protection scope of the present invention.
Claims (3)
1. The data assimilation method based on the adaptive covariance expansion is characterized by comprising the following steps of:
obtaining atmospheric observation values based ont-1, performing mode integration on the analysis value to obtain analysis timetThe forecast field of (1);
according to the forecast field set, estimating to obtaintEnsemble prediction error covariance matrix P for a time instanttEach element of the expansion factor vector is an expansion factor which is a multiplication coefficient in a multiplicative expansion method of the ensemble Kalman filtering method and is used for determining the adjustment degree of variance expansion;
updating the set members of the set Kalman filtering in the initial assimilation process by the expansion factor vectors to increase the variance of a new set Kalman filtering to form new set members, and then performing iterative updating on the new set members by using a set Kalman filtering method to obtain final analysis set members;
using the analysis set member as an initial field to perform pattern prediction;
wherein the analysis time is calculated according to the following formulatThe forecast field of (1):
wherein,is composed oft-1 at the first momenti The number of members of an analysis set,is composed oftAt the first momentiThe members of each of the forecast sets,mas to the number of members of the set,M t-1is thatt-a mode operator at time 1;
calculating an ensemble prediction error covariance matrix according to the following formulaP t :
2. The adaptive covariance based expansion data assimilation method of claim 1, wherein the expansion factor vector is calculated according to the following formula:
constructing an objective function according to the following formula, and solving the minimum value of the objective function to finally obtain an expansion factor vector:
wherein HTIs the transpose of H, dTIs the transpose of d, Tr is the trace of the matrix, R is the diagonal,the background error covariance matrix after dilation:
,for a given statistical prediction error covariance matrix at time t,d is an innovation quantity which represents the deviation between the background field and the observation field in the assimilation process, and H is a tangential operator of a nonlinear observation operator, namely for any state vector x0, the innovation quantity comprises the following components:
3. The adaptive covariance expansion-based data assimilation method of claim 2, wherein the expansion factor vector is usedSubstituting the set analysis member updating process of the set Kalman filtering based on random disturbance to solve each analysis set member:
Priority Applications (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110911181.3A CN113360854B (en) | 2021-08-10 | 2021-08-10 | Data assimilation method based on adaptive covariance expansion |
Applications Claiming Priority (1)
Application Number | Priority Date | Filing Date | Title |
---|---|---|---|
CN202110911181.3A CN113360854B (en) | 2021-08-10 | 2021-08-10 | Data assimilation method based on adaptive covariance expansion |
Publications (2)
Publication Number | Publication Date |
---|---|
CN113360854A CN113360854A (en) | 2021-09-07 |
CN113360854B true CN113360854B (en) | 2021-11-05 |
Family
ID=77540920
Family Applications (1)
Application Number | Title | Priority Date | Filing Date |
---|---|---|---|
CN202110911181.3A Active CN113360854B (en) | 2021-08-10 | 2021-08-10 | Data assimilation method based on adaptive covariance expansion |
Country Status (1)
Country | Link |
---|---|
CN (1) | CN113360854B (en) |
Families Citing this family (4)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN114047563A (en) * | 2021-10-18 | 2022-02-15 | 中国人民解放军国防科技大学 | All-weather assimilation method of infrared hyperspectrum |
CN114048433B (en) * | 2021-10-26 | 2022-06-21 | 南京大学 | Mixed assimilation system and method based on ensemble Kalman filtering framework |
CN114070262B (en) * | 2021-10-26 | 2022-06-21 | 南京大学 | Additional disturbance integrated hybrid set Kalman filtering weather forecast assimilation method and device |
CN116429721A (en) * | 2023-03-15 | 2023-07-14 | 东北农业大学 | Pine nut near infrared spectrum band selection method based on MPA framework |
Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106909983A (en) * | 2017-01-03 | 2017-06-30 | 北京国能日新***控制技术有限公司 | Based on many meteorological sources ultra-short term wind speed forecasting methods of Kalman filter and device |
CN108646319A (en) * | 2018-08-09 | 2018-10-12 | 深圳市雅码科技有限公司 | A kind of heavy showers forecasting procedure and system in short-term |
CN109033029A (en) * | 2018-06-29 | 2018-12-18 | 中国人民解放军国防科技大学 | Method and system for generating deterministic analysis set based on reverse localization |
CN110020462A (en) * | 2019-03-07 | 2019-07-16 | 江苏无线电厂有限公司 | The method that a kind of pair of meteorological data carries out fusion treatment and generate numerical weather forecast |
CN111948736A (en) * | 2019-05-14 | 2020-11-17 | 中国电力科学研究院有限公司 | High-dimensional weather forecast data dimension reduction method based on big data platform |
CN112558188A (en) * | 2021-01-22 | 2021-03-26 | 兰州大学 | Method for improving strong convection forecast by assimilating lightning data |
Family Cites Families (3)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN109001845B (en) * | 2016-07-20 | 2019-12-10 | 中国水利水电科学研究院 | rainfall forecasting method |
EP3438700B1 (en) * | 2017-08-02 | 2020-11-18 | Thales Management & Services Deutschland GmbH | Method for operating a gnss receiver, with gnss signal deselection |
KR101898629B1 (en) * | 2017-09-12 | 2018-09-13 | 서강대학교산학협력단 | Method and device for post-correction of forecasted meteorological parameters by using a H-Infinity filter |
-
2021
- 2021-08-10 CN CN202110911181.3A patent/CN113360854B/en active Active
Patent Citations (6)
Publication number | Priority date | Publication date | Assignee | Title |
---|---|---|---|---|
CN106909983A (en) * | 2017-01-03 | 2017-06-30 | 北京国能日新***控制技术有限公司 | Based on many meteorological sources ultra-short term wind speed forecasting methods of Kalman filter and device |
CN109033029A (en) * | 2018-06-29 | 2018-12-18 | 中国人民解放军国防科技大学 | Method and system for generating deterministic analysis set based on reverse localization |
CN108646319A (en) * | 2018-08-09 | 2018-10-12 | 深圳市雅码科技有限公司 | A kind of heavy showers forecasting procedure and system in short-term |
CN110020462A (en) * | 2019-03-07 | 2019-07-16 | 江苏无线电厂有限公司 | The method that a kind of pair of meteorological data carries out fusion treatment and generate numerical weather forecast |
CN111948736A (en) * | 2019-05-14 | 2020-11-17 | 中国电力科学研究院有限公司 | High-dimensional weather forecast data dimension reduction method based on big data platform |
CN112558188A (en) * | 2021-01-22 | 2021-03-26 | 兰州大学 | Method for improving strong convection forecast by assimilating lightning data |
Non-Patent Citations (1)
Title |
---|
集合资料同化方法在强雷暴天气预报中的应用;张云济,张福青;《气象科技进展》;20180615;全文 * |
Also Published As
Publication number | Publication date |
---|---|
CN113360854A (en) | 2021-09-07 |
Similar Documents
Publication | Publication Date | Title |
---|---|---|
CN113360854B (en) | Data assimilation method based on adaptive covariance expansion | |
Jordan et al. | Evaluating probabilistic forecasts with scoringRules | |
Berkes et al. | Detecting changes in the mean of functional observations | |
CN108304679A (en) | A kind of adaptive reliability analysis method | |
CN112085254B (en) | Prediction method and model based on multi-fractal cooperative measurement gating circulation unit | |
CN110442911B (en) | High-dimensional complex system uncertainty analysis method based on statistical machine learning | |
CN111931983A (en) | Precipitation prediction method and system | |
CN108073782A (en) | A kind of data assimilation method based on the equal weight particle filter of observation window | |
Khalili et al. | An efficient statistical approach to multi-site downscaling of daily precipitation series in the context of climate change | |
CN113139605A (en) | Power load prediction method based on principal component analysis and LSTM neural network | |
Ikbal et al. | Estimating Weibull parameters using maximum likelihood estimation and ordinary least squares: Simulation study and application on meteorological data | |
CN114117852B (en) | Regional heat load rolling prediction method based on finite difference working domain division | |
CN106528679A (en) | Time series analysis method based on multilinear autoregression model | |
Nguyen | Continuous observation hidden Markov model | |
CN111950198A (en) | Ground clutter simulation method based on neural network | |
Zhang et al. | Reference priors for the generalized extreme value distribution | |
Tripathi et al. | Does the Type of Records Affect the Estimates of the Parameters? | |
Li et al. | Bayesian Filtering for High-Dimensional State-Space Models With State Partition and Error Compensation | |
Hao et al. | The Improved Adaptive LASSO Method and Applications in Autoregressive Model | |
CN113810024B (en) | Cost reference particle filtering method based on mixed probability selection operator | |
Ma et al. | Analysis and simulation of the influencing factors on regional water use based on information entropy | |
Poonia et al. | On bivariate modelling using a new statistical distribution with homogeneous marginals | |
Kumar et al. | Data Assimilation applied to a numerical 2D cylinder wake flow | |
Seki et al. | Regression analysis on high-dimensional, block diagonal structure data with focus on latent variables | |
CN114417708A (en) | Slope monitoring design optimization method |
Legal Events
Date | Code | Title | Description |
---|---|---|---|
PB01 | Publication | ||
PB01 | Publication | ||
SE01 | Entry into force of request for substantive examination | ||
SE01 | Entry into force of request for substantive examination | ||
GR01 | Patent grant | ||
GR01 | Patent grant |