CN113568067B - Numerical weather forecasting method and device, computer storage medium and electronic equipment - Google Patents

Numerical weather forecasting method and device, computer storage medium and electronic equipment Download PDF

Info

Publication number
CN113568067B
CN113568067B CN202110810865.4A CN202110810865A CN113568067B CN 113568067 B CN113568067 B CN 113568067B CN 202110810865 A CN202110810865 A CN 202110810865A CN 113568067 B CN113568067 B CN 113568067B
Authority
CN
China
Prior art keywords
observation
sample
simulation
error correction
assimilation
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Active
Application number
CN202110810865.4A
Other languages
Chinese (zh)
Other versions
CN113568067A (en
Inventor
张洪芹
田向军
金哲
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Institute of Atmospheric Physics of CAS
Institute of Tibetan Plateau Research of CAS
Original Assignee
Institute of Atmospheric Physics of CAS
Institute of Tibetan Plateau Research of CAS
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Institute of Atmospheric Physics of CAS, Institute of Tibetan Plateau Research of CAS filed Critical Institute of Atmospheric Physics of CAS
Priority to CN202110810865.4A priority Critical patent/CN113568067B/en
Publication of CN113568067A publication Critical patent/CN113568067A/en
Application granted granted Critical
Publication of CN113568067B publication Critical patent/CN113568067B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01WMETEOROLOGY
    • G01W1/00Meteorology
    • G01W1/10Devices for predicting weather conditions

Abstract

The invention discloses a numerical weather forecasting method and device, a computer storage medium and electronic equipment, and belongs to the field of weather forecasting. Firstly, obtaining a prior value of a state vector by inputting meteorological data through a numerical forecasting mode, and simultaneously constructing a group of state vector collection samples; inputting the prior state vector and the set sample into an observation operator to obtain prior value simulation observation and sample simulation observation; and inputting the obtained information and observation data into an i4DVar assimilation solving module, and obtaining an integral error correction term by minimizing a cost function. The invention improves the classic 4DVar method into the invention of the integral correction four-dimensional variable assimilation algorithm (i4DVar), adds an integral error correction term at the initial time of each sub-time window, and integrally corrects the initial error and the mode error, thereby not only ensuring the assimilation precision, but also reducing the calculation amount.

Description

Numerical weather forecasting method and device, computer storage medium and electronic equipment
Technical Field
The present invention relates to the field of weather forecasting, and in particular, to a method and an apparatus for numerical weather forecasting, a computer storage medium, and an electronic device.
Background
Numerical weather prediction (also called numerical prediction) refers to a method of performing numerical calculation by a large-scale computer under certain initial value and side value conditions according to the actual conditions of the atmosphere, solving a fluid mechanics and thermodynamic equation set describing the weather evolution process, and predicting the atmospheric motion state and the weather phenomenon in a certain time period.
Data assimilation is a data processing technology for providing an initial field for numerical weather forecast, and is widely applied to the field of atmospheric oceans at present. Data assimilation is the process of using optimization methods and theories to fuse observed data with the simulation results of numerical patterns sufficiently in order to obtain the best analytical field.
The four-dimensional variational assimilation method (4DVar) is one of the most important methods in the field of data assimilation, and 4DVar is roughly classified into two types according to a processing method for a mode error: one is a strongly constrained 4DVar method, and the other is a weakly constrained 4DVar method. Both the strong constraint 4DVar method and the weak constraint 4DVar method belong to nonlinear optimization problems essentially, and the optimization algorithm based on tangent and adjoint modes is often adopted for solving the optimization problems, so that the tangency and adjoint modes of numerical modes need to be given, and the difficulty is very high. The strongly constrained 4DVar method usually assumes that a numerical mode is perfect and has no mode error, and all prediction errors of the mode come from initial field errors; and the weak constraint 4DVar method introduces a mode error in the optimization process and optimizes the initial error and the mode error respectively.
The existing strong constraint 4DVar method and the weak constraint 4DVar method have respective disadvantages for processing the mode error and the initial error:
the strong constraint 4DVar method ignores the existence of mode errors, the final analysis increment is an error correction term placed at the initial moment of the assimilation window, and the mode errors of the initial moment and all the assimilation windows are integrally corrected through the error correction term. While in practice mode errors are inevitable, this method of ignoring mode errors is bound to lose assimilation accuracy. Meanwhile, the mode of ignoring the mode error and only adding the correction term at the initial moment of the assimilation window is only suitable for the situation that the influence of the mode error is far less than that of the initial error, and in order to ensure that the initial error is dominant compared with the mode error, the length of the assimilation window is often required to be shortened, so that the strong constraint 4DVar method is difficult to absorb more observation data, and the assimilation precision of the 4DVar is also influenced to a great extent.
The weak constraint 4DVar tries to introduce an error model in the optimization process, the respective optimization of the initial error and the mode error is realized through modeling of the error model, the degree of freedom and the calculation amount of the 4DVar optimization problem are increased suddenly, and further great uncertainty is brought to the final optimization result. And because the mode error source is complex and accurate modeling is almost impossible, great uncertainty exists in the establishment of the error model, and the uncertainty also influences the assimilation precision of the final 4DVar method.
Disclosure of Invention
In order to solve the problems in the prior art, the invention provides a method and a device for numerical weather forecast, a computer storage medium and electronic equipment.
The technical scheme provided by the invention is as follows:
in a first aspect, the present invention provides a method of numerical weather forecasting, the method comprising:
acquiring a prior state vector and a set sample of a current assimilation window;
the assimilation window is divided into a plurality of sub-time windows according to a certain time interval, and an integral error correction term is added to the prior state vector and the set sample at the initial moment of each sub-time window;
respectively inputting the prior state vector and the set sample into an observation operator to obtain a prior simulation observation result and a sample simulation observation result;
constructing a cost function of an overall correction four-dimensional variational assimilation algorithm, initializing an overall error correction term, taking the prior state vector, the set sample, the prior simulation observation result, the sample simulation observation result and the obtained observation data as known information, and obtaining a final overall error correction term by minimizing the cost function;
and adding the final integral error correction term to the prior state vector, and performing mode integration to obtain a weather forecast result of the current assimilation window.
Further, the cost function l of the overall corrected four-dimensional variable assimilation algorithm is as follows:
Figure BDA0003168138100000031
wherein x' is an integral error correction term, QbIs a covariance matrix of x',
Figure BDA0003168138100000032
in order to be a vector of a-priori states,
Figure BDA0003168138100000033
is a prior state vector
Figure BDA0003168138100000034
Input to observation operator Hk(. A) derived prior simulated observation, yobs,kFor observation of data, RkIs a covariance matrix of the observed errors.
Further, obtaining a prior state vector of a current assimilation window by the following method:
will assimilate window [ t0,tS]Divided into (t) by time interval τs-t0) τ sub-time windows;
interpolating from the re-analyzed meteorological data covering the area to be studied to obtain t0Background field of time xb,0(ii) a Or, for t0T is obtained by integrating the forecast model in the background field a period of time before the moment0Background field of time xb,0
Using numerical prediction model Mk(·)、t0Background field of time
Figure BDA0003168138100000035
And an overall error correction term
Figure BDA0003168138100000036
The following background simulation was performed to obtain an assimilation window [ t ]0,tS]Internal state variable
Figure BDA0003168138100000037
The background simulation result of (1);
Figure BDA0003168138100000038
wherein M isk(. is) from tk-1To tkThe numerical value prediction model of (1) is,
Figure BDA0003168138100000039
model M for numerical forecastingk(ii) a state variable of (c),
Figure BDA00031681381000000310
when k is t0Time, right side of formula
Figure BDA00031681381000000311
All state variables
Figure BDA00031681381000000312
The background simulation results of (a) constitute the prior state vector.
Further, acquiring a set sample of the current assimilation window by the following method:
obtaining a historical sample, and obtaining a group of set sample disturbances P by adopting a disturbance state variable method on the historical samplex=(x'1,…,x'N);
Using numerical prediction model Mk(·)、t0Background field of time
Figure BDA00031681381000000313
And aggregate sample perturbation
Figure BDA00031681381000000314
An integrated simulation was performed to obtain an assimilation window [ t ]0,tS]Inner set sample xj,kThe simulation result of (2);
Figure BDA0003168138100000041
wherein j is 1,2, …, N, when k is t0Time, right side of formula
Figure BDA0003168138100000042
All set samples xj,kThe simulation results of (a) make up the aggregate sample.
Further, the obtaining of the final overall error correction term by minimizing the cost function includes:
covariance matrix Q of initialized integral error correction term x' 0, xbPrior state vector
Figure BDA0003168138100000043
A priori simulation of the observation
Figure BDA0003168138100000044
Observation data yobs,kAnd covariance matrix R of observation errorskInputting the cost function l;
expressing the integral error correction term x' as set sample disturbance PxLinear combination of (a) x ═ Pxβ, wherein Px=(x'1,…,x'N),β=(β1,…,βN) Is a coefficient vector;
using sample covariance matrices
Figure BDA0003168138100000045
Covariance matrix Q instead of xb
Changing x' to PxBeta and
Figure BDA0003168138100000046
and substituting the cost function l, and performing equivalent deformation to obtain the following cost function taking beta as a new variable:
Figure BDA0003168138100000047
wherein the content of the first and second substances,
Figure BDA0003168138100000048
when x' is 0
Figure BDA0003168138100000049
Mathematically transforming the cost function with beta as a new variable into a nonlinear least squares form, andsolving by adopting a Gauss-Newton iterative algorithm to obtain beta, and obtaining beta according to the x ═ Pxβ yields the final overall error correction term.
In a second aspect, the present invention provides a numerical weather forecasting apparatus, the apparatus comprising:
the acquisition module is used for acquiring a prior state vector and a set sample of a current assimilation window;
the assimilation window is divided into a plurality of sub-time windows according to a certain time interval, and an integral error correction term is added to the prior state vector and the set sample at the initial moment of each sub-time window;
the simulation module is used for respectively inputting the prior state vector and the set sample into an observation operator to obtain a prior simulation observation result and a sample simulation observation result;
the solving module is used for constructing a cost function of an overall correction four-dimensional variational assimilation algorithm, initializing an overall error correction term, and obtaining a final overall error correction term by minimizing the cost function by taking the prior state vector, the set sample, the prior simulation observation result, the sample simulation observation result and the obtained observation data as known information;
and the forecasting module is used for adding the final integral error correction term to the prior state vector and carrying out mode integration to obtain a weather forecasting result of the current assimilation window.
Further, the cost function l of the overall corrected four-dimensional variable assimilation algorithm is as follows:
Figure BDA0003168138100000051
wherein x' is an integral error correction term, QbIs a covariance matrix of x',
Figure BDA0003168138100000052
in order to be a vector of a-priori states,
Figure BDA0003168138100000053
is a prior state vector
Figure BDA0003168138100000054
Input to observation operator Hk(. A) derived prior simulated observation, yobs,kFor observation of data, RkIs a covariance matrix of the observed errors.
Further, obtaining a prior state vector of the current assimilation window through the following units:
a dividing unit for dividing the assimilation window [ t ]0,tS]Divided into (t) by time interval τs-t0) τ sub-time windows;
an initial background field determination unit for interpolating t from the re-analyzed image data covering the area to be studied0Background field of time xb,0(ii) a Or, for t0T is obtained by integrating the forecast model in the background field a period of time before the moment0Background field of time xb,0
Other background field determining unit for using numerical prediction model Mk(·)、t0Background field of time
Figure BDA0003168138100000055
And an overall error correction term
Figure BDA0003168138100000056
The following background simulation was performed to obtain an assimilation window [ t ]0,tS]Internal state variable
Figure BDA0003168138100000057
The background simulation result of (1);
Figure BDA0003168138100000058
wherein M isk(. is) from tk-1To tkThe numerical value prediction model of (1) is,
Figure BDA0003168138100000059
model M for numerical forecastingk(ii) a state variable of (c),
Figure BDA0003168138100000061
when k is t0Time, right side of formula
Figure BDA0003168138100000062
A prior state vector determination unit for all state variables
Figure BDA0003168138100000063
The background simulation results of (a) constitute the prior state vector.
Further, a sample set of the current assimilation window is obtained by the following units:
an initial set sample determining unit for obtaining a history sample and obtaining a set of set sample disturbances P by using a disturbance state variable method for the history samplex=(x'1,…,x'N);
Other aggregate sample determination units for using the numerical prediction model Mk(·)、t0Background field of time
Figure BDA0003168138100000064
And aggregate sample perturbation
Figure BDA0003168138100000065
An integrated simulation was performed to obtain an assimilation window [ t ]0,tS]Inner set sample xj,kThe simulation result of (2);
Figure BDA0003168138100000066
wherein, j is 1,2, …, N, when k is t0Time, right side of formula
Figure BDA0003168138100000067
A set sample determination unit for all set samples xj,kThe simulation results of (a) make up the aggregate sample.
Further, the obtaining of the final overall error correction term by minimizing the cost function includes:
an input unit for setting the initialized covariance matrix Q of the integral error correction term x' to 0 and xbPrior state vector
Figure BDA0003168138100000068
A priori simulation of the observation
Figure BDA0003168138100000069
Observation data yobs,kAnd covariance matrix R of observation errorskInputting the cost function l;
a first substitution unit for expressing the overall error correction term x' as a set sample perturbation PxLinear combination of (a) x ═ Pxβ, wherein Px=(x'1,…,x'N),β=(β1,…,βN) Is a coefficient vector;
a second substitution unit for using the sample covariance matrix
Figure BDA00031681381000000610
Covariance matrix Q instead of xb
An equivalence transformation unit for changing x' to PxBeta and
Figure BDA00031681381000000611
and substituting the cost function l, and performing equivalent deformation to obtain the following cost function taking beta as a new variable:
Figure BDA0003168138100000071
wherein the content of the first and second substances,
Figure BDA0003168138100000072
is x ═At 0 time
Figure BDA0003168138100000073
A solving unit, which is used for carrying out mathematical transformation on the cost function taking the beta as a new variable into a nonlinear least square form, solving by adopting a Gaussian-Newton iterative algorithm to obtain the beta, and obtaining the beta according to the x ═ Pxβ yields the final overall error correction term.
In a third aspect, the present invention provides a computer storage medium for numerical weather forecasting, comprising a memory for storing processor-executable instructions that, when executed by the processor, implement steps comprising the numerical weather forecasting method of the first aspect.
In a fourth aspect, the present invention provides an electronic device for numerical weather forecast, comprising at least one processor and a memory storing computer-executable instructions, which when executed by the processor, implement the steps of the numerical weather forecast method of the first aspect.
The invention has the following beneficial effects:
firstly, obtaining a prior value (background field) of a state vector by inputting meteorological data in a numerical forecasting mode, and simultaneously constructing a set sample of the state vector; inputting the prior state vector and the set sample into an observation operator to obtain prior value simulation observation and sample simulation observation; and inputting the obtained information and observation data into an i4DVar assimilation solving module, and obtaining an integral error correction term by minimizing a cost function.
The invention divides an assimilation window into a plurality of sub-time windows, and adds an integral error correction term at the initial time of each sub-time window, thereby improving the classic 4DVar method into the integral correction four-dimensional variation assimilation algorithm (i4DVar) of the invention.
Drawings
FIG. 1 is a flow chart of a numerical weather forecasting method of the present invention;
FIG. 2 is a schematic diagram of a numerical weather forecasting apparatus of the present invention;
FIG. 3 is a graph of the cumulative precipitation space distribution at 20/2016 to 03/2016 at 24 hours at 7/21/2016;
FIG. 4 is a plot of the score values for ETS and Far (under-reporting) for 24 hours of cumulative precipitation at 2016 (7/20/03) to 2016 (7/21/03) at different thresholds;
FIG. 5 is a graph comparing the vertical profile distribution of mean root mean square error for 6 hour numerical predictions from CTRL, 4DVar, i4DVar with conventional observations;
FIG. 6 is a graph comparing the vertical profile distribution of the mean root mean square error of 24 hour numerical predictions from CTRL, 4DVar, i4DVar with conventional observations.
Detailed Description
In order to make the technical problems, technical solutions and advantages of the present invention more apparent, the technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings and specific embodiments. It is to be understood that the described embodiments are merely exemplary of the invention, and not restrictive of the full scope of the invention. The components of embodiments of the present invention generally described and illustrated in the figures herein may be arranged and designed in a wide variety of different configurations. Thus, the following detailed description of the embodiments of the present invention, presented in the figures, is not intended to limit the scope of the invention, as claimed, but is merely representative of selected embodiments of the invention. All other embodiments, which can be derived by a person skilled in the art from the embodiments of the present invention without making any creative effort, shall fall within the protection scope of the present invention.
Example 1:
the embodiment of the invention provides a numerical weather forecasting method, as shown in fig. 1, the method comprises the following steps:
s100: and acquiring a prior state vector and a set sample of the current assimilation window.
The method comprises the steps of acquiring meteorological data covering an area to be researched, and acquiring a prior state vector and a set sample of a current assimilation window according to the acquired meteorological data.
The meteorological data is data generated by a global/regional atmospheric numerical prediction mode, mainly comprises elements such as three-dimensional wind speed, atmospheric pressure, atmospheric temperature, atmospheric humidity and the like and corresponding geographical position information thereof, and all the information can be interpolated into a four-dimensional variational assimilation algorithm which is applied by the invention and used for constructing overall correction through the regional or global atmospheric numerical prediction mode.
Variables to be assimilated form a state vector, the state vector comprises a latitudinal wind component U, a radial wind component V, a disturbance temperature T, a disturbance air pressure P, a water-vapor mixing ratio QVAPOR and the like, and the number of the variables in the state vector can be increased or reduced according to actual assimilation problems.
The meteorological data can be input into a data preprocessing WPS module of a numerical forecasting model WRF (or any other atmospheric numerical forecasting mode), and the following operations are performed: (1) determining a region of interest; (2) interpolating topographic data; (3) and interpolating the meteorological data to the forecast area. Data were prepared for numerical simulations using the WRF mode.
When constructing the prior state vector and the set sample, the assimilation window is divided into a plurality of sub-time windows according to a certain time interval, and an integral error correction term is added to the prior state vector and the set sample at the initial moment of each sub-time window.
One specific method for obtaining the prior state vector of the current assimilation window is as follows:
s101: will assimilate window [ t0,tS]Divided into (t) by time interval τs-t0) τ sub-time windows.
In general, the assimilation window t0To tsThe time interval of (a) is 6 hours (h), and the assimilation window can be divided into 6 sub-time windows according to the time interval of each hour.
S102: interpolating from the re-analyzed meteorological data covering the area to be studied to obtain t0Background field of time xb,0(ii) a Or, for t0Ambient field passing a period of time before a momentIntegrating the report model to obtain t0Background field of time xb,0
This step is used to generate an assimilation window [ t ]0,tS]At an initial time t0The method may be: (1) directly interpolated from the global re-analysis image data; (2) several hours in advance, integrated using the numerical mode WRF.
S103: using numerical prediction model Mk(·)、t0Background field of time
Figure BDA0003168138100000091
And an overall error correction term
Figure BDA0003168138100000092
The following background simulation was performed to obtain an assimilation window [ t ]0,tS]Internal state variable
Figure BDA0003168138100000093
The background simulation results of (1).
Figure BDA0003168138100000094
Wherein M isk(. is) from tk-1To tkThe numerical value prediction model of (1) is,
Figure BDA0003168138100000095
model M for numerical forecastingk(ii) a state variable of (c),
Figure BDA0003168138100000096
when k is t0Time, right side of formula
Figure BDA0003168138100000097
All of which are nx
This step is used to perform the construction of the background field at other times in the assimilation window, and the above formula (1) shows that, during the background simulation process, the background field at other times in the assimilation window is constructed
Figure BDA0003168138100000101
The time instant (i.e., the initial time instant of each sub-time window) is increased by the overall error correction term
Figure BDA0003168138100000102
Finally, the whole assimilation window [ t ] is obtained0,tS]Internal state variable
Figure BDA0003168138100000103
The background simulation results of (1).
Therefore, in the present invention, the background field at the current moment is the background field at the previous moment, and is integrated by the numerical prediction model (i.e. the background field at the current moment is integrated by the numerical prediction model)
Figure BDA0003168138100000104
) And if the current time is
Figure BDA0003168138100000105
Then, an integral error correction term is added to the background field of the current moment
Figure BDA0003168138100000106
And then, integrating through a numerical forecasting model to obtain a background field at the next moment, and so on.
S104: all state variables
Figure BDA0003168138100000107
The background simulation results of (a) constitute a prior state vector.
One specific implementation method for obtaining the current assimilation window set sample is as follows:
s105: obtaining a historical sample, and obtaining a group of set sample disturbances P by adopting a disturbance state variable method on the historical samplex=(x'1,…,x'N)。
In the present invention, in order to construct a better quality set of historical samples, a window is assimilated [ t [ [ t ]0,tS]Initial time t0The aggregate sample is generated from historical large aggregate samples to obtain sufficiently large samplesAnd (4) characteristics.
The historical large set sample can be constructed by two methods: (1) setting initial mode integration time by using a numerical prediction mode WRF, and performing long-time mode integration to obtain a historical large set sample; (2) and (4) using a numerical mode forecast result generated in the early stage as a historical large set sample.
Then, a group of set sample disturbances P are obtained by using a large historical set sample and adopting a disturbance state variable methodx,Px=(x'1,…,x'N) (ii) a Will PxAdded to the initial time t0Background field x ofb,0To obtain an initial time t0Aggregate sample of
Figure BDA00031681381000001010
Wherein (j ═ 1,2, …, N).
S106: using numerical prediction model Mk(·)、t0Background field of time
Figure BDA0003168138100000108
And aggregate sample perturbation
Figure BDA0003168138100000109
An integrated simulation was performed to obtain an assimilation window [ t ]0,tS]Inner set sample xj,kThe simulation result of (1).
Figure BDA0003168138100000111
Wherein M isk(. is) from tk-1To tkThe numerical value prediction model of (1) is,
Figure BDA0003168138100000112
model M for numerical forecastingk(ii) a state variable of (c),
Figure BDA0003168138100000113
when k is t0Time, right side of formula
Figure BDA0003168138100000114
This step is used to assimilate the aggregate samples x at other times within the windowj,kThe above equation (2) shows that, for each perturbation sample
Figure BDA0003168138100000115
In particular, in
Figure BDA0003168138100000116
The integral disturbance terms are added at all times
Figure BDA0003168138100000117
Finally, the whole assimilation window [ t ] is obtained0,tS]Inner set sample xj,kThe simulation result of (1).
S107: all set samples xj,kThe simulation results of (a) constitute aggregate samples.
S200: respectively inputting the prior state vector and the set sample into an observation operator to obtain a prior simulation observation result
Figure BDA0003168138100000118
And sample simulation observations
Figure BDA0003168138100000119
Hk(. is) tkObservation operator of a moment, the observation operator being to vector the prior state
Figure BDA00031681381000001110
Spatial mapping from a Pattern State variable to an observed variable yk(having a vector dimension of
Figure BDA00031681381000001111
)。
The invention can adopt Community Radial Transfer Model (CRTM) as an observation operator to obtain a simulated observation result. The CRTM is a fast radiation transmission mode developed by the united states satellite data assimilation center and is highly modular in that it assimilates the microwave and infrared observations of most polar and geostationary satellites.
S300: the obtained multi-time background field
Figure BDA00031681381000001112
(i.e., a prior state vector) and aggregate samples
Figure BDA00031681381000001113
Inputting the observation operator CRTM to obtain the observation of the brightness and temperature of the simulated satellite
Figure BDA00031681381000001114
And
Figure BDA00031681381000001115
(i.e., a priori simulated observations and sample simulated observations), i.e.
Figure BDA00031681381000001116
While
Figure BDA00031681381000001117
From this, a simulated observation disturbance y 'was obtained'k,j=yk,j-yb,kThereby obtaining Py,k=(y'k,1,y'k,2,…,y'k,N) And y'obs,k=yobs,k-yb,k
S400: constructing a cost function of an overall correction four-dimensional variational assimilation algorithm (abbreviated as i4DVar algorithm), initializing an overall error correction term, taking a prior state vector, a set sample, a prior simulation observation result, a sample simulation observation result and obtained observation data as known information, and obtaining a final overall error correction term by minimizing the cost function.
The i4DVar method is a method of simultaneously performing initial error and mode error revision, and requires error revision at multiple times of an assimilation window (i.e., the initial time of each sub-window), not just the initial time of the entire assimilation window. The cost function l of the overall corrected four-dimensional variational assimilation algorithm (i4DVar) is therefore:
Figure BDA0003168138100000121
wherein x 'is an integral error correction term, x' follows a normal distribution, QbA covariance matrix of x' (also called background error covariance matrix),
Figure BDA0003168138100000122
in order to be a vector of a-priori states,
Figure BDA0003168138100000123
is a prior state vector
Figure BDA0003168138100000124
Input to observation operator Hk(. A) derived prior simulated observation, yobs,kAs observation data (i.e. t)kObservation vector of time), RkCovariance matrix (i.e. t) as observation errorkThe observed error covariance matrix at the time of day).
yobs,kThe observation data may be satellite observation data or other observation data, which is downloaded from a data website and mainly includes conventional data, radar data, satellite data, and the like, and the observation error is obtained from an observation data provider or estimated by using a statistical method.
The i4DVar method of the invention is used for observing data yobs,kAnd the observation error covariance matrix RkCovariance matrix Q with background errorbUnder the common constraint of (a) to make an optimal estimate of the overall correction term x'. The input quantity of the i4DVar data assimilation algorithm comprises a priori state vector, an integral error correction term x' is initialized to be 0, and a numerical prediction mode M is adoptedk(ii) performing simulated prediction
Figure BDA0003168138100000125
Observation vector yobs,kAnd its observation error covariance matrix RkAnd a background error covariance matrix B. Inputting the above dataAnd carrying out assimilation analysis on the i4DVar by using a cost function of the overall corrected four-dimensional variable assimilation algorithm.
In general, solving the above i4DVar algorithm requires the observation operator HkTangent operator of
Figure BDA0003168138100000126
With its accompanying operator
Figure BDA0003168138100000127
And numerical prediction mode MkTangential mode of (c)'kAnd companion mode
Figure BDA0003168138100000128
The development difficulty is extremely large. Therefore, the i4DVar idea is coupled with a nonlinear least square set four-dimensional variation assimilation method (NLS-4DVar) in the prior art, and a Gaussian-Newton iterative algorithm is adopted for solving.
The specific solving process comprises the following steps:
s401: covariance matrix Q of initialized integral error correction term x' 0, xbPrior state vector
Figure BDA0003168138100000131
A priori simulation of the observation
Figure BDA0003168138100000132
Observation data yobs,kAnd covariance matrix R of observation errorskThe cost function l is input.
S402: expressing the integral error correction term x' as the set sample disturbance PxLinear combination of (a) x ═ Pxβ, wherein Px=(x'1,…,x'N) Is a perturbation matrix of pattern space set samples, beta ═ beta1,…,βN) Is a coefficient vector.
S403: using sample covariance matrices
Figure BDA0003168138100000133
Covariance matrix Q instead of xb
S404: changing x' to PxBeta and
Figure BDA0003168138100000134
and (3) carrying out equivalent transformation on the cost function l of the formula (3) to obtain the following cost function taking beta as a new variable:
Figure BDA0003168138100000135
wherein the content of the first and second substances,
Figure BDA0003168138100000136
when x' is 0
Figure BDA0003168138100000137
S405: mathematically transforming the cost function (formula (4)) taking the beta as a new variable into a nonlinear least square form, and solving by adopting a Gaussian-Newton iterative algorithm to obtain the beta;
Figure BDA0003168138100000138
where i is the number of iterations (i ═ 1,2, … i)max),βiIs the value of the ith iteration of β.
And according to x' ═ Pxβ yields the final overall error correction term.
In practice, in order to alleviate the problem of sampling error caused by the limited number of the sample sets, a localization scheme is adopted in the solution process, as follows:
Figure BDA0003168138100000139
Figure BDA0003168138100000141
Figure BDA0003168138100000142
Figure BDA0003168138100000143
wherein x is'’iIs the value of the ith iteration of x'. rhoxRevising the matrix, ρ, for the pattern space localizationy,kThe matrix is revised for localization of the observation space. RhoxConstructed from the extent of the region of interest and the number of grid points, py,kUsing a spatial three-dimensional coordinate pair rhoxAnd (4) obtaining by interpolation.
When i ═ imaxThen, the cycle is exited to obtain
Figure BDA0003168138100000144
S500: and adding the final integral error correction term to the prior state vector, and performing mode integration to obtain a weather forecast result of the current assimilation window.
Adding the obtained integral error correction term x' to the prior state vector
Figure BDA0003168138100000145
Mode integration is carried out, so that a better wind field, a better temperature field and a better humidity field can be obtained, and precipitation can be predicted more accurately.
Firstly, obtaining a prior value (background field) of a state vector by inputting meteorological data in a numerical forecasting mode, and simultaneously constructing a set sample of the state vector; inputting the prior state vector and the set sample into an observation operator to obtain prior value simulation observation and sample simulation observation; and inputting the obtained information and observation data into an i4DVar assimilation solving module, and obtaining an integral error correction term by minimizing a cost function.
The invention divides an assimilation window into a plurality of sub-time windows, and adds an integral error correction term at the initial time of each sub-time window, thereby improving the classic 4DVar method into the integral correction four-dimensional variation assimilation algorithm (i4DVar) of the invention.
And the set simulation is introduced into the solution of the i4DVar, the cost function of the i4DVar is minimized by using a Gaussian-Newton iteration method, the calculation of a tangential model and a adjoint model is well avoided through simple mathematical derivation, and the calculation is simpler.
The steps are an iterative solution process in each assimilation window, and the cyclic assimilation process also needs to transmit information between assimilation windows, including the updating of a set sample and a background field, and the specific steps are as follows:
1. updating of next synchronization window background field
Initial time x of next synchronization windowb,0There are two methods of generation: (1) directly interpolated from the global re-analysis image data; (2) obtaining optimal integral correction terms x' and x in current assimilation windowb,0Simulation was performed as follows
Figure BDA0003168138100000151
Then the next synchronization window
Figure BDA0003168138100000152
And constructing a background field corresponding to multiple moments according to the method in the assimilation window in S103.
2. Update of next synchronization window set sample perturbation
Aggregate sample perturbation with current assimilation window
Figure BDA0003168138100000153
Py,kAnd waiting for information, and adopting the following set sample perturbation scheme for preserving dispersion:
Figure BDA0003168138100000154
Figure BDA0003168138100000155
X2=Λ-1/2ZTPy
Figure BDA0003168138100000156
Figure BDA0003168138100000157
Figure BDA0003168138100000161
Figure BDA0003168138100000162
Λ=[∑2+(N-1)I]
Φ is a random orthogonal matrix.
Obtaining the set sample perturbation of the initial time of the next synchronization window, adding to the x of the next synchronization windowb,0An initial set of samples is obtained. Then, a multi-time corresponding set sample is constructed according to a method (S106) in the assimilation window.
And after the information transmission between the assimilation windows is finished, calculating the next assimilation window. This is repeated.
The effect of the invention is explained in detail below with a specific test example:
the assimilation performance of the i4DVar method was evaluated by selecting a one-week assimilation test from 2016, 7, 16, and 03 hours, to 2016, 7, 23, and 03 days. This time period includes a major storm in the North China (35-43N, 113-122E) that occurs between 18 and 21 days 7-2016. And there are two heavy rainfall centers, the first one in taihang mountain and the second one in central south of beijing.
The assimilation performance of i4DVar was demonstrated by comparing the Control Test (CTRL) with the prediction results of 4 DVar. The first is the verification of the precipitation forecast result. The evaluation indicators are the spatial distribution of precipitation, the Root Mean Square Error (RMSE) and correlation (CC) of precipitation forecasts and observations, and the ETS, Far scores.
Fig. 3 is a graph of the cumulative precipitation space distribution in mm from 20/03/2016 to 03/2016/7/21/2016 for 24 hours. Wherein, (a) is the precipitation condition of the hourly national level observation station; (b) a WRF mode prediction result which is a CTRL test of unassimilated satellite data and takes ERA-Interim of ECMWF as an initial value; (c) the prediction result of the WRF mode taking the analysis field of the satellite data of 6 hours from 2016, 7, 19, 21 to 2016, 7, 20, 03 as an initial value is assimilated by the 4DVar method; (d) the prediction result of WRF mode using the i4DVar method to assimilate the analysis field of the same satellite data as in (c) as the initial value.
It can be seen from (b) that in this rainstorm event, the numerical model WRF is weak in precipitation simulation in the south of beijing, strong in precipitation simulation at the junction of inner mongolia, north river and liaoning, and the distribution of precipitation in liaoning is far from the observation of (a). From (c), the numerical prediction result of the analysis field obtained by 4DVar optimizes the precipitation zone of Liaoning. From (d), the precipitation forecast results of the analysis field obtained by the i4DVar method are closer to the actual condition (a). Relative to CTRL, the i4DVar results enhanced precipitation in the south of Beijing, and reduced precipitation intensity at the boundaries of inner Mongolia, Hebei, Liaoning. Meanwhile, false precipitation in the Shandong is weakened. The results show that the i4DVar method can well absorb satellite data and improve a rainfall forecast field. Comparison with 4DVar shows that the i4DVar method has better precipitation improvement capability than the 4DVar method in the precipitation event.
Fig. 4 shows the scoring values of the precipitation of the corresponding precipitation example of fig. 3 at different threshold values. For precipitation forecast, the bigger the ETS (equivalent thread Score) Score value is, the smaller the Far Score value is, indicating that the system has stronger forecast capacity for precipitation. The thresholds of 24 hours of accumulated precipitation are 25mm,50mm,100mm,150mm and 200mm respectively. As can be seen from the results of ETS, the results of ETS scoring for i4DVar are superior to those for CTRL and 4DVar at thresholds of 100mm,150mm, and 200mm, respectively. For the false negative score, the FAR value for i4DVar is slightly greater than the results for CTRL and 4DVar only at the 50mm threshold. The results of ETS and FAR can further show that the i4DVar method has better assimilation performance on satellite data, so that the forecast of strong rainfall is improved.
Table 1 shows the Root Mean Square Error (RMSE) and Correlation Coefficient (CC) of the cumulative precipitation forecast and precipitation observation from 2016, 7, 20, 03 to 2016, 7, 21, 03, and it can be seen from table 1 that i4DVar has a smaller root mean square error and a larger correlation coefficient. Further indicating that the assimilation performance of the i4DVar method is superior to that of the 4DVar method in this case.
TABLE 1
Figure BDA0003168138100000181
Fig. 5 and 6 are graphs comparing vertical profile distributions of mean RMSE from CTRL, 4DVar, i4DVar at 6 hours (from 2016 (7/20/03) to 2016 (7/20/09) and 24 hours (from 2016 (7/20/21) to 2016 (7/21/03)) numerical predictions and 24 hours (from 2016 (7/20/21) respectively) versus conventional observations. In fig. 5 and 6, (a), (b), (c), and (d) correspond to u wind, v wind, air temperature, and humidity variables, respectively.
As can be seen from FIG. 5, the i4DVar method improves significantly better than 4DVar for u/v wind variables, especially at the higher (below 400hPa) and lower (above 700hPa) levels of the mode. For air temperature variables, the conclusion is similar to that of a wind field, and the assimilation advantage of the i4DVar method is obvious. For humidity variables, above 800hPa, the improvement in 4DVar is better than i4 DVar.
As can be seen from FIG. 6, the conclusion for the u/v wind variable is the same as that of FIG. 5, except that the improvement of the region SNAP is more obvious. For the air temperature variable, the improvement of 4DVar is better than that of i4DVar only at 900hPa, and the i4DVar has a more obvious advantage (smaller root mean square error) at other mode layers. For the Humidity variable, the advantage of i4DVar is evident for most of the mode layers. At 700hPa, the root mean square error of i4DVar is large. In summary, in this case, the i4DVar method can better improve the pattern prediction capability for most of the pattern layers of all variables.
Example 2:
an embodiment of the present invention provides a digital weather forecast apparatus, as shown in fig. 2, the retrofit apparatus includes:
the acquisition module 1 is configured to acquire a prior state vector and a set sample of a current assimilation window.
The assimilation window is divided into a plurality of sub-time windows according to a certain time interval, and an integral error correction term is added to the prior state vector and the set sample at the initial moment of each sub-time window.
And the simulation module 2 is used for respectively inputting the prior state vector and the set sample into the observation operator to obtain a prior simulation observation result and a sample simulation observation result.
And the solving module 3 is used for constructing a cost function of the overall correction four-dimensional variational assimilation algorithm, initializing an overall error correction term, and obtaining a final overall error correction term by minimizing the cost function by taking the prior state vector, the set sample, the prior simulation observation result, the sample simulation observation result and the obtained observation data as known information.
And the forecasting module 4 is used for adding the final integral error correction term to the prior state vector, and performing mode integration to obtain a weather forecasting result of the current assimilation window.
The cost function l of the above-mentioned overall corrected four-dimensional variable assimilation algorithm is:
Figure BDA0003168138100000191
wherein x' is an integral error correction term, QbIs a covariance matrix of x',
Figure BDA0003168138100000192
is a priori state directionThe amount of the compound (A) is,
Figure BDA0003168138100000193
is a prior state vector
Figure BDA0003168138100000194
Input to observation operator Hk(. A) derived prior simulated observation, yobs,kFor observation of data, RkIs a covariance matrix of the observed errors.
The prior state vector of the current assimilation window is obtained through the following units:
a dividing unit for dividing the assimilation window [ t ]0,tS]Divided into (t) by time interval τs-t0) τ sub-time windows.
An initial background field determination unit for interpolating t from the re-analyzed image data covering the area to be studied0Background field of time xb,0(ii) a Or, for t0T is obtained by integrating the forecast model in the background field a period of time before the moment0Background field of time xb,0
Other background field determining unit for using numerical prediction model Mk(·)、t0Background field of time
Figure BDA0003168138100000195
And an overall error correction term
Figure BDA0003168138100000196
The following background simulation was performed to obtain an assimilation window [ t ]0,tS]Internal state variable
Figure BDA0003168138100000197
The background simulation results of (1).
Figure BDA0003168138100000198
Wherein M isk(. is) from tk-1To tkNumerical forecasting ofThe model is a model of a human body,
Figure BDA0003168138100000199
model M for numerical forecastingk(ii) a state variable of (c),
Figure BDA00031681381000001910
when k is t0Time, right side of formula
Figure BDA00031681381000001911
A prior state vector determination unit for all state variables
Figure BDA00031681381000001912
The background simulation results of (a) constitute a prior state vector.
Acquiring a set sample of a current assimilation window by:
an initial set sample determining unit for obtaining a history sample and obtaining a set of set sample disturbances P by using a disturbance state variable method for the history samplex=(x'1,…,x'N)。
Other aggregate sample determination units for using the numerical prediction model Mk(·)、t0Background field of time
Figure BDA0003168138100000201
And aggregate sample perturbation
Figure BDA0003168138100000202
An integrated simulation was performed to obtain an assimilation window [ t ]0,tS]Inner set sample xj,kThe simulation result of (1).
Figure BDA0003168138100000203
Wherein j is 1,2, …, N, when k is t0Time, right side of formula
Figure BDA0003168138100000204
A set sample determination unit for all set samples xj,kThe simulation results of (a) constitute aggregate samples.
Further, the obtaining of the final overall error correction term by minimizing the cost function in the present invention includes:
an input unit for setting the initialized covariance matrix Q of the integral error correction term x' to 0 and xbPrior state vector
Figure BDA0003168138100000205
A priori simulation of the observation
Figure BDA0003168138100000206
Observation data yobs,kAnd covariance matrix R of observation errorskThe cost function l is input.
A first substitution unit for expressing the overall error correction term x' as a set sample perturbation PxLinear combination of (a) x ═ Pxβ, wherein Px=(x'1,…,x'N),β=(β1,…,βN) Is a coefficient vector.
A second substitution unit for using the sample covariance matrix
Figure BDA0003168138100000207
Covariance matrix Q instead of xb
An equivalence transformation unit for changing x' to PxBeta and
Figure BDA0003168138100000208
and (3) carrying in a cost function l, and carrying out equivalent transformation to obtain the following cost function taking beta as a new variable:
Figure BDA0003168138100000209
wherein the content of the first and second substances,
Figure BDA00031681381000002010
when x' is 0
Figure BDA00031681381000002011
A solving unit, which is used for carrying out mathematical transformation on the cost function taking the beta as a new variable into a nonlinear least square form, solving by adopting a Gaussian-Newton iterative algorithm to obtain the beta, and obtaining the beta according to the x ═ Pxβ yields the final overall error correction term.
The device provided by the embodiment of the present invention has the same implementation principle and technical effect as the method embodiment 1, and for the sake of brief description, reference may be made to the corresponding content in the method embodiment 1 for the part where the embodiment of the device is not mentioned. It can be clearly understood by those skilled in the art that, for convenience and brevity of description, the specific working processes of the apparatus and the unit described above may all refer to the corresponding processes in the above method embodiment 1, and are not described herein again.
Example 3:
the method of the embodiment 1 provided by the present invention can implement the service logic through a computer program and record the service logic on a storage medium, and the storage medium can be read and executed by a computer, so as to implement the effect of the solution described in the embodiment 1 of the present specification. Accordingly, the present invention also provides a computer storage medium for numerical weather forecasting, comprising a memory for storing processor-executable instructions that, when executed by a processor, implement steps comprising the numerical weather forecasting method of embodiment 1.
The storage medium may include a physical device for storing information, and typically, the information is digitized and then stored using an electrical, magnetic, or optical media. The storage medium may include: devices that store information using electrical energy, such as various types of memory, e.g., RAM, ROM, etc.; devices that store information using magnetic energy, such as hard disks, floppy disks, tapes, core memories, bubble memories, and usb disks; devices that store information optically, such as CDs or DVDs. Of course, there are other ways of storing media that can be read, such as quantum memory, graphene memory, and so forth.
The above description of the storage medium according to method embodiment 1 may also include other implementations. The specific implementation manner may refer to the description of the related method embodiment 1, and is not described in detail here.
Example 4:
the invention also provides an electronic device for numerical weather forecast, which can be a single computer, and can also comprise an actual operation device and the like using one or more methods or one or more embodiment devices of the specification. The electronic device for numerical weather forecast may comprise at least one processor and a memory storing computer-executable instructions, which when executed by the processor, implement the steps of the numerical weather forecast method described in any one or more of the above embodiments 1.
The description of the electronic device according to the method or apparatus embodiment may also include other implementation manners, and a specific implementation manner may refer to the description of related method embodiment 1, which is not described herein in detail.
It should be noted that, the above-mentioned apparatus or system in this specification may also include other implementation manners according to the description of the related method embodiment, and a specific implementation manner may refer to the description of the method embodiment, which is not described herein in detail. The embodiments in the present specification are described in a progressive manner, and the same and similar parts among the embodiments are referred to each other, and each embodiment focuses on the differences from the other embodiments. In particular, for the hardware + program class, storage medium + program embodiment, since it is basically similar to the method embodiment, the description is relatively simple, and for the relevant points, refer to the partial description of the method embodiment.
The foregoing description has been directed to specific embodiments of this disclosure. Other embodiments are within the scope of the following claims. In some cases, the actions or steps recited in the claims may be performed in a different order than in the embodiments and still achieve desirable results. In addition, the processes depicted in the accompanying figures do not necessarily require the particular order shown, or sequential order, to achieve desirable results. In some embodiments, multitasking and parallel processing may also be possible or may be advantageous.
The systems, devices, modules or units illustrated in the above embodiments may be implemented by a computer chip or an entity, or by a product with certain functions. One typical implementing electronic device is a computer. In particular, the computer may be, for example, a personal computer, a laptop computer, a vehicle-mounted human-computer interaction device, a cellular telephone, a camera phone, a smart phone, a personal digital assistant, a media player, a navigation device, an email device, a game console, a tablet computer, a wearable device, or a combination of any of these devices.
For convenience of description, the above devices are described as being divided into various modules by functions, and are described separately. Of course, when implementing one or more of the present description, the functions of each module may be implemented in one or more software and/or hardware, or a module implementing the same function may be implemented by a combination of multiple sub-modules or sub-units, etc. The above-described embodiments of the apparatus are merely illustrative, and for example, the division of the units is only one logical division, and other divisions may be realized in practice, for example, a plurality of units or components may be combined or integrated into another system, or some features may be omitted, or not executed. In addition, the shown or discussed mutual coupling or direct coupling or communication connection may be an indirect coupling or communication connection through some interfaces, devices or units, and may be in an electrical, mechanical or other form.
Those skilled in the art will also appreciate that, in addition to implementing the controller as pure computer readable program code, the same functionality can be implemented by logically programming method steps such that the controller is in the form of logic gates, switches, application specific integrated circuits, programmable logic controllers, embedded microcontrollers and the like. Such a controller may therefore be considered as a hardware component, and the means included therein for performing the various functions may also be considered as a structure within the hardware component. Or even means for performing the functions may be regarded as being both a software module for performing the method and a structure within a hardware component.
The present invention is described with reference to flowchart illustrations and/or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of the invention. It will be understood that each flow and/or block of the flow diagrams and/or block diagrams, and combinations of flows and/or blocks in the flow diagrams and/or block diagrams, can be implemented by computer program instructions. These computer program instructions may be provided to a processor of a general purpose computer, special purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, create means for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.
These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable data processing apparatus to function in a particular manner, such that the instructions stored in the computer-readable memory produce an article of manufacture including instruction means which implement the function specified in the flowchart flow or flows and/or block diagram block or blocks.
These computer program instructions may also be loaded onto a computer or other programmable data processing apparatus to cause a series of operational steps to be performed on the computer or other programmable apparatus to produce a computer implemented process such that the instructions which execute on the computer or other programmable apparatus provide steps for implementing the functions specified in the flowchart flow or flows and/or block diagram block or blocks.
In a typical configuration, a computing device includes one or more processors (CPUs), input/output interfaces, network interfaces, and memory.
It should also be noted that the terms "comprises," "comprising," or any other variation thereof, are intended to cover a non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements does not include only those elements but may include other elements not expressly listed or inherent to such process, method, article, or apparatus. Without further limitation, an element defined by the phrase "comprising an … …" does not exclude the presence of other like elements in a process, method or apparatus that comprises the element.
As will be appreciated by one skilled in the art, one or more embodiments of the present description may be provided as a method, system, or computer program product. Accordingly, one or more embodiments of the present description may take the form of an entirely hardware embodiment, an entirely software embodiment or an embodiment combining software and hardware aspects. Furthermore, one or more embodiments of the present description may take the form of a computer program product embodied on one or more computer-usable storage media (including, but not limited to, disk storage, CD-ROM, optical storage, and the like) having computer-usable program code embodied therein.
One or more embodiments of the present description may be described in the general context of computer-executable instructions, such as program modules, being executed by a computer. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform particular tasks or implement particular abstract data types. One or more embodiments of the present specification can also be practiced in distributed computing environments where tasks are performed by remote processing devices that are linked through a communications network. In a distributed computing environment, program modules may be located in both local and remote computer storage media including memory storage devices.
The embodiments in the present specification are described in a progressive manner, and the same and similar parts among the embodiments are referred to each other, and each embodiment focuses on the differences from the other embodiments. In particular, for the system embodiment, since it is substantially similar to the method embodiment, the description is simple, and for the relevant points, reference may be made to the partial description of the method embodiment. In the description of the specification, reference to the description of the term "one embodiment," "some embodiments," "an example," "a specific example," or "some examples," etc., means that a particular feature, structure, material, or characteristic described in connection with the embodiment or example is included in at least one embodiment or example of the specification. In this specification, the schematic representations of the terms used above are not necessarily intended to refer to the same embodiment or example. Furthermore, the particular features, structures, materials, or characteristics described may be combined in any suitable manner in any one or more embodiments or examples. Furthermore, various embodiments or examples and features of different embodiments or examples described in this specification can be combined and combined by one skilled in the art without contradiction.
Finally, it should be noted that: the above-mentioned embodiments are only specific embodiments of the present invention, which are used for illustrating the technical solutions of the present invention and not for limiting the same, and the protection scope of the present invention is not limited thereto, although the present invention is described in detail with reference to the foregoing embodiments, those skilled in the art should understand that: any person skilled in the art can modify or easily conceive the technical solutions described in the foregoing embodiments or equivalent substitutes for some technical features within the technical scope of the present disclosure; such modifications, changes or substitutions do not depart from the spirit and scope of the present invention in its spirit and scope. Are intended to be covered by the scope of the present invention. Therefore, the protection scope of the present invention shall be subject to the protection scope of the appended claims.

Claims (4)

1. A method of numerical weather forecasting, the method comprising:
acquiring a prior state vector and a set sample of a current assimilation window;
the assimilation window is divided into a plurality of sub-time windows according to a certain time interval, and an integral error correction term is added to the prior state vector and the set sample at the initial moment of each sub-time window;
respectively inputting the prior state vector and the set sample into an observation operator to obtain a prior simulation observation result and a sample simulation observation result;
constructing a cost function of an overall correction four-dimensional variational assimilation algorithm, initializing an overall error correction term, taking the prior state vector, the set sample, the prior simulation observation result, the sample simulation observation result and the obtained observation data as known information, and obtaining a final overall error correction term by minimizing the cost function;
adding the final integral error correction term to the prior state vector, and performing mode integration to obtain a weather forecast result of the current assimilation window;
wherein, the cost function l of the overall corrected four-dimensional variable assimilation algorithm is as follows:
Figure FDA0003517568330000011
wherein x' is an integral error correction term, QbIs a covariance matrix of x',
Figure FDA0003517568330000012
in order to be a vector of a-priori states,
Figure FDA0003517568330000013
is a prior state vector
Figure FDA0003517568330000014
Input to observation operator Hk(. A) derived prior simulated observation, yobs,kFor observation of data, RkA covariance matrix that is an observation error;
obtaining a prior state vector of a current assimilation window by the following method:
will assimilate window [ t0,tS]Divided into (t) by time interval τs-t0) τ sub-time windows;
interpolating from the re-analyzed meteorological data covering the area to be studied to obtain t0Background field of time xb,0(ii) a Or, for t0T is obtained by integrating the forecast model in the background field a period of time before the moment0Background field of time xb,0
Using numerical prediction model Mk(·)、t0Background field of time
Figure FDA0003517568330000015
And an overall error correction term
Figure FDA0003517568330000016
The following background simulation was performed to obtain an assimilation window [ t ]0,tS]Internal state variable
Figure FDA0003517568330000017
The background simulation result of (1);
Figure FDA0003517568330000021
wherein M isk(. is) from tk-1To tkThe numerical value prediction model of (1) is,
Figure FDA0003517568330000022
model M for numerical forecastingk(ii) a state variable of (c),
Figure FDA0003517568330000023
when k is t0Time, right side of formula
Figure FDA0003517568330000024
All state variables
Figure FDA0003517568330000025
The background simulation result of (a) constitutes the prior state directionAn amount;
acquiring a set sample of a current assimilation window by the following method:
obtaining a historical sample, and obtaining a group of set sample disturbances P by adopting a disturbance state variable method on the historical samplex=(x′1,…,x′N);
Using numerical prediction model Mk(·)、t0Background field of time
Figure FDA0003517568330000026
And aggregate sample perturbation
Figure FDA0003517568330000027
An integrated simulation was performed to obtain an assimilation window [ t ]0,tS]Inner set sample xj,kThe simulation result of (2);
Figure FDA0003517568330000028
wherein j is 1,2, …, N, when k is t0Time, right side of formula
Figure FDA0003517568330000029
All set samples xj,kThe simulation results of (a) constitute the collective sample;
the step of obtaining a final overall error correction term by minimizing the cost function includes:
covariance matrix Q of initialized integral error correction term x' 0, xbPrior state vector
Figure FDA00035175683300000210
A priori simulation of the observation
Figure FDA00035175683300000211
Observation data yobs,kAnd covariance matrix R of observation errorskInputting the cost function l;
expressing the integral error correction term x' as set sample disturbance PxLinear combination of (a) x ═ Pxβ, wherein Px=(x′1,…,x′N),β=(β1,…,βN) Is a coefficient vector;
using sample covariance matrices
Figure FDA00035175683300000212
Covariance matrix Q instead of xb
Changing x' to PxBeta and
Figure FDA00035175683300000213
and substituting the cost function l, and performing equivalent deformation to obtain the following cost function taking beta as a new variable:
Figure FDA0003517568330000031
wherein the content of the first and second substances,
Figure FDA0003517568330000032
Figure FDA0003517568330000033
when x' is 0
Figure FDA0003517568330000034
Mathematically transforming the cost function taking the beta as a new variable into a nonlinear least square form, solving by adopting a Gaussian-Newton iterative algorithm to obtain the beta, and obtaining the beta according to the x ═ Pxβ yields the final overall error correction term.
2. A numerical weather forecast apparatus, characterized in that said apparatus comprises:
the acquisition module is used for acquiring a prior state vector and a set sample of a current assimilation window;
the assimilation window is divided into a plurality of sub-time windows according to a certain time interval, and an integral error correction term is added to the prior state vector and the set sample at the initial moment of each sub-time window;
the simulation module is used for respectively inputting the prior state vector and the set sample into an observation operator to obtain a prior simulation observation result and a sample simulation observation result;
the solving module is used for constructing a cost function of an overall correction four-dimensional variational assimilation algorithm, initializing an overall error correction term, and obtaining a final overall error correction term by minimizing the cost function by taking the prior state vector, the set sample, the prior simulation observation result, the sample simulation observation result and the obtained observation data as known information;
the forecasting module is used for adding the final integral error correction term to the prior state vector and carrying out mode integration to obtain a weather forecasting result of the current assimilation window;
wherein, the cost function l of the overall corrected four-dimensional variable assimilation algorithm is as follows:
Figure FDA0003517568330000035
wherein x' is an integral error correction term, QbIs a covariance matrix of x',
Figure FDA0003517568330000036
in order to be a vector of a-priori states,
Figure FDA0003517568330000037
is a prior state vector
Figure FDA0003517568330000038
Input to observation operator Hk(. A) derived prior simulated observation, yobs,kTo count the observationsAccording to RkA covariance matrix that is an observation error;
obtaining a prior state vector of a current assimilation window through the following units:
a dividing unit for dividing the assimilation window [ t ]0,tS]Divided into (t) by time interval τs-t0) τ sub-time windows;
an initial background field determination unit for interpolating t from the re-analyzed image data covering the area to be studied0Background field of time xb,0(ii) a Or, for t0T is obtained by integrating the forecast model in the background field a period of time before the moment0Background field of time xb,0
Other background field determining unit for using numerical prediction model Mk(·)、t0Background field of time
Figure FDA0003517568330000041
And an overall error correction term
Figure FDA0003517568330000042
The following background simulation was performed to obtain an assimilation window [ t ]0,tS]Internal state variable
Figure FDA0003517568330000043
The background simulation result of (1);
Figure FDA0003517568330000044
wherein M isk(. is) from tk-1To tkThe numerical value prediction model of (1) is,
Figure FDA0003517568330000045
model M for numerical forecastingk(ii) a state variable of (c),
Figure FDA0003517568330000046
when k is t0Time, right side of formula
Figure FDA0003517568330000047
A prior state vector determination unit for all state variables
Figure FDA0003517568330000048
The background simulation result of (a) constitutes the prior state vector;
acquiring a set sample of a current assimilation window by:
an initial set sample determining unit for obtaining a history sample and obtaining a set of set sample disturbances P by using a disturbance state variable method for the history samplex=(x′1,…,x′N);
Other aggregate sample determination units for using the numerical prediction model Mk(·)、t0Background field of time
Figure FDA0003517568330000049
And aggregate sample perturbation
Figure FDA00035175683300000410
An integrated simulation was performed to obtain an assimilation window [ t ]0,tS]Inner set sample xj,kThe simulation result of (2);
Figure FDA00035175683300000411
wherein j is 1,2, …, N, when k is t0Time, right side of formula
Figure FDA00035175683300000412
A set sample determination unit for all set samples xj,kThe simulation results of (a) constitute the collective sample;
the step of obtaining a final overall error correction term by minimizing the cost function includes:
an input unit for setting the initialized covariance matrix Q of the integral error correction term x' to 0 and xbPrior state vector
Figure FDA0003517568330000051
A priori simulation of the observation
Figure FDA0003517568330000052
Observation data yobs,kAnd covariance matrix R of observation errorskInputting the cost function l;
a first substitution unit for expressing the overall error correction term x' as a set sample perturbation PxLinear combination of (a) x ═ Pxβ, wherein Px=(x′1,…,x′N),β=(β1,…,βN) Is a coefficient vector;
a second substitution unit for using the sample covariance matrix
Figure FDA0003517568330000053
Covariance matrix Q instead of xb
An equivalence transformation unit for changing x' to PxBeta and
Figure FDA0003517568330000054
and substituting the cost function l, and performing equivalent deformation to obtain the following cost function taking beta as a new variable:
Figure FDA0003517568330000055
wherein the content of the first and second substances,
Figure FDA0003517568330000056
Figure FDA0003517568330000057
when x' is 0
Figure FDA0003517568330000058
A solving unit, which is used for carrying out mathematical transformation on the cost function taking the beta as a new variable into a nonlinear least square form, solving by adopting a Gaussian-Newton iterative algorithm to obtain the beta, and obtaining the beta according to the x ═ Pxβ yields the final overall error correction term.
3. A computer storage medium for numerical weather forecasting, comprising a memory for storing processor-executable instructions that, when executed by the processor, implement steps comprising the numerical weather forecasting method of claim 1.
4. An electronic device for numerical weather forecasting, comprising at least one processor and a memory storing computer-executable instructions that, when executed by the processor, implement the steps of the numerical weather forecasting method of claim 1.
CN202110810865.4A 2021-07-19 2021-07-19 Numerical weather forecasting method and device, computer storage medium and electronic equipment Active CN113568067B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110810865.4A CN113568067B (en) 2021-07-19 2021-07-19 Numerical weather forecasting method and device, computer storage medium and electronic equipment

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110810865.4A CN113568067B (en) 2021-07-19 2021-07-19 Numerical weather forecasting method and device, computer storage medium and electronic equipment

Publications (2)

Publication Number Publication Date
CN113568067A CN113568067A (en) 2021-10-29
CN113568067B true CN113568067B (en) 2022-04-05

Family

ID=78165385

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110810865.4A Active CN113568067B (en) 2021-07-19 2021-07-19 Numerical weather forecasting method and device, computer storage medium and electronic equipment

Country Status (1)

Country Link
CN (1) CN113568067B (en)

Families Citing this family (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN115511192B (en) * 2022-09-30 2023-07-14 中国科学院西北生态环境资源研究院 Precipitation prediction method and system based on lightning data assimilation

Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110045439A (en) * 2018-01-15 2019-07-23 钱维宏 Numerical weather prediction model system based on atmospheric variable Transient Eddy equation group

Family Cites Families (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN104992071B (en) * 2015-07-17 2018-07-13 南京信息工程大学 A kind of initial disturbance method based on pooling information assimilation technique
US20180275621A1 (en) * 2017-03-24 2018-09-27 Mitsubishi Electric Research Laboratories, Inc. Model Predictive Control with Uncertainties
CN109782374B (en) * 2019-01-15 2021-10-01 北京华云星地通科技有限公司 Method and device for optimizing numerical weather forecast through assimilation and inversion of water vapor content
CN110502849B (en) * 2019-08-27 2023-05-30 中国气象局广州热带海洋气象研究所(广东省气象科学研究所) Disturbance mode construction method applied to four-dimensional variation assimilation system
CN111859249B (en) * 2020-06-08 2022-06-14 天津大学 Ocean numerical forecasting method based on analytical four-dimensional set variation

Patent Citations (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN110045439A (en) * 2018-01-15 2019-07-23 钱维宏 Numerical weather prediction model system based on atmospheric variable Transient Eddy equation group

Also Published As

Publication number Publication date
CN113568067A (en) 2021-10-29

Similar Documents

Publication Publication Date Title
Liu et al. Comparison of support vector machine and copula-based nonlinear quantile regression for estimating the daily diffuse solar radiation: A case study in China
TWI745907B (en) A power generation amount prediction device, a power generation amount prediction method and program
Zhang et al. Advances in surrogate modeling for storm surge prediction: storm selection and addressing characteristics related to climate change
CN110232471B (en) Rainfall sensor network node layout optimization method and device
CN107563554B (en) Screening method of statistical downscaling model forecasting factors
CN113379107B (en) Regional ionosphere TEC forecasting method based on LSTM and GCN
CN110597873A (en) Precipitation data estimation method, precipitation data estimation device, precipitation data estimation equipment and storage medium
CN112381337B (en) Multi-source meteorological data fusion processing method, system, terminal and medium
CN103413036A (en) Continuous forest fire weather level forecasting model and application thereof
CN113568067B (en) Numerical weather forecasting method and device, computer storage medium and electronic equipment
Dong et al. Impact of climate change on flood frequency of the Trian Reservoir in Vietnam using RCMs
CN112100922A (en) Wind resource prediction method based on WRF and CNN convolutional neural network
CN112580844A (en) Meteorological data processing method, device, equipment and computer readable storage medium
CN113139327B (en) Ionized layer TEC single-point prediction method and system based on GRU network model
CN110968929A (en) Wind power plant wind speed prediction method and device and electronic equipment
Kang et al. Development of a deep neural network model to estimate solar radiation using temperature and precipitation
CN113435630B (en) Basin hydrological forecasting method and system with self-adaptive runoff yield mode
Han et al. FengWu-GHR: Learning the Kilometer-scale Medium-range Global Weather Forecasting
CN114529035A (en) CART-based wind speed forecasting method of multi-mode integrated model
CN108563733B (en) A kind of land use data assimilation method based on Bayesian inference
Zhang Estimation of Probabilistic Extreme Wind Load Effect with Consideration of Directionality and Uncertainty
CN115759483B (en) Photovoltaic electric field solar irradiance prediction method, electronic equipment and storage medium
CN116306038B (en) Vertical layer matching method and device for meteorological mode and assimilation system
CN113723006B (en) LS-SVM (least squares-support vector machine) -based single-station earth change magnetic field modeling prediction method and system
CN111400844A (en) Parameter scheme set generation method and wind speed forecasting method of meteorological model

Legal Events

Date Code Title Description
PB01 Publication
PB01 Publication
SE01 Entry into force of request for substantive examination
SE01 Entry into force of request for substantive examination
GR01 Patent grant
GR01 Patent grant