阅读说明：本技术 一种局地集合资料同化方法 (Method for assimilating local collection data ) 是由 王世璋 曾明剑 于 2021-09-27 设计创作，主要内容包括：本发明公开了一种局地集合资料同化方法,包括：选定分析时刻并生成分析时刻的模式背景场X～(f)和预报误差样本集合X'；读取模式背景场X～(f)、预报误差样本集合X',以及同化分析需要且能够获得的观测y～(o),并根据经纬度对模式预报区域和观测进行分区,并标记存储观测周围的模式格点的位置信息；根据计算条件和实际需求选定单点、单柱、或多柱组合的分析方式；对背景场X～(f)涉及的所有点、或单柱、或多柱组合进行同化分析；输出分析场。本发明的方法实现了模式变量背景误差协方差在同化过程中的直接计算,避免了直接使用集合样本而导致的分析精度下降的问题(次优解),同时解决了在同化垂直积分型观测时难以进行局地化的问题。(The invention discloses a method for assimilating locally-collected data, which comprises the following steps: mode background field X for selecting analysis time and generating analysis time f And a prediction error sample set X'; read mode background field X f A set of forecast error samples X', and observations y needed and available for assimilation analysis o Partitioning the mode forecasting area and observation according to the longitude and latitude, and marking and storing the position information of mode lattice points around the observation; selecting a single-point, single-column or multi-column combined analysis mode according to the calculation conditions and the actual requirements; to the background field X f All points involved, or single columns, or multi-column combinations, are subjected to assimilation analysis; and outputting an analysis field. The method of the invention realizes the direct calculation of the covariance of the background error of the mode variable in the assimilation process, and avoids the analysis precision caused by directly using the set sampleThe problem of degradation of degree (suboptimal solution) and the problem of difficulty in localization in assimilation vertical integration observation are solved.)

1. Local data collecting systemThe chemical conversion method is characterized by comprising the following steps: mode background field X for selecting analysis time and generating analysis time^fAnd a prediction error sample set X'; read mode background field X^fA set of forecast error samples X', and observations y needed and available for assimilation analysis^oPartitioning the mode forecasting area and observation according to the longitude and latitude, and marking and storing the position information of mode lattice points around the observation; selecting a single-point, single-column or multi-column combined analysis mode according to the calculation conditions and the actual requirements; to the background field X^fAll points involved, or single columns, or multi-column combinations, are subjected to assimilation analysis; and outputting an analysis field.

2. The method of assimilating locally aggregated data as claimed in claim 1, wherein partitioning the model forecast area and observations according to longitude and latitude includes: setting the number of blocks of the partition, pnx blocks in the east-west direction and pny blocks in the north-south direction according to calculation conditions and requirements; traversing each lattice point of the first layer of the background field, and marking the position of the lattice point in the partition according to the longitude and latitude of the lattice point; traversing each observation, and marking the position of each observation in the subarea according to the longitude and latitude of each observation; and traversing each observation, searching the nearest n pattern lattice points around each observation, and recording the number of the pattern lattice points.

3. The method of claim 1, wherein said background field X is selected from a group consisting of a background field X, and a background field X^fAll points involved, or single column, or multi-column combinations, for assimilation analysis include the following steps: according to the selected observation y^oDetermining a mode variable related to a background error covariance matrix, and determining an observation subset influencing a current point, or a single column or a multi-column combination according to partition information and a preset localization distance; interpolating the determined mode variables and the disturbances thereof to the positions corresponding to the determined observation subsets, and respectively storing the positions in H_iX^fAnd H_iX'，H_iIs an interpolation operator from a pattern grid point to an observation grid point;

calculate H_iEach mode variable in X' is associated with H_iCorrelation coefficient of other mode variables in XObtaining the standard deviation S of each variable and the correlation coefficient H of the mode variable on the observation grid_iC_correlIs marked as；

ComputingThe localization coefficient between each element obtains a localization coefficient matrix rho, and is compared with the localization coefficient matrix rhoMaking element-element product (Schur product) to obtain localized product(ii) a If it isIs related to a different mode variable, then multiply that element by ρ_m，ρ_mThe range of (1) is 0 to 1.0;

calculating an adjustment factorλAccording toWhereintrace() Traces representing a computational matrix; computing；

Computing the Innovation (d = y) corresponding to the observation^o-H _oX^f) WhereinH _oIs an observation operator, which projects the mode variable onto the observation variable, may be a non-linear operator,H _othe content of (a) varies depending on the particular observation;

computingAndin which H is_oIs thatH _oThe tangent operator of (1); solving equations according to a conjugate gradient methodObtaining a control variable v; calculating a variable x at selected pattern grid points_mDisturbance X of_m' andthe covariance and localization coefficient of each variable in the system are calculated by element-element product to obtain localized row vector C_m，C_mMultiplying by v and then multiplying by the adjustment factorλTo obtain x_mIncrement of (2)(ii) a ComputingTo obtain x_mAnalysis of (2).

4. The method of claim 1, wherein said background field X is selected from a group consisting of a background field X, and a background field X^fAll points involved, or single column, or multi-column combinations for assimilation analysis also include: before the loop code, adding an OpenMP guide statement, setting the number of each point, or single column or multi-column combination as a private variable, setting the number of observations, observation subsets, error covariance, control variables and the like used for assimilation analysis of each point, or single column or multi-column combination as private variables, and setting the rest as common variables; and after the loop code is finished, adding an OpenMP guide statement and finishing OpenMP parallelism.

Technical Field

The invention relates to a weather forecast technology, in particular to a method for assimilating local collection data.

Background

At present, the level of modern weather forecast is closely related to the accuracy of numerical weather forecast, and the effect of numerical weather forecast depends on the accuracy of initial value. Currently, the main way to improve the accuracy of the initial value is data assimilation. The european mid-range weather forecast center (ECMWF) is a considerable contribution to its advanced data assimilation system in the forefront of numerical weather forecasts.

The method aims at assimilating atmospheric data by improving the accuracy of the mode initial value, and the method is to fit between observation and the mode initial value by using a mathematical method according to the background error covariance so that an analysis field obtained by fitting is closer to the true value of the atmosphere theoretically. Therefore, the development of data assimilation itself is also largely centered on how to perform the fitting, how to construct and use the background error covariance expansion.

After decades of development, the fitting method of data assimilation is developed from early gradual correction and optimal interpolation to several current mainstream methods: variational assimilation and ensemble kalman filter assimilation. The starting point of the data assimilation fitting method can be least square or Bayesian estimation, but the starting point can be finally expressed as a cost function which contains observation and background information and takes the observation error covariance and the background error covariance as weights. The analysis field of data assimilation is the solution that makes this function reach a minimum. Background error covariance for data assimilation developed from early distance correlation to current thermodynamic statistical correlations involving multivariate equilibrium. Compared to earlier, current mainstream data assimilation methods are capable of directly assimilating high spatio-temporal resolution non-conventional observations, can take into account cross-error covariance between pattern prediction variables, with analytical fields with higher accuracy (mainly from contributions from high resolution non-conventional observation data) and higher thermodynamic harmony (contributions from thermodynamic equilibrium equation constraints and flow-dependent background error covariance). In recent years, one of the hot spots of data assimilation development is hybrid assimilation, which is the combination of variational assimilation and ensemble kalman filter assimilation based on given weight coefficients, and is characterized by taking the stable effect (static background error covariance) of variational assimilation and the tracking ability of the ensemble kalman filter flow dependent background error covariance on rapid change errors into consideration.

The data assimilation methods used by the current global major business departments are basically variational assimilation, Kalman filtering assimilation and mixed assimilation, and comprise a European middle-term forecast center, a American environmental climate forecast center, a Japanese weather bureau, a British weather bureau, a French weather bureau, a Canada weather bureau, a China weather bureau and the like. The European middle-term forecasting center, the Japan weather bureau, the British weather bureau, the China weather bureau and the like mainly use a variation method, and the Canada weather bureau uses an ensemble Kalman filtering method. Hybrid assimilation is a technology developed by the main meteorological service departments, and the combination of static background error covariance and flow-dependent background error covariance is considered although the implementation approaches are different.

Another driving force for data assimilation development is closely related to the development of numerical calculation conditions, numerical prediction modes and observation means. When the calculation conditions are insufficient, the data assimilation method needs to make many assumptions and simplifications to meet the calculation condition limit, and even the advanced assimilation method cannot obtain ideal results; when the calculation conditions are sufficient, assumption and simplification are not needed any more, and the assimilation effect can be ensured. Currently, with the continuous improvement of the calculation conditions, on one hand, the numerical mode uses higher and higher resolutions to improve the calculation accuracy, and the number of grids is also increased; on the other hand, in order to adapt to different application scenes and requirements, the types of mode mesh structures are more and more abundant, and the mode mesh structures are not limited to regular quadrilateral meshes, triangular meshes, hexagonal meshes, non-structural variable resolution meshes and other structures and are applied to a numerical prediction mode. Therefore, adapting to the scale of the ultra-large grid number and not being limited by the grid structure becomes the trend of the future data assimilation scheme. In addition, atmospheric observation and observation means are increasing, and a large amount of observation data with different spatial resolutions also put forward the requirement of simultaneous assimilation of multi-source multi-scale observation on the assimilation scheme.

In the face of the development trend of adapting to the scale of the number of the ultra-large grids and being not limited by the grid structure and the requirement of simultaneous assimilation of multi-source multi-scale observation, the existing assimilation method has some defects and can be more ideally solved through the method provided by the invention.

a. Prior art general expression forms

The existing assimilation method can use a cost functionJRepresentation (refer to Gao et al, 2004 version):

， （1）

where v is the control variable, C is the square root of the background error covariance, d is the difference between the observed and the initial value of the pattern projected into the observed variable space, H is the tangential operator of the above projection, and R is the observed error covariance matrix. The goal of data assimilation is to solve so thatJIs equal to the solution of the control variable v of 0. When in useJWhen the gradient of (d) is equal to 0, an expression is obtained

， （2）

Where I is the identity matrix.

The main difference between the variational method, the ensemble kalman filtering method, and the hybrid method is how to solve equation (2).

b. Traditional variational method

In mainstream variational and assimilative systems (e.g., WRFDA (Barker et al, 2012), GSI (Hu et al, 2017), ARPS3Dvar (Gao et al, 2004), etc.), equation (2) is solved by a descent algorithm, such as a conjugate gradient method. First the square root C of the background error covariance matrix needs to be reconstructed from the statistical information. The covariance matrix is very large due to the background error of the complete mode space (10)⁷×10⁷) Direct statistical calculations are not possible, so the C matrix is usually constructed using some function. C is usually expressed as the product of three components: c_p、C_v、C_hIn which C is_pIs a physical transformation between a mode variable and a control variable, C_vIs a vertical error covariance matrix, C_hIs a horizontal error matrix. C_pBecause it is often the case that a difference equation is involved,and is therefore closely related to the pattern mesh structure. Different grid resolutions and different grid structures need to be reconstructed C_p。C_vThe calculation is typically accomplished by eigenvalue decomposition. C_hAnd the calculation is realized in the mainstream variational assimilation system by a recursive filtering method.

Of the three C components, C is closely related to the present invention_hThe corresponding recursive filtering method will be described in detail below. The other two components are not described in detail because the specific implementation forms are various and the relevance of the two components to the present invention is weak. According to the definition of the recursive filter (Gao et al, 2004):

, （3）

whereinψ _i Is the firstiInitial value of a grid point (e.g., v at the second in a certain iteration)iThe value of the point),ϕ _i is fromiTonThe value after the filtering of (a) is,χ _i is the value after one complete recursive filtering,βthe decorrelation coefficients of the filtering are determined:

, （4）

whereinLIs at a decorrelation scalesIs the space grid distance (e.g., Δ) in any directionsIn the east and west directions, ΔxIn the north-south direction, Δy) The number of passes of the recursive filtering isN _pass. The recursive filter is usually calculated sequentially from west to east in the east-west direction and then from south to north in the north-south direction of the pattern prediction region. It is therefore necessary to use a regular quadrilateral mesh.

After the calculation scheme for C is determined, the variational assimilation system uses a descent algorithm to solve for v. Here, the conjugate gradient method (Shewchuk 1994) is taken as an example. First, the variables are defined:

， （5）

then calculate in an iterative manner

。 （6）

When in useLess than a certain threshold, e.g. 10^-6When the calculation is stopped, v at that time is output_i+1The value is obtained. Incremental passing of analysis of the pattern variable spaceAnd (4) calculating.

It should be noted that the variational approach is usually to solve over the entire modal prediction region.

c. Ensemble Kalman filtering method

There are many variations of the ensemble kalman filtering method, among which the method with efficient parallel computing power and the ability to assimilate multiple observations simultaneously is the LETKF method (hunt. et al, 2007). This method is also most closely related to the present invention and is described in detail herein. Other variants such as EnSRF, EAKF, ETKF, etc. are not described herein. The LETKF method estimates the background error covariance in an aggregate manner, so that C of equation (2) no longer represents the square root of the background error covariance matrix in LETKF, but rather the aggregate membership X is averaged with the aggregateDifference of (2)。

First, define the variables

。 （7）

Then toThe transpose is subjected to SVD singular vector decomposition to obtainWhere E is the left singular vector matrix, Σ is the singular value matrix, and D is the right singular vector matrix. Substituting the SVD decomposition result into the formula (2) to obtain:

。 （8）

where I is the identity matrix. Incremental passing of analysis of the pattern variable spaceAnd (4) calculating.

The LETKF updates the schema variables on each grid point in turn in the serial version, but the computation of each grid point is independent of each other, so that the parallelization is high. Furthermore, since it is a single point calculation, its localization scheme can be implemented by adjusting the observation error covariance R. With respect to the grid points being analyzed by the LETKF, distant observations give larger error values and close observations give relatively smaller error values. The localization adjustment coefficients are determined from the distance correlation function. The value of the function in the area close to the grid point of the current analysis mode is close to 1.0, the value of the function in the area far away from the grid point of the current analysis mode is close to 0.0, and the distribution form of the function is close to Gaussian distribution.

d. Mixed assimilation method

The core of the existing hybrid assimilation approach is the combination of the static background error covariance matrix and the flow-dependent background error covariance matrix (Wang et al, 2007). There are two forms of mixing:

， （9）

and ， （10）

wherein B represents a background error covariance matrix, subscripts hybrid, s, d represent a mixed background error covariance, a static (static) background error covariance, and a dynamic (dynamic) background error covariance, respectively, x represents a vector of a mode variable space,is the mixing weight coefficient. In the aspect of mixing and assimilating, the invention directly adopts the form of the formula (9) without any innovation.

One of the main reasons for the deficiencies of the above-mentioned prior art is the construction of the background error covariance. The traditional variational method uses recursive filtering to construct horizontal error covariance, so that the method is difficult to adapt to grids with different structures, and LETKF is minimized in an aggregation space, so that the calculation accuracy is relatively low, aggregation must be relied on, and path integral type observation is difficult to process. Therefore, how to construct the background error covariance matrix becomes a key to solve the above problem.

Disclosure of Invention

The invention mainly aims to provide a method for assimilating locally-collected data, which is not influenced by a mode grid structure and can assimilate multisource multiscale observation data simultaneously compared with a variational assimilation method; compared with an LETKF method, the method can better process path integral type observation, has higher calculation efficiency and calculation precision, and is easy to realize the combination of static and dynamic background error covariance. The method is compatible with various grid types for large-scale computing requirements, has the capabilities of multi-source and multi-scale observation and simultaneous assimilation, and gives consideration to the computational accuracy of assimilation to a certain extent.

The technical scheme adopted by the invention is as follows: a method for assimilating locally aggregated data, comprising: mode background field X for selecting analysis time and generating analysis time^fAnd a prediction error sample set X'; read mode background field X^fA set of forecast error samples X', and observations y needed and available for assimilation analysis^oRoot of Chinese angelicaPartitioning the mode forecasting area and observation according to the longitude and latitude, and marking and storing the position information of mode grid points around the observation; selecting a single-point, single-column or multi-column combined analysis mode according to the calculation conditions and the actual requirements; to the background field X^fAll points involved, or single columns, or multi-column combinations, are subjected to assimilation analysis; and outputting an analysis field.

Further, the partitioning the mode forecast area and the observation according to the longitude and latitude includes: setting the number of blocks of the partition, pnx blocks in the east-west direction and pny blocks in the north-south direction according to calculation conditions and requirements; traversing each lattice point of the first layer of the background field, and marking the position of the lattice point in the partition according to the longitude and latitude of the lattice point; traversing each observation, and marking the position of each observation in the subarea according to the longitude and latitude of each observation; and traversing each observation, searching the nearest n pattern lattice points around each observation, and recording the number of the pattern lattice points.

Further, the pair of background fields X^fAll points involved, or single column, or multi-column combinations, for assimilation analysis include the following steps: according to the selected observation y^oDetermining a mode variable related to a background error covariance matrix, and determining an observation subset influencing a current point, or a single column or a multi-column combination according to partition information and a preset localization distance; interpolating the determined mode variables and the disturbances thereof to the positions corresponding to the determined observation subsets, and respectively storing the positions in H_iX^fAnd H_iX'，H_iIs an interpolation operator from a pattern grid point to an observation grid point; calculate H_iEach mode variable in X' is associated with H_iThe correlation coefficients of other mode variables in X' are obtained to obtain the standard deviation S of each variable and the correlation coefficient H of the mode variable on the observation grid_iC_correlIs marked as(ii) a ComputingThe localization coefficient between each element obtains a localization coefficient matrix rho, and is compared with the localization coefficient matrix rhoMaking element-element product (Schur product) to obtain localized product(ii) a If it isIs related to a different mode variable, then multiply that element by ρ_m，ρ_mThe range of (1) is 0 to 1.0; calculating an adjustment factorλAccording toWhereintrace() Traces representing a computational matrix; computing(ii) a Computing the Innovation (d = y) corresponding to the observation^o-H _oX^f) WhereinH _oIs an observation operator, which projects the mode variable onto the observation variable, may be a non-linear operator,H _othe content of (a) varies depending on the particular observation; computingAndin which H is_oIs thatH _oThe tangent operator of (1); solving equations according to a conjugate gradient methodObtaining a control variable v; calculating a variable x at selected pattern grid points_mDisturbance X of_m' andthe covariance and localization coefficient of each variable in the system are calculated by element-element product to obtain localized row vector C_m，C_mMultiplying by v and then multiplying by the adjustment factorλTo obtain x_mIncrement of (2)(ii) a ComputingTo obtain x_mAnalysis of (2).

Further, the pair of background fields X^fAll points involved, or single column, or multi-column combinations for assimilation analysis also include: before the loop code, adding an OpenMP guide statement, setting the number of each point, or single column or multi-column combination as a private variable, setting the number of observations, observation subsets, error covariance, control variables and the like used for assimilation analysis of each point, or single column or multi-column combination as private variables, and setting the rest as common variables; and after the loop code is finished, adding an OpenMP guide statement and finishing OpenMP parallelism.

The invention has the advantages that: the present invention provides a method for assimilating locally integrated data. The method constructs a local background error covariance of the mode variable projected to an observation grid, realizes the direct calculation of the mode variable background error covariance in the assimilation process, avoids the problem of analysis precision reduction (suboptimal solution) caused by the direct use of an aggregate sample by the traditional local aggregation assimilation method, and solves the problem that the local aggregation assimilation method is difficult to perform local aggregation in the assimilation vertical integral observation.

Compared with the traditional local area collection assimilation method, the local area collection data assimilation method improves the calculation precision, introduces a conjugate gradient descent algorithm and a single-column and multi-column combined analysis method, can greatly reduce repeated calculation of observation overlapping regions, remarkably accelerates the speed of assimilation analysis, and provides calculation efficiency guarantee for future service or actual use of the method.

Compared with the traditional local area collection assimilation method, the local area collection data assimilation method can realize mixed assimilation by combining static background errors, and is beneficial to improving the assimilation analysis effect when the number of forecast error samples is insufficient.

Compared with the assimilation method of the traditional variational framework using the recursive filtering technology, the assimilation method of the local collection data is not limited by a mode prediction grid structure, can adapt to more application scenes, and is easy to transplant.

In addition to the objects, features and advantages described above, other objects, features and advantages of the present invention are also provided. The present invention will be described in further detail below with reference to the drawings.

Drawings

The accompanying drawings, which are incorporated in and constitute a part of this application, illustrate embodiments of the invention and, together with the description, serve to explain the invention and not to limit the invention.

FIG. 1 is a graph comparing the computational efficiency of the local congregation assimilation method of the invention with LETKF;

FIG. 2 is a mode variable (a) of the present inventionu,(b)v,(c)θ,(d)q _vA background error map of (1);

FIG. 3 is a model variable (a) of the assimilation observation of the present invention as observation type 2u,(b)v,(c)θ,(d)q _vA background error map of (1);

FIG. 4 is a model variable (a) of the assimilation observation of the present invention as observation type 3u,(b)v,(c)θ,(d)q _vBackground error map of (1).

Detailed Description

In order to make the objects, technical solutions and advantages of the present invention more apparent, the present invention is described in further detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are merely illustrative of the invention and are not intended to limit the invention.

FIG. 1 is a graph comparing the calculated efficiencies of the LETKF and the local area aggregation assimilation method (the black solid line is a trend line of linear regression fitting), the vertical axis is the ratio of the wall clock time required for the local area aggregation assimilation method to complete all the lattice point analyses to the time used by the LETKF, and the horizontal axis is the total number observed;

FIG. 2 mode variables (a)u,(b)v,(c)θ,(d)q _vBackground error ofVertical distribution of difference plots (cross-labeled solid line), local ensemble assimilation method (LEDA) analysis error (solid line), and LETKF analysis error (dashed line). The vertical axis is the mode vertical layer, and the horizontal axis is the error magnitude. Assimilation observations are observation type 1;

FIG. 3 is a model variable (a) of the assimilation observation of the present invention as observation type 2u,(b)v,(c)θ,(d)q _vA background error map of (1);

FIG. 4 is a model variable (a) of the assimilation observation of the present invention as observation type 3u,(b)v,(c)θ,(d)q _vBackground error map of (1).

The invention relates to a method for assimilating locally-collected data, which comprises the following steps:

I) mode background field X for selecting analysis time and generating analysis time^fAnd a prediction error sample set X';

II) read mode background field X^fA set of forecast error samples X', and observations y needed and available for assimilation analysis^oPartitioning the mode forecasting area and observation according to the longitude and latitude, and marking and storing the position information of mode lattice points around the observation;

III) selecting a single-point, single-column or multi-column combined analysis mode according to the calculation conditions and the actual requirements (the multi-column combined analysis consumes more memory, but has faster assimilation analysis speed);

IV) selection based on III), according to assimilation formulas (1), (2), (5) and (6), and using the formula of covariance of local background errorTo the background field X^fAll points involved, or single columns, or multi-column combinations, are subjected to assimilation analysis;

v) outputting an analysis field after the step IV) is finished.

In the step II), the pattern prediction area and the observation partition comprise the following specific steps:

(1) setting the number of blocks of the partition, pnx blocks in the east-west direction and pny blocks in the north-south direction according to calculation conditions and requirements;

(2) traversing each grid point at the first layer of the background field, marking the position of the grid point in the subarea according to the longitude and latitude of the grid point, wherein pnx and pny are both 10, and the subarea numbers of the point at the middle most of the background field and the adjacent point are both (5, 5);

(3) traversing each observation, and marking the position of each observation in the subarea according to the longitude and latitude of each observation;

(4) based on the results of (2) and (3), each observation is traversed, the nearest n pattern grid points (n is given according to the needs of users and is usually 5-9) around each observation are searched, and the number of the pattern grid points is recorded.

In the step IV), the assimilation analysis of each point, or single column, or multi-column combination comprises the following specific steps:

(1) according to the observation y selected in II)^oDeterminingInThe mode variables (including the number of vertical layers involved) to which the matrix relates; determining an observation subset influencing the current point, or a single column or a multi-column combination according to the partition information in the step II) and a preset localization distance;

(2) interpolating the mode variable (background field or set average) and the disturbance (difference between set member and set average) determined in (1) to the positions corresponding to the observation subsets determined in (1), and respectively storing the positions in H_iX^fAnd H_iX'（H_iThe method is characterized in that the method is an interpolation operator from a mode grid point to an observation grid point, an interpolation algorithm is determined by a user, and distance weight interpolation is taken as an example in the invention); the position of the mode variable used for interpolation is determined by step II);

(3) calculate H_iEach mode variable (i.e., mode variable corresponding to each row) in X' is associated with H_iThe correlation coefficients of other mode variables (other rows) in X' are obtained to obtain the standard deviation S and H of each variable_iC_correl(is described as）；

(4) ComputingThe localization coefficient between each element obtains a localization coefficient matrix rho, and is compared with the localization coefficient matrix rhoMaking element-element product (Schur product) to obtain localized product(ii) a If it isIs related to a different mode variable, then multiply that element by ρ_m，ρ_mThe range of (1) is 0 to 1.0; preferably,. rho_m=0.5；

(5) Calculating an adjustment factorλAccording toWhereintrace() Traces representing a computational matrix;

(6) computing

(7) Computing the Innovation (d = y) corresponding to the observation^o-H _oX^f) WhereinH _oIs an observation operator, which projects the mode variable onto the observation variable, may be a non-linear operator,H _othe content of (a) varies depending on the particular observation;

(8) computingAndin which H is_oIs thatH _oThe tangent operator of (1);

(9) calculating a control variable v according to the formula (5) and the formula (6) of the conjugate gradient method;

(10) calculating a variable x at selected pattern grid points_mDisturbance X of_m' andthe covariance and localization coefficient of each variable in the system are used to obtain a localized row vector C by using element-element product (Schur product)_m，C_mDot-multiplied with v and then multiplied by an adjustment coefficientλTo obtain x_mIncrement of (2)；

(11) ComputingTo obtain x_mAnalysis of (2).

In the step IV), assimilation analysis is performed on all points, or a single column, or a combination of multiple columns, and parallel computation is performed using methods such as OpenMP, so as to improve computation efficiency, including the following steps:

(1) before the loop code of the step IV), adding an OpenMP guide statement, setting the number of each point, or single column or multi-column combination in the step IV) as a private variable, setting the observation quantity, observation subset, error covariance, control variable and the like used for assimilation analysis of each point, or single column or multi-column combination as private variables, and setting the rest as common variables;

(2) step IV) is performed normally;

(3) and after the loop code in the step IV) is finished, adding an OpenMP guide statement, and finishing OpenMP parallelism.

A specific embodiment is given below, which does not need local numerical prediction, uses simulated observation, and facilitates the user to test the assimilation algorithm under the condition of limited calculation conditions and no actual observation:

(I.1) from the NCDC official network https:// www.ncdc. noaa. gov/data-access/model-data downloads GFS data with a resolution of 0.5 °, where a 36 hour forecast of 12UTC is downloaded as a mode background field X on any day^f. The GFS prediction data predicted to that time, 5 days prior to that time, is downloaded for use in generating the set of prediction error samples X' (the time lag generating set member). The GFS analysis field at that time is downloaded as a true value (for verifying the analysis results). The WRF v3.9.1 mode and the WPS v4.1 are used for converting the downloaded GFS data into regional data, and the grid resolution and the regional range can be set by a user. Here, 15km grid resolution and 150 × 150 grid points are taken as an example.

(I.2) generating a simulated observation:

simulation observation type 1: conventional sounding observations, say H_o=H _o= I, where I is the identity matrix. The observed variable being two components of the horizontal wind fielduAndvtemperature of the bodyθRatio of water to steamq _vAnd ground pressureps. The density of the type 1 simulation observation in the east-west direction and the south-north direction is one for every 7 grid points, and the observation error isuAndv0.5 m/s, a bit temperatureθ0.5K, water-vapor mixing ratioq _v0.0005 kg/kg, and ground pressurepsIs 10 Pa. Except thatpsAnd the other observation variables are taken one for each 2 mode layers in the vertical direction.

Simulation observation type 2: the whole layer can reduce the water vapor content (pwv). H_o=H _o=1/gΣ（q_vDp), Σ denotes the integral from the first layer of the pattern to the highest layer, dp denotes the difference in gas pressure between each layer, and qv is the water vapor mixing ratio of each layer. The observation density is one for every 4 grid points in the east-west direction and the south-north direction, and the magnitude of the observation error is 1.0.

Simulated observation type 3: wind profile observation, containing two variables: wind direction and wind speed. Of wind directionH _o=arctan(u/v) + n pi, where n is an integer, such thatH _oBetween plus and minus 0.5 pi; wind direction H_o=u/sδv-v/sδu. Of wind speedH _o=（u²+v²）^0.5(ii) a Wind direction H_o=u/sδu+v/sδvWherein s = (u)²+v²）^0.5. The observation density is that every 3 grid points in the east-west direction and the south-north direction are one, the vertical direction comprises all mode layers, and the sizes of observation errors of the wind direction and the wind speed are respectively 0.5 m/s and 0.1 radian.

And (I.3) setting operation parameters of the assimilation system, including localization radiuses, search distances and the like of the three types of observation.

(I.4) reading mode background field X^fAnd a set of prediction error samples X', and (i.2) the three observations.

(I.5) partitioning the mode data and observed longitude and latitude, and marking mode grid points around the observation:

the east-west direction and the south-north direction of the forecast area are respectively divided into 10 blocks, namely pnx and pny are equal to 10; traversing each lattice point of the first layer of the background field, and marking the position of the lattice point in the partition according to the longitude and latitude of the lattice point; traversing each observation, and marking the position of each observation in the subarea according to the longitude and latitude of each observation; and traversing each observation, searching the nearest 5 pattern lattice points around each observation, recording the number of the pattern lattice points, and only searching in the adjacent partitions.

(i.6) if the multi-column combination analysis method is used, the number of groups and the pattern grid points corresponding to each group are recorded (in this example, 9-column combination analysis is used).

(I.7) assimilation analysis was performed for each multi-column combination:

determining an observation subset influencing the current multi-column combination according to the partition information in the step II) and a preset localization distance; interpolating the mode variable (background field or set average) and its disturbance (difference between set member and set average) to the position of observation grid in observation subset, and storing them in H_iX^fAnd H_iX'（H_iIs the distance weight interpolation operator of the mode grid points to the observation grid points, where the mode variable positions participating in the interpolation are determined by (i.5); calculate H_iEach mode variable (i.e., mode variable corresponding to each row) in X' is associated with H_iThe correlation coefficients of other mode variables (other rows) in X' are obtained to obtain the standard deviation S and H of each variable_iC_correl(is described as) (ii) a ComputingThe localization coefficient among each element obtains a localization coefficient matrix rho, and

andmaking element-element product (Schur product) to obtain localized product(ii) a If it isIs related to a different mode variable, then multiply that element by ρ_m(= 0.5); calculating an adjustment factorλAccording to(ii) a ComputingComputing the Innovation (d = y) corresponding to the observation^o-H _oX^f) WhereinH _oIs an observation operator that projects the mode variable onto the observed variable, which may be a non-linear operator, in this exampleH _oThe content of (1) is determined according to the formula (I.2); computingAndin which H is_oIs thatH _oThe tangent operator of (a), determined by (i.2); calculating a control variable v according to the formula (5) and the formula (6) of the conjugate gradient method; calculating variables at each pattern grid point contained in the current multi-column combinationx_mDisturbance X of_m' andcovariance of each variable in (a) and localization coefficients (including p)_m) The locally formed row vector C is obtained as an element-element product (Schur product)_m，C_mDot-multiplied with v and then multiplied by an adjustment coefficientλTo obtain x_mIncrement of (2)(ii) a ComputingTo obtain x_mAnalysis of (2); and entering the next multi-column combination (if the OpenMP parallel mode is started, determining the number n of the multi-column combinations which are calculated simultaneously according to the number of used CPU cores, and updating the n multi-column combinations simultaneously in I.7).

And (I.8) completing assimilation analysis of all mode grid points and mode variables, outputting analysis results, and ending the assimilation operation process.

The key technical problem details solved by the local collection data assimilation scheme designed by the invention are as follows:

1. by spatial projection, a direct calculation of the local background error covariance matrix (square root) is achieved.

In order to solve the problem that the covariance matrix of the background error is difficult to be directly calculated, the invention redesigns the square root C matrix of the covariance of the background error in the formula (2). In general, the C matrix in equation (2) involves all the mode variables of the complete prediction grid region, and thus the degree of freedom is very high. And the invention projects the C matrix. (2) HC in the formula can be decomposed into two parts H_oAnd H_iC, wherein H_oAnd H_iRespectively the transformation of the mode variables into the observation variables and the projection (interpolation) of the mode mesh into the observation mesh. When the observation involves an integral stratification (e.g. pwv observation), H_iA whole-slice projection (interpolation) containing multiple gas columns. This approach solves the problem of the path integral that the LETKF is difficult to processThe problem of type observation. Use of H in the invention_iC is the square root of the covariance matrix of the background error and is noted as. As can be seen by comparison with C,is greatly reduced because ofInvolving only calculation of H_oThe required mode variables. Meanwhile, the concept of localization is performed with reference to LETKF: observations far from the current analysis grid can be disregarded, and for any analysis grid,involving only one locally-wide mode variable, and thus allowing explicit computation to be performed directly. Here will beReferred to as local background error covariance. When estimating C using set-wise approach, one could theoretically first calculate B = XX locally^TAnd applying a localization coefficient matrix to B through Schur product, and obtaining localized coefficients through SVD decomposition. However, the method has large calculation amount and low calculation efficiency, so that the invention designs a direct estimation methodThe matrix method comprises the following steps:

， （11）

where ρ is a matrix of localized coefficients, C_correlIs a matrix of correlation coefficients for background errors, S is a diagonal matrix, the elements are the standard deviations of the mode variables,λis a magnitude adjustment factor such thatTrace of (D) and (C)_correlEqual (i.e., the two matrices are of comparable magnitude). If the adjustment coefficient is not multiplied, the eigenvalue of the background error covariance matrix estimated by equation (11) is the original background error covariance matrix B (B = SC)_correlS^T) The square of the eigenvalue is doubled, and the background error covariance is too large, so that the assimilation analysis result is unreasonable.λIs calculated asWhereinNIs C_correlThe rank (or number of rows),trace() The traces of the computation matrix are represented.

Here pairThe size of the matrix is described. Suppose H_oIs onenobs×nhThe matrix of (a) is,nobsis the number of observations within a certain radius around a pattern lattice point,nhis the number of mode variables (the product of the number of lattice points and the number of types of mode variables) involved in the observation, thenIs of a size ofnh×nh. In the actual calculation, the calculation results are,nhof the order of about 10⁰~10³Much smaller than 10 for the complete C matrix, depending on the number of local observations and the number of mode variables involved in the observations⁷And the magnitude is completely within the range of the current single-core computing capability. Due to the fact thatThe matrix is locally processed, soIs a diagonal-element-based matrix with a condition number (ratio of maximum eigenvalue to minimum eigenvalue) much smaller than X of LETKF and therefore less ill-conditioned. Furthermore, since the projection process does not make any assumptions, the rank of the background error covariance is not reduced.

In addition, due toThe elements can be directly obtained by collecting sample statistics, and localization can also be directly calculated according to longitude and latitude differences among observations, so that the elements are usedThe local data collection assimilation method can avoid the limitation of the grid structure of the recursive filtering method.

Another advantage of directly calculating the background error covariance is that the method offers the possibility to adjust the covariance size between the mode variables. Because numerical patterns are highly non-linear, the covariance between the pattern variables tends to be difficult to accurately characterize the trend of change between the two pattern variables. If the covariance is incorrect, or the magnitude is too large, an error or too large increment is generated, so that the assimilation analysis cannot achieve the effect of reducing the error. The size of the covariance among the variables is properly adjusted, so that overlarge and unreasonable analysis increment is favorably reduced, and the assimilation analysis effect is favorably improved. The specific mode is as follows: record C_correlIf the two variable types are different, multiplying the two variable types by a preset coefficient rho_m。ρ_mThe size of (a) is between 0 and 1. In the present invention, 0.5 is suggested.

2. The method for simultaneously analyzing by using the conjugate gradient method and the multiple gas columns realizes the high-efficiency and quick calculation of assimilation analysis.

The formulae (5) and (6) used in the conjugate gradient method requireThe matrix is symmetric and positive definite. It was observed that when estimating C using the pattern of LETKF (pooled samples),the condition number of (2) is very high, and the solution of the conjugate gradient method cannot be converged. Therefore, the solution can be achieved only by the SVD method. Used in the present inventionThe matrix (11 type) is constructed by greatly reducing the condition numberThe matrix enables convergence of the solution process of the conjugate gradient method. The conjugate gradient method is mainly based on the dot product of the vector, so the calculation speed is far faster than the SVD decomposition speed, especially in the case of large background error covariance matrix.

Due to the fact thatThe method is equivalent to the square root of a background error covariance matrix thinned in the horizontal direction, so that the localization process does not aim at a certain pattern lattice point, and a single point, a single column (one gas column) and a plurality of columns (a plurality of gas columns) can be selected for simultaneous analysis, thereby reducing the repeated calculation of observation overlapping regions, and accelerating the speed and the operating efficiency of assimilation analysis.

3. Facilitating assimilation schemes for mixing different types of covariance

Compared with the LETKF method which cannot use the static background error covariance, the local set data assimilation method can conveniently realize the mixing of the static background error covariance and the dynamic background error covariance according to the formula (9). When the static background error covariance used is estimated as a distance correlation function,may be represented as:

， （12）

wherein the subscripti,jTo represent the first of the matrixiGo to the firstjColumn, superscript Ens denotes aggregate samples,ρthe distance-related coefficient is represented by,Sthe standard deviation of the mode variable corresponding to the ith row is shown. When the weight coefficient y is reduced to 0,becomes a complete static covariance.

Examples

(1) Compared with LETKF using matrix decomposition (eigenvalue decomposition or SVD decomposition), the method has higher calculation speed

A set of test comparisons is presented here, the tests taking into account different observed densities. It can be seen from fig. 1 that the local congregation assimilation method of the present invention is much smaller in calculation time than the lettf method, and the time required for assimilation of about 15 ten thousand observations is only 20% or less of the lettf. This result is due, on the one hand, to the reduction of the iterative calculations and the application of the conjugate gradient method by the local collective assimilation method according to the invention, and, on the other hand, to the reduction of the number of assimilating units by the multi-column combination (e.g. the 9-column combination here). In LETKF, the number of units is the total number of lattice points; in the local area aggregation assimilation method, the number of the units is the total number of the lattice points divided by the number of the lattice points contained in the multi-air column combination. The more gas columns within each cell, the fewer the total number of cells. But the observations involved in each cell also increase, so the number of gas columns per cell cannot be increased wirelessly. The calculation speed can be obviously improved by using about 9 air columns.

(2) The analysis error of the local area collective assimilation method of the present invention is comparable or smaller compared to LETKF

A set of test comparisons is given here, the tests taking into account different observation types and combinations (see D2.2 embodiment (i.2) to generate simulated observations). As can be seen from FIG. 2, when conventional sounding observations (type 1) are assimilatedThe error of u and v analyzed by the local collective assimilation method is smaller than that of LETKF at the lower layer. The errors of the analyzed theta in the middle and lower layers are smaller than those of the LETKF. Two methods are in water vaporq _vIs roughly equivalent analytically. As can be seen from fig. 3, when assimilating observation type 2 (reducible water vapor content pwv), the locally integrated assimilation method of the present invention was approximately equivalent to the analysis effect of LETKF on u and v, but did not cause an increase in the low layer θ error as with LETKF. The adjustment level of the two to the water vapor is basically equivalent. It can be seen from fig. 4 that when assimilating nonlinear observations (wind direction, wind speed — observation type 3), both do not analyze as well as direct observations (type 1), but the wind field components u and v of the local collective assimilation method of the present invention are slightly smaller than lettf at the lower and upper levels, and slightly larger than the latter only at the middle level. Meanwhile, the local area collection assimilation method avoids the increase of errors of non-observed variables (theta and qv) caused by LETKF. Overall, the local area aggregation assimilation method provided by the invention has at least LETKF equivalent in analysis effect, has better analysis capability on direct observation, and has less negative influence (error increase) on indirect observation.

(3) The local area set assimilation method has higher expandability and easy portability

Compared with a variation assimilation method using a recursive filtering technology, the local set assimilation method is not limited by the conditions of grid structures, continuous change of resolution ratios and the like, and can be applied to different types of numerical modes; compared with the LETKF method, the local ensemble assimilation method can use mixed background error covariance, even completely static background error covariance, and can complete assimilation analysis under the condition of limited ensemble number or no ensemble member.

The above description is only for the purpose of illustrating the preferred embodiments of the present invention and is not to be construed as limiting the invention, and any modifications, equivalents, improvements and the like that fall within the spirit and principle of the present invention are intended to be included therein.