A three-dimensional radar echo extrapolation method based on deformation convolution dynamic field constraint

Abstract

The application relates to a three-dimensional radar echo extrapolation method based on a deformation convolution dynamic field constraint, which comprises the following steps: S1, continuously collecting radar reflectivity and radial velocity data within 2 hours; S2, converting the radar reflectivity data into three-dimensional data in a Cartesian rectangular coordinate system, mapping the radial velocity data into a radial rectangular coordinate system and constructing a sample set; S3, acquiring middle-layer radial convergence feature information in a radial velocity graph; S4, building a radar echo extrapolation model; S5, introducing deformation convolution to automatically extract the dynamic field information of a convection system; S6, dividing the sample set into a training set, a verification set and a test set; and S7, obtaining a reflectivity image in the next 1 hour. The application fully considers the three-dimensional structural features and dynamic characteristics of a weather system, simultaneously combines multi-source physical constraint information, and realizes accurate simulation and prediction of a complex evolution process of the weather system.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of weather radar, and relates to technologies related to radar signal processing, echo extrapolation and radar data analysis, in particular to a three-dimensional radar echo extrapolation method based on deformation convolution dynamic field constraint. BACKGROUND

[0002] The development and evolution of weather systems are accompanied by complex atmospheric activities, in which horizontal atmospheric flow dominates the advection motion, and vertical convection plays a key role in the life cycle of the weather system. The existing radar echo extrapolation technology has deficiencies in simulating the complex changes such as generation and disappearance of weather systems. The traditional method relies on radar reflectivity to estimate the cloud cluster motion vector, which is only suitable for simple linear mode, while deep learning can capture non-linear motion, but still needs to be improved in reflecting echo generation and disappearance.

[0003] In view of this, in order to improve the accuracy and reliability of weather radar echo extrapolation, it is necessary to analyze the weather system as a three-dimensional whole, so that the model can deeply learn the spatio-temporal correlation between different height layers. At present, some scholars have carried out research on this, but most of the models are only based on historical radar echo data, without fully considering the environmental thermal instability, vertical wind shear, water vapor vertical distribution and other physical constraint information, while these factors are crucial to the development and evolution of severe convective weather. Recent studies have begun to incorporate multi-source data, but there is still a lack of related research on integrating the strong downdraft within the convective system, i.e. the mid-level radial convergence (MARC) feature, into the extrapolation model. MARC, as an important indicator of the dynamics of convective systems, is of great significance to predicting the trajectory and speed of cloud clusters, and its inclusion can improve the accuracy of radar extrapolation technology in predicting high-intensity reflection areas.

[0004] In order to further improve the accuracy of radar extrapolation technology in predicting high-intensity reflection areas, a three-dimensional radar echo extrapolation method that integrates radial velocity information is urgently needed. SUMMARY

[0005] The purpose of the present application is to overcome the shortcomings of the prior art and provide a three-dimensional radar echo extrapolation method based on deformation convolution dynamic field constraint, which fully considers the three-dimensional structural features and dynamic characteristics of weather systems, and combines multi-source physical constraint information to achieve accurate simulation and prediction of the complex evolution process of weather systems.

[0006] The technical problem of the present application is solved by the following technical solution:

[0007] A three-dimensional radar echo extrapolation method based on deformation convolution dynamic field constraint, the steps of the method are:

[0008] S1, the radar operates in VCP21 body scan mode, that is, 9 scans of different elevations are completed every 6 minutes, and radar reflectivity and radial velocity data are collected continuously for 2 hours;

[0009] S2, for the data of each radar scan time: convert the radar reflectivity data from polar coordinates into three-dimensional data in the Cartesian rectangular coordinate system; perform de-aliasing on the radial velocity data, and map the radial velocity data into the radial rectangular coordinate system; use the reflectivity and radial velocity images in the first hour as the input of the extrapolation model, and use the reflectivity image in the second hour as the label to build a sample set;

[0010] S3, use the middle layer radial convergence automatic identification algorithm to obtain the middle layer radial convergence feature information in the radial velocity image;

[0011] S4, build a radar echo extrapolation model taking three-dimensional radar reflectivity and radial velocity images in one hour as input and taking reflectivity images in the subsequent 1 hour as output;

[0012] S5, introduce deformation convolution in the radar echo extrapolation model to automatically extract the dynamic field information of the convective system;

[0013] S6, divide the sample set into a training set, a validation set and a test set, train the radar echo extrapolation model on the training set, optimize the hyperparameters on the validation set, and test the radar echo extrapolation model on the test set;

[0014] S7, use the radar reflectivity data in any one hour, convert the radar reflectivity data at each radar scan time from polar coordinates to three-dimensional data in the Cartesian rectangular coordinate system using the steps of S2, perform de-aliasing on the radial velocity data at each elevation, map the radial velocity data to the radial rectangular coordinate system, and use the above two kinds of data in this hour as the input of the model, that is, the reflectivity image in the subsequent 1 hour can be obtained.

[0015] Moreover, the S2 converts the radar reflectivity data at each radar scan time from polar coordinates to three-dimensional data in the Cartesian rectangular coordinate system by first converting the radar base data from polar coordinates to Cartesian coordinates using bilinear interpolation in the horizontal dimension, and then interpolating in the vertical dimension to obtain three-dimensional equidistant grid point field data;

[0016] The S2 performs de-aliasing on the radial velocity data at each elevation refers to solving the velocity aliasing problem using a de-velocity aliasing algorithm based on minimizing the velocity difference between regions, and mapping the de-aliasing processed radial velocity data to the radial rectangular coordinate system.

[0017] Moreover, the S4 adopts an "independent encoding and fusion" strategy to process ten consecutive three-dimensional radar echo images within the past one hour, and the reflectivity data and the mid-level radial convergence structure information are respectively mapped to features by respective encoders, and the two types of feature information are integrated by a decoder to generate an extrapolated echo sequence in the future, and the radar echo extrapolation model completes the effective extrapolation of the three-dimensional radar echo data by jointly using a ConvLSTM network based on 3D convolution and a convolutional neural network CNN;

[0018] The ConvLSTM network based on 3D convolution replaces the two-dimensional convolution operation of the ConvLSTM unit with a three-dimensional convolution operation to process the three-dimensional radar data corresponding to each radar scanning time, and uses the cyclic structure unique to the ConvLSTM unit to extract and synthesize the radar three-dimensional data at different times within an hour; each ConvLSTM unit contains three key components, a forget gate, an input gate and an output gate, which receives the data input X t at the current time point, the state output H t-1 at the previous time point and the long-term memory unit C t-1 as inputs, and generates the updated long-term memory unit C t and the state output H t at the current time point as outputs, and the operations performed by the ConvLSTM unit and the calculation of the output results are as follows:

[0019] i t =σ(W xi *X t +W hi *H t-1 +W ci *C t-1 +b i )

[0020] f t =σ(W xf *X t +(W hf *H t-1 +W cf *C t-1 +b f )

[0021] o t =σ(W xo *X t +W ho *H t-1 +W co *C t-1 +b o )

[0022]

[0023] where * represents convolution operation, represents Hadamard product, σ is activation function, W and b represent corresponding convolution kernel weight and bias term respectively, tanh is hyperbolic tangent function, i t , f t and o t represent input gate, forget gate and output gate respectively.

[0024] Moreover, the deformed convolution of S5 refers to adjusting the convolution kernel form according to actual needs, further expanding the extraction and representation of effective speed characteristics and its associated information, and the two ways of deformed convolution are:

[0025] (1) Based on the consideration of the maximum statistical distance between positive and negative speed points, the height K H of the convolution kernel in each row is fixedly set to 15, and the width K W of the convolution kernel is dynamically calculated and determined according to its specific radial position in the image by the following formula;

[0026]

[0027] where r represents the distance of the center of the convolution kernel from the top of the radial velocity image;

[0028] In order to realize the quantification of the feature parameters including the maximum speed difference dv max and the average speed difference dv mean , the weight distribution in the convolution kernel is determined by the following formula;

[0029]

[0030] where i represents the distance of the pixel point from the top of the template, j represents the distance of the pixel point from the left side of the template, (i,j) is the coordinate position of the pixel point in the template, Ω is the relevant area of feature calculation, K l (i,j) is the weight size at (i,j) coordinate;

[0031] (2) The template matching method is adopted to implement the partitioned speed matching and calculation, limit the distance between the speed extreme points, and evaluate the relative position between the center of the maximum inflow speed and the center of the outflow speed, accurately identify the key speed characteristics in the airflow field, and the airflow field speed characteristic identification template is a regular hexagon, which is specially used to identify the 4 types of airflow field structures of convergence, divergence, cyclone and anticyclone;

[0032] Based on the consideration of statistical distance, the deformed convolution for determining the rotation type is fixedly set to 11 in the height K H of the convolution kernel in each row, and the width K W of the convolution kernel is dynamically calculated and determined according to its specific position in the image by the following formula;

[0033]

[0034] where h e {-5, -4, -3,..., 3, 4, 5}.

[0035] Further, the S6 trains the radar echo extrapolation model using a weighted mean squared error (WMSE) as a model loss function to guide the model to pay more attention to the prediction accuracy of the strong echo area, and the WMSE is defined as follows:

[0036]

[0037] where I h,x,y represents the reflectivity size at the coordinate (h, x, y) in the label image, represents the reflectivity size at (h, x, y) in the predicted image, Weight h,x,y represents the weight value of the position, and different weights can be set according to the different echo intensities, so that the model pays more attention to the strong echo area.

[0038] The positive effects that can be produced by the present application are:

[0039] The present application analyzes the weather system as a three-dimensional whole, and incorporates the mid-level radial convergence feature in the convective system to improve the accuracy of predicting the high-intensity reflection area, uses a 3D convolution ConvLSTM network and a convolutional neural network CNN to build a radar echo extrapolation model, adopts an independent coding and fusion strategy to code the reflectivity data and the MARC structure information respectively and then integrates them by a decoder, and introduces a deformation convolution method to expand the extraction and representation of the effective velocity feature and its associated information. Through the convergence feature effectiveness test and the deformation convolution test verification, the present application can significantly improve the radar echo extrapolation prediction accuracy by incorporating the dynamic field structure constraint into the extrapolation model and applying the deformation convolution. BRIEF DESCRIPTION OF DRAWINGS

[0040] Figure 1 is a flowchart of the present application;

[0041] Figure 2 is a schematic diagram of the three-dimensional radar echo extrapolation model based on the deformation convolution dynamic field constraint of the present application;

[0042] Figure 3 is a region internal velocity fading flowchart of the single-elevation radial velocity map of the present application;

[0043] Figure 4 is an inter-region velocity fading flowchart of the single-elevation radial velocity map of the present application;

[0044] Figure 5 A schematic diagram of a template defined region in the polar coordinate system of the present application;

[0045] Figure 6 A schematic diagram of a template defined region in the radial rectangular coordinate system of the present application;

[0046] Figure 7 A schematic diagram of a three-dimensional echo extrapolation model under the constraint of the present application MARC;

[0047] Figure 8 A schematic diagram of a gas flow field velocity feature identification template of the present application. DETAILED DESCRIPTION

[0048] The present application will be further described in detail by specific embodiments, the following examples are only descriptive, not limiting, and cannot limit the protection scope of the present application.

[0049] A three-dimensional radar echo extrapolation method based on deformation convolution dynamic field constraint, the innovation of which is that the steps of the method are:

[0050] (1) Collect S-band Doppler weather radar data from May 2015 to October 2016 in units of hours, involving 12 radar sites at 15 different locations in North China and East China. In VCP21 body scanning mode, the radar will complete conical scanning at 9 preset elevation angles in 6 minutes, so as to continuously collect radar reflectivity and radial velocity data within 2 hours.

[0051] (2) Convert the radar reflectivity data at each radar scanning time from the polar coordinate system to the three-dimensional data in the Cartesian rectangular coordinate system, perform velocity de-aliasing on the radial velocity data at each elevation angle, and map the radial velocity data to the radial rectangular coordinate system. The radar image data of 20 consecutive time points is defined as a sample. Among them, the first 10 time points cover 1 hour (6minx10), and the corresponding radar image data includes both reflectivity and radial velocity, which is used as the input of the radar echo extrapolation model, aiming to predict the next 10 consecutive radar images within 1 hour. The last 10 radar images in the sample are only composed of reflectivity, which is used as the label data for supervised model training and performance evaluation. This embodiment contains 9540 samples, and there is no time overlap between the samples, which together constitute a sample set.

[0052] The radial velocity data at each elevation angle is de-aliased by using the velocity de-aliasing algorithm MVDR based on the minimization of the velocity difference between regions. Velocity aliasing refers to the situation where the true radial velocity exceeds the non-aliasing interval [-V N ,V NWhen the velocity value is measured, the velocity value will be "folded" and deviated, which is shown as a false jump of positive and negative maximum velocity value in the radial velocity map. The velocity ambiguity resolution algorithm MVDR includes two parts: intra-region velocity ambiguity resolution and inter-region velocity ambiguity resolution.

[0053] The intra-region velocity ambiguity resolution first calculates the velocity difference between each point and its adjacent points in the connected region to identify the possible velocity ambiguity points. Based on these points, the region is further divided into multiple sub-regions. Next, according to the principle of minimizing the velocity difference between sub-regions, the ambiguity resolution is performed. Let the number of isolated regions be n r , the intra-region velocity ambiguity resolution process of the monopulse radial velocity map is shown in Figure 3 .

[0054] The inter-region velocity ambiguity resolution is a method to handle the ambiguity of velocity values in isolated velocity regions far from the radar center. The specific processing steps are as follows: First, start from the region close to the radar and gradually advance the ambiguity resolution operation to the region far from the radar. Second, give greater weight to regions close to the radar and with higher reliability to ensure the accuracy of the reference velocity. Let the number of regions to be executed be n a , these regions are sorted in ascending order of distance from the radar, and the inter-region velocity ambiguity resolution process of the monopulse radial velocity map is shown in Figure 4 .

[0055] After the velocity ambiguity resolution, nine high-resolution radial velocity information at each elevation angle in the polar coordinate system is obtained. On this basis, the radial velocity data is mapped to the radial rectangular coordinate system, thereby overcoming the problems of polar coordinate representation being not conducive to computer processing and the loss of azimuth information in Cartesian coordinate system. The specific definition of this coordinate system is as follows: In the plane rectangular coordinate system, the horizontal axis is positive to the right, representing the angle information in the polar coordinate system, where 0° at the origin represents the north direction, and the angle resolution is 1°; the vertical axis is positive downward, representing the radial distance of each point to the radar, and the distance resolution is 1 kilometer. In order to prevent the convergence zone near the north direction from being cut off, the image is expanded by 20°.

[0056] (3) When the MARC automatic recognition algorithm extracts the MARC feature information in the radial velocity map, there are mainly two steps. First, identify the convergence points in the radial velocity map at each elevation angle, and calculate their convergence strength and other characteristic parameters through template calculation to determine the strong convergence region at each elevation angle. Then, perform longitudinal matching on the strong convergence region to identify the strong convergence region that exists simultaneously at two or more elevation angles as a strong convergence field, and give its characteristic vector to obtain its three-dimensional structure and position information. The specific steps are as follows:

[0057] 1. Convergence point identification and convergence strength calculation

[0058] In the reflectivity map under the radial rectangular coordinate system, the area with reflectivity value higher than 35 dBZ is extracted as a mask, and only the foreground area in the mask is reserved as the target research area in the subsequent process.

[0059] For the processed single elevation radial velocity map, the boundaries of all positive velocity areas (v > 0) and negative velocity areas (v≤0) are extracted, the adjacent boundary points of the positive and negative velocity areas are determined, and the boundary points that meet the condition that the upper and lower velocity values change from positive to negative are identified as the convergence points. For each amplitude convergence point p, the feature vector is calculated in its related area Ω, where each parameter is respectively the maximum velocity difference dv max , the average velocity difference dv mean , and the coordinates of the positive and negative velocity extreme points and The calculation process of each parameter is as follows:

[0060]

[0061] In the formula, n + , and represent the number of positive velocity pixel points, the total sum of positive velocity, and the maximum positive velocity, respectively, n - , and represent the number of negative velocity pixel points, the total sum of negative velocity, and the maximum negative velocity, respectively. Ω represents the positive and negative velocity areas that meet the specific spatial and scale constraints with point p as the center. Here, a template with point p as the center is used to select the related area. Specifically, the spatial constraint requires that the deviation of the line connecting the maximum positive velocity point and the maximum negative velocity point from the radial direction should not be too large; the scale constraint specifies that the radial distance between the positive and negative velocity points should not be too far. Based on this, the template is designed as two isosceles triangles with a common vertex, and the area covered by the template is Ω. The template limited area in the polar coordinate system is shown in Figure 5 .

[0062] In the polar coordinate system, the size of the two isosceles triangle templates is fixed, but it needs to be rotated with the change of the azimuth angle to keep consistent with the radial direction, increasing the calculation complexity. In the radial rectangular coordinate system, the template can be directly aligned with the radial direction, but coordinate system mapping is required. When the radial distance r of the template center point from the radar increases, the central angle corresponding to the bottom edge of the template becomes smaller. The template limited area in the radial rectangular coordinate system is shown in Figure 6 . If the coordinates of the point to be measured p are (θ, r), then the arc length covered at the position extending a distance h from point p in the radial direction can be obtained by the arc length calculation formula, which is approximately equal to (4 / 7)|h|×2, and the corresponding central angle range can also be determined as [θ-k h , θ+kh ], where:

[0063]

[0064] 2. Single elevation strong convergence region positioning

[0065] When the characteristic parameter of a convergence point satisfies dv max ≥ 25 m / s and dv mean ≥ 12 m / s at the same time, it is recorded as a "strong convergence point". The judgment of the mean velocity difference aims to prevent misleading high values caused by individual noise points, and the threshold of 12 m / s is derived from the statistical analysis results of the sample. For each strong convergence region, the size of its minimum bounding rectangle ω, the center point coordinates, and the characteristic vector of the convergence point are recorded. Finally, the representative point p0 of each strong convergence region is determined by traversing and comparing the characteristic vectors of all points in the strong convergence region.

[0066] 3. Longitudinal matching of strong convergence regions

[0067] The method of searching from low elevation to high elevation layer by layer is adopted to match the strong convergence regions between adjacent elevation layers. When the center distance of the bounding rectangles of the convergence regions in adjacent high and low elevation layers is less than 15 kilometers, it is considered that the two convergence regions are related. In this process, there are two special cases that need to be explained:

[0068] First, if there is a "one-to-many" or "many-to-one" convergence region association, this study selects the one with the largest dv max to retain, and the others are excluded.

[0069] Second, for strong convergence regions that cannot find matching objects in adjacent elevations, especially when these unmatched convergence regions appear in higher elevation layers, they may represent newly emerging strong convergence regions, so they should be retained and continue to try to match with convergence regions in higher elevation layers.

[0070] Once a strong convergence region can obtain at least two associated matches in the vertical direction at different elevations, it can be judged as a significant MARC.

[0071] (4) Building a radar echo extrapolation model that takes three-dimensional radar reflectivity and radial velocity images within one hour as input and the reflectivity image of the next hour as output. This study adopts a dual-encoder input mechanism, which processes data according to the physical properties and semantic differences of the input data using the "independent encoding and fusion" strategy. Specifically, the reflectivity data and MARC structure information are processed through their respective encoders to obtain feature mappings, and the decoder integrates these two types of feature information to generate the extrapolated echo sequence in the future.

[0072] The model core construction module is ConvLSTM and CNN. In view of the characteristics of three-dimensional radar data, the ConvLSTM unit is adaptively improved, that is, the two-dimensional convolution is replaced by three-dimensional convolution to form a 3DConvLSTM unit, which significantly improves the ability of the model to process three-dimensional data. Finally, through the joint use of 3DConvLSTM and CNN units, the effective extrapolation of three-dimensional radar echo data is realized. The structure of the three-dimensional echo extrapolation model under the MARC constraint is as shown in Figure 7

[0073] The reflectivity encoder receives three-dimensional echo data at each time point as input, and the CNN module maps the original image data to the feature space, and then fuses the echo features from different height levels to express the redistribution of echo intensity caused by spatiotemporal atmospheric motion. The reflectivity encoder includes three three-dimensional convolution layers, wherein the first convolution layer keeps the spatial dimension of the input image unchanged, promoting effective conversion from image space to feature space; the subsequent two convolution layers realize deep fusion of cross-level features through specific parameter configuration. The reflectivity encoder parameter details are shown in Table 1.

[0074] Table 1 Reflectivity encoder parameter information table

[0075]

[0076] The dynamic field encoder processes the radial velocity data at the same time point, extracts the dynamic properties of the convection system and the dynamic behavior of the cloud cluster, including the three-dimensional structure and position information of the MARC features. The encoder is composed of a combination of 3DConvLSTM modules and CNN modules, which deeply calculate spatiotemporal sequence data to express potential data patterns and dependencies, thereby providing accurate constraint conditions for radar echoes. The dynamic field encoder parameter details are shown in Table 2.

[0077] Table 2 Dynamic field encoder parameter information table

[0078]

[0079] The decoder integrates the hidden representations output by the reflectivity encoder and the dynamic field encoder, and through a recurrent decoding operation, reconstructs the morphology, motion and evolution law of the echo. The decoder extracts echo features at future time points from the hidden state output by the 3DConvLSTM, and uses deconvolution technology to reversely map these features from the abstract feature space back to the original image space, completing the upsampling processing of the data, and finally generating the predicted radar echo sequence. The decoder parameter details are shown in Table 3.

[0080] Table 3 Decoder parameter information table

[0081] ​

[0082] (5) This embodiment adopts a deformed convolution method to improve the radar echo extrapolation model. The convolution kernel form is adjusted according to actual needs, the calculation process is simplified, and the extraction and representation of the effective velocity feature and its associated information are further expanded, so as to realize the rapid and automatic identification of MARC. The radar echo extrapolation model improved based on deformed convolution is as shown in Figure 2 . Specifically, the deformed convolution mainly has the following two methods:

[0083] 1. Based on the consideration of the maximum statistical distance between positive and negative velocity points, the height K H of the convolution kernel in each row is fixedly set to 15, and the width K W of the convolution kernel is dynamically calculated and determined according to the specific radial position thereof in the image by the following formula.

[0084]

[0085] In the formula, r represents the distance of the center of the convolution kernel from the top of the radial velocity image.

[0086] In order to realize the quantification of the feature parameters including the maximum velocity difference dv max and the average velocity difference dv mean , the weight distribution in the convolution kernel is determined by the following formula.

[0087]

[0088] In the formula, i represents the distance of the pixel point from the top of the template, j represents the distance of the pixel point from the left side of the template, (i,j) is the coordinate position of the pixel point in the template, Ω is the relevant area for feature calculation, and the definition process thereof is described in step (3). K l (i,j) is the weight size at the (i,j) coordinate.

[0089] Through the above weight setting, the deformed convolution can perform convolution operation along the specified direction, process the input radial velocity data, and finally generate the corresponding feature vector.

[0090] 2. In the complex airflow environment of the convection system, in addition to the airflow convergence structure, the airflow divergence, rotation and other structural features presented by the relative position of the maximum inflow and outflow velocity centers are important indication information for the potential occurrence of severe convective weather. The present application adopts a template matching method to implement zoned velocity matching and calculation, limits the distance between the velocity extreme points, and evaluates the relative position between the maximum inflow velocity center and the outflow velocity center, so as to accurately identify the key velocity features in the airflow field. The airflow field velocity feature identification template is as shown in Figure 8 .

[0091] The feature recognition template is a regular hexagon, which is specially used for identifying four types of airflow field structures, i.e. convergence, divergence, cyclone and anticyclone. In the template, 1) isosceles triangles a and b with the same size and a common vertex are responsible for capturing the two types of features, i.e. convergence and divergence, the height of which is set according to the maximum statistical distance between positive and negative velocity points in the convergence and divergence field, and the length of the bottom edge is determined according to the maximum allowable value of Δθ; 2) rhombuses c and d with the same size are used to capture the two types of features, i.e. cyclone and anticyclone, and the spatial constraint condition is that the offset angle of the line connecting the maximum positive velocity point and the maximum negative velocity point with the radial direction needs to be large enough, and the radial distance between the positive / negative velocity center points should not be too large. The above triangle group and rhombus group form a complementary regular hexagonal region in spatial distribution, and the center point of the region is denoted as p, and the convergence, divergence, cyclone and anticyclone properties of the point are obtained from the a, b region data and the c, d region data in the template.

[0092] Based on the consideration of statistical distance, the deformation convolution for determining the rotation type is performed on each row of the convolution kernel with a height K H Fixedly set to 11, the width K of the convolution kernel W Then, the height of the convolution kernel is dynamically calculated according to the specific position of the feature in the image by the following formula.

[0093]

[0094] In the formula, h∈{-5,-4,-3,…,3,4,5}.

[0095] (6) Dividing the sample set into a training set, a validation set and a test set, training the radar echo extrapolation model on the training set, performing hyperparameter tuning on the validation set, and performing model testing on the test set.

[0096] Each radar completes a volume scan every 6 minutes on average, and 10 three-dimensional radar reflectivity and radial velocity images in succession within one hour are used as the input of the radar echo extrapolation model, and then 10 reflectivity images in succession within the next hour are used as the label, a total of 9540 samples, and the samples are randomly allocated to the training set, the validation set and the test set according to a ratio of 4:1:1. Finally, 6360 samples are obtained in the training set, and 1590 samples are obtained in the validation set and the test set.

[0097] The radar echo extrapolation model is trained, and the weighted mean square error WMSE is used as the model loss function to guide the model to pay more attention to the prediction accuracy of the strong echo area. The WMSE is defined as follows:

[0098]

[0099] In the formula, I h,x,y represents the reflectivity size at the coordinates (h, x, y) in the label image, Weight represents the reflectance at (h,x,y) in the predicted image. h,x,y The weight value representing this location can be set differently based on the echo intensity, allowing the model to focus more on areas with strong echoes. In this embodiment, Weight... h,x,y The settings are as follows:

[0100]

[0101] (7) Using radar reflectivity data and radial velocity data within any hour, the radar reflectivity data at each radar scan moment is converted from polar coordinate system to three-dimensional data in Cartesian rectangular coordinate system. The radial velocity data at each elevation angle is deblurred and mapped to radial rectangular coordinate system. The above two types of data within that hour are used as input to the model to obtain the reflectivity image for the following hour.

[0102] This embodiment uses the Critical Success Index (CSI), Heidke Skill Score (HSS), and Fair Skill Score (ETS) to measure the consistency between the estimated proportion of reflectance exceeding a certain reflectance threshold and the actual proportion of reflectance exceeding that threshold. Higher values ​​for CSI, ETS, and HSS indicate better model performance. The definitions of CSI, HSS, and ETS are as follows:

[0103]

[0104] The formula for calculating REF is as follows:

[0105]

[0106] TP, FN, FP, and TN are defined by the following confusion matrix:

[0107] Table 4 Confusion Matrix

[0108]

[0109] The experimental results and analysis of this embodiment include three parts. The first two parts are the effectiveness verification of the model architecture, including the determination of different deformable convolution strategies and the effectiveness of incorporating dynamic field convergence features. The third part is the comparative experiment of the model of this embodiment with other typical models.

[0110] 1. Comparative testing of different deformable convolution strategies

[0111] Significant MARC requires the simultaneous fulfillment of the maximum velocity difference dv. max ≥25m / s and average velocity difference dv mean≥ 12 m / s, where the former is referred to as condition 1 and the latter is referred to as condition 2. Among them, condition 1 is the primary condition and is also an important criterion for determining whether the MARC reaches a "significant level"; condition 2 is a supplementary condition for excluding false "significant MARC" caused by noise. Statistical analysis of all actual significant MARC results in the modeling data shows that the proportion of MARC that can simultaneously meet the "significant MARC" standard (condition 1) among the MARC that can meet condition 2 is more than 95%. Therefore, the embodiment designs a comparative test scheme of "different deformation convolution strategies", that is, based on Figure 2 The following four models are generated in sequence and the effect comparison test of the four models is organized:

[0112] Model 1, let Figure 2 The "deformation convolution" in the above formula adopts a convolution kernel calculation method based on the maximum statistical distance between positive and negative velocity points, and the convergence-divergence type dynamic field three-dimensional distribution is obtained by the encoder using condition 1 for identifying MARC. Other modules of Figure 2 are trained and generated by the training sample set and the verification sample set.

[0113] Model 2, only replace the "condition 1" for identifying MARC in model 1 with "condition 2", and generate it using the same sample set.

[0114] Model 3, jointly use condition 1 of model 1 and condition 2 of model 2, fuse the two convergence-divergence type dynamic field three-dimensional distributions, and generate it using the same sample set with other modules of Figure 2 .

[0115] Model 4, select the "MARC condition" used in the model with the best test quality among model 1 to model 3 to obtain the convergence-divergence dynamic field three-dimensional distribution by the encoder; let Figure 2 The "deformation convolution" in the above formula adopts a convolution kernel calculation method based on the maximum statistical distance between positive and negative velocity points, and the rotation type dynamic field three-dimensional distribution is obtained by the encoder; then fuse the above two types of three-dimensional dynamic distribution fields, and generate it using the same sample set with other modules of Figure 2 .

[0116] Using the test sample set not participating in the model training to organize the test of the above four models, calculating the average of CSI, HSS and ETS three indicators based on the prediction results given by the model extrapolation for 1 hour, and obtaining the evaluation results based on all test samples as shown in Table 5.

[0117] ①Comparing model 2 and model 4 under each threshold, it can be seen that in the problem of establishing a strong convective reflectivity echo extrapolation model, it is not the more the better to integrate the type of dynamic field information, and the integration of convergence-divergence type dynamic field information has the largest improvement on the radar echo extrapolation effect;

[0118] The model 2 performs best. This result can be explained. In the process of obtaining radial velocity data, various errors are inevitable due to various interferences and limitations of the equipment itself. Even if the velocity ambiguity work has been performed in the embodiment, it is still impossible to completely avoid misleading high values caused by individual noise points. The maximum velocity difference is susceptible to noise, resulting in a less effective extrapolation than the extrapolation based on the average velocity difference.

[0119] In particular, the model 2 that wins under the different deformation convolution strategies is recorded as 3DConvLSTM-Dconv.

[0120] Table 5: Test scores obtained by models under four deformation convolution strategies

[0121]

[0122] 2. Effectiveness of the integrated convergence feature and comparison test of different integration schemes

[0123] To verify the effectiveness of the three-dimensional radar echo extrapolation model under the constraint of the dynamic field proposed in the embodiment, the test sample set not participating in the model training is used to test the 3DConvLSTM-MARC, 3DConvLSTM-Dconv and 3DConvLSTM models. Figure 7 The test results are shown in Table 6.

[0124] Table 6: Score results of three extrapolation models with and without integrated dynamic field information

[0125]

[0126] Under the three comprehensive evaluation indexes of CSI, HSS and ETS, the 3DConvLSTM model without integrated dynamic field information has the lowest score. After introducing the MARC distribution information, the CSI score is improved by 0.99% to 1.53%, and the HSS score is improved by 0.77% to 1.88%, of which the improvement rates under the 40dBZ threshold reach 5.8% and 4.5%, respectively. The improved extrapolation model 3DConvLSTM-Dconv based on deformation convolution further improves the 3DConvLSTM-MARC model. Compared with the 3DConvLSTM model, the CSI score is improved by 1.42% to 1.98%, and the HSS score is improved by 1.2% to 2.43%, of which the improvement rates under the 40dBZ threshold are 7.4% and 5.8%, respectively.

[0127] It can be seen that the power field structure constraint based on radial velocity data is integrated into the radar echo extrapolation model, which can significantly improve the accuracy of the model in the strong reflectivity area. Compared with the 3DConvLSTM-MARC model, the 3DConvLSTM-Dconv extrapolation model can more comprehensively capture the shape and change rule of the convective system, especially when the prediction time is longer, the extrapolation results provided by the model show higher accuracy.

[0128] 3. Comparison test with typical extrapolation model based on reflectivity data only

[0129] In order to verify the effectiveness of the extrapolation model in the embodiment, other typical extrapolation models based on reflectivity data only are selected for comparison test, and the comparison models include constant extrapolation Persistence, optical flow and ConvLSTM. The CSI, HSS and ETS scores of the above three comparison models and the model in the embodiment under different reflectivity thresholds for 2km height reflectivity echo are shown in Table 7.

[0130] Table 7 Score results of different models

[0131]

[0132] Under various reflectivity thresholds, the CSI, HSS and ETS scores of the model in the embodiment are the highest. Taking the CSI score as an example, compared with the traditional optical flow algorithm, the result of the deep learning method ConvLSTM is improved by 2.22%~4.57% compared with it, and the 3DConvLSTM-Dconv model in the embodiment is further improved by 4.74%~6.47% compared with the ConvLSTM method, especially in the strong reflectivity (40dBZ) area, the 3DConvLSTM-Dconv model improves the CSI by 26%, 56% and 179% compared with ConvLSTM, optical flow and constant extrapolation, respectively.

[0133] In general, the scores of the extrapolation model based on deep learning are better than those of the traditional extrapolation algorithm, which embodies the advantage of the model in the nonlinear process modeling. The 3DConvLSTM-Dconv model further improves the prediction accuracy on the basis of the ConvLSTM method, which proves that the modeling scheme of increasing the data dimension from 2D to 3D and then integrating the power field constraint based on radial velocity data is effective and feasible, and the improvement of the model quality is significant.

[0134] Although the embodiments of the present application and the drawings are disclosed, those skilled in the art can understand that various substitutions, changes and modifications are possible without departing from the spirit and scope of the present application and the appended claims, therefore, the scope of the present application is not limited to the disclosed content of the embodiments and the drawings.​

Claims

1. A three-dimensional radar echo extrapolation method based on deformation convolution dynamic field constraint, characterized in that: Steps of the method are: S1, the radar operates in VCP21 volume scan mode, that is, 9 scans of different elevation angles are completed every 6 minutes, and radar reflectivity and radial velocity data within 2 hours are continuously collected; S2, for the data of each radar volume scan time: convert the radar reflectivity data from polar coordinates into three-dimensional data in the Cartesian rectangular coordinate system; perform de-aliasing processing on the radial velocity data, and map the radial velocity data into the radial rectangular coordinate system; use the reflectivity and radial velocity images within the first hour as the input of the extrapolation model, and use the reflectivity image within the second hour as the label to construct a sample set; S3, use the middle layer radial convergence automatic identification algorithm to obtain the middle layer radial convergence feature information in the radial velocity image; S4, build a radar echo extrapolation model with three-dimensional radar reflectivity and radial velocity images within one hour as input and reflectivity images within the next 1 hour as output; S5, introduce deformation convolution in the radar echo extrapolation model to automatically extract the dynamic field information of the convection system; The deformation convolution refers to adjusting the shape of the convolution kernel according to actual needs, further expanding the extraction and representation of effective velocity features and their associated information, and the two ways of deformation convolution are: (1) Based on the maximum statistical distance between positive and negative velocity points, the deformation convolution is in the convolution kernel height K of each row H Fixed to 15, the width K of the convolution kernel W Then, according to its specific radial position in the image, it is dynamically calculated and determined by the following formula; In the formula, r represents the distance between the center of the convolution kernel and the top of the radial velocity image; To realize the quantification of the feature parameters including the maximum speed difference dv max and the average speed difference dv mean The weight distribution within the convolution kernel is determined by the following formula; In the formula, i represents the distance of a pixel point from the top of the template, j represents the distance of a pixel point from the left side of the template, (i,j) is the coordinate position of the pixel point in the template, Ω is the relevant region for feature calculation, K l (i,j) is the weight size at the (i,j) coordinate. (2) Use a template matching method to implement partitioned velocity matching and calculation, limit the distance between velocity extreme points, and evaluate the relative position between the center of the maximum inflow velocity and the center of the outflow velocity to accurately identify key velocity features in the airflow field. The airflow field velocity feature identification template is a regular hexagon, which is specially used to identify 4 types of airflow field structures: convergence, divergence, cyclone and anticyclone. Based on the statistical distance of the consideration, for the determination of the rotation type of the deformation convolution in each row of the convolution kernel height K H Fixed set to 11, the width of the convolution kernel K w Then according to its specific position in the image is dynamically calculated by the following formula; In the formula, h∈{-5,-4,-3,…,3,4,5}; S6, divide the sample set into a training set, a validation set and a test set, train the radar echo extrapolation model on the training set, optimize the hyperparameters on the validation set, and test the radar echo extrapolation model on the test set; S7, use radar reflectivity data within any hour, convert radar reflectivity data at each radar scan time from polar coordinates to three-dimensional data in the Cartesian rectangular coordinate system using the steps of S2, perform de-aliasing processing on the radial velocity data at each elevation angle, map the radial velocity data to the radial rectangular coordinate system, and use the above two types of data within the hour as the input of the model, that is, the reflectivity image within the next 1 hour can be obtained.

2. The method of claim 1, wherein: In the formula, r represents the distance between the center of the convolution kernel and the top of the radial velocity image; The S2 converts the radar reflectivity data at each radar scan time from polar coordinates to three-dimensional data in the Cartesian rectangular coordinate system by first converting the radar base data from polar coordinates to Cartesian coordinates using the bilinear interpolation method in the horizontal dimension, and then interpolating in the vertical dimension to obtain three-dimensional equidistant grid field data; The S2 performs de-aliasing processing on the radial velocity data at each elevation angle refers to using a de-velocity aliasing algorithm based on minimizing the velocity difference between regions to solve the velocity aliasing problem, and mapping the de-aliasing processed radial velocity data to the radial rectangular coordinate system.

3. The method of claim 1, wherein: The S4 takes an "independent coding and fusion" strategy to process three-dimensional radar echo images of ten continuous time points in the past one hour, and the reflectivity data and the mid-level radial convergence structure information are respectively obtained by the respective encoder to obtain the feature mapping, the decoder integrates the two types of feature information, and decodes to generate the extrapolated echo sequence in the future, and the radar echo extrapolation model is completed by jointly using the ConvLSTM network based on 3D convolution and the convolutional neural network CNN to effectively extrapolate the three-dimensional radar echo data; The ConvLSTM network based on 3D convolution is to replace the two-dimensional convolution operation of the ConvLSTM unit with a three-dimensional convolution operation to process the radar three-dimensional data corresponding to each radar scanning time, and to use the specific cycle structure of the ConvLSTM unit to extract and synthesize the radar three-dimensional data at different times within one hour; Each ConvLSTM unit contains three key components: a forget gate, an input gate, and an output gate. This unit receives the data input X at the current time point. t The state output H at the previous time point t-1 and long-term memory unit C t-1 As input, and to produce an updated long-term memory unit C t and the current state output H t As output, the operations performed by the ConvLSTM unit and the calculation of its output results are shown below: i t = σ(W xi * X t + W hi * H t-1 + W ci * C t-1 + b i ) f t = σ(W xf * X t + W hf * H t-1 + W cf * C t-1 + b f ) o t = σ(W xo * X t + W ho * H t-1 + W co * C t-1 + b o ) where * represents convolution operation, represents Hadamard product, σ is an activation function, W and b represent the corresponding convolution kernel weight and bias term respectively, tanh is hyperbolic tangent function, i t , f t and o t represent input gate, forget gate and output gate respectively.

4. The method of claim 1, wherein: The S6 trains the radar echo extrapolation model using a weighted mean squared error (Weighted Mean Squared Error, WMSE) as a model loss function to guide the model to pay more attention to the prediction accuracy of the strong echo area, and the WMSE is defined as follows: In the formula, I h，x，y represents the reflectivity size at the coordinate (h, x, y) in the label image, represents the reflectivity size at (h, x, y) in the predicted image, Weight h，x，y represents the weight value of the position, and different weight values can be set according to different echo intensities, so that the model pays more attention to the strong echo area.

Images