CN108254750B

CN108254750B - Down-blast intelligent identification early warning method based on radar data

Info

Publication number: CN108254750B
Application number: CN201711193420.6A
Authority: CN
Inventors: 卞浩瑄; 王兴; 周显信; 王介君; 朱文镇
Original assignee: Nanjing University of Information Science and Technology
Current assignee: Nanjing Walker Weize Technology Co ltd
Priority date: 2017-11-24
Filing date: 2017-11-24
Publication date: 2021-07-30
Anticipated expiration: 2037-11-24
Also published as: CN108254750A

Abstract

The invention discloses an intelligent downburst identification early warning method based on radar data. Firstly, an optical flow field of a vertical section of the reflectivity factor is analyzed by adopting an optical flow method, and the law of the storm core evolution along with time is obtained. And then, performing function fitting on the descending process of the storm core top height by adopting a Lagrange mechanical model. And then, matching and identifying the positive and negative speed pair images in the radial speed field of the middle layer of the storm core by adopting a histogram and a babbit coefficient statistical analysis method, and finally realizing intelligent early warning of downburst flow by integrating judgment of a series of thresholds. By the method, potential downburst can be identified in the early stage of storm core lifting and rapid sinking, and the sinking rate of the storm core and the time for the storm core to reach the ground can be estimated, so that intelligent identification and early warning of the downburst are realized.

Description

Down-blast intelligent identification early warning method based on radar data

Technical Field

The invention belongs to the field of atmospheric science, and relates to a method for recognizing a radial convergence and divergence image of a radial velocity field of a Doppler weather radar data to realize intelligent prediction and early warning of downburst flow.

Background

Downburst is a strong divergent straight-line airflow formed when a cold outflow formed in a strong convection storm reaches a near-ground layer, and wind direction shear and wind speed mutation caused by the downburst airflow form a serious threat to the flight safety of an aircraft, particularly the stable take-off and landing of the aircraft. Although the great disasters caused by downburst flow cause the formation and development of the disasters frequently occurring in groups and fatalities, the strict monitoring on the formation and development of the disasters has been highly emphasized by aviation and meteorological departments for a long time, and the monitoring, forecasting and early warning are very difficult due to the small spatial scale, the fast development speed of life and consumption and the short duration of life and consumption, so that the disasters are a worldwide problem for a long time.

The Doppler weather radar is used as an important means for monitoring and analyzing strong convection weather, and plays a significant role in monitoring and forecasting downburst. However, as downburst has the characteristics of small spatial dimension, high living and disappearing speed and the like, under the limitation of the existing detection means and technical conditions, a great technical bottleneck still exists in the intelligent identification and early warning of downburst.

Disclosure of Invention

The invention provides a downburst flow identification tracking and early warning algorithm based on Doppler weather radar data, which comprises the steps of firstly obtaining the law of the evolution of a storm core along with time by adopting an optical flow method, then carrying out function fitting on the descending process of the top height of the storm core, then carrying out matching identification on positive and negative speed pair images in a radial velocity field in a middle layer of the storm core by adopting a histogram and babbit coefficient statistical analysis method, and finally realizing the intelligent early warning of the downburst flow by synthesizing a series of threshold judgment.

The technical scheme adopted by the invention is as follows:

step 1: the following analysis is performed on the basic reflectivity factor in the live radar data: when the horizontal space continuous coverage range of the reflectivity value on one elevation surface exceeding the threshold TH1 is larger than the threshold TH2, the following steps are carried out; otherwise, no processing is performed. The calculation method of the step comprises the following steps:

wherein the content of the first and second substances,

representing radar elevation, theta representing body sweep azimuth, r representing range bin number,

and the echo intensity value of any position on a certain elevation surface under polar coordinates is represented. S represents the area of the region where the echo intensity value exceeds TH 2.

Step 2: the radar data at the current time and the previous Tx times are respectively analyzed as follows: the position where the maximum reflectance value occurs is first calculated, and then the vertical profile of the reflectance factor across that position is calculated and plotted as a vertical profile. These vertical profiles use consistent two-dimensional coordinates and scales.

And step 3: and (3) performing flow field analysis on the vertical section maps of all adjacent moments drawn in the step (2) by adopting an optical flow method to obtain the moving trend of the storm core (namely the large-value echo area in the image).

And 4, step 4: and 3, analyzing the optical flow field to only obtain the speed and the direction of the movement of the storm center in the past period of time, and in order to obtain the trend and the speed of the movement of the storm center in the future period of time, carrying out extrapolation prediction by using the calculated optical flow field. The extrapolation prediction needs to be based on a certain mechanical relationship, which can be linear or nonlinear, and the key to the quality of the prediction effect depends on whether the mechanical relationship accurately represents the current rule of storm digestion development. The method adopts a Lagrange mechanical model to fit the relationship of the storm core movement trend evolving along with time. If the predicted storm core rate of descent exceeds a particular threshold TH3, continuing the following steps; otherwise, no processing is performed.

And 5: since the rapid descent of the storm core is only a necessary, but not sufficient, condition for the occurrence of downbursts, it is necessary to further analyze several characteristics of the radar radial velocity field. And respectively calculating histograms of a positive speed area and a negative speed area in a range which takes the storm core position as the center and takes the threshold TH4 as the radius in the current radar radial speed map. In order to obtain the histogram more simply and effectively, the calculation process does not analyze an RGB color image drawn by a radial velocity field, but directly analyzes a positive velocity value and a negative velocity value in the radial velocity field, takes TH5 as the precision and 0-TH 6 as the value range of the abscissa of the histogram, and respectively accumulates and counts the frequency of different velocity values, thereby obtaining a positive velocity distribution (outflow) histogram and a negative velocity (inflow) distribution histogram.

Step 6: matching analysis is carried out on the two histograms by adopting a Babbitt coefficient, and the calculation method comprises the following steps:

wherein P is₁And P₂Which represent the percentage of any color value (i.e., velocity value) in the positive and negative velocity distribution histograms, respectively, to the total. Multiplying the percentage of each same abscissa i, squaring, and solving an accumulated value to obtain a calculated result, namely the similarity value (the Papanicolaou coefficient factor value) of the two images. The value can objectively and quantitatively reflect the matching degree of the positive and negative speed pairs, the value range is 0 to 1, and the larger the value is, the higher the similarity is. When the value is greater than the threshold TH7, it is determined that a downburst is about to occur as a high probability event.

And 7: and (4) calculating the descending trend and speed of the storm core and the time of arriving at the ground by combining the steps and the judgment of a series of characteristic thresholds, and finally realizing intelligent identification tracking and early warning of the downburst.

The values of the thresholds are different according to different regions and different seasons, and a group of reference value ranges are given, as shown in table 1:

TABLE 1 ranges of determination thresholds for downburst

Compared with the prior art, the invention has the beneficial effects that:

the method mainly adopts a series of methods such as optical flow analysis, linear extrapolation, image similarity matching and the like for a plurality of characteristics shown on radar echo images and radial velocity fields in the early and middle periods of occurrence of downburst, realizes intelligent identification tracking and early warning for the occurrence of potential downburst, and provides an improved method and implementation steps for key technical problems.

When the storm core sinks rapidly, the middle layer convergence process cannot present symmetrical positive and negative speed pairs on a radar radial velocity diagram, and great difficulty is brought to image recognition. In order to solve the problem, the invention provides a histogram and a statistical analysis method of the Babbitt coefficient of the positive and negative speed regions, and focuses on carrying out graphic image analysis on the converged speed magnitude and percentage in each direction so as to make accurate judgment.

Drawings

Fig. 1 is a main algorithm flow of the present invention.

FIG. 2 is a vertical cross-section of the evolution of the reflectivity factor kernel over time, with events 2015-06-0113: 01,2015-06-0113: 07,2015-06-0113: 13 from left to right, respectively.

Fig. 3(a) is a reflectance factor core optical flow field generated by the 3 image calculations shown in fig. 2.

Fig. 3(b) is a cross-sectional view (forecast) generated by a 12 minute extrapolation (i.e., 2015-06-0113: 25).

FIG. 3(c) is a 2015-06-0113: 25 reflectance factor kernel vertical cross-section (live).

FIG. 4(a) shows the radar radial velocity field 2015-06-0113: 13.

Fig. 4(b) shows the positive speed region searched.

Fig. 4(c) shows the negative velocity region retrieved.

Fig. 5(a) is a histogram of the positive velocity region.

Fig. 5(b) a histogram of the negative velocity region.

FIG. 6 is a graph of probability density distribution for different velocity values in positive and negative velocity regions.

Detailed Description

The invention is further described below with reference to the accompanying drawings.

In order to examine the practical effect of the method of the present invention, the doppler weather radar-based data of a region in the state of north Hubei province before and after 21 st day 1/6/2015 is selected for analysis in this embodiment. At 21 hours, about 32 minutes, a ship turn-over and sinking event of 'star of east' in China occurs, and the ship sails to the vicinity of a water channel of a great mars in Yangtze river, Guest county, Jingzhou city, encounters sudden downburst and strong wind and rainstorm brought by the sudden downburst, so that the ship turns over and sinks, and 442 people are in distress. Since monitoring of highly convective weather such as downburst is a continuous process, radar-based data of 2015-6-113: 13(UTC, the same below) is used as an example.

Step 1: reading the base data "Z _ RADR _ I _ Z9716_20150601131300_ O _ DOR _ SA _ cap. bin. bz 2", first processing the reflectivity factor at the 1 st elevation angle, and calculating a plurality of regions with reflectivity values greater than TH1(TH1 ═ 50dBZ), wherein the region with the largest area is located at [112.96 °, 29.86 °]Nearby, about 122km²Greater than TH2(TH 2-100 km)²). When the determination condition of the above expression (1) is satisfied, the following steps are continued.

Step 2: and drawing a vertical section view which is defined as VCS _1313 after crossing the area with the maximum reflectivity factor, and then calculating and drawing reflectivity factor vertical section views of radar data of 13:01 and 13:07 on the same day by using the same section angle, position and length, wherein the reflectivity factor vertical section views are respectively defined as VCS _1301 and VCS _ 1307. The radar data of the last 2 moments, i.e., Tx 2, is selected. Vertical cross-sections using identical two-dimensional coordinates and scales, in km, are plotted from the 3 time radar data as shown in fig. 2.

And step 3: the calculation of the optical flow field is performed for VCS _1301 and VCS _1307, VCS _1307 and VCS _1313, respectively, using an optical flow method. Since the optical flow method is applied on the premise that the image object has a small displacement amount and the displacement vector is globally or locally smooth, and the echo image is typically a non-rigid body and does not conform to the premise assumption of the optical flow method, the following improvement is made in the calculation.

(1) Filtering out areas with the reflectivity factors smaller than 45dBZ in VCS _1307 and VCS _1307, only calculating an optical flow field in a core area of the reflectivity factors, and iterating by using 4 layers of pyramid layering in the calculation process to obtain an optical flow vector field which is marked as Vo_p1。

(2) Will Vo_p1Each velocity vector is decomposed into 2 scalars parallel to the x-axis and the y-axis, the median values in the 2 directions are respectively calculated, and then the 2 scalars are combined into 1 velocity vector V_mThis vector characterizes to some extent the general trend of the development of the core shift in the vertical profile.

(3) Using V_mFor the overall position of the echo image in VCS _1307Move with a displacement of 1V_m. Then, the image is subjected to optical flow calculation with VCS _1307 to obtain an optical flow vector field, which is denoted as Vo_p2. Then to Vo_p2And V_mMatrix addition calculation is carried out to obtain an optical flow field which is marked as V_1-2。

(4) The same algorithm is adopted to calculate the optical flow fields of VCS _1307 and VCS _13013, which are marked as V_2-3。

(5) The weighted value of the optical flow field is calculated using equation (9):

V＝a·V_1-2+b·V_2-3 (3)

in this embodiment, a takes the value of 0.375, and b takes the value of 0.625. Thus, an optical flow velocity vector field V is calculated and plotted as an image as shown in fig. 3 (a).

It should be noted that, in the calculation process of the present embodiment, the cross-sectional reflectance factor is used as an input instead of directly using the cross-sectional view in the RGB format as the input data for the optical flow calculation. Since the optofluidic method usually deals with a 0,255 gray-scale image and the reflectance value does not exceed this range, no further conversion of the reflectance value is necessary.

And 4, step 4: and calculating the descending speed of the vertical section core ceiling. This embodiment extracts the velocity value V at the position as in (right) coordinates (35.5,8.2) of fig. 2 from the optical flow velocity vector field V_top，V_topThis value is greater than TH3(TH3 is 0.7) at 0.77(km/6min), indicating that the estimated storm core rate of descent exceeds a certain threshold, and the following steps are continued.

And 5: the radar radial velocity data of 13:13 is read and the velocity field image is shown in fig. 4 (a). With [112.96 °, 29.86 ° ] (storm core position) in step 1 as the center and TH4(TH4 ═ 28km) as the radius, a positive velocity region and a negative velocity region are extracted, respectively, and images drawn by these regions are shown in fig. 4(b) and 4 (c). Then, the statistical distribution of the velocity values (color values of the image) corresponding to each pixel in the positive and negative velocity regions is calculated, and the histograms are drawn as shown in fig. 5(a) and 5 (b). Fig. 5(a) is an outflow histogram, which is generated by calculating each pixel value of the RED channel in fig. 4 (b); fig. 5(b) is an inflow histogram calculated from each pixel value of the GREEN channel in fig. 4 (c). The abscissa is [0,255] and the ordinate is the cumulative sum of the different velocity values (color values of the image) which are directly related to the calculated image size. The more valuable information here is the proportion of the different velocity values to the image, i.e. their probability distribution.

Thus, the histograms (i.e., the mass distribution maps) represented in fig. 5(a) and 5(b) are further converted into probability density distributions of different velocity values in the positive and negative velocity regions, respectively, as shown in fig. 6. In fig. 6, the abscissa represents the radial velocity value, wherein the negative velocity is in units of "-m/s". The probability density distribution data is substituted into equation (8), and the barbituric coefficient ρ of the positive and negative velocity histograms is calculated to be 0.89, which is greater than the threshold TH7(TH7 is 0.75). As such, it may be determined that a downburst is about to occur as a probable event.

Step 6: in order to estimate the time for the downburst to affect the near-surface, a linear extrapolation is performed using the optical flow field calculated in the above step. FIG. 3(b) is a vertical section using a Lagrangian mechanics model, based on graph VCS _13013, and an optical flow velocity vector field V as a motion vector, extrapolated for the next 12 minutes (i.e., 13: 25). In order to reduce the common divergence phenomenon of the extrapolated image, the interpolation filling of blank singular points is also carried out on the image subjected to the optical flow extrapolation. FIG. 3(c) is a 13:26 radar live plot.

Comparing fig. 3(b) and fig. 3(c), it can be seen that the spatial distribution and the moving trend of the echo body in the radar vertical section are substantially consistent with those of the simultaneous scene by using the 12-minute extrapolation prediction made by the above method. However, the echo above 15dBZ in fig. 3(b) is overall weaker than live, and the echo kernel (greater than 50dBZ) sinks at a slower top height than live. This is mainly because the linear mechanical extrapolation model adopted is difficult to fit the actual trend of storm movement, and meanwhile, the process of the storm generation and elimination development is nonlinear and extremely complex.

According to the business program developed by the algorithm and the calculation process, the radar base data information is read and analyzed in real time, and intelligent identification tracking and early warning can be realized on the occurrence of the potential downburst in the early stage of storm core lifting and rapid sinking. From the calculation result of the embodiment, effective early warning can be made by 2-3 individual scanning times in advance.

Claims

1. An intelligent downburst flow identification early warning method based on radar data is characterized by comprising the following steps:

step 1: the following analysis is performed on the basic reflectivity factor in the live radar data: when the horizontal space continuous coverage range of the reflectivity value on one elevation surface exceeding the threshold TH1 is larger than the threshold TH2, the following steps are carried out; otherwise, no treatment is carried out; the calculation method of the step comprises the following steps:

wherein the content of the first and second substances,

representing the echo intensity value of any position on a certain elevation surface under polar coordinates; s represents the area of the region where the echo intensity value exceeds TH 2;

step 2: the radar data at the current time and the previous Tx times are respectively analyzed as follows: firstly, calculating the position where the maximum reflectivity value appears, then calculating the vertical section of the reflectivity factor passing through the position, and drawing a vertical section diagram; these vertical profiles use consistent two-dimensional coordinates and scales;

and step 3: performing flow field analysis on the vertical section maps at each adjacent moment drawn in the step 2 by adopting an optical flow method to obtain a storm core, namely the trend of movement of a large-value echo area in the image;

respectively calculating optical flow fields of VCS _1301, VCS _1307 and VCS _1313 by adopting an optical flow method; since the optical flow method is applied under the precondition that the image target has a small displacement amount and the displacement vector is globally or locally smooth, and the echo image is a typical non-rigid body and does not conform to the precondition assumption of the optical flow method, the following improvement is made during calculation:

(1) filtering out areas with the reflectivity factors smaller than 45dBZ in VCS _1307 and VCS _1307, only calculating an optical flow field in a core area of the reflectivity factors, and iterating by using 4 layers of pyramid layering in the calculation process to obtain an optical flow vector field marked as V_op1；

(2) Will V_op1Each velocity vector is decomposed into 2 scalars parallel to the x-axis and the y-axis, the median values in the 2 directions are respectively calculated, and then the 2 scalars are combined into 1 velocity vector V_mThe vector represents the general trend of the development of the core motion in the vertical profile to a certain extent;

(3) using V_mThe total position of the echo image in VCS _1307 is moved by 1V_m(ii) a Then, the image is subjected to optical flow calculation with VCS _1307 to obtain an optical flow vector field, which is denoted as V_op2Then, V is adjusted_op2And V_mMatrix addition calculation is carried out to obtain an optical flow field which is marked as V_1-2；

(4) The same algorithm is adopted to calculate the optical flow fields of VCS _1307 and VCS _13013, which are marked as V_2-3；

V＝a·V_1-2+b·V_2-3 (3)

and 4, step 4: fitting the relationship of the storm core movement trend evolving along with time by adopting a Lagrange mechanical model; if the predicted storm core rate of descent exceeds a particular threshold TH3, continuing the following steps; otherwise, no treatment is carried out;

and 5: since the rapid descent of the storm core is only a necessary condition for the occurrence of downburst, not a sufficient condition, it is necessary to further analyze several characteristics of the radar radial velocity field; respectively calculating histograms of a positive speed area and a negative speed area in a range which takes a storm core position as a center and a threshold TH4 as a radius in a current radar radial velocity diagram; in order to obtain the histogram more simply and effectively, the calculation process does not analyze an RGB color image drawn by a radial velocity field, but directly analyzes a positive velocity value and a negative velocity value in the radial velocity field, takes TH5 as the precision and 0-TH 6 as the value range of the abscissa of the histogram, and respectively accumulates and counts the frequency of different velocity values, so as to obtain positive velocity distribution, namely an outflow histogram and negative velocity, namely an inflow distribution histogram;

wherein P is₁And P₂Respectively representing the percentage of any color value, namely the speed value, in the positive speed distribution histogram and the negative speed distribution histogram to the total amount; multiplying the percentage of each same abscissa i, squaring, and solving an accumulated value to obtain a calculated result, namely a similarity value of the two images, namely a Papanicolaou coefficient factor value; the value objectively and quantitatively reflects the matching degree of the positive and negative speed pairs, the value range is 0 to 1, and the larger the value is, the higher the similarity is represented; when the value is greater than the threshold TH7, determining that the downburst is about to occur as a high probability event;