Regional convective scale ensemble forecasting method and system

Abstract

The invention discloses a regional convective scale ensemble forecasting method and system. According to the scheme, an initial field is optimized by adopting a background field deviation correction technology; enhancing the disturbance of the tropical region based on the observation disturbance of normal distribution; constructing multi-scale singular vector mixed initial value disturbance; dynamically adjusting a side boundary disturbance coefficient to improve dispersion; optimizing the cut-off wave number of the SPPT scheme to 50 to represent the small-scale uncertainty; and a self-adaptive post-processing technology is developed to guarantee the stability of the system. The implementation effect shows that the 24-hour heavy rainfall FSS score is improved by 9.3%-10.6%, the rainstorm probability forecasting skill is improved by 18%-29%, the calculation efficiency is improved, and the service operation fault rate is reduced to 0. According to the method, the extreme weather prediction capability and the prediction product reliability are remarkably improved.

Description

Technical Field

[0001] The present invention relates to the field of weather forecasting, and particularly to a regional convective scale ensemble forecasting method and system. Background Art

[0002] The accuracy and reliability of weather forecasting are crucial for fields such as disaster warning, agricultural production, and transportation. With the development of numerical weather forecasting technology, the ensemble forecasting method has gradually become an important means to improve forecasting accuracy. Ensemble forecasting generates multiple perturbed members to quantify the uncertainties of initial conditions, model parameters, and boundary conditions, thereby providing probabilistic forecasting products to help decision-makers evaluate the likelihood of weather events occurring.

[0003] Traditional regional ensemble forecasting systems usually rely on a model framework with a relatively low resolution (such as 10 km), making it difficult to accurately capture meso-scale and small-scale weather phenomena (such as convective storms, local heavy precipitation, etc.). In addition, the prior art has the following limitations:

[0004] Large initial condition errors: Traditional initial value perturbation methods (such as Ensemble Transform Kalman Filter ETKF) are difficult to fully characterize multi-scale uncertainties, especially with insufficient perturbations in tropical regions. During the background field assimilation process, there are systematic biases between the 6-hour forecast field and the analysis field of the global model, affecting the quality of the initial field.

[0005] Fixed lateral boundary perturbations with insufficient dispersion: The lateral boundary perturbations of regional ensemble forecasting usually adopt static coefficients and cannot adapt to the uncertainty changes under different weather situations, resulting in a low dispersion in the later stage of forecasting.

[0006] Insufficient expression of model physical process uncertainties: The truncation wave number of traditional Stochastic Physics Parameterization Tendencies (SPPT) is relatively low, making it difficult to reflect the small-scale error growth characteristics in convective scale models.

[0007] Computing resources and stability issues: High-resolution ensemble forecasting has high requirements for computing resources. Existing systems are prone to integration overflow or missing member data during the simulation of extreme weather events (such as typhoons, heavy rains), affecting the operational stability.

[0008] In recent years, convective scale ensemble forecasting (3 km resolution) has become a research hotspot, aiming to improve the forecasting skills of heavy precipitation and severe convective weather. However, how to optimize initial perturbations, improve the characterization of lateral boundary uncertainties, enhance model physical process perturbations, and ensure the stable operation of the system remains a technical problem to be solved urgently.

[0009] Therefore, it is necessary to develop a new regional convective scale ensemble forecasting method to overcome the above technical defects and improve the forecasting accuracy and the prediction ability of extreme weather events. Summary of the Invention

[0010] In view of the above deficiencies in the prior art, the present invention provides a regional convective scale ensemble forecasting method and system.

[0011] In order to achieve the above object of the invention, the technical solution adopted by the present invention is as follows:

[0012] A regional convective scale ensemble forecasting method, comprising the following steps:

[0013] S1. Obtain regional meteorological data, including the global forecast system background field, observation data, ensemble boundary perturbation field, and historical meteorological data;

[0014] S2. Based on the data obtained in S1, optimize the background field of the control forecast member through bias correction technology, and generate the initial conditions of multiple perturbation members through the observation perturbation method;

[0015] S3. According to the initial conditions generated in S2, construct an initial value perturbation field by fusing multi-scale singular vectors and observation perturbations, and generate the lateral boundary conditions of the ensemble members through a dynamically adjusted lateral boundary perturbation scheme;

[0016] S4. Randomly perturb the model parameters using an improved physical process perturbation method, and run a numerical weather prediction model for ensemble forecasting;

[0017] S5. Based on the adaptive post-processing technology, perform data conversion and product generation on the forecast results to ensure output stability.

[0018] Further, the specific steps of the bias correction technology in S2 are as follows:

[0019] S201. Statistically calculate the deviation field between the historical global forecast system background field and the corresponding analysis field. The calculation formula is:

[0020] δx b (t) = x f (t - 6) - x a (t)

[0021] where x f (t - 6) is the 6-hour forecast field of the global forecast system, x a (t) is the assimilation analysis field, t is the current time, and x b (t) is the deviation field;

[0022] S202. Perform time averaging on the historical deviation field to calculate the average deviation, and use the average deviation to correct the current background field, which is expressed as:

[0023]

[0024] where x b-new (t) is the new corrected background field, is the mean deviation.

[0025] Furthermore, the specific steps of the observation perturbation method in S2 are as follows:

[0026] S211. Generate random numbers based on the normal distribution and superimpose them on the unperturbed observation data to generate perturbed observations, expressed as:

[0027] O n = O ctl + R n × ε

[0028] where O n is the perturbed observation, O ctl is the unperturbed observation, R n is the standard normal random distribution number, and ε is the standard deviation of the observation error;

[0029] S212. Set dynamic threshold perturbations for the radar reflectivity data, including a basic threshold, an upward floating threshold, and a downward floating threshold, and allocate them to different ensemble members;

[0030] S213. Use the perturbed observations and the background field for three-dimensional variational assimilation to generate an initial perturbation field containing medium and small scale uncertainties.

[0031] Furthermore, the specific steps of the multi-scale singular vector initial value perturbation method in S3 are as follows:

[0032] S301. Calculate the multi-scale singular vector perturbation field based on different horizontal resolutions and optimized durations;

[0033] S302. Perform optimal truncation on the singular vectors in the spectral space, retaining the low-frequency perturbations of the large-scale singular vectors and the high-frequency perturbations of the small-scale singular vectors;

[0034] S303. Linearly superimpose the multi-scale singular vector perturbations and the observation perturbation analysis field to form a fused initial value perturbation field, expressed as:

[0035] Perb j = SVPerb j + ObPerb j

[0036] where SVperb j is the multi-scale SV optimal spectral truncation mixed initial value perturbation, Obsperb j is the observation perturbation, Perb j is the optimally fused initial value perturbation, and j is the ensemble member.

[0037] Furthermore, the specific steps of the dynamic lateral boundary perturbation scheme in S3 are as follows:

[0038] S311. Extract the difference between the background fields of the global ensemble forecast perturbation members and the control forecast to generate the lateral boundary perturbation field, expressed as:

[0039] bckg_perb i = bckg_GEPS i - bckg_GEPS0

[0040] bckg_CAEPS i = bckg_ctl + bckg_Perb i × r

[0041] In the formula, i is the perturbation member, bckg_GEPS i is the global background field of the perturbation member i, bckg_GEPS0 is the global background field of the control forecast, bckg_perb i is the global background forecast perturbation field, r is the lateral boundary perturbation coefficient, and bckg_ctl is the background field of the control forecast;

[0042] S312. Dynamically calculate the lateral boundary perturbation coefficient according to the ratio of the historical root mean square error to the dispersion, expressed as:

[0043]

[0044] In the formula, RMSE(k,t) is the root mean square error at the k-th layer and the t-th forecast lead time, spread(k,t) is the corresponding ensemble dispersion, nt is the forecast lead time, k is the vertical level, and r(k) is the lateral boundary perturbation coefficient;

[0045] S313. Multiply the perturbation coefficient by the lateral boundary perturbation field to generate the dynamically adjusted lateral boundary condition.

[0046] Furthermore, the specific steps of the physical process perturbation method in S4 are as follows:

[0047] S401. Use the stochastic physical process tendency perturbation method to perform lognormal distribution perturbation on the key parameters of the cloud microphysics and boundary layer parameterization schemes;

[0048] S402. Increase the spectral space truncation wave number to 50 to retain the uncertainty information of smaller scales;

[0049] S403. Generate a three-dimensional spatio-temporal correlated random field based on the first-order autoregressive process. The specific calculation formula is:

[0050]

[0051] In the formula, is a three-dimensional random function of Gaussian distribution, λ j , φ j , tj represent longitude, latitude, and time respectively, and α l,m (t) is the spectral coefficient of the random field at time t, l and m are the total horizontal wave number and the zonal wave number, and L is the horizontal truncation scale of the random field.

[0052] Furthermore, the specific steps of the adaptive post-processing technology in S5 are as follows:

[0053] S501. Detect the integrity and continuity of the output data of the ensemble members, and mark the missing or abnormal members;

[0054] S502. If the number of missing members exceeds the preset threshold, terminate the post-processing process for the current forecast period;

[0055] S503. Interpolate, calculate diagnostic quantities, and convert the format of the valid member data to generate standardized forecast products.

[0056] Furthermore, the diagnostic quantity calculation in S503 includes:

[0057] S5031. Interpolate the model surface data to the isobaric surface, use the cubic spline interpolation method to process continuous variables, and the linear interpolation method to process discontinuous variables;

[0058] S5032. Calculate dynamic factors, thermal factors, and water vapor factors, including the K index, convective available potential energy, and water vapor flux divergence;

[0059] S5033. Generate probabilistic forecast products, including ensemble mean, dispersion, and the probability of extreme weather events.

[0060] A regional convective-scale ensemble forecasting system includes:

[0061] A data processing module for integrating multiple data sources, including ensemble forecast perturbation data, observational data, ensemble boundary perturbation data, and global background field data, and performing preprocessing;

[0062] An ensemble member generation module configured to optimize the initial conditions through perturbation analysis and assimilation analysis to generate multiple perturbed initial fields; the perturbation analysis includes background field continuity and periodic deviation correction technology, specifically dynamically correcting the control forecast background field based on the historical statistical background field deviation;

[0063] A cloud analysis and fusion module for combining cloud physics data and multi-source information to generate a fused cloud physics field;

[0064] A model integration module that runs a numerical weather prediction model based on the initial conditions and lateral boundary conditions to simulate the future weather evolution;

[0065] The forecast output and post - processing module is used to convert the model integration results into data products and visualization graphs, and adopt adaptive post - processing technology to automatically adjust the product generation process according to effective ensemble members;

[0066] The evaluation and verification module is used to evaluate the forecast results through historical data and real - time observations, and provide feedback for system optimization.

[0067] The present invention has the following beneficial effects:

[0068] 1. Significantly improve the forecast accuracy. Through multi - scale initial value perturbation and dynamic lateral boundary perturbation technologies, the uncertainty expression of the initial field and boundary conditions is optimized, making the FSS score of 24 - hour heavy precipitation increase by 9.3%, and the forecast errors of the 850hPa temperature field and near - surface wind field are significantly reduced. The 3 - km high - resolution model can capture the characteristics of small - scale weather systems more accurately, and the forecast accuracy of the position and intensity of the rainstorm center is better than that of the 10 - km system.

[0069] 2. Enhance the prediction ability of extreme weather. The improved initial value perturbation and physical process perturbation schemes (such as optimizing the SPPT truncation wave number to 50) improve the prediction ability of extreme events such as severe convection and typhoons, and the 48 - hour heavy rain probability score increases by 29%. New products such as typhoon track intensity and neighborhood precipitation probability are added to provide more reliable decision - making support for disaster prevention.

[0070] 3. Optimize the calculation efficiency and resource utilization. By adopting parallel computing and dynamic perturbation coefficient adjustment technologies, a 72 - hour forecast is completed under 960 computing nodes (peak value), with a running time of 70 minutes. The resource utilization rate is increased by 100% compared with the 10 - km system, and the storage requirement only increases to 2T.

[0071] 4. Improve the uncertainty representation. The initial value scheme that combines multi - scale singular vectors and observational perturbations makes up for the lack of perturbations in the tropical region; the dynamic lateral boundary perturbation technology improves the matching degree of ensemble dispersion and root - mean - square error by 20%.

[0072] 5. Improve the system stability and service continuity. The adaptive post - processing technology ensures the stability of product output through an automatic fault - tolerance mechanism, with a member overflow rate of 0% during business operation, supporting zero - failure meteorological services for major events such as the Hangzhou Asian Games. Brief Description of the Drawings

[0073] Figure 1 It is a schematic flow chart of the regional convective - scale ensemble forecasting method of the present invention. Detailed Embodiments

[0074] The specific embodiments of the present invention will be described below to facilitate those skilled in the art of the present technology to understand the present invention. However, it should be clear that the present invention is not limited to the scope of the specific embodiments. For those of ordinary skill in the art of the present technology, as long as various changes are within the spirit and scope of the present invention defined and determined by the appended claims, these changes are obvious, and all inventions and creations using the concept of the present invention are within the scope of protection.

[0075] A regional convective scale ensemble forecasting method includes the following steps:

[0076] S1. Obtain regional meteorological data, including the global forecast system background field, observation data, ensemble boundary perturbation field, and historical meteorological data;

[0077] In this embodiment, data from multiple sources are integrated, including ensemble forecast perturbations (GEPS-SVs), observation data, ensemble boundary perturbations (GEPS-Bdy), and the global background field (such as GFS data). These initial fields are perturbed from the observation data to reflect the uncertainty of the atmospheric state. These data are preprocessed to provide basic information for subsequent steps.

[0078] S2. Based on the data obtained in S1, optimize the background field of the control forecast member through bias correction technology, and generate the initial conditions of multiple perturbed members through the observation perturbation method;

[0079] In this embodiment, the background field continuity and periodicity bias correction technology

[0080] The initial value of the control forecast for the 3 km convective scale ensemble forecast is generated by performing three-dimensional variational assimilation on the background field of the global model's 6-hour forecast and the observation data. To address the bias problem of the background field, bias correction technology can provide a more accurate background field for the control forecast of the convective scale ensemble forecast. The background field bias correction algorithm is to statistically analyze the bias of the historical background field in real time, conduct sensitivity tests on the statistical time length, obtain the relatively optimal number of historical days, and then subtract the statistically obtained average background field bias from the original background field to complete the bias correction of the control forecast background field.

[0081] To correct the background field bias, a bias correction algorithm for the control forecast of the convective scale ensemble forecast is constructed using historical background field information.

[0082] This algorithm statistically analyzes the bias of the historical background field in real time, selects the statistically obtained average background field bias based on the optimal time, and corrects the bias of the background field. The specific steps are as follows:

[0083] S201. Statistically analyze the bias field between the historical global forecast system background field and the corresponding analysis field.

[0084] The deviation fields on a daily basis within a period of time before the starting forecast time are statistically obtained from the 6-hour forecast fields of NCEP-GFS on a daily basis and the assimilation analysis fields of NCEP-GFS at the corresponding times. Let x f (t - 6) represent the 6-hour forecast of NCEP-GFS, and x a (t) represent the analysis field of NCEP-GFS at the corresponding time, and δx b (t) represent the background field deviation obtained through statistics. t is the starting forecast time, and the calculation formula for the deviation field is:

[0085] δx b (t) = x f (t - 6) - x a (t)

[0086] S202. Perform time averaging on the historical deviation fields to calculate the average deviation, and use the average deviation to correct the current background field.

[0087] Statistically average the deviation fields on a daily basis within the obtained period. In the formula, represents the background field deviation of the statistical average. n represents the sample size within the statistical period. Since it is a daily statistics, it can also be regarded as the number of days of historical data. For the value of n, corresponding sensitivity tests were carried out, and based on the results, the relatively optimal time was determined, n = 40, that is, the historical data 40 days before the starting forecast time are used as the samples for statistics.

[0088]

[0089] Use the obtained average background field deviation for deviation correction, which is expressed as:

[0090]

[0091] In the formula, x b_new (t) represents the new corrected background field, and x b (t) is the original background field.

[0092] Through the above steps, the deviation correction of the background field for the control forecast member is completed.

[0093] S3. According to the initial conditions generated in S2, construct an initial perturbation field that fuses multi-scale singular vectors and observational perturbations, and generate the lateral boundary conditions of the ensemble members through a dynamically adjusted lateral boundary perturbation scheme;

[0094] In this embodiment, the specific steps of the multi-scale singular vector initial perturbation method are as follows:

[0095] S301. Calculate the multi-scale singular vector perturbation field based on different horizontal resolutions and optimized durations;

[0096] In the process of solving for SVs, the integrated variables of the GRAPES global tangent linear model and adjoint model include the perturbation values of the zonal wind u, meridional wind v, perturbation potential temperature θ′, and perturbation dimensionless pressure ∏′, denoted as x = (u′, v′, θ′″′ T ). The initial value perturbation vector x is integrated backward for a certain period of time through the GRAPES global TLM (denoted as: tangent linear operator L) (the integration time of TLM and ADM is the optimal time interval), and the evolved perturbation vector x t = Lx. The singular vector is the initial value perturbation that grows fastest within the optimal time interval in the tangent linear model. Its calculation can be transformed into maximizing the ratio J(x) of the norm of the evolved perturbation vector to the norm of the initial value perturbation vector:

[0097]

[0098] where ‖·‖ is defined as the Euler inner product (·,·); E is the transformation operator that transforms the physical space perturbation vector x into a dimensionless vector in the Euler space for mathematical calculations.

[0099] The above equation can be transformed and expressed as:

[0100]

[0101] The above equation can be transformed into the eigenvalue decomposition problem of the matrix (ELE -1 ) T (ELE -1 ) = E -1 L T E 2 LE -1 :

[0102]

[0103] where λ is 2 and are the is-th eigenvalue and eigenvector of the matrix E -1 L T E 2 LE -1 respectively, and L T is the GRAPES global adjoint model. Correspondingly, λ is and are the is-th singular value and right singular vector in the Euler space of the matrix ELE -1 respectively.

[0104] Furthermore, a projection operator P is introduced to constrain the calculation region of the singular vector, that is, the target region. Equation (2.7) becomes:

[0105]

[0106] The above-mentioned transformation operator E is the key factor for measuring the magnitude of the perturbation. The formula for calculating the modulus of the total dry energy (TE, total energy) of the GRAPES singular vector is as follows:

[0107]

[0108] In the formula, the sum of the first two terms is the modulus of the perturbation kinetic energy (KE, kinetic energy), and the sum of the last two terms is the modulus of the perturbation potential energy (PE, potentital energy). λ and are the longitude and latitude of the spherical coordinate system of the GRAPES model respectively, is the height terrain-following coordinate, C p is the specific heat at constant pressure of dry air, ρ r , T r , θ r and ∏ r are the reference density, reference temperature, reference potential temperature, and reference dimensionless pressure respectively. The transformation operator E based on the modulus of the total dry energy is:

[0109]

[0110] Table 1 shows the parameter settings of SVs at different scales calculated based on the CMA global tangent linear and adjoint models. In this embodiment, different-scale SV initial value perturbations are constructed using 2.5° SV (optimal time interval is 48 h), 1.5° SV (optimal time interval is 24 h), and 0.5° SV (optimal time interval is 6 h). The model integration time steps are 1200 s, 900 s, and 300 s respectively. The target regions of the 2.5° SV initial value perturbations are the Northern and Southern Hemispheres, and the typhoon SV is calculated. The specific formula is as shown in 2.11; the target regions of the 1.5° SV initial value perturbations are the Northern and Southern Hemispheres, and the formula is as shown in 2.12; the target region of the 0.5° SV initial value perturbations is the eastern region of China (102°E - 135°E, 17°N - 50°N), and the formula is as shown in 2.13. The linearized physical processes of 2.5° SV and 1.5° SV are dry physical processes, including the linearized planetary boundary layer (PBL) scheme and the linearized mountain blocking (MB) drag scheme. The linearized physical process of 0.5° SV is a wet physical process, including the linearized PBL and MB schemes, the linearized large-scale condensation and cumulus convection parameterization methods to capture the initial value perturbation structures at smaller scales related to thermal uncertainties. Moreover, the energy modes of different-scale SVs all adopt dry energy modes. The perturbation vertical levels are from layer 4 to layer 60. The number of 2.5° and 1.5° SVs is 30, and the number of 0.5° SVs is 20. It should be noted that although the impact of the linearized large-scale condensation scheme on SVs is mainly reflected in the temperature and specific humidity terms, since we adopt the dry energy mode, in the case study of this embodiment, the impact of the linearized large-scale condensation process is mainly reflected in the temperature term.

[0111] Table 1 Parameter Settings of Different-Scale SVs

[0112]

[0113]

[0114]

[0115]

[0116] Among them, j is the number of initial value perturbations. Here, j is taken as 7. The 7 initial value perturbations are added to the control forecast initial value in the form of positive and negative pairs to obtain 14 perturbed member initial values. They are 2.5° SV with the target region in the Northern Hemisphere, 2.5° SV with the target region in the Southern Hemisphere, and 2.5° SV of tropical cyclones respectively. They are respectively the numbers of and They are 1.5° SVs with the target area being the Northern Hemisphere and 1.5° SVs with the target area being the Southern Hemisphere; and are respectively and the numbers of; is 0.5° SV with the target area being East Asia, is the number of; respectively represent random numbers of 2.5° SV, 1.5° SV, and 0.5° SV; are respectively the initial value perturbations of 2.5° SV, 1.5° SV, and 0.5° SV.

[0117] S302. Optimally truncate the singular vectors in the spectral space, retaining the low-frequency perturbations of large-scale singular vectors and the high-frequency perturbations of small-scale singular vectors;

[0118] A multi-scale singular vector initial value perturbation method for convective-scale ensemble forecasting is constructed for the multi-scale characteristics of the initial value errors of the convective-scale model. That is, based on the tangent linear adjoint model of the CMA global ensemble forecasting, different-scale SVs initial value perturbations are constructed using a 2.5° horizontal resolution (optimization duration 48h), a 1.5° horizontal resolution (optimization duration 24h), and a 0.5° horizontal resolution (optimization duration 6h). The initial value perturbations are integrated by the model to obtain the ensemble member forecast fields at 0h, and they are dynamically downscaled into the convective-scale ensemble forecasting model to obtain the initial value fields formed by different-scale SVs initial value perturbations of the convective-scale ensemble forecasting.

[0119] S303. Linearly superpose the multi-scale singular vector perturbations and the observed perturbation analysis field to form a fused initial value perturbation field.

[0120] Based on the ensemble member analysis field obtained from the observed perturbations, calculate the differences between the perturbed members and the analysis field of the control forecast as the observed initial value perturbations. Linearly combine the observed perturbations and the multi-scale SVs initial value perturbations to construct an optimal fused initial value perturbation method to make up for the problem of insufficient perturbation distribution of the multi-scale SVs initial value perturbations in the southern region, with the expectation of obtaining the initial value perturbation structure in the entire Chinese region.

[0121] In this embodiment, the multi-scale SV optimal spectral truncation hybrid initial value perturbations and the observed perturbations are combined to construct an optimal fused initial value perturbation method. The specific formula is as follows:

[0122] Perb j = SVperb j + Obsperb j

[0123] where j is the ensemble member, SVperb jFor the multi-scale SV optimal spectral truncation hybrid initial value perturbation, Obsperb j For the observation perturbation, Perb j For the initial value perturbation of optimal fusion. It should be noted here that the observations of each member formed by the observation perturbation are assimilated with the NCEP 6-hour background field to obtain the assimilation analysis fields of each member, and a mixing process is carried out between the assimilation analysis field and the NCEP analysis field, so that the finally formed assimilation analysis field combines the large-scale components of the NCEP analysis field and the small-scale components of the assimilation analysis field formed by the observation perturbation, which is consistent with the assimilation process of CMA-MESO.

[0124] The specific steps of the dynamic lateral boundary perturbation scheme are as follows:

[0125] S311. Extract the difference between the global ensemble forecast perturbation members and the background field of the control forecast to generate the lateral boundary perturbation field.

[0126] In the convective scale ensemble forecast system, the hybrid lateral boundary perturbation technology is used, that is, based on the background field driving data of the CMA-GEPS global ensemble forecast system, the background field perturbation values required by the convective scale ensemble forecast system are extracted, the perturbation values are multiplied by the lateral boundary perturbation coefficient, and added to the control forecast background field to form the perturbation member background field, and the lateral boundary fields of the ensemble members are formed through dynamic downscaling.

[0127] Aiming at the problem of insufficient dispersion in the convective scale ensemble forecast with a fixed lateral boundary perturbation coefficient, the hybrid lateral boundary perturbation scheme in the convective scale ensemble forecast is optimized and improved. Based on the relationship between the historical root mean square error and dispersion, a scheme for dynamically adjusting the lateral boundary perturbation coefficient is constructed to obtain a better relationship between the dispersion skills of the convective scale ensemble forecast.

[0128] The specific method of the hybrid lateral boundary perturbation is that the source of the control forecast lateral boundary is the same as that of CMA-MESO, which is NCEP-GFS. To characterize the forecast uncertainty of CMA-GEPS, we adopt the hybrid lateral boundary perturbation method, that is, extract the global background forecast perturbation field of the CMA-GEPS perturbation members relative to the control forecast, multiply it by the lateral boundary perturbation amplitude, and add it to the control forecast global background driving data to construct the global hybrid background field, and obtain the lateral boundary ensemble of each member through the way of dynamic downscaling. The specific formula is shown as follows:

[0129] bckg_Perb i =bckg_GEPS i -bckg_GEPS0

[0130] bckg_CAEPS i =bckg_ctl+bckg_Perb i ×r

[0131] where \(i\) is the perturbation member, and \(bckg\_GEPS\) i is the global background field of the perturbation member \(i\), \(bckg\_GEPS0\) is the global background field of the control forecast, and \(bckg\_perb\) i is the perturbation field of the global background forecast. \(r\) is the lateral boundary perturbation coefficient, and \(bckg\_ctl\) is the background field of the control forecast. The lateral boundary perturbation variables are the \(u\) and \(v\) wind field variables, and the dimensionless air pressure \(p_i\) and potential temperature \(\theta\) fields are constrained according to \(u\) and \(v\). This is mainly based on the geostrophic adjustment theory. The convective scale model focuses on small-scale processes. When the horizontal scale of the ageostrophic perturbation is small, the air pressure field adapts to the wind field. That is, during the adaptation process, the initial ageostrophic state causes horizontal divergence and convergence as well as vertical motion, exciting fast inertia-gravity waves. The fast dispersion of this wave causes the energy of the local ageostrophic perturbation to quickly disperse outward, while establishing the corresponding air pressure field and reaching the quasi-geostrophic balance state.

[0132] S312. Dynamically calculate the lateral boundary perturbation coefficient according to the ratio of the historical root mean square error to the dispersion.

[0133] The main purpose of the lateral boundary perturbation coefficient is to ensure a good fit between the ensemble dispersion and the root mean square error of the convective scale ensemble forecast. For this reason, in addition to the static perturbation amplitude adjustment scheme, that is, the \(u\) and \(v\) perturbation coefficients at different vertical levels are all taken as the same value, we have also established a set of dynamic perturbation amplitude adjustment schemes, that is, the average perturbation coefficient \(r\) is calculated in real time according to the ratio relationship between the root mean square error and the dispersion of the first 10 historical time steps, and is used for the construction of the lateral boundary perturbation for the next time step start. And the coefficients at different vertical levels are different. The specific formula for the lateral boundary perturbation coefficient is as follows:

[0134]

[0135] where \(RMSE(k,t)\) is the root mean square error at the \(k\)-th layer and the \(t\)-th forecast time, \(spread(k,t)\) is the corresponding ensemble dispersion, \(n_t\) is the forecast time, \(k\) is the vertical level, and \(r(k)\) is the lateral boundary perturbation coefficient.

[0136] Calculate the ratio of the root mean square error to the dispersion at different vertical levels and forecast times (excluding the initial time because the initial value perturbation mainly plays a role at the initial time), and average the ratio over the forecast time to obtain the \(u\) and \(v\) perturbation coefficients at different vertical levels.

[0137] S4. Use the improved physical process perturbation method to randomly perturb the model parameters and run the numerical weather prediction model for ensemble forecasting;

[0138] In this embodiment, the specific steps of the physical process perturbation method are as follows:

[0139] S401. Use the stochastically perturbed physical process tendency method to perturb the key parameters of the cloud microphysics and boundary layer parameterization schemes with a lognormal distribution.

[0140] Stochastically Perturbed Parameterization Tendencies (SPPT) is to stochastically perturb the physical tendencies in view of the error effects of sub-grid physical process parameterization schemes (such as cloud physics, convection, turbulent diffusion, etc.). First, a stochastic model needs to be constructed. CMA-REPS V4.0 adopts the three-dimensional spatio-temporal correlation stochastic model generation technology independently developed by us, which has a normal distribution characteristic based on the first-order autoregressive process. The introduction is as follows:

[0141] Define ψ j (λ j , φ j , t j ) as a three-dimensional random function of Gaussian distribution (mean μ = 0, standard deviation σ):

[0142]

[0143] S402. By increasing the spectral space truncation wave number to 50, retain the uncertainty information of smaller scales.

[0144] , by increasing the value of the truncation wave number, retain the random information of higher frequencies in the spectral space, so as to form a stochastic perturbation distribution containing the uncertainty of smaller scales. Based on the design architecture of the SPPT scheme, increasing the truncation wave number will bring more computational consumption. Therefore, we comprehensively considered the system resource consumption and computational duration, and determined through experiments that the new maximum truncation wave number is 50 (the original value is 28). From the new spatial distribution of the stochastic model, the optimized SPPT scheme can represent a smaller spatial scale of uncertainty and better meet the requirements of convective-scale probabilistic forecasting. Taking the heavy rain caused by the northeast cold vortex as an example, an index of the relative importance of initial value perturbation and model perturbation was designed to characterize the influence of initial value perturbation and model perturbation on the ensemble forecast dispersion. The experimental results show that neither adding initial value perturbation nor model perturbation has a reducing effect on the ensemble dispersion, and this reducing effect will weaken with the increase of the forecast lead time. Initial value perturbation plays an important role in both isobaric surface elements and near-surface elements, and model perturbation has an effect on the dispersion of 850 hPa and near-surface elements.

[0145] S403. Generate a three-dimensional spatio-temporal correlated random field based on the first-order autoregressive process. The specific calculation formula is:

[0146]

[0147] In the formula, is a three-dimensional random function of Gaussian distribution, λj , φ j , t j represent longitude, latitude, and time respectively, and α l,m α(t) is the spectral coefficient of the random field at time t, l and m are the total horizontal wavenumber and zonal wavenumber respectively, L is the horizontal truncation scale of the random field, and α l,m (t) is shown as follows:

[0148]

[0149] Based on the random field introduce a stretching function defined as:

[0150]

[0151] β is taken as a constant with a value of -1.27, and ψ max represents the upper boundary of the random field ψ j (λ j , φ j , t j ). Through the construction of the stretching function, the upper and lower boundaries can be given.

[0152] The random perturbation values constructed by the above method vary with space and time, making each key parameter have an independent random perturbation field, and ensuring that each parameter varies within its reasonable statistical experience range. The above formula is a common first-order Markov process (also known as red noise). At the initial time t = 0, R(t) is a random number with a normal distribution N(0, 1), and the time correlation scale of the generated random perturbation field is determined by the parameter τ. The random perturbation field is expanded in the horizontal space through spherical harmonic functions, so that the obtained random perturbation field not only varies with time but also has spatial structure characteristics. In addition, considering the change in the normal distribution characteristic structure of the random perturbation field, the PDF distribution of the random perturbation field is changed through the stretching function.

[0153] S5. Perform data conversion and product generation on the forecast results based on the adaptive post-processing technology to ensure output stability.

[0154] During the operation of the numerical weather prediction (NWP) operational system, situations where the system fails to run successfully occur occasionally due to various reasons. For the ensemble prediction system, multiple ensemble members, initial condition perturbations, and model perturbations will further increase the probability of such situations, which will cause delays in data products and graphical products of the post-processing module of the entire system, thereby affecting the normal provision of probabilistic service products. To solve this problem, a post-processing technology with adaptive ensemble members has been developed and applied to the CMA-REPS V4.0 convective-scale ensemble prediction operational system to ensure the stable and reliable output of ensemble prediction products and improve the stability of system operation. This technology discriminates the ensemble members that have successfully completed the calculation, automatically adjusts the calculation and processing of relevant ensemble probability products in the forecast data of the system, and automatically adjusts the drawing of forecast graphical products.

[0155] The CMA model unified post-processing system is a unified post-processing and product development system that can process the CMA global, regional, typhoon, and ensemble prediction systems. Except for the different file reading, other processes such as interpolation processing, diagnostic quantity calculation, and data output share the same processing method. The processing process of this post-processing system mainly includes the reading of the model's raw output, the reading of the system control file, the encoding of the model surface data, and the calculation and encoding of the pressure surface data and diagnostic quantity data of forecast elements, etc. Specifically, it can be divided into the following steps:

[0156] (1) Read the model's raw output and the control file varname.dat, which specifies which elements and vertical levels are to be processed;

[0157] (2) Calculate temperature and pressure from the model surface pressure function and potential temperature. The CMA model has an equidistant longitude-latitude grid in the horizontal direction and a terrain-following coordinate in the vertical direction, using the Charney-Phillips leapfrog grid. Temperature, humidity, cloud amount, and vertical velocity are placed on the whole model levels, and pressure and horizontal wind speed are calculated on the half levels. Therefore, in the process of post-processing the model from the model surface to the isobaric surface, the model levels need to be processed first. The basic forecast variables representing temperature and pressure in the model system are potential temperature and the EXNER pressure function. In post-processing, the model forecast variables potential temperature and Exner pressure need to be converted into temperature and pressure through calculation formulas first.

[0158] (3) Vertical interpolation of three-dimensional elements from the model surface to the isobaric surface. In the post-processing process of numerical weather prediction models, vertical interpolation is a very important part. Usually, different interpolation schemes are adopted according to the characteristics of model elements. For example, linear interpolation is often used for the vertical interpolation of discontinuous cloud water substances. In the CMA system, linear interpolation on ln(P) is used for humidity and cloud amount at each level, and cubic spline interpolation (abbreviated as Spline interpolation) is used for the vertical interpolation of relatively continuous variables such as height, temperature, and wind.

[0159] (4) Calculate sea-level pressure and other diagnostic quantities. Diagnostic quantities are calculated based on the basic model forecast quantities. Quantitative calculations and analyses are performed on the terms of the thermodynamic and kinetic equations that describe the physical and dynamic processes governing various weather phenomena and weather systems to understand their relative roles. By calculating diagnostic quantities, we can obtain more important weather information from a limited number of conventional meteorological observation data, which helps us make in-depth and quantitative explanations of the dynamic and thermodynamic characteristics of various weather systems and weather processes, and on this basis, establish conceptual models of various weather systems and weather processes to provide reliable information for weather forecasting.

[0160] The diagnostic quantities output by the CMA post-processing system can be classified into dynamic factors, thermodynamic factors, water vapor factors, etc. according to their physical properties, mainly including K-index, lifting condensation level, convective available potential energy, convective inhibition energy, lifting index, temperature advection, 0-degree layer height, total index, water vapor flux, water vapor flux divergence, temperature dew point difference, etc.

[0161] This embodiment also provides a regional convective-scale ensemble forecasting system, including:

[0162] A data processing module, configured to integrate multiple data sources, including ensemble forecast perturbation data, observational data, ensemble boundary perturbation data, and global background field data, and perform preprocessing;

[0163] An ensemble member generation module, configured to optimize the initial conditions through perturbation analysis and assimilation analysis to generate multiple perturbed initial fields; the perturbation analysis includes background field continuity and periodic deviation correction techniques, specifically dynamically correcting the control forecast background field based on the background field deviation statistically obtained from history.

[0164] A cloud analysis and fusion module, configured to combine cloud physics data and multi-source information to generate a fused cloud physics field;

[0165] A model integration module, which runs a numerical weather prediction model based on the initial conditions and lateral boundary conditions to simulate the future weather evolution;

[0166] A forecast output and post-processing module, configured to convert the model integration results into data products and visualization graphics, and adopt an adaptive post-processing technology to automatically adjust the product generation process according to effective ensemble members;

[0167] An evaluation and verification module, configured to evaluate the forecast results through historical data and real-time observations and provide feedback for system optimization.

[0168] For the member adaptive function of the post-processing of regional convective scale ensemble forecasts and the operation process of the product subsystem, it is designed based on the ECFLOW workflow management and monitoring software to achieve the effect of automatic scheduling and adjustment of the forecast member modules. The Python language and Shell script language are used to develop and run the scheduling of each module in the post-processing of the regional convective scale ensemble forecast and the product subsystem in China, generate a total scheduling file that realizes the node relationship, trigger relationship, and job type of each module, and configure a parameter file for the total environment variables of the subsystem that is adapted to the environment variables of the total scheduling file of the subsystem, and finally realize the member adaptation of the downstream modules of each member integration module.

[0169] Specifically, the function is designed as follows:

[0170] (1) Define a function get_names_of_grib2 in the Grib2Processor class of the grib2_process.py script. The function is to obtain all grib2 file names under a specified path, sort them in ascending order according to the member number and forecast lead time, store them in a list named grb2_files_list, and return the list.

[0171] (2) Define a function get_valid_grib2_files in the Grib2Processor class. The function is to obtain grib2 files with continuous forecast lead times starting from 0. Specifically: First, define an empty dictionary grb2_file_dict, traverse the elements in the grb2_files_list list returned by the get_names_of_grib2 function in (1). For a certain ensemble member, if the file for a certain forecast lead time is missing, then add the ensemble member number to the variable invalid_x_member (i.e., the list of missing member numbers), and directly enter the next member loop; if all 73 grib2 files for this member exist, then store them in the grb2_file_dict dictionary in the form of key (member number, forecast lead time) and value grib2 file name, and enter the next member loop. In this way, the finally obtained grb2_file_dict dictionary contains all grib2 files with continuous forecast lead times starting from 0, that is, the grib2 file information of valid members.

[0172] (3) After obtaining the grib2 data path in the main program main.py, call the functions in (1)-(2) above to obtain the valid member dictionary information grib2_filedict.

[0173] (4) Implement the judgment of the number of valid members in the function (read_grib2_by_*) for reading grib2 data, so as to determine whether the program continues to run or exits. Specifically, pass the grib2_filedict variable into the read_grib2_by_* function. In this function, it will loop through 15 ensemble members. At the beginning of the loop, first judge whether the current member and the current forecast lead time are in the valid member dictionary grib2_filedict. If so, continue to read the grib2 data; otherwise, add the current ensemble member number to the invalid member list invalid_member and directly enter the loop of the next ensemble member. After completing all the ensemble member loops in this way, the final invalid_member is obtained. At this time, if the number of elements in invalid_member exceeds the preset value (set to 5 in this system), it is considered that there are too many missing ensemble members under the current forecast lead time, and the python program is directly exited without performing subsequent calculations and plotting operations; otherwise, the subsequent calculations and plotting are completed normally.

[0174] The present invention is described with reference to the flowcharts and / or block diagrams of methods, apparatuses (systems), and computer program products according to embodiments of the present invention. It should be understood that each process and / or block in the flowchart and / or block diagram, and the combination of processes and / or blocks in the flowchart and / or block diagram, can be implemented by computer program instructions. These computer program instructions can be provided to the processor of a general-purpose computer, a special-purpose computer, an embedded processor, or other programmable data processing devices to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing devices generate a machine for implementing the functions specified in Figure 1 one process or multiple processes and / or blocks Figure 1 the device for the function specified in one block or multiple blocks.

[0175] These computer program instructions can also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to work in a specific manner, so that the instructions stored in the computer-readable memory generate a manufactured article including an instruction device, and the instruction device implements the functions specified in Figure 1 one process or multiple processes and / or blocks Figure 1 the function specified in one block or multiple blocks.

[0176] These computer program instructions can also be loaded onto a computer or other programmable data processing device, so that a series of operation steps are executed on the computer or other programmable device to generate a computer-implemented process. Thus, the instructions executed on the computer or other programmable device provide for implementing the functions specified in Figure 1 one process or multiple processes and / or blocks Figure 1Steps of the functions specified in one or more boxes.

[0177] In the present invention, specific embodiments are used to illustrate the principles and implementation manners of the present invention. The descriptions of the above embodiments are only used to help understand the method of the present invention and its core idea. At the same time, for those of ordinary skill in the art, according to the idea of the present invention, there will be changes in the specific implementation manners and application scopes. In summary, the content of this specification should not be construed as a limitation to the present invention.

[0178] Those of ordinary skill in the art will realize that the embodiments described herein are for helping the reader understand the principles of the present invention, and it should be understood that the protection scope of the present invention is not limited to such specific statements and embodiments. Those of ordinary skill in the art can make various other specific deformations and combinations that do not depart from the essence of the present invention according to the technical revelations disclosed in the present invention, and these deformations and combinations are still within the protection scope of the present invention.

Claims

1. A regional convective scale ensemble forecasting method, characterized in that It includes the following steps: S1. Obtain regional meteorological data, including the Global Forecast System background field, observational data, ensemble boundary perturbation field, and historical meteorological data; S2. Based on the data obtained in S1, optimize the background field of the control forecast member through bias correction technology, and generate the initial conditions of multiple perturbation members through the observational perturbation method; S3. According to the initial conditions generated in S2, construct an initial value perturbation field by fusing multi-scale singular vectors and observational perturbations, and generate the lateral boundary conditions of the ensemble members through a dynamically adjusted lateral boundary perturbation scheme; S4. Use the improved physical process perturbation method to randomly perturb the model parameters, and run the numerical weather prediction model for ensemble forecasting; S5. Based on the adaptive post-processing technology, perform data conversion and product generation on the forecast results to ensure output stability.

2. The regional convective scale ensemble forecasting method according to claim 1, wherein The specific steps of the bias correction technology in S2 are as follows: S201. Statistically calculate the deviation field between the historical Global Forecast System background field and the corresponding analysis field. The calculation formula is: δx b (t) = x f (t - 6) - x a (t) where x f (t - 6) is the 6-hour forecast field of the global forecast system, x a (t) is the assimilation analysis field, t is the current time, x b (t) is the deviation field; S202. Perform time averaging on the historical deviation field to calculate the average deviation, and use the average deviation to correct the current background field, expressed as: where x b-new (t) is the corrected new background field, is the average deviation.

3. The regional convective scale ensemble forecasting method according to claim 1, wherein The specific steps of the observational perturbation method in S2 are as follows: S211. Generate random numbers based on the normal distribution, and superimpose them on the unperturbed observational data to generate perturbed observations, expressed as: O n = O ctl + R n × ε Among them, O n is the disturbance observation, O ctl is the undisturbed observation, R n is the standard normal random distribution number, and ε is the standard deviation of the observation error; S212. Set dynamic threshold perturbations for radar reflectivity data, including a base threshold, an upper threshold, and a lower threshold, and distribute them to different ensemble members; S213. Use the perturbed observations and the background field for three-dimensional variational assimilation to generate an initial perturbation field containing medium and small-scale uncertainties.

4. The regional convective scale ensemble forecasting method according to claim 1, characterized in that The specific steps of the multi-scale singular vector initial value perturbation method in S3 are as follows: S301. Calculate the multi-scale singular vector perturbation field based on different horizontal resolutions and optimization durations; S302. Perform optimal truncation on the singular vectors in the spectral space, and retain the low-frequency perturbations of the large-scale singular vectors and the high-frequency perturbations of the small-scale singular vectors; S303. Linearly superimpose the multi-scale singular vector perturbations and the observational perturbation analysis field to form a fused initial value perturbation field, expressed as: Perb j = SVPerb j + ObPerb j where SVperb j is the multi-scale SV optimal spectral truncation hybrid initial value perturbation, Obsperb j is the observation perturbation, Perb j is the optimal fusion initial value perturbation, and j is the ensemble member.

5. The regional convective scale ensemble forecasting method according to claim 1, characterized in that The specific steps of the dynamic lateral boundary perturbation scheme in S3 are as follows: S311. Extract the difference between the global ensemble forecast perturbation member and the background field of the control forecast to generate a lateral boundary perturbation field, expressed as: bckg_Perb i = bckg_GEPS i - bckg_GEPS0 bckg_CAEPS i = bckg_ctl + bckg_Perb i × r where \(i\) is the perturbed member, and \(bckg\_GEPS\) i is the global background field of the perturbed member \(i\), \(bckg\_GEPS0\) is the global background field of the control forecast, and \(bckg\_perb\) i is the global background forecast perturbation field, \(r\) is the lateral boundary perturbation coefficient, and \(bckg\_ctl\) is the background field of the control forecast; S312. Dynamically calculate the lateral boundary perturbation coefficient according to the ratio of the historical root mean square error to the dispersion, expressed as: In the formula, RMSE(k,t) is the root mean square error at the k-th layer and the t-th forecast time, spread(k,t) is the corresponding ensemble dispersion, nt is the forecast time, k is the vertical level, and r(k) is the lateral boundary perturbation coefficient; S313. Multiply the perturbation coefficient by the lateral boundary perturbation field to generate dynamically adjusted lateral boundary conditions.

6. The regional convective scale ensemble forecasting method according to claim 1, wherein The specific steps of the physical process perturbation method in S4 are as follows: S401. Adopt the random physical process tendency perturbation method to perform lognormal distribution perturbations on the key parameters of the cloud microphysics and boundary layer parameterization schemes; S402. By increasing the spectral space truncation wave number to 50, retain smaller-scale uncertainty information; S403. Generate a three-dimensional spatio-temporal correlated random field based on a first-order autoregressive process. The specific calculation formula is as follows: In the formula, is a three-dimensional random function of Gaussian distribution, λ j , φ j , t j represent longitude, latitude and time respectively, α l,m (t) is the spectral coefficient of the random field at time t, l and m are the total horizontal wave number and zonal wave number, and L is the horizontal truncation scale of the random field.

7. The regional convective scale ensemble forecasting method according to claim 1, characterized in that The specific steps of the adaptive post-processing technology in S5 are as follows: S501. Detect the integrity and continuity of the output data of the ensemble members, and mark the missing or abnormal members. S502. If the number of missing members exceeds a preset threshold, terminate the post-processing process for the current forecast period. S503. Interpolate, calculate diagnostic quantities, and perform format conversion on the valid member data to generate standardized forecast products.

8. The regional convective scale ensemble prediction method according to claim 7, wherein The diagnostic quantity calculation in S503 includes: S5031. Interpolate the model surface data to the isobaric surface, using cubic spline interpolation for continuous variables and linear interpolation for discontinuous variables. S5032. Calculate dynamic factors, thermodynamic factors, and water vapor factors, including the K index, convective available potential energy, and water vapor flux divergence. S5033. Generate probability forecast products, including ensemble mean, dispersion, and the probability of extreme weather events.

9. A system for the regional convective scale ensemble forecasting method according to any one of claims 1-8, characterized in that, Including: A data processing module for integrating multiple data sources, including ensemble forecast perturbation data, observational data, ensemble boundary perturbation data, and global background field data, and performing preprocessing. An ensemble member generation module configured to optimize the initial conditions through perturbation analysis and assimilation analysis to generate multiple perturbed initial fields. The perturbation analysis includes a background field continuity and periodic deviation correction technology, specifically dynamically correcting the control forecast background field based on the historically statistical background field deviation. A cloud analysis and fusion module for combining cloud physics data and multi-source information to generate a fused cloud physics field. A model integration module that runs a numerical weather prediction model based on the initial conditions and lateral boundary conditions to simulate the future weather evolution. A forecast output and post-processing module for converting the model integration results into data products and visualization graphics, and using an adaptive post-processing technology to automatically adjust the product generation process according to the valid ensemble members. An evaluation and verification module for evaluating the forecast results through historical data and real-time observations and providing feedback for system optimization.

Images