A tornado probability prediction method based on ensemble kalman filter and large eddy simulation

Abstract

The application particularly relates to a tornado probability prediction method based on ensemble Kalman filtering and large eddy simulation. The method comprises the following steps: acquiring an ensemble initial field and determining a first target scale region; performing a first data extrapolation operation on the ensemble initial field to obtain a quasi-assimilation ensemble initial field corresponding to the first target scale region; performing an assimilation operation on the quasi-assimilation ensemble initial field to obtain an ensemble final analysis field; performing spatial gridding on the first target scale region to obtain a plurality of grid cell columns; determining, for each member, the helicity information corresponding to the final analysis field of the member, discarding a plurality of members, and determining the effective grid cell columns corresponding to each of the retained members; determining a target grid cell column and a second target scale region based on the effective grid cell columns; and acquiring tornado occurrence probability information corresponding to the second target scale region based on the horizontal wind component data contained in the final analysis field of each of the retained members.

Description

Technical Field

[0001] This invention relates to the field of meteorological early warning technology, and specifically to a tornado probability forecasting method based on ensemble Kalman filtering and large eddy simulation. Background Technology

[0002] Tornadoes are powerful whirlwinds generated by funnel-shaped clouds extending from the base of thunderstorm clouds to the ground. They are characterized by their small spatial and temporal scale, rapid growth and dissipation, and strong destructive power, which can cause significant damage and economic losses to power facilities (such as land and sea wind power plants, photovoltaic power stations, and power grid systems).

[0003] In existing technologies, high-resolution numerical models are commonly used for storm prediction. The essence of numerical models is to solve fluid dynamics and thermodynamic equations to simulate the evolution of variables such as temperature, air pressure, humidity, and wind speed in the atmosphere. Its operation process includes: integrating real-time data from satellite cloud images, radar echoes, ground observation stations, and radiosondes to generate the initial model field of the global three-dimensional atmosphere (i.e., the "current atmospheric state"); calculating nonlinear equations on a supercomputer to extrapolate changes in atmospheric state over the next few hours to days; and converting the calculation results into storm prediction results.

[0004] However, this numerical model has significant shortcomings. First, due to unavoidable errors in the initial field of the global three-dimensional atmosphere, these errors increase considerably with the extension of integration time during supercomputer calculations. Second, the parameterization of physical processes within the model also contains errors. These two types of errors significantly increase the uncertainty of the model results, making it difficult to make accurate forecasts for tornadoes with smaller spatiotemporal scales (sub-kilometer scale). Therefore, tornado forecasts often deviate significantly from actual weather conditions. Summary of the Invention

[0005] To address the aforementioned issues, this invention provides a tornado probability forecasting method based on ensemble Kalman filtering and large eddy simulation. This method can directly and effectively forecast the probability of tornado occurrence in a target area, and it has the significant advantage of high forecast accuracy.

[0006] To achieve the above objectives, the technical solution adopted by the present invention is as follows:

[0007] In a first aspect, the present invention provides a tornado probability forecasting method based on ensemble Kalman filtering and large eddy simulation, the method comprising:

[0008] Obtain a set of initial fields containing the initial fields of each of the multiple members, divide the set of initial fields into regions of different scales, and determine the first target scale region;

[0009] Perform the first data extrapolation operation on the initial field of the set to obtain the quasi-assimilated initial field of the set corresponding to the first target scale region;

[0010] The initial field of the quasi-assimilation set corresponding to the first target scale region is assimilated to obtain the final analysis field of the set containing the final analysis fields of each of the multiple members.

[0011] The first target-scale region is spatially gridded to obtain multiple grid cell columns;

[0012] For each member, the helicity information corresponding to that member is determined based on the final analysis field corresponding to that member. The helicity information includes the storm uplift helicity corresponding to each grid cell column in the first target scale region.

[0013] Based on the helicity information of each member, multiple members are discarded, and the effective grid cell column corresponding to each retained member is determined.

[0014] Based on the effective grid cell columns corresponding to each retained member, the target grid cell column is determined, and the second target scale region is delineated within the first target scale region with the target grid cell column as the center.

[0015] Based on the horizontal wind component data contained in the final analysis field of each retained member, the probability information of tornado occurrence corresponding to the second target scale region is obtained.

[0016] As a preferred embodiment of the present invention, the method for discarding members is specifically as follows:

[0017] Based on the vertical vorticity and vertical wind speed data contained in the final analysis field corresponding to a member, the storm upheurism of each grid cell column in the first target scale region is calculated.

[0018] The storm's upward helicity is compared with a helicity threshold. If the value is greater than the threshold, the first label is assigned to the grid cell column.

[0019] Count the total number of the first marker. If the total number of the first marker is zero, discard that member.

[0020] As a preferred embodiment of the present invention, the method for determining the effective grid cell column corresponding to each retained member is as follows:

[0021] For any retained member, shallow potential assessment is performed on each grid cell column within the first target scale region; the assessment region is formed with the grid cell column as the center during the assessment.

[0022] The total number of grid cell pillars with the first label in the statistical judgment area is counted. If the total number is greater than the shallow potential judgment threshold, the third label is assigned to the grid cell pillar located at the center.

[0023] From all grid cell columns with a third marker, select the grid cell column with the largest storm uphill spiral as the valid grid cell column corresponding to the retained member.

[0024] As a preferred embodiment of the present invention, the method for calculating the storm's upward helicity is as follows:

[0025] Multiple horizontal model layers are divided along the height direction of the grid cell column. The product of vertical vorticity and vertical velocity corresponding to each horizontal model layer is calculated, and all products are integrated to obtain the storm uphelicality corresponding to the grid cell column.

[0026] As a preferred embodiment of the present invention, during cyclic assimilation, linear interpolation is used to preprocess multiple types of observation data with different acquisition time intervals, so that the multiple types of observation data are input into the ensemble Kalman filter assimilation system at the same frequency.

[0027] As a preferred embodiment of the present invention, when dividing the initial field of the set into regions of different scales, a multi-layer nested grid is used for proportional division, and the region with the smallest scale is selected as the first target scale region.

[0028] As a preferred embodiment of the present invention, the specific method for obtaining the tornado occurrence probability information corresponding to the second target scale region is as follows:

[0029] For each grid cell column to be calculated in the second target scale region, obtain the horizontal wind component data contained in the final analysis field of each retained member, and calculate the synthetic wind speed corresponding to each retained member based on the horizontal wind component data.

[0030] Determine whether the composite wind speed is greater than the tornado speed threshold. If it is, assign a fourth label to the grid cell column to be calculated.

[0031] Calculate the ratio of the total number of fourth tags corresponding to each grid cell column to the total number of retained members, and use this ratio as the probability of tornado occurrence for each grid cell.

[0032] The total number of fourth tags corresponding to any grid cell column to be calculated is less than or equal to the total number of retained members;

[0033] In a second aspect, the present invention also provides an electronic device, including a processor and a memory;

[0034] The processor is connected to the memory;

[0035] Memory, used to store executable program code;

[0036] The processor runs the program corresponding to the executable program code stored in the memory to execute the aforementioned tornado probability prediction method based on ensemble Kalman filtering and large eddy simulation.

[0037] Thirdly, the present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the aforementioned method for tornado probability forecasting based on ensemble Kalman filtering and large eddy simulation.

[0038] In summary, the present invention has the following beneficial effects:

[0039] This invention expands the membership of the original ensemble initial field by using the ensemble Kalman filtering method, resulting in an ensemble initial field with a larger number of members, thereby reducing the impact of uncertainty in model forecasting caused by initial field errors on the final forecast results.

[0040] By combining multiple sets of multi-source observation data to perform cyclic assimilation of the initial field of the alignment assimilation set, the matching degree between the initial field and the actual atmospheric state is greatly enhanced, providing a prerequisite for subsequent large eddy simulation calculations. Compared with the existing technology that usually uses tornado diagnostic parameters, such as storm uplift helicity, to predict tornadoes, this invention innovatively proposes a shallow potential set diagnostic algorithm for tornado generation to determine nested regions. By further improving the resolution through model nesting, the method of calculating the probability of tornado occurrence in the region corresponding to each grid cell column in the nested region makes up for the shortcomings of the existing technology. It can directly calculate the probability of tornado occurrence in the target region instead of indirect prediction through tornado diagnostic parameters. It has the significant advantage of high prediction accuracy and has high practical value. Attached Figure Description

[0041] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0042] Figure 1 This is a flowchart of the method of the present invention;

[0043] Figure 2 This is a schematic diagram of a judgment region centered at (i, j) in one embodiment. Detailed Implementation

[0044] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed merely to enable those skilled in the art to better understand and implement the subject matter described herein, and are not intended to limit the scope, applicability, or examples set forth in the claims. The function and arrangement of the elements discussed may be changed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the various examples. For example, the described methods may be performed in a different order than described, and steps may be added, omitted, or combined. Furthermore, features described in some examples may be combined in other examples.

[0045] like Figure 1 As shown, this invention provides a tornado probability forecasting method based on ensemble Kalman filtering and large eddy simulation, the method comprising:

[0046] Obtain a set of initial fields containing the initial fields of each of the multiple members, divide the set of initial fields into regions of different scales, and determine the first target scale region;

[0047] An ensemble initial field is a set of multiple initial fields, which refers to a set of initial atmospheric states containing multiple members. Each member represents a possible initial condition and corresponds to an initial field. Each initial condition contains various meteorological parameters such as temperature, humidity, and wind field. The set of multiple members is used to quantify the uncertainty of the initial conditions.

[0048] In this embodiment, the ensemble initial field adopts the ensemble forecast field of a global model. Specifically, it can be selected from the GFS analysis field and forecast field of the global weather forecast model developed by the world's leading meteorological agency, NCEP; or the next-generation reanalysis dataset ERA5 developed by ECMWF. The respective ensemble forecast versions of the two, GEFS and ERA5 Ensemble Forecasts, contain 21 and 51 members, respectively.

[0049] Since the directly acquired initial ensemble field corresponds to a global scale, while the target area for forecasting is often focused on a scale of tens of kilometers, it is necessary to divide the scale region corresponding to the initial ensemble field (e.g., mesoscale region division in existing technologies) to reduce the computational load on supercomputers. Furthermore, after region division, the resolution represented by each individual region unit gradually increases, thereby improving the computational accuracy of supercomputers. Specifically, in this embodiment, WRF-ARW is used to downscale the initial field of the global model. WRF-ARW is a numerical weather prediction system provided by world-renowned meteorological research institutions, primarily used for short-term weather forecasting simulation, atmospheric process simulation, and long-term climate simulation. It features high accuracy and high spatiotemporal resolution, making it an important tool for meteorological and atmospheric research.

[0050] When dividing the mesoscale regions, WRF-ARW was used to create a three-layer nested network for regional division and ensemble forecasting, resulting in three regions with a resolution ratio of 9:3:1: region F1 (first layer), region F2 (second layer), and region F3 (third layer). Each region is divided into multiple grid cells, with side lengths of 4050m, 1350m, and 450m respectively, showing a progressively increasing resolution across the three layers. Taking my country as an example, in terms of scope, region F1 (first layer) focuses on the Asian tectonic plate level, allowing the selection of the entire Asian region as the target; region F2 (second layer) focuses on the national land area level, allowing the selection of specific regions as targets (e.g., Northeast, Northwest, North China, Southwest, East China, or South China, depending on the location of the tornado prediction target); and region F3 (third layer) focuses on the provincial level, allowing the selection of the province where the target is predicted.

[0051] In this embodiment, the third layer region F3 (resolution 450 m) with the smallest scale is selected to correspond to the first target scale region D1. The subsequent initialization of the ensemble initial field mode and forward integration are all performed within the range of the first target scale region D1. Therefore, only the data of the corresponding region and scale in each member are extracted.

[0052] The first target scale region D1 is divided into grid cells, specifically in the form of grid cell columns. This is because the data corresponding to each member is mapped to three-dimensional space on a spatial scale. Therefore, at the vertical height above the horizontal bottom surface of the grid cell, multiple horizontal model layers are further divided according to different isobaric surfaces. This results in the grid cell column being divided into multiple approximately cubic data space structures at the vertical height. The data of the members is arranged on the edges, faces, and points of the data space structure, and is calculated by a supercomputer. In this embodiment, for the height space at a specific pressure surface height (e.g., horizontal altitude 500m), the number of horizontal model layers is artificially increased to make the horizontal model layers at that pressure surface height denser. This will better characterize the development process of vortices near the ground, especially the generation process of frictional vorticity, which is very important for tornado formation.

[0053] It should also be noted that although the Earth is spherical and the areas of the grid regions divided based on latitude and longitude are not equal, the data provided by the members of the initial field of the ensemble obtained in this invention has already undergone digital terrain transformation. Therefore, in subsequent calculations, the bottom edge of the grid cell column can be set as a square to participate in the calculation.

[0054] The first data extrapolation operation is performed on the initial field of the set, specifically by initializing the mode and performing forward extrapolation (integration) to obtain the quasi-assimilated initial field of the set corresponding to the first target scale region D1. The number of members included in the quasi-assimilated initial field of the set has been expanded, as detailed below.

[0055] Since this invention provides a weather phenomenon forecasting system, it is necessary to forecast a certain time in advance (using supercomputer integral extrapolation). Furthermore, the members in the acquired initial ensemble field are not significantly different. To improve prediction accuracy, it is necessary to simulate and extrapolate the possible development forms of the weather system, thereby improving the accuracy of subsequent tornado probability predictions. This necessitates expanding the number of members to better quantify the uncertainty of the model results. In this embodiment, ten different physical parameterization schemes are designed to specifically configure the 20 members to achieve the expansion purpose, as shown in Table 1, where GEFSmem represents different members.

[0056] Table 1 Scheme Configuration Table

[0057] Group Member Group Microphysics scheme Road surface mode Near-surface layer scheme Boundary layer scheme Longwave and shortwave radiation schemes Configuration 1 GEFS mem01-mem04 Milbrandt–Yau Pleim–Xiu LSM Pleim–Xiu ACM2 (Pleim) CAM Configuration 2 GEFS mem01-mem04 Milbrandt–Yau Pleim–Xiu LSM Pleim–Xiu ACM2 (Pleim) GFDL Configuration 3 GEFS mem05-mem08 Milbrandt–Yau Pleim–Xiu LSM Monin–Obukhov MYJ GFDL Configuration 4 GEFS mem05-mem08 Morrison Noah LSM MM5 Monin–Obukhov YSU RRTM Long Wave Dudhia Short Wave Configuration 5 GEFS mem09-mem12 Morrison Pleim–Xiu LSM Pleim–Xiu ACM2 (Pleim) RRTM Long Wave Dudhia Short Wave Configuration 6 GEFS mem09-mem12 Thompson Noah LSM MM5 Monin–Obukhov YSU Goddard Configuration 7 GEFS mem13-mem16 Morrison Noah LSM QNSE QNSE CAM Configuration 8 GEFS mem13-mem16 WSM 6-class Noah LSM MYNN MYNN 2.5 level TKE Goddard Configuration 9 GEFS mem17-mem20 WDM 6-class Pleim–Xiu LSM Pleim–Xiu ACM2 (Pleim) Goddard Configuration 10 GEFS mem17-mem20 WDM 6-class RUC LSM MYNN MYNN 2.5 level TKE GFDL

[0058] As can be seen, in this embodiment, the initial number of members obtained based on GEFS is 20. Four members are grouped together, resulting in five member groups. Each member group is assigned two physical process parameterization schemes, forming ten configuration groups in total. This expands the number of members from 20 to 40. Due to the existence of initial field errors and model parameterization scheme errors, insufficient members in the ensemble forecast increase the uncertainty in the final forecast results. This embodiment expands the number of members through different physical parameterization scheme configurations and multiple initial conditions, effectively reducing the uncertainty of the final forecast results caused by these two types of errors and improving prediction accuracy.

[0059] After determining the ensemble forecast model configuration and expanding the membership, model initialization can be performed. In this embodiment, model initialization is set to occur 6 hours before the start of assimilation, i.e., 6-hour forward forecast calculations are performed on 40 members. This is because the spatiotemporal scale of the input multi-source observation data is relatively small in the subsequent Kalman filter cyclic assimilation step. Therefore, during the 6-hour model initialization process, it is necessary to generate corresponding small- and medium-scale information in region F3 through multi-layer nested downscaling, and finally obtain the quasi-assimilation ensemble initial field under the first target scale region D1 corresponding to the third layer region F3. The members in this quasi-assimilation ensemble initial field can be used for Kalman filter cyclic assimilation with the multi-source observation data with small spatiotemporal scales.

[0060] Next, the initial field of the quasi-assimilation set corresponding to the first target scale region is assimilated to obtain the final analysis field of the set containing the final analysis fields of each of the multiple members.

[0061] The purpose of this step is to obtain the final analysis field of the ensemble, which matches the actual atmospheric state significantly better than the initial field of the quasi-assimilated ensemble. This can narrow the divergence of ensemble members in the forecast calculation process. In meteorology, ensemble Kalman filtering (EnKF) is mainly applied to large-scale data assimilation problems, where observational data can be treated as random variables affected by disturbances. Generally, EnKF-based correlation methods aim to simplify or realize the calculation of error statistics related to flow dependence. Unlike traditional methods, the EnKF method does not solve the time evolution equation of the model state probability density function, but instead uses the Monte Carlo method to estimate the error statistics in the forecast process. During assimilation, the EnKF method uses the dynamic equation to perform time integration on a large ensemble of model states, and then calculates the probability density function matrix at different times based on this ensemble. Compared to the four-dimensional variational assimilation (4DVar) method, the EnKF data assimilation method does not require an adjoint model in highly nonlinear microphysical processes at the convective scale, and is therefore more suitable.

[0062] This embodiment is based on the ARPS EnKF data assimilation system, using multi-source observation data to perform cyclic assimilation on the initial field of the alignment assimilation set in the first target-scale region D1. Before assimilation, the model variables in each member need to be projected into the observation variable space using the forward observation operator for direct assimilation. The assimilated observation data mainly includes radar observation data, surface meteorological observation data, satellite observation data, and radiosonde observation data. The radar observation data specifically is Doppler radar data, which includes radar reflectivity factor and radar radial wind; the surface meteorological observation data specifically is encrypted ground station data, which includes variables such as wind direction, wind speed, temperature, dew point temperature, and air pressure. After passing through the quality control module of the EnKF data assimilation system, these data, along with the model variables transformed by the observation operator, enter the EnKF analysis module for assimilation analysis.

[0063] During assimilation, different observational data will employ different assimilation time intervals based on their varying spatiotemporal scales, resolutions, and characteristics. Due to the extremely small spatiotemporal scale and rapid formation and dissipation of tornadoes, radar data is the most crucial element in the assimilation process. The entire assimilation framework should be designed around radar data assimilation. For example, a 6-minute cyclic interval for radar data assimilation ensures that storm information in the model is updated in real-time and fully constructed. Considering the potential for missing lower-level storm data when the storm is far from the radar, encrypted ground meteorological observation data is needed to supplement this information. Typically, if the temporal resolution of ground meteorological observation data differs from that of radar observation data (e.g., a 5-minute interval), this embodiment also employs linear interpolation to input ground meteorological observation data and radar observation data at the same frequency for rapid 6-minute cyclic assimilation. For satellite observation data, the temporal resolution is usually lower; for example, the highest temporal resolution of observation data from my country's FY4 series geostationary satellites is a 1-hour interval. Therefore, when using different assimilation time intervals, a 1-hour cyclic input can be used. Since radiosonde data and wind profiler radar data mainly reflect weather-scale information, they are also input cyclically at one-hour intervals here.

[0064] Cyclic assimilation is performed on the first target scale region D1. After the first assimilation of the quasi-assimilation ensemble initial field, the system will perform at least 6 minutes of near-term ensemble forecasting according to the configuration in Table 1 above. The quasi-assimilation ensemble initial field obtained after the first assimilation will be used as the "input" for the next assimilation process, thus becoming the ensemble initial field for the next assimilation process. This process is repeated multiple times, assimilating the data with multiple sets of multi-source observation data at different times. This makes the data in the quasi-assimilation ensemble initial field increasingly closer to the near-real atmospheric state reflected by the multi-source observation data. The final assimilation result is the final ensemble analysis field. Because it has undergone multiple cyclic assimilation processes, the final ensemble analysis field has a much higher degree of matching with the real atmospheric state than the quasi-assimilation ensemble initial field. Furthermore, the final ensemble analysis field can maintain usable accuracy for at least 3 hours, providing data support for subsequent tornado eddy simulations. For example, the final ensemble analysis field obtained at 17:00 after cyclic assimilation can be used for the probability forecast of tornadoes within the next 3 hours.

[0065] The first target scale region is spatially gridded to obtain multiple grid cell columns; the horizontal base of each grid cell column is a square with equal side length.

[0066] Then, using the 40 members in the final analysis field of the ensemble, the probability of tornado occurrence in the target area is predicted. The target area is located in the first target scale region D1. However, since the resolution of the first target scale region D1 is only 450m, and a reliable tornado forecast requires a regional resolution of at least 50m, it is necessary to further improve the regional resolution to determine the location of the second target scale region D2.

[0067] Determining the location of the second target-scale region D2 requires the following steps.

[0068] First, for each member, the helicity information corresponding to that member is determined based on the final analysis field corresponding to that member. The helicity information includes the storm uplift helicity corresponding to each grid cell column in the first target scale region.

[0069] This embodiment provides a tornado generation shallow potential ensemble diagnostic algorithm based on 0-6km storm updraft helicity to determine the location of nested simulation regions. Research shows that the magnitude of the maximum 0-6km storm updraft helicity is highly correlated with the occurrence of strong convective processes. The formula for calculating the 0-6km storm updraft helicity UH is:

[0070]

[0071] Wherein, UH represents the storm's upward helicity from 0 to 6 km, with the integral level ranging from 0 km to 6 km. It is vertical vorticity. Vertical velocity.

[0072] The final analysis field comprises 40 members, each containing vertical vorticity. and vertical velocity Therefore, for any grid cell column, the storm uphill spiral degree UH corresponding to each of the 40 members can be calculated.

[0073] In summary, the storm uplift helicity UH is calculated for all grid cell columns in the first target scale region D1. Specifically, this is based on the vertical vorticity provided by individual members in the final analysis field. and vertical wind speed Calculate the product, then integrate sequentially to calculate the storm uplift helicity UH of all grid cell columns in the first target scale region D1, and then compare the storm uplift helicity UH with the helicity threshold F. uh The comparison is performed. If the value is greater than the value, the first label Z1 is assigned to the grid cell column; otherwise, the second label Z2 is assigned. The helicity threshold F is used. uh A comparison threshold set manually can be obtained based on publicly available data from the World Meteorological Organization.

[0074] It is important to note that among the 40 members mentioned above, some exhibit significant deviations and therefore need to be discarded to reduce interference with subsequent judgments and forecasts, and to decrease computational load. The specific principle for discarding these members is as follows: for each member, after calculating the storm upheurism UH for all grid cell columns, obtain the total number of first markers. If the total number of first markers is zero, then that member is discarded. In other words, based on the information provided by this member, the storm upheurism UH for all grid cell columns is less than the helicity threshold F. uh In this case, the member can be considered to have no predictive significance, and the member to be retained can be obtained. The number of retained members is set to N, where N≤40.

[0075] Next, for any retained member, shallow potential assessment is performed on each grid cell column within the first target scale region D1; during the assessment, a assessment region is formed centered on the grid cell column; the total number of grid cell columns with the first label Z1 in the assessment region is counted, and if the total number is greater than the shallow potential assessment threshold F... T Then, assign the third marker Z3 to the grid cell column located at the center; from all grid cell columns with the third marker Z3, select the grid cell column with the largest storm upheaval UH as the valid grid cell column corresponding to the reserved member.

[0076] The specific implementation method of the above steps is given below:

[0077] For each grid cell column, a shallow potential assessment is performed. The assessment region is centered on the grid cell column, and a ring-shaped region of equal units can be captured around it. The total number of grid cell columns with the first marker within the assessment region is then counted. If the total number exceeds the shallow potential assessment threshold F... T Then, the third label Z3 is assigned to the central grid cell column.

[0078]

[0079] Where, if represents an "if" logical judgment; the value of the first flag Z1 is considered to be 1, and the value of the second flag Z2 is considered to be 0; k is the current member; This represents the sum of the marker values ​​in the judgment region centered at (i, j); for example... Figure 2 As shown, (i, j) represents the horizontal position of the grid cell column; (a, b) is the horizontal position range of the judgment region, where R represents the horizontal base side length of the grid cell column; F uh The helicity threshold; This represents the storm uplift spiral of each individual grid cell column in the judgment region; the total number of grid cell columns in the judgment region is (2R+1). 2 .

[0080] For example, when k=1, R=1, F uh =200m 2 s -2 hour, The square judgment region consists of 9 grid cells centered at the horizontal position (i, j). The storm uphill spiral is calculated independently for each grid cell. and compared it with 200m 2 s -2 A comparison is performed, assigning the first label Z1 to the grid cell cylinders whose comparison result is greater than or equal to the value, and assigning the second label Z2 to the others. Then, the values ​​of all labels are summed to obtain the result. ,Will Shallow potential judgment threshold F T If the former is greater than the latter, then assign the third label Z3 to the grid cell column of (i, j).

[0081] For example, given R=1, when it is necessary to determine whether the grid cell column with a horizontal position of (22, 22) corresponding to the first target scale region D1 and the 16th member has a shallow potential for tornado generation (i.e., it is necessary to determine whether the grid cell column has a third marker), it is only necessary to calculate And compare it with the shallow potential judgment threshold F T (Assuming = 5) Compare the results. If, among the 9 grid cell columns in the judgment region, 6 grid cell columns possess the first marker Z1, then... =6, which is greater than the shallow potential judgment threshold F. T =5, at this time the grid cell column of (22, 22) is divided into the third mark Z3.

[0082] Next, from all grid cell columns with the third marker Z3 (i.e. all grid cell columns with shallow tornado generation potential), the grid cell column with the largest storm upheurism UH is selected as the valid grid cell column corresponding to the retained member.

[0083] Based on the above calculations, for a selected retained member, after traversing and calculating all grid cell pillars, only some grid cell pillars are assigned the third marker Z3. Then, from all grid cell pillars possessing the third marker Z3, the grid cell pillar with the largest storm uphill helix UH is selected. This grid cell pillar is defined as the effective grid cell pillar corresponding to the retained member. Therefore, each member that has not been discarded corresponds to one effective grid cell pillar. The specific calculation method is as follows:

[0084]

[0085] Among them, (i k j k) represents the horizontal position of the effective grid cell column; Max indicates finding the maximum storm upheaven spiral UH; the where function indicates finding the storm upheaven spiral UH of the grid cell column with shallow tornado generation potential in member k; Loc indicates the positioning action.

[0086] Then, based on the effective grid cell columns corresponding to each retained member, the target grid cell column is determined, and the second target area is delineated within the first target area with the target grid cell column as the center;

[0087] The specific method is to obtain the effective grid cell pillars corresponding to all retained members in the final analysis field of the set, perform set averaging calculation, find the target grid cell pillar, and the target grid cell pillar is the nested positioning center of the second target scale region.

[0088]

[0089] Among them, (i c j c ) represents the horizontal position of the target mesh cell column; N represents the number of members to be retained.

[0090] Determine the target mesh cell column (i c j c After that, the second target scale region D2 is delineated with the target grid unit column as the center. Specifically, since the target grid unit column is a grid unit column in the first target scale region D1, and its horizontal base side length is 450m, which does not meet the accuracy of the subsequent tornado occurrence probability forecast, the second target scale region D2 is divided into multiple grid unit columns with a horizontal base side length of 50m at a resolution of 50m. The area of ​​the second target scale region D2 is determined according to the number of newly divided grid unit columns in the second target scale region D2, for example, 120×120km.

[0091] Because the second target scale region D2 fully considers the shallow potential for tornado formation, subsequent calculations only need to be performed on each grid cell column in the second target scale region D2 independently, without needing to calculate other areas outside the second target scale region D2. This can greatly save supercomputer computing power and improve forecast timeliness.

[0092] Finally, based on the horizontal wind component data contained in the final field of each of the retained members, the probability information of tornado occurrence corresponding to the second target area is obtained.

[0093] Specifically, this embodiment provides a method for probabilistically predicting tornadoes based on the wind speed prediction results of members, combined with a grid-by-grid probabilistic prediction method. The specific method is as follows: for each grid cell column to be calculated in the second target area, obtain the horizontal wind component data contained in the final analysis field corresponding to all retained members, and calculate the synthetic wind speed corresponding to each retained member based on the horizontal wind component data; determine whether the synthetic wind speed is greater than the tornado wind speed judgment threshold. If it is greater, assign a fourth label to the grid cell column to be calculated; calculate the ratio of the total number of fourth labels corresponding to each grid cell column to be calculated to the total number of retained members, as the tornado occurrence probability corresponding to each grid cell to be calculated.

[0094] The following examples will provide further explanation.

[0095] Select the column of grid cells to be calculated in the second target scale region (i q i e The composite wind speed is calculated sequentially based on the horizontal wind speed components provided by all retained members, and the composite wind speed is compared with the tornado speed by the threshold EF0. min The magnitude of the combined wind speed is determined if it exceeds the tornado speed threshold EF0. min (For example, set to 29 ms) -1 If the value of the grid cell column to be calculated is less than 1, then a fourth label Z4 is assigned to the column, indicating that the column has reached the minimum tornado speed level with the data of the retained members. Otherwise, a fifth label Z5 is assigned.

[0096]

[0097] Here, if indicates a conditional logical judgment; the value of the fourth flag Z4 is considered 1, and the value of the fifth flag Z5 is considered 0; k is the current member; This represents the label value of the column (q, e) of the mesh element to be computed under the k-th member; EF0 min Determine the threshold for tornado speed. Iterate through all retained members and count the total number of fourth markers Z4 obtained for the grid cell column to be calculated.

[0098] Next, based on all retained members, probabilistic calculations are performed on all the grid cell cylinders to be calculated in the second target scale region D2, specifically as follows:

[0099]

[0100] in, This represents the probability of a tornado occurring in the region corresponding to the grid cell column with horizontal position (q, e). This represents the label value of the grid cell column (q, e) to be calculated under the k-th member.

[0101] The above formula means that the probability of a tornado occurring in the area corresponding to a grid cell column to be calculated is numerically equal to the probability of the combined wind speed being greater than the tornado speed threshold EF0 among all retained members N. min The number of members is the proportion of all retained members N. The higher this proportion, the more retained members in this region support their prediction of a tornado occurring. Since the resolution of the grid cell columns to be calculated has reached the 50 m level, the computational requirements for explicit tornado prediction are met.

[0102] In another possible embodiment, before the second target scale region D2, there is a step of determining the quasi-second target scale region D2'. The quasi-second target scale region D2' and the second target scale region D2 are established based on the same center, with the former having a larger area and a smaller resolution than the latter. The quasi-second target scale region D2' is gridded using the same method, with the grid side length set to 150m. Before making a probabilistic forecast for the second target scale region D2, a probabilistic forecast can be made for each grid cell column in the quasi-second target scale region D2', and the forecast time interval between the two regions is set to 5min-15min. This allows sufficient development time for near-surface vortices in the quasi-second target scale region D2', which is beneficial for providing more scientific and complete tornado forecast information.

[0103] In another possible embodiment, the present invention also provides an electronic device, including a processor and a memory;

[0104] The processor is connected to the memory;

[0105] Memory, used to store executable program code;

[0106] The processor runs the program corresponding to the executable program code stored in the memory to execute the aforementioned tornado probability prediction method based on ensemble Kalman filtering and large eddy simulation.

[0107] In another possible embodiment, the present invention also provides a computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the above-described method for tornado probability forecasting based on ensemble Kalman filtering and large eddy simulation.

[0108] Several embodiments of this disclosure have been described above. These descriptions are exemplary and not exhaustive, and are not limited to the disclosed embodiments. Many modifications and variations will be apparent to those skilled in the art without departing from the scope and spirit of the described embodiments. The terminology used herein is chosen to best explain the principles, practical applications, or technological improvements to the embodiments in the market, or to enable others skilled in the art to understand the embodiments disclosed herein.

Claims

1. A tornado probability prediction method based on ensemble Kalman filter and large eddy simulation, characterized in that, The method comprises: acquiring a set initial field comprising a plurality of initial fields corresponding to a plurality of members respectively, and performing regional division of different scales on the set initial field, and determining a first target scale region; performing a first data deduction operation on the set initial field to obtain a quasi-assimilation set initial field corresponding to the first target scale region; performing an assimilation operation on the quasi-assimilation set initial field corresponding to the first target scale region to obtain a set final field comprising a plurality of final fields corresponding to the plurality of members respectively; performing spatial gridding on the first target scale region to obtain a plurality of grid cell columns; for each member, determining, based on the final field corresponding to the member, a helicity information corresponding to the member, the helicity information comprising storm updraft helicity corresponding to each grid cell column in the first target scale region; performing a discard processing on the plurality of members based on the helicity information corresponding to each member respectively, and determining an effective grid cell column corresponding to each retained member respectively; determining a target grid cell column based on the effective grid cell column corresponding to each retained member respectively, and demarcating a second target region in the first target region with the target grid cell column as the center; acquiring tornado occurrence probability information corresponding to the second target region based on horizontal wind component data contained in the final field corresponding to each retained member respectively.

2. The tornado probability prediction method based on ensemble Kalman filter and large eddy simulation according to claim 1, characterized in that, The discard processing on the member is specifically: calculating storm updraft helicity of each grid cell column in the first target scale region based on vertical vorticity data and vertical wind speed data contained in the final field corresponding to a member; comparing the storm updraft helicity with a helicity threshold value, if greater, assigning a first mark to the grid cell column; statistically counting the total number of the first marks, if the total number of the first marks is zero, discarding the member.

3. The method of claim 2, wherein the method is based on an ensemble Kalman filter and large eddy simulation. The method for determining the effective grid cell column corresponding to each retained member is specifically: for any retained member, performing a shallow potential judgment on each grid cell column in the first target scale region; when judging, forming a judgment region with the grid cell column as the center; statistically counting the total number of grid cell columns with the first mark in the judgment region, if the total number is greater than a shallow potential judgment threshold value, assigning a third mark to the grid cell column located at the center; selecting, from all grid cell columns with the third mark, a grid cell column with the maximum storm updraft helicity as the effective grid cell column corresponding to the retained member.

4. The tornado probability prediction method based on ensemble Kalman filter and large eddy simulation according to claim 3, characterized in that, The calculation method of the storm updraft helicity is specifically: dividing a plurality of horizontal pattern layers in the height direction of the grid cell column, calculating the product of the vertical vorticity and the vertical speed corresponding to each horizontal pattern layer, and performing integral calculation on all products to obtain the storm updraft helicity corresponding to the grid cell column.

5. The tornado probability prediction method based on ensemble Kalman filter and large eddy simulation according to claim 1, characterized in that, When performing regional division of different scales on the set initial field, a multi-layer nested network is used for equal proportion division, and the smallest scale region is selected as the first target scale region.

6. The tornado probability prediction method based on ensemble Kalman filter and large eddy simulation according to claim 1, characterized in that, The specific method for acquiring tornado occurrence probability information corresponding to the second target region is: For each to-be-calculated grid cell column in the second target area, horizontal wind component data contained in a final field corresponding to each reserved member is obtained, and a synthetic wind speed corresponding to each reserved member is calculated based on the horizontal wind component data; It is judged whether the synthetic wind speed is greater than a tornado speed judgment threshold, and if so, a fourth mark is assigned to the to-be-calculated grid cell column; A ratio of a total number of fourth marks corresponding to each to-be-calculated grid cell column to a total number of reserved members is calculated as a tornado occurrence probability corresponding to each to-be-calculated grid cell.

7. The tornado probability prediction method based on ensemble Kalman filter and large eddy simulation according to claim 1, characterized in that, The first data deduction operation comprises a mode initialization step and a forward deduction calculation step.

8. An electronic device, comprising: The processor and the memory are connected. The processor and the memory are connected. The memory is used for storing executable program codes. The processor runs a program corresponding to the executable program codes by reading the executable program codes stored in the memory, so as to execute the tornado probability prediction method based on the ensemble Kalman filtering and large eddy simulation according to any one of claims 1-7.

9. A computer-readable storage medium, characterized in that, The computer program is stored on the memory and is executed by the processor to implement the tornado probability prediction method based on the ensemble Kalman filtering and large eddy simulation according to any one of claims 1-7.

Images