A multivariable frequency domain assimilation method and system

Abstract

The present invention relates to a multivariate frequency domain assimilation method and system, comprising the following steps: 1) determining a study area, an assimilation window length, correction vectors, and assimilation vectors, and preparing observation products of a land surface process model and its driving data, an observation operator and its parameter data, and assimilation vectors; 2) generating a set of correction vectors for the land surface process model; 3) forward integrating the land surface process model and the observation operator using the set of correction vectors to obtain a set of forecast correction vectors and a set of forecast assimilation vectors for each integration moment within the assimilation window; 4) assimilating the correction vectors in the frequency domain using the forecast set of correction vectors and the forecast set of assimilation vectors with the observation products of the assimilated vectors; and 5) performing assimilation at subsequent moments in time according to steps 2-4 based on the assimilated correction vectors. To overcome the time domain limitations of sequential and continuous assimilation methods, the present invention utilizes multiple surface state variable observation products and the land surface process model to enable frequency domain assimilation.

Description

Technical Field

[0001] The present invention relates to a multivariate data assimilation method and system, and in particular to a multivariate frequency domain assimilation method and system. Background Art

[0002] Relying on its inherent physical processes and dynamic mechanisms, land surface process models can continuously simulate surface state variables related to water and energy cycles in time and space, such as surface temperature, soil moisture, and soil temperature at the surface-atmosphere interface. They are an important method for studying large-scale and basin-scale water and energy cycles. However, the quality of input data and the incompleteness of model parameterization schemes of land surface process models result in low simulation accuracy.

[0003] Ground-based observations can measure surface variables continuously over time, but the number of observation points is typically small and spatially representative, making it difficult to dynamically monitor surface state variables over large areas. Land surface data assimilation methods, by organically combining the strengths of ground-based observation products and land surface model simulations, can effectively improve the accuracy of land surface process models in simulating surface state variables, providing a new approach for accurately and continuously simulating regional water and heat fluxes.

[0004] Currently, two main assimilation methods are commonly used, both domestically and internationally: sequential assimilation and continuous assimilation. Sequential assimilation uses observational products to update the state variables of the land surface process model based on the uncertainty of the observational products and the land surface process model simulation at the integration moment, thereby sequentially obtaining the optimal a posteriori estimate of the land surface process model state variables. Continuous assimilation constructs an objective function within the assimilation window, calculating the optimal state variables and the deviation between the observed variables and the initial state variables, and then uses an iterative optimization algorithm to solve for the optimal state variables.

[0005] Both sequential and continuous assimilation seek the optimal solution for the state variables of land surface process models within the time domain. However, the evolution of surface process states typically involves physical information across multiple timescales. State variables with high spatial uniformity, such as wind speed, vary slowly, while state variables with low spatial uniformity, such as surface temperature, vary rapidly. This results in the time series information of different surface state variables containing different frequency components. Sequential and continuous assimilation cannot account for these differences in the properties of state variables, and assimilation efficiency urgently needs to be improved.

[0006] To overcome the limitations of current sequential assimilation and continuous assimilation methods in the time domain, a multivariable frequency domain assimilation method and system was invented using observation products of multiple surface state variables and land surface process models. Summary of the Invention

[0007] The technical problem to be solved by the present invention is to address the deficiencies of the existing technology and provide a multivariable frequency domain assimilation method and system.

[0008] The technical solutions of the present invention are as follows:

[0009] A multivariate frequency domain assimilation method comprises the following steps:

[0010] Step 1: Determine the study area, assimilation window length, correction vector and assimilation vector, and prepare observation products including land surface process model and its driving data, observation operator and its parameter data, and assimilation vector;

[0011] Step 2: Generate a set of correction vectors for the land surface process model;

[0012] Step 3: Use the correction vector set to forward integrate the land surface process model and the observation operator to obtain the forecast set of the correction vector and the forecast set of the assimilation vector at each integration moment in the assimilation window;

[0013] Step 4: Use the forecast ensemble of the correction vector and the forecast ensemble of the assimilated vector to assimilate the correction vector in the frequency domain with the aid of the observation product of the assimilated vector;

[0014] Step 5: Based on the assimilated correction vector, perform assimilation at subsequent moments according to Steps 2 to 4.

[0015] In the multivariate frequency domain assimilation method, step 2 generates a set of correction vectors using the historical mean vector and historical variance vector corresponding to the assimilation time of the correction vector. The specific steps are:

[0016] Step 2011: Generate N×M random variables u1 uniformly distributed on (0, 1) nm and u2 nm , 1≤n≤N, 1≤m≤M, where n, m represent the subscript variables of the random variables, N is the dimension of the correction vector, and M is the size of the correction vector set;

[0017] Step 2012: Using N×M u1 nm and u2 nm Generate N×M random variables U that obey the standard normal distribution nm , the formula is as follows:

[0018]

[0019] Step 2013: Using N×M U nm Generate a random correction vector set {X Nm} 1≤m≤M :

[0020]

[0021] Where ′ represents the vector or matrix transpose, [μ1,…,μ N]′ and [σ1,…,σ N ]′ are the correction vectors X N =[x1,…,x n ,…,x N ]′ corresponds to the historical mean vector and historical variance vector at the assimilation moment.

[0022] The multivariate frequency domain assimilation method, wherein the specific steps of step 3 are:

[0023] Step 3011: Set {X Nm} 1≤m≤M Each member X Nm Input the land surface process model and forward integrate under the driving data to generate the forecast results of the correction vector at each integration time in the assimilation window. The integration formula is:

[0024]

[0025] Where t represents the current assimilation time, w represents the number of integration steps in the assimilation window, W represents the length of the assimilation window, and F t , F t+1 ,…,F t+w denote the driving data at the integration time t, t+1, ..., t+w, М denotes the land surface process model, M(X Nm ,F t ) means to convert X Nm Enter М and press F t Driven by one step forward, X t+wNm =[x t+w1m ,…,x t+wnm ,…,x t+wNm ]′ means to convert X Nm Enter М and press F t , F t+1 ,…, F t+w The prediction result of forward integration w steps under driving, x t+w1m ,…,x t+wnm ,…,x t+wNm Represents X t+wNm The forecast results of the 1st…, nth,…, Nth dimensional correction variables.

[0026] Step 3012: X in step 3011 t+wNm The forecast set X that constitutes the correction vector TNM :

[0027]

[0028] Step 3013: X TNM Each member in is input into the observation operator and combined with the observation operator parameter data to obtain the forecast set of the assimilated vector. The formula is:

[0029]

[0030] Where H represents the observation operator, L represents the dimension of the assimilation vector, and P t+w represents the parameter data of the observation operator at time t+w, H(X t+wNm ,P t+w ) means to convert X t+wNm Enter H and combine with P t+w Get the forecast result Y of the assimilation vector at time t+w t+wLm =[y t+w1m ,…,y t+wlm ,…,y t+wLm ]′,y t+w1m ,…,y t+wlm ,…,y t+wLm represents the forecast results of the 1st, lth, ..., Lth dimensional assimilated variables, Γ t+wLM Y t+wL1 ,…,Y t+wLM The forecast ensemble composed of TLM Because Γ tLM ,…,Γ t+WLM The set of assimilated vector forecast ensembles.

[0031] The multivariate frequency domain assimilation method, wherein the specific steps of step 4 are:

[0032] Step 4011: Generate an assimilated vector observation product set using the historical mean vector and historical variance vector of the assimilated vector observation product at each integration moment in the corresponding assimilation window. The specific generation process is as follows:

[0033] a) Generate (W+1)×L×M random variables v1 uniformly distributed on (0, 1) wlm and v2 wlm , w, l, m represent the subscript variables of random variables;

[0034] b) Using v1 wlm and v2 wlm Generate (W+1)×L×M random variables V that obey the standard normal distribution wlm , the formula is as follows:

[0035]

[0036] c) Using V wlm Generate a set O of W+1 assimilated vector observation products of size M that obey the L-dimensional normal distribution TLM :

[0037]

[0038] Where, δ t+wl , ν t+wl are the historical variance and historical mean of the l-th dimension assimilated variable in the assimilation vector at time t+w, respectively. t+wlm is the mth random variable of the generated lth dimension assimilated variable;

[0039] Step 4012: Calculate the Kalman gain matrix at each integration moment in the assimilation window using the correction vector forecast set, the assimilated vector forecast set, and the observation product set. The specific calculation process is as follows:

[0040] Utilize X TNM and Y TLM Calculate the covariance matrix between the correction vector and the assimilation vector at each integration moment in the assimilation window. The calculation formula is:

[0041]

[0042] Where, P t+wn H t+wi is the forecast covariance between the n-th dimension correction variable and the i-th dimension assimilation variable at time t+w, is the forecast covariance matrix between the N-dimensional correction vector and the L-dimensional assimilation vector at time t+w, P TN H TL is the set of forecast covariance matrices between the N-dimensional correction vector and the L-dimensional assimilation vector;

[0043] Take advantage of Y TLM Calculate the covariance matrix of the assimilation vector forecast results at each integration moment. The calculation formula is:

[0044]

[0045] Where H t+wi H t+wl is the forecast covariance between the assimilated variables of the i-th and l-th dimensions at time t+w, is the forecast covariance matrix between L assimilated variables at time t+w, H TL H TL is the set of forecast covariance matrices between assimilated variables;

[0046] Using O TLM Calculate the covariance matrix of the assimilation vector observation product at each integration moment in the assimilation window. The calculation formula is:

[0047]

[0048] Where, O t+wi O t+wl is the observed covariance between the assimilated variables of the i-th dimension and the l-th dimension at time t+w, is the observation covariance matrix between the assimilated variables at time t+w, O TL O TL is the set of observation covariance matrices;

[0049] d) Using P TN H TL 、H TL H TL and O TL O TL Calculate the Kalman gain matrix at each integration moment in the assimilation window The calculation formula is as follows:

[0050]

[0051] Step 4013: Concatenate the mean vector of the assimilated vector forecast set at each integration moment in the assimilation window and the assimilated vector observation product into a time series signal. The concatenation process is as follows:

[0052] a) Set Γ tLM ,…,Γ t+WLM The forecast mean vector of is connected in series to form a time series signal. The connection method is:

[0053]

[0054] Where Y tl ,…,Y t+wl ,…,Y t+Wl are the predicted mean values ​​of the l-th dimension assimilated variables at each integration moment, Y Wl is the time series signal of the predicted mean value of the l-th dimension assimilation variable in the assimilation window;

[0055] b) Assimilate the vector observation product O at each integration moment tL =[o t1 ,…,o tl ,…,o tL ]′,…,O t+wL =[o t+w1 ,…,o t+wl ,…,o t+wL ]′,…,O t+WL =[o t+W1 ,…,o t+Wl ,…,o t+WL ]′ are connected in series to form a timing signal. The connection method is:

[0056]

[0057] In the formula, o Wl is the time series signal composed of the observed values ​​of the l-th dimension assimilated variables, O Wl is the time series signal matrix of the assimilated vector observation product;

[0058] Step 4014: concatenate the Kalman gain matrices at each integration moment within the assimilation window into a Kalman gain matrix time series signal. The concatenation method is:

[0059]

[0060] Where g Wnl The Kalman gain matrix G t+wnl , a timing signal consisting of the elements in the nth row and lth column of the set t≤t+w≤t+W;

[0061] Step 4015: Perform frequency decomposition on the time series signal of the mean vector of the assimilated vector forecast set and the time series signal of the assimilated variable observation product to obtain their respective frequency components. The frequency decomposition process is as follows:

[0062] a) Change Y Wl Perform frequency decomposition to obtain frequency components. The frequency decomposition method is as follows:

[0063]

[0064] Where Ry Wl [j] is Y Wl The cosine harmonic coefficient value with intermediate frequency j, Iy Wl [j] is Y Wl The sinusoidal harmonic coefficient value with intermediate frequency j;

[0065] b) Change o Wl Perform frequency decomposition to obtain frequency components. The frequency decomposition method is as follows:

[0066]

[0067] Where, Ro Wl [j] is o Wl The cosine harmonic coefficient value with intermediate frequency j, Io Wl [j] is o Wl The sinusoidal harmonic coefficient value with intermediate frequency j;

[0068] Step 4016: g Wnl Perform frequency decomposition to obtain the frequency components of the time series signals of each matrix element of the Kalman gain matrix. The frequency decomposition method is as follows:

[0069]

[0070] Where Rg Wnl [j] is g Wnl The cosine harmonic coefficient value with intermediate frequency j, Ig Wnl [j] is g Wnl The sinusoidal harmonic coefficient values ​​with intermediate frequency j.

[0071] Step 4017: Calculate the correction amount of the correction vector at different frequency components, and synthesize the correction amount of the correction vector at all frequency components into the correction amount of the correction vector at each integration time in the assimilation window. The synthesis process is as follows:

[0072] a) Calculate the correction amount of each dimension of the correction vector at different frequency components. The calculation method is:

[0073]

[0074] Where Rκ Wn [j] is the correction amount of the cosine harmonic with frequency j in the n-th dimension correction variable time series signal, Iκ Wn [j] is the correction amount of the sinusoidal harmonic with frequency j in the n-th dimension correction variable time series signal;

[0075] b) The correction amount of each dimensional correction variable on all frequency components 0≤j≤W / 2 is synthesized into the correction amount of each dimensional correction variable at each integration time in the assimilation window. The synthesis formula is as follows:

[0076]

[0077] Where, ξ t+wn is the correction amount of the n-th dimension correction variable at time t+w;

[0078] Step 4018: Calculate the mean vector of the correction vector at each integration moment in the assimilation window. The calculation method is:

[0079]

[0080] Where, ξ n is the mean value of the correction amount of the n-th dimension correction variable at all integration moments, N is the mean vector of the correction amounts of all N-dimensional correction variables;

[0081] Step 4019: Correct the correction vector using the mean vector of the correction amount of the correction vector, specifically:

[0082] X tN =X tN +Ξ N

[0083] Where, X tN is the value of the correction vector at assimilation time t and the assimilated value after correction.

[0084] The technical solution of the present invention to solve the above technical problems is as follows: a system for a multivariable frequency domain assimilation method, comprising an assimilation preparation module, a correction vector set module, an assimilation vector forecast module, and an assimilation module;

[0085] The assimilation preparation module is used to determine the study area, assimilation window length, correction vector and assimilation vector, and prepare observation products of land surface process model and its driving data, observation operator and its parameter data, and assimilation vector;

[0086] The correction vector set module is used to generate a correction vector set of the land surface process model;

[0087] The assimilation vector forecast module is used to use the correction vector set to forward integrate the land surface process model and the observation operator to obtain the forecast set of the correction vector and the forecast set of the assimilation vector at each integration moment in the assimilation window;

[0088] The assimilation module is used to realize assimilation of the correction vector in the frequency domain by using the forecast set of the correction vector and the forecast set of the assimilation vector with the aid of the observation product of the assimilation vector.

[0089] The system, wherein the correction vector set module includes a uniform distribution random variable generation unit, a standard normal distribution random variable generation unit and a normal distribution random variable generation unit;

[0090] The uniformly distributed random variable generating unit is used to generate N×M random variables u1 uniformly distributed on (0, 1) nm and u2 nm , 1≤n≤N, 1≤m≤M, where n, m represent the subscript variables of the random variables, N is the dimension of the correction vector, and M is the size of the correction vector set;

[0091] The standard normal distribution random variable generating unit is used to use N×M u1 nm and u2 nm Generate N×M random variables U that obey the standard normal distribution nm , the formula is as follows:

[0092]

[0093] The normal distribution random variable generating unit is used to utilize N×M U nm Generate a random correction vector set {X Nm} 1≤m≤M :

[0094]

[0095] Where ′ represents the vector or matrix transpose, [μ1,…,μ N ]′ and [σ1,…,σ N ]′ are the correction vectors X N =[x1,…,x n,…,x N ]′ corresponds to the historical mean vector and historical variance vector at the assimilation moment.

[0096] The system, wherein the assimilation vector prediction module includes a correction vector set integration unit, a correction vector prediction set unit and an assimilation vector prediction set unit;

[0097] The correction vector set integration unit is used to convert {X Nm} 1≤m≤M Each member X Nm Input the land surface process model and forward integrate under the driving data to generate the forecast results of the correction vector at each integration time in the assimilation window. The integration formula is:

[0098]

[0099] Where t represents the current assimilation time, w represents the number of integration steps in the assimilation window, W represents the length of the assimilation window, and F t , F t+1 ,…,F t+w denote the driving data at the integration time t, t+1, ..., t+w, М denotes the land surface process model, M(X Nm ,F t ) means to convert X Nm Enter М and press F t Driven by one step forward, X t+wNm =[x t+w1m ,…,x t+wnm ,…,x t+wNm ]′ means to convert X Nm Enter М and press F t , F t+1 ,…, F t+w The prediction result of forward integration w steps under driving, x t+w1m ,…,x t+wnm ,…,x t+wNm Represents X t+wNm The forecast results of the 1st…, nth,…, Nth dimensional correction variables.

[0100] The correction vector prediction set unit is used to generate the X t+wNm Construct the correction vector forecast set X TNM :

[0101]

[0102] The assimilation vector forecasting set unit is used to convert X TNM Each member in is input into the observation operator and combined with the observation operator parameter data to obtain the forecast set of the assimilated vector. The formula is:

[0103]

[0104] Where H represents the observation operator, L represents the dimension of the assimilation vector, and P t+w represents the parameter data of the observation operator at time t+w, H(X t+wNm ,P t+w ) means to convert X t+wNm Enter H and combine with P t+w Get the forecast result Y of the assimilation vector at time t+w t+wLm =[y t+w1m ,…,y t+wlm ,…,y t+wLm ]′,y t+w1m ,…,y t+wlm ,…,y t+wLm represents the forecast results of the 1st, lth, ..., Lth dimensional assimilated variables, Γ t+wLM Y t+wL1 ,…,Y t+wLM The forecast ensemble composed of TLM Because Γ tLM ,…,Γ t+WLM The set of assimilated vector forecast ensembles.

[0105] The system, wherein the assimilation module includes an observation product set generation unit, a Kalman gain matrix calculation unit, a vector signal series connection unit, a matrix signal series connection unit, a vector signal frequency decomposition unit, a matrix signal frequency decomposition unit, a correction amount calculation unit, a correction amount mean calculation unit, and a correction unit;

[0106] The observation product set generation unit is used to generate an assimilated vector observation product set using the historical mean vector and historical variance vector of the assimilated vector observation product at each integration moment in the corresponding assimilation window. The specific steps are as follows:

[0107] a) Generate (W+1)×L×M random variables v1 uniformly distributed on (0, 1) wlm and v2 wlm , w, l, m represent the subscript variables of random variables;

[0108] b) Using v1 wlm and v2 wlm Generate (W+1)×L×M random variables V that obey the standard normal distribution wlm , the formula is as follows:

[0109]

[0110] c) Using V wlm Generate a set O of W+1 assimilated vector observation products of size M that obey the L-dimensional normal distribution TLM :

[0111]

[0112] Where, δ t+wl , ν t+wl are the historical variance and historical mean of the l-th dimension assimilated variable in the assimilation vector at time t+w, respectively. t+wlm is the mth random variable of the generated lth dimension assimilated variable;

[0113] The Kalman gain matrix calculation unit is used to calculate the Kalman gain matrix at each integration moment in the assimilation window using the correction vector forecast set, the assimilation vector forecast set, and the observation product set. The specific steps are as follows:

[0114] a) Using X TNM and Y TLM Calculate the covariance matrix between the correction vector and the assimilation vector at each integration moment in the assimilation window. The calculation formula is:

[0115]

[0116] Where, P t+wn H t+wi is the forecast covariance between the n-th dimension correction variable and the i-th dimension assimilation variable at time t+w, is the forecast covariance matrix between the N-dimensional correction vector and the L-dimensional assimilation vector at time t+w, P TN H TL is the set of forecast covariance matrices between the N-dimensional correction vector and the L-dimensional assimilation vector;

[0117] b) Using Y TLM Calculate the covariance matrix of the assimilation vector forecast results at each integration moment. The calculation formula is:

[0118]

[0119] Where H t+wi H t+wl is the forecast covariance between the assimilated variables of the i-th and l-th dimensions at time t+w, is the forecast covariance matrix between L assimilated variables at time t+w, H TL H TL is the set of forecast covariance matrices between assimilated variables;

[0120] c) Using O TLM Calculate the covariance matrix of the assimilation vector observation product at each integration moment in the assimilation window. The calculation formula is:

[0121]

[0122] Where, O t+wi O t+wl is the observed covariance between the assimilated variables of the i-th dimension and the l-th dimension at time t+w, is the observation covariance matrix between the assimilated variables at time t+w, O TL O TL is the set of observation covariance matrices;

[0123] d) Using P TN H TL 、H TL H TL and O TL O TL Calculate the Kalman gain matrix at each integration moment in the assimilation window The calculation formula is as follows:

[0124]

[0125] The vector signal concatenation unit is used to concatenate the mean vector of the assimilated vector forecast set at each integration moment in the assimilation window and the assimilated vector observation product into a time series signal. The concatenation process is as follows:

[0126] a) Set Γ tLM ,…,Γ t+WLM The forecast mean vector of is connected in series to form a time series signal. The connection method is:

[0127]

[0128] Where Y tl ,…,Y t+wl ,…,Y t+Wl are the predicted mean values ​​of the l-th dimension assimilated variables at each integration moment, Y Wl is the time series signal of the predicted mean value of the l-th dimension assimilation variable in the assimilation window;

[0129] b) Assimilate the vector observation product O at each integration moment tL =[o t1 ,…,o tl ,…,o tL ]′,…,O t+wL =[o t+w1 ,…,o t+wl ,…,o t+wL ]′,…,O t+WL =[o t+W1 ,…,o t+Wl ,…,o t+WL ]′ are connected in series to form a timing signal. The connection method is:

[0130]

[0131] In the formula, oWl is the time series signal composed of the observed values ​​of the l-th dimension assimilated variables, O Wl is the time series signal matrix of the assimilated vector observation product;

[0132] The matrix signal series connection unit is used to connect the Kalman gain matrices at each integration moment in the assimilation window in series to form a Kalman gain matrix time series signal. The series connection method is:

[0133] g Wnl =[G tnl ,...,G t+wnl ,...,G t+Wnl ] 1≤n≤N 1≤l≤L

[0134] Where g Wnl The Kalman gain matrix G t+wnl , a timing signal consisting of the elements in the nth row and lth column of the set t≤t+w≤t+W;

[0135] The vector signal frequency decomposition unit is used to perform frequency decomposition on the time series signal of the mean vector of the assimilated vector forecast set and the time series signal of the assimilated variable observation product to obtain their respective frequency components. The frequency decomposition process is as follows:

[0136] a) Change Y Wl Perform frequency decomposition to obtain frequency components. The frequency decomposition method is as follows:

[0137]

[0138] Where Ry Wl [j] is Y Wl The cosine harmonic coefficient value with intermediate frequency j, Iy Wl [j] is Y Wl The sinusoidal harmonic coefficient value with intermediate frequency j;

[0139] b) Change o Wl Perform frequency decomposition to obtain frequency components. The frequency decomposition method is as follows:

[0140]

[0141] Where, Ro Wl [j] is o Wl The cosine harmonic coefficient value with intermediate frequency j, Io Wl [j] is o Wl The sinusoidal harmonic coefficient value with intermediate frequency j;

[0142] The matrix signal frequency decomposition unit is used to decompose g WnlPerform frequency decomposition to obtain the frequency components of the time series signals of each matrix element of the Kalman gain matrix. The frequency decomposition method is as follows:

[0143]

[0144] Where Rg Wnl [j] is g Wnl The cosine harmonic coefficient value with intermediate frequency j, Ig Wnl [j] is g Wnl The sinusoidal harmonic coefficient value with intermediate frequency j;

[0145] The correction amount calculation unit is used to calculate the correction amount of the correction vector on different frequency components and synthesize the correction amounts of the correction vector on all frequency components into the correction amount of the correction vector at each integration time in the assimilation window. The synthesis process is as follows:

[0146] a) Calculate the correction amount of each dimension of the correction vector at different frequency components. The calculation method is:

[0147]

[0148] Where Rκ Wn [j] is the correction amount of the cosine harmonic with frequency j in the n-th dimension correction variable time series signal, Iκ Wn [j] is the correction amount of the sinusoidal harmonic with frequency j in the n-th dimension correction variable time series signal;

[0149] The correction amount of each dimensional correction variable on all frequency components 0≤j≤W / 2 is synthesized into the correction amount of each dimensional correction variable at each integration time in the assimilation window. The synthesis formula is as follows:

[0150]

[0151] Where, ξ t+wn is the correction amount of the n-th dimension correction variable at time t+w;

[0152] The correction amount mean calculation unit is used to calculate the mean vector of the correction amount of the correction vector at each integration moment in the assimilation window, and the calculation method is:

[0153]

[0154] Where, ξ n is the mean value of the correction amount of the n-th dimension correction variable at all integration moments, N is the mean vector of the correction amounts of all N-dimensional correction variables;

[0155] The correction unit is used to correct the correction vector using the mean vector of the correction amount of the correction vector, and the correction method is:

[0156] X tN =X tN +Ξ N

[0157] Where, X tN is the value of the correction vector at assimilation time t and the assimilated value after correction.

[0158] The beneficial effects of the present invention are as follows: compared with traditional time domain data assimilation methods, the present invention strengthens the contribution of different frequency components in the assimilated variable observation products to the optimization of land surface process model state variables through frequency decomposition of correction variables and assimilation variable time series signals, improves the assimilation effect of multivariable assimilation due to the use of a single correction factor, and has broad application prospects in improving the simulation effects of ecological models, hydrological models, atmospheric models, etc. BRIEF DESCRIPTION OF THE DRAWINGS

[0159] Figure 1 This is a flow chart of a multivariate frequency domain assimilation method according to the present invention;

[0160] Figure 2 This is a flowchart for the specific implementation of step 2 of the present invention;

[0161] Figure 3 This is a flowchart for the specific implementation of step 3 of the present invention;

[0162] Figure 4 This is a flowchart for the specific implementation of step 4 of the present invention;

[0163] Figure 5 This is a flowchart for the specific implementation of the multivariate frequency domain assimilation method of the present invention;

[0164] Figure 6 This is a system block diagram of a multivariable frequency domain assimilation method according to the present invention;

[0165] Figure 7 This is a structural block diagram of the correction vector set module of the present invention;

[0166] Figure 8 This is a structural block diagram of the assimilation vector prediction module of the present invention;

[0167] Figure 9 This is a structural block diagram of the assimilation module of the present invention;

[0168] In the accompanying drawings, the components represented by the reference numerals are as follows:

[0169] 1. Assimilation preparation module, 2. Correction vector set module, 3. Assimilation vector forecast module, 4. Assimilation module; 201. Uniform distribution random variable generation unit, 202. Standard normal distribution random variable generation unit, 203. Normal distribution random variable generation unit; 301. Correction vector set integration unit, 302. Correction vector forecast set unit, 303. Assimilation vector forecast set unit; 401. Observation product set generation unit, 402. Kalman gain matrix calculation unit, 403. Vector signal concatenation unit, 404. Matrix signal concatenation unit, 405. Vector signal frequency decomposition unit, 406. Matrix signal frequency decomposition unit, 407. Correction amount calculation unit, 408. Correction amount mean calculation unit, 409. Correction unit. DETAILED DESCRIPTION

[0170] The principles and features of the present invention are described in detail below with reference to the accompanying drawings. The examples given are only used to explain the present invention and are not used to limit the scope of the present invention.

[0171] like Figure 1 As shown in FIG, a multivariate frequency domain assimilation method includes the following steps:

[0172] Step 1: Determine the study area, assimilation window length, correction vector and assimilation vector, and prepare observation products including land surface process model and its driving data, observation operator and its parameter data, and assimilation vector;

[0173] Step 2: Generate a set of correction vectors for the land surface process model;

[0174] Step 3: Use the correction vector set to forward integrate the land surface process model and the observation operator to obtain the forecast set of the correction vector and the forecast set of the assimilation vector at each integration time in the assimilation window;

[0175] Step 4: Use the forecast ensemble of the correction vector and the forecast ensemble of the assimilated vector to assimilate the correction vector in the frequency domain with the aid of the observation product of the assimilated vector;

[0176] Step 5: Based on the assimilated correction vector, perform assimilation at subsequent moments according to Steps 2 to 4.

[0177] The correction vector of the land surface process model should be such that the improvement of its simulation accuracy can improve the prediction accuracy of the correction vector of the land surface process model at subsequent moments. The land surface process model and observation operator describe the physical connection between the correction vector and the assimilation vector from different perspectives. The clm2, hyssib, mosaic, Commonly used land surface process models such as noahmp.4.0.1 can be used. Global driving data required for land surface process model integration can be downloaded from http: / / hydrology.princeton.edu / data / pgf / v2 / 0.5deg / 3hourly / , or from http: / / www.ncdc.ac.cn / portal / . Observational products of assimilated variables can be downloaded from https: / / www.crensed.ac.cn / auth / login or https: / / fluxnet.fluxdata.org / . The observation operator is typically a forward model in which the assimilation vector is derived from the correction vector.

[0178] like Figure 2 As shown, step 2 uses the historical mean vector and historical variance vector of the correction vector corresponding to the assimilation time to generate a correction vector set. The specific steps are:

[0179] Step 2011: Generate N×M random variables u1 uniformly distributed on (0, 1) nm and u2 nm , 1≤n≤N, 1≤m≤M, where n, m represent the subscript variables of the random variables, N is the dimension of the correction vector, and M is the size of the correction vector set;

[0180] Step 2012: Using N×M u1 nm and u2 nm Generate N×M random variables U that obey the standard normal distribution nm , the formula is as follows:

[0181]

[0182] Step 2013: Using N×M U nm Generate a random correction vector set {X Nm} 1≤m≤M :

[0183]

[0184] Where ′ represents the vector or matrix transpose, [μ1,…,μ N ]′ and [σ1,…,σ N ]′ are the correction vectors X N =[x1,…,x n ,…,x N ]′ corresponds to the historical mean vector and historical variance vector at the assimilation moment.

[0185] like Figure 3 As shown, the specific steps of step 3 are:

[0186] Step 3011: Set {X Nm} 1≤m≤M Each member X Nm Input the land surface process model and forward integrate under the driving data to generate the forecast results of the correction vector at each integration time in the assimilation window. The integration formula is:

[0187]

[0188] Where t represents the current assimilation time, w represents the number of integration steps in the assimilation window, W represents the length of the assimilation window, and F t , F t+1 ,…,F t+w denote the driving data at the integration time t, t+1, ..., t+w, М denotes the land surface process model, M(X Nm ,F t ) means to convert X Nm Enter М and press F t Driven by one step forward, X t+wNm =[x t+w1m ,…,x t+wnm ,…,x t+wNm ]′ means to convert X Nm Enter М and press F t , F t+1 ,…, F t+w The prediction result of forward integration w steps under driving, x t+w1m ,…,x t+wnm ,…,x t+wNm Represents X t+wNm The forecast results of the 1st…, nth,…, Nth dimensional correction variables.

[0189] Step 3012: X in step 3011 t+wNm Construct the correction vector forecast set X TNM :

[0190]

[0191] Step 3013: X TNMEach member in is input into the observation operator and combined with the observation operator parameter data to obtain the forecast set of the assimilated vector. The formula is:

[0192]

[0193] Where H represents the observation operator, L represents the dimension of the assimilation vector, and P t+w represents the parameter data of the observation operator at time t+w, H(X t+wNm ,P t+w ) means to convert X t+wNm Enter H and combine with P t+w Get the forecast result Y of the assimilation vector at time t+w t+wLm =[y t+w1m ,…,y t+wlm ,…,y t+wLm ]′,y t+w1m ,…,y t+wlm ,…,y t+wLm represents the forecast results of the 1st, lth, ..., Lth dimensional assimilated variables, Γ t+wLM Y t+wL1 ,…,Y t+wLM The forecast ensemble composed of TLM Because Γ tLM ,…,Γ t+WLM The set of assimilated vector forecast ensembles.

[0194] like Figure 4 As shown, the specific steps of step 4 are:

[0195] Step 4011: Generate an assimilated vector observation product set using the historical mean vector and historical variance vector of the assimilated vector observation product at each integration moment in the corresponding assimilation window. The specific generation process is as follows:

[0196] a) Generate (W+1)×L×M random variables v1 uniformly distributed on (0, 1) wlm and v2 wlm , w, l, m represent the subscript variables of random variables;

[0197] b) Using v1 wlm and v2 wlm Generate (W+1)×L×M random variables V that obey the standard normal distribution wlm , the formula is as follows:

[0198]

[0199] c) Using V wlm Generate a set O of W+1 assimilated vector observation products of size M that obey the L-dimensional normal distribution TLM :

[0200]

[0201] Where, δ t+wl , ν t+wl are the historical variance and historical mean of the l-th dimension assimilated variable in the assimilation vector at time t+w, respectively. t+wlm is the mth random variable of the generated lth dimension assimilated variable;

[0202] Step 4012: Calculate the Kalman gain matrix at each integration moment in the assimilation window using the correction vector forecast set, the assimilated vector forecast set, and the observation product set. The specific calculation process is as follows:

[0203] a) Using X TNM and Y TLM Calculate the covariance matrix between the correction vector and the assimilation vector at each integration moment in the assimilation window. The calculation formula is:

[0204]

[0205] Where, P t+wn H t+wi is the forecast covariance between the n-th dimension correction variable and the i-th dimension assimilation variable at time t+w, is the forecast covariance matrix between the N-dimensional correction vector and the L-dimensional assimilation vector at time t+w, P TN H TL is the set of forecast covariance matrices between the N-dimensional correction vector and the L-dimensional assimilation vector;

[0206] b) Using Y TLM Calculate the covariance matrix of the assimilation vector forecast results at each integration moment. The calculation formula is:

[0207]

[0208] Where H t+wi H t+wl is the forecast covariance between the assimilated variables of the i-th and l-th dimensions at time t+w, is the forecast covariance matrix between L assimilated variables at time t+w, H TL H TL is the forecast ensemble of the forecast covariance matrix between the assimilated variables;

[0209] Using O TLM Calculate the covariance matrix of the assimilation vector observation product at each integration moment in the assimilation window. The calculation formula is:

[0210]

[0211] Where, O t+wi Ot+wl is the observed covariance between the assimilated variables of the i-th dimension and the l-th dimension at time t+w, is the observation covariance matrix between the assimilated variables at time t+w, O TL O TL is the set of observation covariance matrices;

[0212] d) Using P TN H TL 、H TL H TL and O TL O TL Calculate the Kalman gain matrix at each integration moment in the assimilation window The calculation formula is as follows:

[0213]

[0214] Step 4013: Concatenate the mean vector of the assimilated vector forecast set at each integration moment in the assimilation window and the assimilated vector observation product into a time series signal. The concatenation process is as follows:

[0215] a) Set Γ tLM ,…,Γ t+WLM The forecast mean vector of is connected in series to form a time series signal. The connection method is:

[0216]

[0217] Where Y tl ,…,Y t+wl ,…,Y t+Wl are the predicted mean values ​​of the l-th dimension assimilated variables at each integration moment, Y Wl is the time series signal of the predicted mean value of the l-th dimension assimilation variable in the assimilation window;

[0218] b) Assimilate the vector observation product O at each integration moment tL =[o t1 ,…,o tl ,…,o tL ]′,…,O t+wL =[o t+w1 ,…,o t+wl ,…,o t+wL ]′,…,O t+WL =[o t+W1 ,…,o t+Wl ,…,o t+WL ]′ are connected in series to form a timing signal. The connection method is:

[0219]

[0220] In the formula, o Wlis the time series signal composed of the observed values ​​of the l-th dimension assimilated variables, O Wl is the time series signal matrix of the assimilated vector observation product;

[0221] Step 4014: concatenate the Kalman gain matrices at each integration moment within the assimilation window into a Kalman gain matrix time series signal. The concatenation method is:

[0222] g Wnl =[G tnl ,...,G t+wnl ,...,G t+Wnl ]1≤n≤N 1≤l≤L

[0223] Where g Wnl The Kalman gain matrix G t+wnl , a timing signal consisting of the elements in the nth row and lth column of the set t≤t+w≤t+W;

[0224] Step 4015 is as follows Figure 5 As shown in Figure 2, the time series signal of the mean vector of the assimilated vector forecast ensemble and the time series signal of the assimilated variable observation product are frequency decomposed to obtain their respective frequency components. The frequency decomposition process is as follows:

[0225] a) Change Y Wl Perform frequency decomposition to obtain frequency components. The frequency decomposition method is as follows:

[0226]

[0227] Where Ry Wl [j] is Y Wl The cosine harmonic coefficient value with intermediate frequency j, Iy Wl [j] is Y Wl The sinusoidal harmonic coefficient value with intermediate frequency j;

[0228] b) Change o Wl Perform frequency decomposition to obtain frequency components. The frequency decomposition method is as follows:

[0229]

[0230] Where, Ro Wl [j] is o Wl The cosine harmonic coefficient value with intermediate frequency j, Io Wl [j] is o Wl The sinusoidal harmonic coefficient value with intermediate frequency j;

[0231] Step 4016 Figure 5 As shown: g Wnl Perform frequency decomposition to obtain the frequency components of the time series signals of each matrix element of the Kalman gain matrix. The frequency decomposition method is as follows:

[0232]

[0233] Where Rg Wnl [j] is g Wnl The cosine harmonic coefficient value with intermediate frequency j, Ig Wnl [j] is g Wnl The sinusoidal harmonic coefficient values ​​with intermediate frequency j.

[0234] Step 4017 Figure 5 As shown in the figure: the correction amount of the correction vector is calculated on different frequency components, and the correction amount of the correction vector on all frequency components is synthesized into the correction amount of the correction vector at each integration time in the assimilation window. The synthesis process is as follows:

[0235] a) Calculate the correction amount of each dimension of the correction vector at different frequency components. The calculation method is:

[0236]

[0237] Where Rκ Wn [j] is the correction amount of the cosine harmonic with frequency j in the n-th dimension correction variable time series signal, Iκ Wn [j] is the correction amount of the sinusoidal harmonic with frequency j in the n-th dimension correction variable time series signal;

[0238] b) The correction amount of each dimensional correction variable on all frequency components 0≤j≤W / 2 is synthesized into the correction amount of each dimensional correction variable at each integration time in the assimilation window. The synthesis formula is as follows:

[0239]

[0240] Where, ξ t+wn is the correction amount of the n-th dimension correction variable at time t+w;

[0241] Step 4018 Figure 5 As shown in Figure 2, the mean vector of the correction vector at each integration moment in the assimilation window is calculated using the following method:

[0242]

[0243] Where, ξ n is the mean value of the correction amount of the n-th dimension correction variable at all integration moments, N is the mean vector of the correction amounts of all N-dimensional correction variables;

[0244] Step 4019 Figure 5 As shown: the correction vector is corrected using the mean vector of the correction amount of the correction vector, specifically:

[0245] XtN =X tN +Ξ N

[0246] Where, X tN is the value of the correction vector at assimilation time t and the assimilated value after correction.

[0247] like Figure 6 As shown, a system for implementing a multivariable frequency domain assimilation method includes an assimilation preparation module 1, a correction vector set module 2, an assimilation vector prediction module 3, and an assimilation module 4;

[0248] The assimilation preparation module 1 is used to determine the study area, assimilation window length, correction vector and assimilation vector, and prepare observation products of land surface process model and its driving data, observation operator and its parameter data, and assimilation vector;

[0249] The correction vector set module 2 is used to generate a correction vector set of the land surface process model;

[0250] The assimilation vector forecast module 3 is used to use the correction vector set to forward integrate the land surface process model and the observation operator to obtain the forecast set of the correction vector and the forecast set of the assimilation vector at each integration moment in the assimilation window;

[0251] The assimilation module 4 is configured to assimilate the correction vector in the frequency domain using the forecast set of the correction vector and the forecast set of the assimilated vector with the aid of the observation product of the assimilated vector.

[0252] like Figure 7 As shown, the correction vector set module includes a uniform distribution random variable generation unit 201, a standard normal distribution random variable generation unit 202 and a normal distribution random variable generation unit 203;

[0253] The uniformly distributed random variable generating unit 201 is used to generate N×M random variables u1 uniformly distributed on (0, 1) nm and u2 nm , 1≤n≤N, 1≤m≤M, where n, m represent the subscript variables of the random variables, N is the dimension of the correction vector, and M is the size of the correction vector set;

[0254] The standard normal distribution random variable generating unit 202 is used to use N×M u1 nm and u2 nm Generate N×M random variables U that obey the standard normal distribution nm , the formula is as follows:

[0255]

[0256] The normal distribution random variable generating unit 203 is used to use N×M U nm Generate a random correction vector set {X Nm} 1≤m≤M :

[0257]

[0258] Where ′ represents the vector or matrix transpose, [μ1,…,μ N ]′ and [σ1,…,σ N ]′ are the correction vectors X N =[x1,…,x n ,…,x N ]′ corresponds to the historical mean vector and historical variance vector at the assimilation moment.

[0259] like Figure 8 As shown, the assimilation vector prediction module includes a correction vector set integration unit 301, a correction vector prediction set unit 302 and an assimilation vector prediction set unit 303;

[0260] The correction vector set integration unit 301 is used to convert {X Nm} 1≤m≤M Each member X Nm Input the land surface process model and forward integrate under the driving data to generate the forecast results of the correction vector at each integration time in the assimilation window. The integration formula is:

[0261]

[0262] Where t represents the current assimilation time, w represents the number of integration steps in the assimilation window, W represents the length of the assimilation window, and F t , F t+1 ,…,F t+w denote the driving data at the integration time t, t+1, ..., t+w, М denotes the land surface process model, M(X Nm ,F t ) means to convert X Nm Enter М and press F t Driven by one step forward, X t+wNm =[x t+w1m ,…,x t+wnm ,…,x t+wNm ]′ means to convert X Nm Enter М and press F t , F t+1 ,…, F t+w The prediction result of forward integration w steps under driving, x t+w1m ,…,x t+wnm ,…,xt+wNm Represents X t+wNm The forecast results of the 1st…, nth,…, Nth dimensional correction variables.

[0263] The correction vector prediction collection unit 302 is used to t+wNm Construct the correction vector forecast set X TNM :

[0264]

[0265] The assimilation vector forecast collection unit 303 is used to convert X TNM Each member in is input into the observation operator and combined with the observation operator parameter data to obtain the forecast set of the assimilated vector. The formula is:

[0266]

[0267] Where H represents the observation operator, L represents the dimension of the assimilation vector, and P t+w represents the parameter data of the observation operator at time t+w, H(X t+wNm ,P t+w ) means to convert X t+wNm Enter H and combine with P t+w Get the forecast result Y of the assimilation vector at time t+w t+wLm =[y t+w1m ,…,y t+wlm ,…,y t+wLm ]′,y t+w1m ,…,y t+wlm ,…,y t+wLm represents the forecast results of the 1st, lth, ..., Lth dimensional assimilated variables, Γ t+wLM Y t+wL1 ,…,Y t+wLM The forecast ensemble composed of TLM Because Γ tLM ,…,Γ t+WLM The set of assimilated vector forecast ensembles.

[0268] like Figure 9 As shown, the assimilation module includes an observation product set generation unit 401, a Kalman gain matrix calculation unit 402, a vector signal series connection unit 403, a matrix signal series connection unit 404, a vector signal frequency decomposition unit 405, a matrix signal frequency decomposition unit 406, a correction amount calculation unit 407, a correction amount mean calculation unit 408, and a correction unit 409;

[0269] The observation product set generation unit 401 is used to generate an assimilated vector observation product set using the historical mean vector and historical variance vector of the assimilated vector observation products at each integration moment in the corresponding assimilation window. The specific steps are as follows:

[0270] a) Generate (W+1)×L×M random variables v1 uniformly distributed on (0, 1) wlm and v2 wlm , w, l, m represent the subscript variables of random variables;

[0271] b) Using v1 wlm and v2 wlm Generate (W+1)×L×M random variables V that obey the standard normal distribution wlm , the formula is as follows:

[0272]

[0273] c) Using V wlm Generate a set O of W+1 assimilated vector observation products of size M that obey the L-dimensional normal distribution TLM :

[0274]

[0275] Where, δ t+wl , ν t+wl are the historical variance and historical mean of the l-th dimension assimilated variable in the assimilation vector at time t+w, respectively. t+wlm is the mth random variable of the generated lth dimension assimilated variable;

[0276] The Kalman gain matrix calculation unit 402 is used to calculate the Kalman gain matrix at each integration time in the assimilation window using the correction vector forecast set, the assimilation vector forecast set, and the observation product set. The specific steps are as follows:

[0277] a) Using X TNM and Y TLM Calculate the covariance matrix between the correction vector and the assimilation vector at each integration moment in the assimilation window. The calculation formula is:

[0278]

[0279] Where, P t+wn H t+wi is the forecast covariance between the n-th dimension correction variable and the i-th dimension assimilation variable at time t+w, is the forecast covariance matrix between the N-dimensional correction vector and the L-dimensional assimilation vector at time t+w, P TN H TL is the set of forecast covariance matrices between the N-dimensional correction vector and the L-dimensional assimilation vector;

[0280] b) Using Y TLM Calculate the covariance matrix of the assimilation vector forecast results at each integration moment. The calculation formula is:

[0281]

[0282] Where H t+wi H t+wl is the forecast covariance between the assimilated variables of the i-th and l-th dimensions at time t+w, is the forecast covariance matrix between L assimilated variables at time t+w, H TL H TL is the forecast ensemble of the forecast covariance matrix between the assimilated variables;

[0283] c) Using O TLM Calculate the covariance matrix of the assimilation vector observation product at each integration moment in the assimilation window. The calculation formula is:

[0284]

[0285] Where, O t+wi O t+wl is the observed covariance between the assimilated variables of the i-th dimension and the l-th dimension at time t+w, is the observation covariance matrix between the assimilated variables at time t+w, O TL O TL is the set of observation covariance matrices;

[0286] d) Using P TN H TL 、H TL H TL and O TL O TL Calculate the Kalman gain matrix at each integration moment in the assimilation window The calculation formula is as follows:

[0287]

[0288] The vector signal concatenation unit 403 is used to concatenate the mean vector of the assimilated vector forecast set at each integration moment in the assimilation window and the assimilated vector observation product into a time series signal. The concatenation process is as follows:

[0289] a) Set Γ tLM ,…,Γ t+WLM The forecast mean vector of is connected in series to form a time series signal. The connection method is:

[0290]

[0291] Where Y tl ,…,Y t+wl ,…,Y t+Wl are the predicted mean values ​​of the l-th dimension assimilated variables at each integration moment, Y Wlis the time series signal of the predicted mean value of the l-th dimension assimilation variable in the assimilation window;

[0292] b) Assimilate the vector observation product O at each integration moment tL =[o t1 ,…,o tl ,…,o tL ]′,…,O t+wL =[o t+w1 ,…,o t+wl ,…,o t+wL ]′,…,O t+WL =[o t+W1 ,…,o t+Wl ,…,o t+WL ]′ are connected in series to form a timing signal. The connection method is:

[0293]

[0294] In the formula, o Wl is the time series signal composed of the observed values ​​of the l-th dimension assimilated variables, O Wl is the time series signal matrix of the assimilated vector observation product;

[0295] The matrix signal series connection unit 404 is used to connect the Kalman gain matrices at each integration moment in the assimilation window in series to form a Kalman gain matrix time series signal. The connection method is:

[0296] g Wnl =[G tnl ,...,G t+wnl ,...,G t+Wnl ]1≤n≤N 1≤l≤L

[0297] Where g Wnl The Kalman gain matrix G t+wnl , a timing signal consisting of the elements in the nth row and lth column of the set t≤t+w≤t+W;

[0298] The vector signal frequency decomposition unit 405 is used to perform frequency decomposition on the time series signal of the mean vector of the assimilated vector forecast set and the time series signal of the assimilated variable observation product to obtain their respective frequency components. The frequency decomposition process is as follows:

[0299] a) Change Y Wl Perform frequency decomposition to obtain frequency components. The frequency decomposition method is as follows:

[0300]

[0301] Where Ry Wl [j] is Y Wl The cosine harmonic coefficient value with intermediate frequency j, Iy Wl[j] is Y Wl The sinusoidal harmonic coefficient value with intermediate frequency j;

[0302] b) Change o Wl Perform frequency decomposition to obtain frequency components. The frequency decomposition method is as follows:

[0303]

[0304] Where, Ro Wl [j] is o Wl The cosine harmonic coefficient value with intermediate frequency j, Io Wl [j] is o Wl The sinusoidal harmonic coefficient value with intermediate frequency j;

[0305] The matrix signal frequency decomposition unit 406 is used to decompose g Wnl Perform frequency decomposition to obtain the frequency components of the time series signals of each matrix element of the Kalman gain matrix. The frequency decomposition method is as follows:

[0306]

[0307] Where Rg Wnl [j] is g Wnl The cosine harmonic coefficient value with intermediate frequency j, Ig Wnl [j] is g Wnl The sinusoidal harmonic coefficient value with intermediate frequency j;

[0308] The correction amount calculation unit 407 is used to calculate the correction amount of the correction vector at different frequency components and synthesize the correction amounts of the correction vector at all frequency components into the correction amount of the correction vector at each integration time in the assimilation window. The synthesis process is as follows:

[0309] a) Calculate the correction amount of each dimension of the correction vector at different frequency components. The calculation method is:

[0310]

[0311] Where Rκ Wn [j] is the correction amount of the cosine harmonic with frequency j in the n-th dimension correction variable time series signal, Iκ Wn [j] is the correction amount of the sinusoidal harmonic with frequency j in the n-th dimension correction variable time series signal;

[0312] b) The correction amount of each dimensional correction variable on all frequency components 0≤j≤W / 2 is synthesized into the correction amount of each dimensional correction variable at each integration time in the assimilation window. The synthesis formula is as follows:

[0313]

[0314] Where, ξ t+wn is the correction amount of the n-th dimension correction variable at time t+w;

[0315] The correction amount mean calculation unit 408 is used to calculate the mean vector of the correction amount of the correction vector at each integration moment in the assimilation window. The calculation method is:

[0316]

[0317] Where, ξ n is the mean value of the correction amount of the n-th dimension correction variable at all integration moments, N is the mean vector of the correction amounts of all N-dimensional correction variables;

[0318] The correction unit 409 is used to correct the correction vector using the mean vector of the correction amount of the correction vector, and the correction method is:

[0319] X tN =X tN +Ξ N

[0320] Where, X tN is the value of the correction vector at assimilation time t and the assimilated value after correction.

[0321] The above embodiments are only used to illustrate the present invention rather than to limit the present invention. Any supplement, modification or equivalent replacement of the technical solution of the present invention by those skilled in the art shall fall within the scope of protection of the appended claims of the present invention.

Claims

1. A multivariate frequency domain assimilation method, characterized in that: The following steps are involved: Step 1: Determine the study area, assimilation window length, correction vector and assimilation vector, and prepare observation products including land surface process model and its driving data, observation operator and its parameter data, and assimilation vector; Step 2: Generate a set of correction vectors for the land surface process model; Step 3: Use the correction vector set to forward integrate the land surface process model and observation operator to obtain the forecast set of correction vectors and the forecast set of assimilation vectors at each integration time in the assimilation window, specifically: Step 3011: Set {X Nm } 1≤m≤M Each member X Nm Input the land surface process model and forward integrate under the driving data to generate the forecast results of the correction vector at each integration time in the assimilation window. The integration formula is: Where t represents the current assimilation time, w represents the number of integration steps in the assimilation window, W represents the length of the assimilation window, and F t , F t+1 ,…,F t+w represent the driving data at the integration time t, t+1, ..., t+w, M represents the land surface process model, M(X Nm ,F t ) means to convert X Nm Enter M and enter F t Driven by one step forward, X t+wNm =[x t+w1m ,…,x t+wnm ,…,x t+wNm ]′ means to convert X Nm Enter M and enter F t , F t+1 ,…,F t+w The prediction result of forward integration w steps under driving, x t+w1m ,…,x t+wnm ,…,x t+wNm Represents X t+wNm The forecast results of the 1st…, nth,…, Nth dimensional correction variables; Step 3012: X in step 3011 t+wNm Construct the correction vector forecast set X TNM : x t+wNM ={X t+wN1 ,…,X t+wNm ,…,X t+wNM } X TNM ={χ tNM ,…,x t+wNM …,x t+WNM } Step 3013: X TNM Each member in is input into the observation operator and combined with the observation operator parameter data to obtain the forecast set of assimilation vectors. The formula is: Y t+wLm =H(X t+wNm ,P t+w )0≤w≤W 1≤m≤M Γ t+wLM ={And t+wL1 ,…,AND t+wLM } Y TLM ={Γ tLM ,…,G t+WLM } Where H represents the observation operator, L represents the dimension of the assimilation vector, and P t+w represents the parameter data of the observation operator at time t+w, H(X t+wNm ,P t+w ) means to convert X t+wNm Enter H and combine with P t+w Get the forecast result Y of the assimilation vector at time t+w t+wLm =[y t+w1m ,…,y t+wlm ,…,y t+wLm ]′,y t+w1m ,…,y t+wlm ,…,y t+wLm represents the forecast results of the 1st, lth, ..., Lth dimensional assimilated variables, Γ t+wLM Y t+wL1 ,…,Y t+wLM The forecast ensemble composed of TLM Because Γ tLM ,…,Γ t+WLM The set of assimilated vector forecast ensembles composed of; Step 4: Use the forecast ensemble of the correction vector and the forecast ensemble of the assimilated vector to assimilate the correction vector in the frequency domain with the aid of the observation product of the assimilated vector; Step 5: Based on the assimilated correction vector, perform assimilation at subsequent moments according to Steps 2 to 4.

2. A multivariate frequency domain assimilation method according to claim 1, characterized in that: The step 2 generates a set of correction vectors using the historical mean vector and historical variance vector corresponding to the assimilation time of the correction vector. The specific steps are: Step 2011: Generate N×M random variables u1 uniformly distributed on (0, 1) nm and u2 nm , 1≤n≤N, 1≤m≤M, where n, m represent the subscript variables of the random variables, N is the dimension of the correction vector, and M is the size of the correction vector set; Step 2012: Using N×M u1 nm and u2 nm Generate N×M random variables U that obey the standard normal distribution nm , the formula is as follows: Step 2013: Using N×M U nm Generate a random correction vector set {X Nm } 1≤m≤M : X Nm =[x 1m ,…,x Nm ]′ x 1m =U 1m ×σ1+μ1,…,x Nm =U Nm ×s N +m N Where ′ represents the vector or matrix transpose, [μ1,…,μ N ]′ and [σ1,…,σ N ]′ are the correction vectors X N =[x1,…,x n ,…,x N ]′ corresponds to the historical mean vector and historical variance vector at the assimilation moment.

3. The multivariate frequency domain assimilation method according to claim 1, wherein: The specific steps of step 4 are: Step 4011: Generate an assimilated vector observation product set using the historical mean vector and historical variance vector of the assimilated vector observation product at each integration moment in the corresponding assimilation window. The specific generation process is as follows: a) Generate (W+1)×L×M random variables v1 uniformly distributed on (0, 1) wlm and v2 wlm , w, l, m represent the subscript variables of random variables; b) Using v1 wlm and v2 wlm Generate (W+1)×L×M random variables V that obey the standard normal distribution wlm , the formula is as follows: c) Using V wlm Generate a set O of W+1 assimilated vector observation products of size M that obey the L-dimensional normal distribution TLM : about t+wlm =V wlm ×δ t+wl +v t+wl Where, δ t+wl , ν t+wl are the historical variance and historical mean of the l-th dimension assimilated variable in the assimilation vector at time t+w, respectively. t+wlm is the mth random variable of the generated lth dimension assimilated variable; Step 4012: Calculate the Kalman gain matrix at each integration moment in the assimilation window using the correction vector forecast set, the assimilated vector forecast set, and the observation product set. The specific calculation process is as follows: a) Using X TNM and Y TLM Calculate the covariance matrix between the correction vector and the assimilation vector at each integration moment in the assimilation window. The calculation formula is: Where, P t+wn H t+wi is the forecast covariance between the n-th dimension correction variable and the i-th dimension assimilation variable at time t+w, is the forecast covariance matrix between the N-dimensional correction vector and the L-dimensional assimilation vector at time t+w, P TN H TL is the set of forecast covariance matrices between the N-dimensional correction vector and the L-dimensional assimilation vector; b) Using Y TLM Calculate the covariance matrix of the assimilation vector forecast results at each integration moment. The calculation formula is: Where H t+wi H t+wl is the forecast covariance between the assimilated variables of the i-th and l-th dimensions at time t+w, is the forecast covariance matrix between L assimilated variables at time t+w, H TL H TL is the forecast ensemble of the forecast covariance matrix between the assimilated variables; c) Using O TLM Calculate the covariance matrix of the assimilation vector observation product at each integration moment in the assimilation window. The calculation formula is: Where, O t+wi O t+wl is the observed covariance between the assimilated variables of the i-th dimension and the l-th dimension at time t+w, is the observation covariance matrix between the assimilated variables at time t+w, O TL O TL is the set of observation covariance matrices; d) Using P TN H TL 、H TL H TL and O TL O TL Calculate the Kalman gain matrix G at each integration moment in the assimilation window t+wNL , the calculation formula is as follows: O t+wi H t+wl =O t+wi O t+wl +H t+wi H t+wl Step 4013: Concatenate the mean vector of the assimilated vector forecast set at each integration moment in the assimilation window and the assimilated vector observation product into a time series signal. The concatenation process is as follows: a) Set Γ tLM ,…,Γ t+WLM The forecast mean vector of is connected in series to form a time series signal. The connection method is: AND Wl =[And tl ,...,AND t+wl ,...,AND t+Wl ] Where Y tl ,…,Y t+wl ,…,Y t+Wl are the predicted mean values ​​of the l-th dimension assimilated variables at each integration moment, Y Wl is the time series signal of the predicted mean value of the l-th dimension assimilation variable in the assimilation window; b) Assimilate the vector observation product O at each integration moment tL =[o t1 ,…,o tl ,…,o tL ]′,…,O t+wL =[o t+w1 ,…,o t+wl ,…,o t+wL ]′,…,O t+WL =[o t+W1 ,…,o t+Wl ,…,o t+WL ]′ are connected in series to form a timing signal. The connection method is: the Wl =[o tl ,...,o t+wl ,...,o t+Wl ]1≤l≤L The WL =[o W1 ,...,o Wl ,...,o WL ]′ In the formula, o Wl is the time series signal composed of the observed values ​​of the l-th dimension assimilated variables, O Wl is the time series signal matrix of the assimilated vector observation product; Step 4014: concatenate the Kalman gain matrices at each integration moment within the assimilation window into a Kalman gain matrix time series signal. The concatenation method is: g Wnl =[G tnl ,…,G t+wnl ,…,G t+Wnl ]1≤n≤N 1≤l≤L Where g Wnl The Kalman gain matrix G t+wnl , a timing signal consisting of the elements in the nth row and lth column of the set t≤t+w≤t+W; Step 4015: Perform frequency decomposition on the time series signal of the mean vector of the assimilated vector forecast set and the time series signal of the assimilated variable observation product to obtain their respective frequency components. The frequency decomposition process is as follows: a) Change Y Wl Perform frequency decomposition to obtain frequency components. The frequency decomposition method is as follows: Where Ry Wl [j] is Y Wl The cosine harmonic coefficient value with intermediate frequency j, Iy Wl [j] is Y Wl The sinusoidal harmonic coefficient value with intermediate frequency j; b) Change o Wl Perform frequency decomposition to obtain frequency components. The frequency decomposition method is as follows: Where, Ro Wl [j] is o Wl The cosine harmonic coefficient value with intermediate frequency j, Io Wl [j] is o Wl The sinusoidal harmonic coefficient value with intermediate frequency j; Step 4016: g Wnl Perform frequency decomposition to obtain the frequency components of the time series signals of each matrix element of the Kalman gain matrix. The frequency decomposition method is as follows: Where Rg Wnl [j] is g Wnl The cosine harmonic coefficient value with intermediate frequency j, Ig Wnl [j] is g Wnl The sinusoidal harmonic coefficient value with intermediate frequency j; Step 4017: Calculate the correction amount of the correction vector at different frequency components, and synthesize the correction amount of the correction vector at all frequency components into the correction amount of the correction vector at each integration time in the assimilation window. The synthesis process is as follows: a) Calculate the correction amount of each dimension of the correction vector at different frequency components. The calculation method is: Where Rκ Wn [j] is the correction amount of the cosine harmonic with frequency j in the n-th dimension correction variable time series signal, Iκ Wn [j] is the correction amount of the sinusoidal harmonic with frequency j in the n-th dimension correction variable time series signal; b) The correction amount of each dimensional correction variable on all frequency components 0≤j≤W / 2 is synthesized into the correction amount of each dimensional correction variable at each integration time in the assimilation window. The synthesis formula is as follows: Where, ξ t+wn is the correction amount of the n-th dimension correction variable at time t+w; Step 4018: Calculate the mean vector of the correction vector at each integration moment in the assimilation window. The calculation method is: X N =[ξ1,…,ξ n ,...,x N ]′ Where, ξ n is the mean value of the correction amount of the n-th dimension correction variable at all integration moments, N is the mean vector of the correction amounts of all N-dimensional correction variables; Step 4019: Correct the correction vector using the mean vector of the correction amount of the correction vector, specifically: X tN =X tN +Ξ N Where, X tN is the value of the correction vector at assimilation time t and the assimilated value after correction.

4. A system for implementing a multivariate frequency domain assimilation method according to any one of claims 1 to 3, characterized in that: It includes assimilation preparation module, correction vector collection module, assimilation vector forecast module and assimilation module; The assimilation preparation module is used to determine the study area, assimilation window length, correction vector and assimilation vector, and prepare observation products of land surface process model and its driving data, observation operator and its parameter data, and assimilation vector; The correction vector set module is used to generate a correction vector set of the land surface process model; The assimilation vector forecast module is used to use the correction vector set to forward integrate the land surface process model and the observation operator to obtain the forecast set of the correction vector and the forecast set of the assimilation vector at each integration moment in the assimilation window; The assimilation module is used to realize assimilation of the correction vector in the frequency domain by using the forecast set of the correction vector and the forecast set of the assimilation vector with the aid of the observation product of the assimilation vector.

5. The system according to claim 4, wherein: The correction vector set module includes a uniform distribution random variable generation unit, a standard normal distribution random variable generation unit and a normal distribution random variable generation unit; The uniformly distributed random variable generating unit is used to generate N×M random variables u1 uniformly distributed on (0, 1) nm and u2 nm , 1≤n≤N, 1≤m≤M, where n, m represent the subscript variables of the random variables, N is the dimension of the correction vector, and M is the size of the correction vector set; The standard normal distribution random variable generating unit is used to use N×M u1 nm and u2 nm Generate N×M random variables U that obey the standard normal distribution nm , the formula is as follows: The normal distribution random variable generating unit is used to utilize N×M U nm Generate a random correction vector set {X Nm } 1≤m≤M : X Nm =[x 1m ,…,x Nm ]′ x 1m =U 1m ×σ1+μ1,…,x Nm =U Nm ×s N +m N Where ′ represents the vector or matrix transpose, [μ1,…,μ N ]′ and [σ1,…,σ N ]′ are the correction vectors X N =[x1,…,x n ,…,x N ]′ corresponds to the historical mean vector and historical variance vector at the assimilation moment.

6. The system according to claim 4, wherein: The assimilation vector prediction module includes a correction vector set integration unit, a correction vector prediction set unit and an assimilation vector prediction set unit; The correction vector set integration unit is used to convert {X Nm } 1≤m≤M Each member X Nm Input the land surface process model and forward integrate under the driving data to generate the forecast results of the correction vector at each integration time in the assimilation window. The integration formula is: Where t represents the current assimilation time, w represents the number of integration steps in the assimilation window, W represents the length of the assimilation window, and F t , F t+1 ,…,F t+w represent the driving data at the integration time t, t+1, ..., t+w, M represents the land surface process model, M(X Nm ,F t ) means to convert X Nm Enter M and enter F t Driven by one step forward, X t+wNm =[x t+w1m ,…,x t+wnm ,…,x t+wNm ]′ means to convert X Nm Enter M and enter F t , F t+1 ,…,F t+w The prediction result of forward integration w steps under driving, x t+w1m ,…,x t+wnm ,…,x t+wNm Represents X t+wNm The forecast results of the 1st…, nth,…, Nth dimensional correction variables; The correction vector prediction set unit is used to generate the X t+wNm Construct the correction vector forecast set X TNM : x t+wNM ={X t+wN1 ,…,X t+wNm ,…,X t+wNM } X TNM ={χ tNM ,…,x t+wNM …,x t+WNM } The assimilation vector forecasting set unit is used to convert X TNM Each member in is input into the observation operator and combined with the observation operator parameter data to obtain the forecast set of the assimilated vector. The formula is: Y t+wLm =H(X t+wNm ,P t+w )0≤w≤W 1≤m≤M Γ t+wLM ={And t+wL1 ,…,AND t+wLM } Y TLM ={Γ tLM ,…,G t+WLM } Where H represents the observation operator, L represents the dimension of the assimilation vector, and P t+w represents the parameter data of the observation operator at time t+w, H(X t+wNm ,P t+w ) means to convert X t+wNm Enter H and combine with P t+w Get the forecast result Y of the assimilation vector at time t+w t+wLm =[y t+w1m ,…,y t+wlm ,…,y t+wLm ]′,y t+w1m ,…,y t+wlm ,…,y t+wLm represents the forecast results of the 1st, lth, ..., Lth dimensional assimilated variables, Γ t+wLM Y t+wL1 ,…,Y t+wLM The forecast ensemble composed of TLM Because Γ tLM ,…,Γ t+WLM The set of assimilated vector forecast ensembles.

7. The system according to claim 4, wherein: The assimilation module includes an observation product set generation unit, a Kalman gain matrix calculation unit, a vector signal series connection unit, a matrix signal series connection unit, a vector signal frequency decomposition unit, a matrix signal frequency decomposition unit, a correction amount calculation unit, a correction amount mean calculation unit, and a correction unit; The observation product set generation unit is used to generate an assimilated vector observation product set using the historical mean vector and historical variance vector of the assimilated vector observation product at each integration moment in the corresponding assimilation window. The specific steps are as follows: a) Generate (W+1)×L×M random variables v1 uniformly distributed on (0, 1) wlm and v2 wlm , w, l, m represent the subscript variables of random variables; b) Using v1 wlm and v2 wlm Generate (W+1)×L×M random variables V that obey the standard normal distribution wlm , the formula is as follows: c) Using V wlm Generate a set O of W+1 assimilated vector observation products of size M that obey the L-dimensional normal distribution TLM : the t+wlm =V wlm ×d t+wl +n t+wl Where, δ t+wl , ν t+wl are the historical variance and historical mean of the l-th dimension assimilated variable in the assimilation vector at time t+w, respectively. t+wlm is the mth random variable of the generated lth dimension assimilated variable; The Kalman gain matrix calculation unit is used to calculate the Kalman gain matrix at each integration moment in the assimilation window using the correction vector forecast set, the assimilation vector forecast set, and the observation product set. The specific steps are as follows: a) Using X TNM and Y TLM Calculate the covariance matrix between the correction vector and the assimilation vector at each integration moment in the assimilation window. The calculation formula is: Where, P t+wn H t+wi is the forecast covariance between the n-th dimension correction variable and the i-th dimension assimilation variable at time t+w, is the forecast covariance matrix between the N-dimensional correction vector and the L-dimensional assimilation vector at time t+w, P TN H TL is the set of forecast covariance matrices between the N-dimensional correction vector and the L-dimensional assimilation vector; b) Using Y TLM Calculate the covariance matrix of the assimilation vector forecast results at each integration moment. The calculation formula is: Where H t+wi H t+wl is the forecast covariance between the assimilated variables of the i-th and l-th dimensions at time t+w, is the forecast covariance matrix between L assimilated variables at time t+w, H TL H TL is the forecast ensemble of the forecast covariance matrix between the assimilated variables; c) Using O TLM Calculate the covariance matrix of the assimilation vector observation product at each integration moment in the assimilation window. The calculation formula is: Where, O t+wi O t+wl is the observed covariance between the assimilated variables of the i-th dimension and the l-th dimension at time t+w, is the observation covariance matrix between the assimilated variables at time t+w, O TL O TL is the set of observation covariance matrices; d) Using P TN H TL 、H TL H TL and O TL O TL Calculate the Kalman gain matrix at each integration moment in the assimilation window The calculation formula is as follows: O t+wi H t+wl =O t+wi O t+wl +H t+wi H t+wl The vector signal concatenation unit is used to concatenate the mean vector of the assimilated vector forecast set at each integration moment in the assimilation window and the assimilated vector observation product into a time series signal. The concatenation process is as follows: a) Set Γ tLM ,…,Γ t+WLM The forecast mean vector of is connected in series to form a time series signal. The connection method is: AND Wl =[And tl ,…,AND t+wl ,...,AND t+Wl ] Where Y tl ,…,Y t+wl ,…,Y t+Wl are the predicted mean values ​​of the l-th dimension assimilated variables at each integration moment, Y Wl is the time series signal of the predicted mean value of the l-th dimension assimilation variable in the assimilation window; b) Assimilate the vector observation product O at each integration moment tL =[o t1 ,…,o tl ,…,o tL ]′,…,O t+wL =[o t+w1 ,…,o t+wl ,…,o t+wL ]′,…,O t+WL =[o t+W1 ,…,o t+Wl ,…,o t+WL ]′ are connected in series to form a timing signal. The connection method is: the Wl =[o tl ,...,o t+wl ,…,o t+Wl ]1≤l≤L The WL =[o W1 ,...,o Wl ,...,o WL ]′ In the formula, o Wl is the time series signal composed of the observed values ​​of the l-th dimension assimilated variables, O Wl is the time series signal matrix of the assimilated vector observation product; The matrix signal series connection unit is used to connect the Kalman gain matrices at each integration moment in the assimilation window in series to form a Kalman gain matrix time series signal. The series connection method is: g Wnl =[G tnl ,...,G t+wnl ,...,G t+Wnl ]1≤n≤N 1≤l≤L Where g Wnl The Kalman gain matrix G t+wnl , a timing signal consisting of the elements in the nth row and lth column of the set t≤t+w≤t+W; The vector signal frequency decomposition unit is used to perform frequency decomposition on the time series signal of the mean vector of the assimilated vector forecast set and the time series signal of the assimilated variable observation product to obtain their respective frequency components. The frequency decomposition process is as follows: a) Change Y Wl Perform frequency decomposition to obtain frequency components. The frequency decomposition method is as follows: Where Ry Wl [j] is Y Wl The cosine harmonic coefficient value with intermediate frequency j, Iy Wl [j] is Y Wl The sinusoidal harmonic coefficient value with intermediate frequency j; b) Change o Wl Perform frequency decomposition to obtain frequency components. The frequency decomposition method is as follows: Where, Ro Wl [j] is o Wl The cosine harmonic coefficient value with intermediate frequency j, Io Wl [j] is o Wl The sinusoidal harmonic coefficient value with intermediate frequency j; The matrix signal frequency decomposition unit is used to decompose g Wnl Perform frequency decomposition to obtain the frequency components of the time series signals of each matrix element of the Kalman gain matrix. The frequency decomposition method is as follows: Where Rg Wnl [j] is g Wnl The cosine harmonic coefficient value with intermediate frequency j, Ig Wnl [j] is g Wnl The sinusoidal harmonic coefficient value with intermediate frequency j; The correction amount calculation unit is used to calculate the correction amount of the correction vector on different frequency components and synthesize the correction amounts of the correction vector on all frequency components into the correction amount of the correction vector at each integration time in the assimilation window. The synthesis process is as follows: a) Calculate the correction amount of each dimension of the correction vector at different frequency components. The calculation method is: Where Rκ Wn [j] is the correction amount of the cosine harmonic with frequency j in the n-th dimension correction variable time series signal, Iκ Wn [j] is the correction amount of the sinusoidal harmonic with frequency j in the n-th dimension correction variable time series signal; b) The correction amount of each dimensional correction variable on all frequency components 0≤j≤W / 2 is synthesized into the correction amount of each dimensional correction variable at each integration time in the assimilation window. The synthesis formula is as follows: Where, ξ t+wn is the correction amount of the n-th dimension correction variable at time t+w; The correction amount mean calculation unit is used to calculate the mean vector of the correction amount of the correction vector at each integration moment in the assimilation window, and the calculation method is: X N =[ξ1,…,ξ n ,…,x N ]′ Where, ξ n is the mean value of the correction amount of the n-th dimension correction variable at all integration moments, N is the mean vector of the correction amounts of all N-dimensional correction variables; The correction unit is used to correct the correction vector using the mean vector of the correction amount of the correction vector, and the correction method is: X tN =X tN +Ξ N Where, X tN is the value of the correction vector at assimilation time t and the assimilated value after correction.

Images