CN110988883A

CN110988883A - Intelligent squall line characteristic identification early warning method in radar echo image

Info

Publication number: CN110988883A
Application number: CN201911309889.0A
Authority: CN
Inventors: 卞浩瑄; 王兴; 陆冰鉴; 王丽娟; 苗春生; 陈苏婷; 朱文镇; 张海阳
Original assignee: Nanjing University of Information Science and Technology
Current assignee: Nanjing Walker Weize Technology Co ltd
Priority date: 2019-12-18
Filing date: 2019-12-18
Publication date: 2020-04-10
Anticipated expiration: 2039-12-18
Also published as: CN110988883B

Abstract

The invention provides an intelligent squall line characteristic identification early warning method in a radar echo image, which takes Doppler weather radar detection data as a main data source, analyzes the spatial distribution and the intensity of basic reflectivity detected by a radar, and realizes intelligent squall line characteristic identification early warning on the radar through a series of steps of numerical value preprocessing, filtering, image characteristic extraction, central axis analysis on a target object, squall line shape analysis and the like. The application of the method can automate and objectify the work of subjective analysis and interpretation of radar echo images by weather professionals in the past, and improve the accuracy and timeliness of the related services of squall line identification and strong convection weather warning.

Description

Intelligent squall line characteristic identification early warning method in radar echo image

The technical field is as follows:

the invention belongs to the field of geoscience, and relates to a method for intelligently identifying and early warning a squall line phenomenon represented in a radar echo image, in particular to a method for intelligently identifying and early warning squall line characteristics in a radar echo image.

Background art:

the squall line is an obvious strip-shaped or arc-like image appearing on a radar echo image when a weather radar detects, and the echo intensity of the shape area is obviously higher than that of the peripheral area outside the area. A meteorological "squall line" refers to a narrow, highly convective air band formed when meteorological elements, such as atmospheric pressure, temperature, or wind, change abruptly. Generally, it is a cloud and rain belt composed of a plurality of thunderstorm monomers or super monomers arranged approximately in a straight line or arc, the width of which is about tens to a hundred kilometers, the length of which is about tens to hundreds of kilometers, and the duration from the initial formation to the basic dissipation of which can reach hours to tens of hours.

The squall line formation, development process, and squall line size and spatial distribution may be clearly detected by a weather radar. In the forming stage, the radar echo image generally shows a plurality of relatively independent convection monomers, the outline of each monomer is clear, the echo intensity of the central area of each monomer exceeds 40dBZ, the horizontal scale of each monomer is about 15-20 kilometers, and the spatial distribution of each monomer is relatively dispersed. In the process of continuous development, the volumes of the convection monomers are gradually increased, the heights of the monomers are gradually raised to be more than 10 kilometers, and the convection monomers are gradually gathered and combined to form banded spatial distribution. In the vigorous stage, an obvious strip-shaped or long arc-shaped high-value echo area appears on a radar echo image, the echo intensity reaches 50 or even more than 60dbZ, and the echo top height rises to more than 15 kilometers. The outer edge echo intensity gradient of the high-value echo area is large and is represented as a strong super monomer or a strong convection monomer group. In this case, the squall line characteristics may be displayed on the radar echo image.

When the squall line occurs, strong convective weather phenomena such as steep rise of air pressure, sudden drop of temperature, sudden change of wind direction, strengthened wind power and the like often occur, and meanwhile, disastrous weather such as thunderstorms, lightning, hail or tornadoes can be accompanied. Although the squall line is not the direct cause of the disastrous weather, the occurrence of the squall line often predicts the occurrence of devastating winds, tornadoes, hail, and brief strong precipitation, and once it occurs, its destructive power is very great. For example, in 2009, 6, 3, 20 hours and 30 minutes, the Shanghai region in Henan is suddenly attacked by strong convection weather, the wind power reaches 11 grades at most, rainstorm and hail are accompanied, and the whole process lasts for about 50 minutes. According to statistics after disasters, 137 villages and towns in the whole city are in disaster at different degrees, more than 4800 houses collapse, 769.3 ten thousand trees are damaged, 309.84 ten thousand mu of lodging wheat is obtained, 20 people lose lives, and the direct economic loss is up to 15 hundred million yuan.

For experienced weather professionals, the squall line formation and development process can be accurately interpreted through subjective analysis of radar echo images. For a computer program, how to accurately identify squall line characteristics in a radar echo image and avoid missing judgment and misjudgment as much as possible is a complex algorithm problem. In fact, with the rapid development of computer technology, the continuous innovation of weather detection equipment and technology is driven, weather observation and detection data which can be obtained by weather service personnel in real time in the last 5-10 years are increased explosively, and the information content of the mass data which is pushed and updated in real time exceeds the capability of people for subjective analysis and interpretation in time. Therefore, it is imperative that a computer be used to perform automated and intelligent identification and early warning on the squall line in the radar echo image, and that the requirement for high accuracy of identification be met.

The invention content is as follows:

in order to realize objective quantitative identification of squall line characteristics which may cause thunderstorms, lightning, hails or tornadoes and other disastrous weather in weather, the invention provides a method for intelligently identifying and early warning squall line characteristics in radar echo images, and the specific technical scheme adopted by the invention is as follows:

an intelligent squall line characteristic identification early warning method in a radar echo image comprises the following specific steps:

step 1, defining the elevation angle number of radar volume sweep in a radar base data file as α, wherein the value range of the elevation angle number is [0, max _ α ], and max _ α represents the maximum elevation angle number of the radar volume sweep, the initial value of α is 0, namely α is 0, defining a high-pass filtering threshold value as Filt _ A, the value range of Filt _ A is a positive integer from 0 to 100, and the initial value of Filt _ A is any positive integer between [60 and 100 ];

step 2, reading a radar base data file to be analyzed, and extracting basic reflectivity data R (α) with an elevation number of α from the file;

step 3, converting R (α) from a polar coordinate form to a plane rectangular coordinate form, and marking as RP (α, x, y), wherein α is an elevation number, x and y are respectively an abscissa and an ordinate of a radar detection point, and the radar detection point is hereinafter referred to as a pixel;

step 4, carrying out high-pass filtering processing on the basic reflectivity data in the RP (α, x, y), wherein the high-pass filtering threshold value is Filt _ A, and recording the basic reflectivity data subjected to the high-pass filtering processing as RPF (α, x, y);

step 5, carrying out binarization processing on the RPF (α, x, y), setting the pixel value subjected to high-pass filtering in the step 4 as 0, setting the pixel value not subjected to high-pass filtering in the step 4 as 1, and marking the pixel value as RPFB (α, x, y);

6, carrying out noise reduction treatment on the RPFB (α, x, y), eliminating isolated clutter points and forming isolated points after high-pass filtering;

and 7, identifying the central axis of the target object from the RPFBD (α, x, y), and specifically comprising the following steps:

step 7.1, sequentially analyzing each pixel value in RPFBD (α, x, y), if the current pixel value is 0, making the central axis identification value of the pixel be 0, and continuously analyzing the next pixel, if the current pixel value is 1, counting the number of pixels C1 with the pixel value being 1 in 8 pixels immediately adjacent to the current pixel value, if C1 is 3, making the central axis identification value of the pixel be 0, and continuously analyzing the next pixel value, if C1 is not 3, making the central axis identification value of the pixel be 1, and continuously analyzing the next pixel value, and after the analysis of each pixel value in RPFBD (α, x, y) is finished, replacing the central axis identification value of each pixel with the pixel value, and outputting a result of RPFBD _ k1(α, x, y);

and 7.2, sequentially analyzing each pixel value in RPFBD _ k1(α, x, y), if the current pixel value is 0, making the central axis identification value of the pixel be 0, continuously analyzing the next pixel, if the current pixel value is 1, counting the number of pixels C2 with the pixel value of 1 in 8 pixels immediately adjacent to the current pixel value, and if C2 belongs to [3, 4] s]If so, the central axis identification value of the pixel is made to be 0, and the next pixel value is continuously analyzed; if it is not

When each pixel value in RPFBD _ k1(α, x, y) is analyzed, replacing the pixel value with the central axis identification value of each pixel, and outputting a result which is recorded as RPFBD _ k2(α, x, y);

and 7.3, sequentially analyzing each pixel value in RPFBD _ k2(α, x, y), if the current pixel value is 0, making the central axis identification value of the pixel be 0, continuously analyzing the next pixel, if the current pixel value is 1, counting the number of pixels C3 with the pixel value of 1 in 8 pixels immediately adjacent to the current pixel value, and if C3 belongs to [3,4, 5] s]If so, the central axis identification value of the pixel is made to be 0, and the next pixel value is continuously analyzed; if it is not

When each pixel value in RPFBD _ k2(α, x, y) is analyzed, replacing the pixel value with the central axis identification value of each pixel, and outputting a result which is recorded as RPFBD _ k3(α, x, y);

and 7.4, sequentially analyzing each pixel value in RPFBD _ k3(α, x, y), if the current pixel value is 0, making the central axis identification value of the pixel be 0, continuously analyzing the next pixel, if the current pixel value is 1, counting the number of pixels C4 with the pixel value of 1 in 8 pixels immediately adjacent to the current pixel value, and if C4 belongs to [3,4,5, 6]]If so, the central axis identification value of the pixel is made to be 0, and the next pixel value is continuously analyzed; if it is not

When each pixel value in RPFBD _ k3(α, x, y) is analyzed, replacing the pixel value with the central axis identification value of each pixel, and outputting a result which is recorded as RPFBD _ k4(α, x, y);

and 7.5, sequentially analyzing each pixel value in RPFBD _ k4(α, x, y), if the current pixel value is 0, making the central axis identification value of the pixel be 0, continuously analyzing the next pixel, if the current pixel value is 1, counting the number of pixels C5 with the pixel value of 1 in 8 pixels immediately adjacent to the current pixel value, and if C5 belongs to [3,4,5,6, 7] for]If so, the central axis identification value of the pixel is made to be 0, and the next pixel value is continuously analyzed; if it is not

When each pixel value in RPFBD _ k4(α, x, y) is analyzed, replacing the pixel value with the central axis identification value of each pixel, and outputting a result which is recorded as RPFBD _ k5(α, x, y);

step 7.6, performing multiple rounds of iterative processing on the RPFBD _ k5(α, x, y) according to the method from step 7.1 to step 7.5, taking the RPFBD _ k5(α, x, y) of the previous round as the RPFBD (α, x, y) of step 6.1 in the next round, stopping iteration until no new pixel value is set to be 0, and recording the output result as RPFBD _ kn (α, x, y);

step 8, determining whether the squall line exists in the RPFBD _ kn (α, x, y);

and step 9: defining a data set Rlt and a counter FN, wherein the initial value of FN is 0;

if the RPFBD _ kn (α, x, y) in the step 8 exists in the squall line, recording a data group [ α, Filt _ A, FN ] into a data set Rlt, and entering the step 10, if the squall line does not exist in the step 8, firstly making a counter FN +1 and a high-pass filtering threshold Filt _ A as Filt _ A-FN × Filt _ B, wherein Filt _ B is a fixed parameter and takes a positive integer from 1 to 10, then judging the sizes of Filt _ A and the fixed parameter Filt _ N, if Filt _ A is not less than Filt _ N, returning to the step 4, otherwise, if Filt _ A is less than Filt _ N, entering the step 10, wherein the value range of Filt _ N is 0< Filt _ N < Filt _ A;

step 10, resetting the elevation angle number α to α +1 and the high-pass filtering threshold Filt _ a to an initial value, returning to step 2, continuing processing according to the method from step 2 to step 9 until α is max _ α, and entering step 10;

step 11, analyzing Rlt the data set, if Rlt is empty, the squall line was not identified by the radar-based data file to be analyzed in step 2, if Rlt is not empty, the squall line was identified by the radar-based data file to be analyzed in step 2, and the squall line characteristics became stronger the greater the number of data sets recorded in Rlt, and the squall line characteristics became stronger the smaller the FN value in the data sets [ α, Filt _ A, FN ] recorded in Rlt, indicating the squall line characteristics became stronger in the elevation plane, otherwise, the squall line characteristics became weaker.

Preferably, the specific step method of the noise reduction processing in step 6 is as follows:

step 6.1, starting from RPFB (α,0,0), defining a square area with the side length of SL, respectively counting the number of pixels with the pixel value of 0 and the pixel value of 1 in the square area, and respectively recording as Stat _0 and Stat _1, wherein SL is a positive integer and SL belongs to [3,50 ];

step 6.2: if Stat _0> kxStat _1, setting all pixel values in the square area to be 0; otherwise, setting all pixel values in the square area as 1, wherein k is an empirical threshold, and belongs to [0.1,10 ];

and 6.3, sliding the square area to the right or downwards by taking SL as a step size, and repeating the processing procedures of the steps 5.1 and 5.2 until all areas in the RPFB (α, x, y) are processed by the steps 5.1 and 5.2, wherein the processed result is recorded as RPFBD (α, x, y).

Preferably, the specific steps for determining whether a squall line exists in step 8 are as follows:

step 8.1, a straight line r (X ', y', omega) is defined in a normal line mode, wherein r represents the shortest distance from the coordinate origin to the straight line, the absolute value of the maximum value of r is an integer closest to the diagonal length of the image corresponding to RPFBD _ kn (α, X, y), omega represents the included angle of the straight line and the positive direction of the X axis, omega belongs to [0 DEG ], 180 DEG and (X ', y') represents a point in a plane rectangular coordinate system;

step 8.2, constructing a two-dimensional matrix MAT [ p, q ] by taking omega _ i as a horizontal unit step length and r _ i as a vertical unit step length, wherein p represents a row, p belongs to [ -n.r _ i, - (n-1). r _ i,..,. 2. r _ i,0, r _ i,2r _ i,.., (n-1). r _ i, n.r _ i ], q represents a column, q belongs to [0, omega _ i, 2. omega _ i, 3. omega _ i,. 180. omega _ i ],. all units in MAT [ p, q ] have initial values of 0, omega _ i is divided by the included angle between a straight line and an X axis, omega _ i takes a positive integer which can be divided by 180, r _ i is a granularity divided by the shortest distance between the straight line and an origin, n.r _ i and n.r _ i represent the granularities divided by the positive directions of the straight line and the positive direction of the X axis, and the maximum rpr _ i and the diagonal (r _ i) is larger than 5 and the maximum rpr _ i) is equal to the corresponding index of a scale, and is equal to the shortest distance, and the maximum rpr _ i, and the opposite to the maximum rpr _ i, and the opposite scale, and the opposite to the opposite scale;

step 8.3, sequentially traversing each pixel in RPFBD _ kn (α, x, y), skipping if the current pixel value is 0, and continuing traversing the next pixel, if the current pixel value is 1, setting the coordinate of the current pixel as (i, j), namely x ═ i, y ═ j, substituting i into x ', j into y' according to the linear expression in step 8.1, and then substituting q, q ∈ [0, ω _ i,2 ω _ i,3 ω _ i, ·,180- ω _ i ] in the MAT transverse direction into ω in the linear expression in step 8.1 one by one, thereby obtaining a series of r (i, j, q);

step 8.4: analysis of MAT [ p, q ]]The numerical value of each unit in the table is screened out

The squall line is absent in the RPFBD _ kn (α, x, y) if the number of screened units is 0, the squall line is present in the RPFBD _ kn (α, x, y) if the number of screened units is more than 0, and the number of screened units is the number of the central axis of the squall line, wherein Max (MAT) represents MAT [ p, q ] and]the maximum value of all unit values in the table; trd is an empirical parameter, Trd ∈ (0, 1); the smaller the Trd value is, the lower the missing report rate of the squall line is, but the false report rate is increased; conversely, the larger the Trd value is, the higher the missing report rate of squall line identification is, but the false report rate is reduced.

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

the invention provides an intelligent squall line characteristic identification early warning method in a radar echo image, which takes Doppler weather radar detection data as a main data source, analyzes the spatial distribution and the intensity of basic reflectivity detected by a radar, and realizes intelligent squall line characteristic identification early warning on the radar through a series of steps of numerical value preprocessing, filtering, image characteristic extraction, central axis analysis on a target object, squall line shape analysis and the like. The application of the method can automate and objectify the work of subjective analysis and interpretation of radar echo images by weather professionals in the past, and improves the accuracy and timeliness of the related services of squall line identification and strong convection weather warning.

In practical business, a meteorological professional cannot accurately determine whether a squall line exists only by studying and reading a radar echo image of a certain elevation plane. The intelligent identification early warning method can comprehensively analyze the basic reflectivity data of different elevation surfaces, and further improves the accuracy and reliability of automatically identifying the characteristics of the squall line.

Description of the drawings:

FIG. 1 is a control flow chart in the embodiment;

FIG. 2 is a schematic structural diagram of a matrix MAT according to an embodiment;

FIG. 3 is a basic reflectivity image plotted from the radar data R (0) in the example;

FIG. 4 is an image drawn by RP (0, x, y) in the example;

FIG. 5 is an image rendered by RPF (0, x, y) in an embodiment;

FIG. 6 is an image rendered by RPFB (0, x, y) in an embodiment;

FIG. 7 is an image rendered by RPFBD (0, x, y) in an embodiment;

fig. 8 is an image drawn by RPFBD _ kn (0, x, y) in the embodiment;

FIG. 9 is a schematic line representation illustrating squall line locations in an embodiment.

The specific implementation mode is as follows:

the invention is described in detail below with reference to the figures and the specific embodiments.

The first embodiment is as follows:

the embodiment of the invention provides an intelligent squall line characteristic identification and early warning method in a radar echo image, which comprises the following specific steps:

step 2, reading a radar base data file to be analyzed, and extracting basic reflectivity data with an elevation angle number of α from the file, wherein the basic reflectivity data is marked as R (α);

step 3, converting R (α) from a polar coordinate form to a plane rectangular coordinate form, and marking as RP (α, x, y), wherein α is an elevation number, x and y respectively represent the abscissa and ordinate of radar detection points, each radar detection point is determined by (α, x, y), and the radar detection points are referred to as pixels;

and 5, carrying out binarization processing on the RPF (α, x, y), setting the pixel value subjected to high-pass filtering in the step 3 to be 0, and setting the pixel value not subjected to filtering to be 1, thereby obtaining the RPF (α, x, y) after binarization, and marking the RPF as RPFB (α, x, y).

6, carrying out noise reduction treatment on the RPFB (α, x, y), eliminating isolated clutter points and isolated points formed after high-pass filtering, wherein the noise reduction treatment specifically comprises the following steps:

step 6.2: if Stat _0> kxStat _1, setting all pixel values in the square area to be 0; otherwise, all pixel values in the square area are set to 1, where k is an empirical threshold, and k ∈ [0.1,10 ].

And 6.3, sliding the square area to the right or downwards by taking SL as a step size, and repeating the processing procedures of the steps 6.1 and 6.2 until all areas in the RPFB (α, x, y) are processed by the steps 6.1 and 6.2, wherein the processed result is marked as RPFBD (α, x, y), for example, square areas RPFB (α,0,0) -RPFB (α 1, SL-1, SL-1), RPFB (α,0, SL) -RPFB (α, SL-1, 2SL-1), RPFB (α,0, 2SL) -RPFB (α -1, 3SL-1),. 9., RPFB (α -1, 0) -RPFB (α, 2SL-1, SL-1), RPFB (α -1), SL-FB (α, 2-1), RPFB (3662-1, 2 SL-631), SL-FB (α, α), SL-FB (α), SL-1, SL-2, and so on.

It should be added that, if the maximum value of x or y in the RPFB (α, x, y) is not an integral multiple of SL, there is a special case that when a square region with a side length of SL slides to the edge of the RPFB (α, x, y) image, i.e. the rightmost side or the bottommost side of the RPFB (α, x, y), a side length overflow occurs, and at this time, the pixel with the square region overflow is set to 0 by default, and then the processing is performed according to the calculation procedures of steps 6.1 and 6.2.

Step 7, identifying a central axis of the target object from the RPFBD (α, x, y), in order to express more clearly, regarding the RPFBD (α, x, y) in a two-dimensional array form as a binary image, regarding a point with a pixel value of 1 in the RPFBD (α, x, y) as a radar basic reflectivity high value area in the image, namely the target object to be identified, and regarding a point with a pixel value of 0 as a background image, wherein the target object is represented as a plurality of color blocks with uneven size on the image, and the purpose of the step is to extract the central axis of the color block from the color blocks, and the specific steps are as follows:

step 7.1, analyzing each pixel value in RPFBD (α, x, y) in sequence from a fixed direction, such as from top to bottom and from left to right, if the current pixel value is 0, making the central axis identification value of the pixel be 0, continuing to analyze the next pixel, if the current pixel value is 1, counting the number of pixels C1 with the pixel value of 1 in 8 pixels immediately adjacent to the current pixel, if C1 is 3, namely, there are 3 pixels with high radar basic reflectivity around the current pixel, making the central axis identification value of the pixel be 0, continuing to analyze the next pixel value, if C1 is not equal to 3, making the central axis identification value of the pixel be 1, continuing to analyze the next pixel value, and after the analysis of each pixel value in RPFBD (α, x, y) is completed, replacing the central axis identification value of each pixel with the pixel value, outputting the result of RPFBD _ k2, x _ k, y (RPFBD, rpy), and traversing all pixels in the steps of RPFBD (α, x, y).

And 7.2, analyzing each pixel value in RPFBD _ k1(α, x, y) from a fixed direction, such as from top to bottom and from left to right in turn, if the current pixel value is 0, making the central axis identification value of the pixel be 0, continuing to analyze the next pixel, if the current pixel value is 1, counting the number of pixels C2 with the pixel value of 1 in 8 pixels immediately adjacent to the current pixel value, and if C2 belongs to [3, 4] E]If 3 or 4 pixels with high radar basic reflectivity exist around the current pixel, the central axis identification value of the pixel is 0, and the next pixel value is continuously analyzed; if it is not

The central axis identification value of the pixel is set to be 1, the next pixel value is analyzed continuously, after each pixel value in RPFBD _ k1(α, x, y) is analyzed, the central axis identification value of each pixel replaces the pixel value, the output result is recorded as RPFBD _ k2(α, x, y), and the step completes the traversal of all pixels in RPFBD _ k1(α, x, y).

And 7.3, analyzing each pixel value in RPFBD _ k2(α, x, y) from a fixed direction, such as from top to bottom and from left to right in turn, if the current pixel value is 0, making the central axis identification value of the pixel be 0, continuing to analyze the next pixel, if the current pixel value is 1, counting the number of pixels C3 with the pixel value of 1 in 8 pixels immediately adjacent to the current pixel value, and if C3 belongs to [3,4, 5] E]If 3-5 pixels with high radar basic reflectivity exist around the current pixel, the central axis identification value of the pixel is 0, and the next pixel value is continuously analyzed; if it is not

The central axis identification value of the pixel is set to be 1, the next pixel value is analyzed continuously, after each pixel value in RPFBD _ k2(α, x, y) is analyzed, the central axis identification value of each pixel replaces the pixel value, the output result is recorded as RPFBD _ k3(α, x, y), and the step completes the traversal of all pixels in RPFBD _ k2(α, x, y).

And 7.4, analyzing each pixel value in RPFBD _ k3(α, x, y) from a fixed direction, such as from top to bottom and from left to right in turn, if the current pixel value is 0, making the central axis identification value of the pixel be 0, continuing to analyze the next pixel, if the current pixel value is 1, counting the number of pixels C4 with the pixel value of 1 in 8 pixels immediately adjacent to the current pixel value, and if C4 belongs to [3,4,5, 6] E]If 3-6 pixels with high radar basic reflectivity exist around the current pixel, the central axis identification value of the pixel is 0, and the next pixel value is continuously analyzed; if it is not

Make the pixel inThe axis identification value is 1, the next pixel value is analyzed continuously, when each pixel value in RPFBD _ k3(α, x, y) is analyzed, the axis identification value of each pixel replaces the pixel value, the output result is recorded as RPFBD _ k4(α, x, y), and the step completes the traversal of all pixels in RPFBD _ k3(α, x, y).

And 7.5, analyzing each pixel value in RPFBD _ k4(α, x, y) from a fixed direction, such as from top to bottom and from left to right in turn, if the current pixel value is 0, making the central axis identification value of the pixel be 0, continuing to analyze the next pixel, if the current pixel value is 1, counting the number of pixels C5 with the pixel value of 1 in 8 pixels immediately adjacent to the current pixel value, and if C5 belongs to [3,4,5,6, 7] for the next pixel]If 3-7 pixels with high radar basic reflectivity exist around the current pixel, the central axis identification value of the pixel is 0, and the next pixel value is continuously analyzed; if it is not

And when each pixel value in the RPFBD _ k4(α, x, y) is analyzed, replacing the central axis identification value of each pixel with the pixel value, and outputting the result as RPFBD _ k5(α, x, y).

And 7.6, shrinking each color block represented by the radar basic reflectivity high value area to the center position thereof through the calculation of the steps 7.1 to 7.5, wherein the position of each color block is not changed but the area is reduced visually from the image, in order to obtain the central axis of each color block, performing multiple rounds of iterative processing according to the methods in the steps 7.1 to 7.5, wherein the RPFBD _ k5(α, x and y) of the previous round is used as the RPFBD (α, x and y) of the step 7.1 in the next round, and the iterative calculation is stopped until no new pixel value is set to be 0, and the output result is recorded as RPFBD _ kn (α, x and y), and the image of the RPFBD _ kn (α, x and y) presents a large number of disorderly, free and inconsistent line segments.

In other words, if a squall line exists, one or more straight lines should be found from RPFBD _ kn (α, x, y) such that a large number of disorderly, chapter-free, and short line segments in RPFBD _ kn (α, x, y) are generally distributed around the straight lines, in order to determine whether the central axis of the squall line exists in the RPFBD _ kn (α, x, y), the following steps are employed for analysis:

step 8.1, defining a straight line r (X ', y', omega) by adopting a normal line, wherein r represents the shortest distance from the origin of coordinates to the straight line, and theoretically r belongs to [ - ∞, ∞ ], in the invention, the absolute value of the maximum value of r is an integer closest to the diagonal length of an image corresponding to RPFBD _ kn (α, X, y), [ omega ] represents the included angle of the straight line and the positive direction of the X axis, and belongs to [ 0], 180 DEG ], and (X ', y') represents a point in a plane rectangular coordinate system, wherein (r, omega) is fixed and unchanged for any straight line in the plane rectangular coordinate system, and different (r, omega) form different straight lines;

step 8.2, constructing a two-dimensional matrix MAT [ p, q ] with ω _ i as a horizontal unit step and r _ i as a vertical unit step, wherein p represents a row, p [ -n · r _ i, - (n-1) · r _ i, · 2 · r _ i, -r _ i,0, r _ i,2r _ i,. q [ -1) · r _ i, n · r _ i ], q represents a column, q ∈ [0, ω _ i,2 ω _ i,3 ω _ i,. and 180 ω _ i ], (n-1) · r _ i, n · r _ i ], initial values of all cells in MAT [ p, q ] are 0, ω _ i is used to define a granularity of division of a straight line with a positive direction of the X-axis, the smaller the value of ω _ i, the number of scales on the horizontal (column), the finer the included angle is the finer, the smaller the number of scales, the finer the included angle is the finer the smaller the number of scales, the smaller the included angle is the smaller the number of scales, the smaller the number of the smaller the included angle is the smaller the linear scale is the smaller the linear r _ i is the smaller the linear r is the longer the linear r is the.

Step 8.3, sequentially traversing each pixel in RPFBD _ kn (α, x, y), if the current pixel value is 0, skipping, and continuing traversing the next pixel, if the current pixel value is 1, setting the coordinate of the current pixel as (i, j), namely x ═ i, y ═ j, substituting i into x ', j into y ', and further substituting q, q ∈ [0, ω _ i,2 ω _ i,3 ω _ i, ·,180- ω _ i ] in the MAT transverse direction into ω in the linear expression in step 8.1 one by one according to the linear expression in step 8.1, thereby obtaining a series of r (i, j, q), for example, r (i, j, 0), r (i, j, ω _ i), r (i, j, 2 ω _ i),. r (i, j, 180- ω _ i), for each r (i, j, q), finding any one of the absolute difference between p and p [ +1, p ', p ];

If the number of screened units is greater than 0, the squall line exists in RPFBD _ kn (α, x, y);

the screened units are the number of central axes of squall lines identified by the method, the number may be 0 or an integer greater than 0, the (p, q) corresponding to the units respectively represent a straight line, and as can be seen from the image, a large number of disorderly, seal-free and long and short line segments in RPFBD _ kn (α, x, y) are distributed just around the straight lines, wherein Max (MAT) represents the maximum value of all unit values in MAT [ p, q ], Trd is an empirical parameter, Trd belongs to (0,1), the value of the Trd belongs to the meteorological factors such as the number of pixels in a radar basic reflectivity image, seasons in which strong convection weather easily occurs, geographic spaces in which radar is detected, and the like, the smaller the Trd value is, the lower the missing rate of identifying squall lines is, but the false alarm rate is increased, otherwise, the higher the missing rate of identifying squall lines is higher, but the false alarm rate is reduced, and when the identification is actually used, the weather phenomenon occurring is strong convection, the larger the missing rate of identifying squall lines is analyzed, and the success rate of identifying the historical squall lines is improved.

if the RPFBD _ kn (α, x, y) exists in step 8, the data set [ α, Filt _ a, FN ] is recorded in the dataset Rlt and step 10 is entered, if no squall line exists in step 8, the counter FN +1 is first set and the high pass filter threshold Filt _ a is set to Filt _ a-FN × Filt _ B, where Filt _ B is a constant parameter and takes a positive integer from 1 to 10, the magnitudes of Filt _ a and Filt _ N are determined, if Filt _ a ≧ Filt _ N, step 4 is returned, otherwise, if Filt _ a < Filt _ N, step 10 is entered, where Filt _ N is a constant parameter and 0< Filt _ N < Filt _ a, and the physical meaning of Filt _ N is that a basic reflectivity smaller than Filt _ N is unlikely to occur.

And step 10, resetting the elevation angle number α to α +1, resetting the mid-high pass filter threshold to an initial value, returning to step 2, continuing the processing from step 2 to step 9 until α is max _ α, and entering step 11.

Step 11, analyzing Rlt a data set, wherein the data set may be empty or may be composed of 1 or more [ α, Filt _ a, FN ] records, if Rlt a data set is empty, the radar-based data file to be analyzed in step 2 does not identify a squall line, if Rlt a data set is not empty, the radar-based data file to be analyzed in step 2 identifies a squall line, and the greater the number of data sets recorded in Rlt indicates that a squall line is identified in each of a plurality of radar detection elevation planes, the stronger the squall line characteristic is, the smaller the value of FN in [ α, Filt _ a, FN ] the data set recorded in Rlt starts from 0, the smaller the value indicates that the squall line characteristic is stronger on the elevation plane, otherwise, the squall line characteristic is relatively weaker.

The information recorded by the data set Rlt may be used to provide objective and quantitative squall line evaluation indexes directly to weather professionals, or may be integrated into a relevant service system to achieve the functional goals of squall line identification and early warning in a richer presentation form.

The first application embodiment:

in the application embodiment, Doppler weather radar detection data of squall line weather phenomenon is selected as analysis data of the application embodiment of the invention. The values of various parameters related to the invention in the application embodiment are as follows:

serial number	Parameter name	Parameter value
				1	max_α	2
2	Filt_A	60dBZ
				3	SL	5
4	k	0.25
			5	r_i	20
6	ω_i	20
			7	n	57
8	n·r_i	1140
			9	Trd	0.75
10	Filt_B	5
			11	Filt_N	50

As shown in fig. 1, the embodiment of the present application includes the following specific steps:

step 1, defining a variable α ∈ [0, max _ α ], α represents the elevation angle number of the radar body sweep, max _ α represents the maximum elevation angle number of the radar body sweep, max _ α is determined by the hardware model of the radar and the body sweep parameters, α has an initial value of 0, namely α ═ 0, defining a variable Filt _ a, which is used for recording a threshold value of high-pass filtering, the initial value of Filt _ a usually takes a positive integer between [60 and 100], in the embodiment of the application, max _ α ═ 2, and Filt _ a ═ 60 dBZ.

Step 2, reading a radar base data file to be analyzed, and extracting α basic reflectivity data on the elevation surface from the file, wherein the basic reflectivity data is marked as R (α). in the embodiment of the application, a basic reflectivity image drawn by R (0) is shown in FIG. 3.

Step 3, converting R (α) from a polar coordinate form into data in a plane rectangular coordinate form^[1，2]The average of the values denoted as RP (α, x,y). An image drawn by RP (0, x, y) is shown in fig. 4.

And 4, carrying out high-pass filtering on the basic reflectivity value in the RP (α, x, y), wherein the filtering threshold value is Filt _ A, the high-pass filtered basic reflectivity value is marked as RPF (α, x, y), wherein α is an elevation angle number, x and y respectively represent the abscissa and the ordinate of radar detection points, each radar detection point is determined by (α, x, y), and is hereinafter referred to as a pixel, and an image drawn by the RPF (0, x, y) is shown in FIG. 5.

And step 5, carrying out binarization processing on the RPF (α, x, y), setting all filtered pixel values in the step 4 to be 0 and all unfiltered pixel values to be 1, thereby obtaining the RPF (α, x, y) after binarization, and marking the RPF as RPFB (α, x, y), wherein an image drawn by the RPFB (0, x, y) is shown in FIG. 6.

And 6, carrying out noise reduction treatment on the RPFB (α, x, y), and simultaneously removing isolated clutter points and isolated points filtered by high-pass filtering, wherein the specific method comprises the following steps:

step 6.1, starting from RPFB (α,0,0), a square with a side length SL is defined, the number of pixels with pixel values of 0 and 1 in the square area is counted respectively and is recorded as Stat _0 and Stat _ l, where SL is a positive integer and SL belongs to [3,50], and in this application example, SL is 5.

Step 6.2: if Stat _0> kxStat _1, setting the values of all pixels in the square area to 0; conversely, the values of all pixels in the square region are set to 1, where k is an empirical threshold, and k ∈ [0.1,10], and in this embodiment, k is 0.25.

Step 6.3, sliding the square area to the right or downward by using SL as a step size, and repeating the calculation processes of steps 6.1 and 6.2, for example, square areas RPFB (α,0,0) -RPFB (α,4, 4), RPFB (

α

0,0, 5) -RPFB (α,4, 9), RPFB (α,0, 10) -RPFB (α,4, 14),.. once, RPFB (α,4, 0) -RPFB (α, 9, 4), RPFB (7, 4, 5) -RPFB (α, 9, 9), RPFB (α,4, 10 RPFB (α, 9, 14),. once, and so on, until all areas in RPFB (α, x, y) are subjected to calculation once for 6.1 and 6.2, RPFB (α, x, y) is recorded as fbd (fbd) by step 6, fby), and then drawing an image as shown in fig. 7, 7.

It should be added that, if the maximum value of x or y in RPFB (α, x, y) is not an integral multiple of SL, there is a special case that when a square with a side length of SL slides to the rightmost side or the bottommost side of RPFB (α, x, y) with SL as a sliding step, side length overflow occurs, and at this time, the pixel with square overflow is set to 0 by default, and then the processing is performed according to the calculation procedures of steps 6.1 and 6.2.

Step 7, identifying a central axis of a target object from RPFBD (α, x, y), in order to express more clearly, the RPFBD (α, x, y) in a two-dimensional array form can be regarded as a binary image as shown in FIG. 7, wherein a point with a pixel value of 1 corresponds to a radar basic reflectivity high value area in the image, namely the target object to be identified, a point with a pixel value of 0 is regarded as a background image, and the point is black in FIG. 7. the target object is represented as a plurality of white blocks with uneven size on the image, the purpose of the step is to extract the central axis of the white blocks from the white blocks, and the specific steps are as follows:

step 7.1, analyzing the value of each pixel point in RPFBD (α, x, y) from a fixed direction, such as from top to bottom and from left to right, skipping if the value is 0, and continuing to analyze the next pixel, and if the value is 1, proceeding to step 7.2.

And 7.2, counting the number of pixels C1 with the pixel value of 1 in 8 pixels immediately adjacent to the coordinate (i, j) by taking the coordinate (i, j) as the center, setting the value of the current position to be 0 if C1 is 3, namely 3 pixels with high radar basic reflectivity exist around RPFBD (α, i, j), namely RPFBD (α, i, j) is 0, completing the traversal of all the pixels according to the steps 7.1 and 7.2, and outputting the result to be marked as RPFBD _ k1(α, x, y).

And 7.3, sequentially analyzing the value of each pixel point in RPFBD _ k1(α, x, y), skipping and continuously analyzing the next pixel if the value is 0, counting the number C2 of the pixel value of 1 in 8 pixels immediately adjacent to the pixel value if the value is 1, and setting the value of the current position to be 0 if C2 belongs to [3, 4], namely 3 or 4 pixels with high radar basic reflectivity exist around the current pixel, finishing the traversal of all the pixel points in the step, and marking the output result as RPFBD _ k2(α, x, y).

And 7.4, similar to the step 7.3, sequentially analyzing the value of each pixel point in the RPFBD _ k2(α, x, y), skipping and continuously analyzing the next pixel if the value is 0, counting the number C3 of the pixel value 1 in 8 pixels immediately adjacent to the pixel value if the value is 1, and setting the value of the current position to be 0 if the value C3 belongs to [3,4,5 ].

And 7.5, similar to the step 7.4, sequentially analyzing the value of each pixel point in the RPFBD _ k3(α, x, y), skipping and continuously analyzing the next pixel if the value is 0, counting the number C4 of the pixel value 1 in 8 pixels immediately adjacent to the pixel value if the value is 1, and setting the value of the current position to be 0 if the value C4 belongs to [3,4,5,6 ].

And 7.6, similar to the step 7.5, sequentially analyzing the value of each pixel point in the RPFBD _ k4(α, x, y), skipping and continuously analyzing the next pixel if the value is 0, counting the number of pixels C5 with the pixel value of 1 in 8 pixels immediately adjacent to the pixel value if the value is 1, and setting the value of the current position to be 0 if C4 belongs to [3,4,5,6,7 ].

In the analysis process of the step 7.1 to the step 7.6, for the pixel point at the image edge, a part of the pixel immediately adjacent to the pixel point exceeds the image edge, at this time, the pixel value exceeding the image edge is set as 0, and then analysis is carried out;

and 7.7, after the calculation of the steps 7.1 to 7.6, shrinking each color block represented by the radar basic reflectivity high value area to the center position, wherein the position of each color block is not changed but the area is reduced visually from the image, in order to obtain the central axis of each color block, multiple rounds of iterative calculation of the steps 7.1 to 7.6 are required, the RPFBD _ k5(α, x, y) of the previous round is used as the RPFBD (α, x, y) of the step 7.1 in the next round, the iterative calculation is stopped until no new pixel value is set to be 0, the output result is recorded as RPFBD _ kn (α, x, y), the image of the RPFBD _ kn (α, x, y) is presented as a large number of disordered and non-uniform line segments, as shown in fig. 8, and fig. 8 shows a pure black background image around for clarity.

step 8.1, a straight line r (X ', y', omega) is defined in a normal line mode, wherein r represents the shortest distance from the coordinate origin to the straight line, and theoretically r belongs to [ - ∞, ∞ ], in the invention, the absolute value of the maximum value of r is an integer closest to the diagonal length of the image corresponding to RPFBD _ kn (α, X, y), omega represents the included angle of the straight line and the positive direction of the X axis, and belongs to [0 DEG ], 180 DEG, and (X ', y') represents a point in a plane rectangular coordinate system, wherein (r, omega) is fixed and is not changed for any straight line in the plane rectangular coordinate system, and different (r, omega) form different straight lines.

Step 8.2: using omega _ i as horizontal unit step length and r _ i as vertical unit step length to construct a two-dimensional matrix MAT [ p, q ]]Where p represents a row, p [ -n-r _ i, - (n-1) r _ i,. -, -2-r _ i,0, r _ i,2r _ i, -, (n-1) r _ i, n-r _ i]Q denotes a column, q ∈ [0, ω _ i,2 ω _ i,3 ω _ i]。MAT[p，q]The initial values of all the cells in the list are 0. The smaller the value of the omega _ i, the more the number of scales on the transverse direction (column) is, the finer the granularity of the angle division is, and the omega _ i is generally taken asSimilarly, r _ i is used to define the granularity divided by the shortest distance between a straight line and the origin, the smaller the value of r _ i, the larger the number of scales in the longitudinal direction (row), the finer the granularity divided by the shortest distance, n · r _ i and-n · r _ i respectively represent the maximum scale and the minimum scale in the longitudinal direction (row) of the MAT, and r _ i is usually 1, 2, 5 or an integer greater than 5, and n · r _ i is closest to the diagonal length of the image corresponding to RPFBD _ kn (α, x, y) and can be divided by r _ i

Where n · r _ i takes the diagonal length closest to the image and can be divided by r _ i, where n · r _ i equals 1140 and n equals 57.

Step 8.3, sequentially traversing each pixel in RPFBD _ kn (α, x, y), if the value of the pixel is 0, skipping no processing, and then traversing the next pixel until all pixels are traversed once, if the value of the pixel is 1, further analyzing if the coordinate of the current pixel is (i, j), i.e., x ═ i, y ═ j, i is substituted into x', j according to the linear expression in step 8.1, and then substituting each q, q ∈ [0, ω _ i,2 ω _ i,3 ω _ i, and.., 180 — ω _ i ] in the linear expression in step 8.1, so as to obtain a series of r (i, j, q) such as r (i, j, 0), r (i, j, ω _ i), r (i, j, r [, i, j, r [ -, r (i, j, ω _ i), r (i, j, r, j, ω _ i, r, i, j, q), r (i, q, i, q, i, q, i, q, i, q, i, q, i, q, i, q, i, q, i, q, i, q, i, q, i.

Step 8.4: analysis of MAT [ p, q ]]The numerical value of each unit is screened out

The (p, q) s corresponding to these cells each represent a straight line around which a number of randomly and randomly sized line segments in RPFBD _ kn (α, x, y) are distributed, where Max (MAT) represents MAT [ p, q ] in which]The maximum value of all unit values in the radar image is Trd, the Trd belongs to an empirical parameter (0,1), and the value of the Trd is closely related to meteorological factors such as the number of pixels in the radar basic reflectivity image, seasons in which strong convection weather easily occurs, geographic space detected by the radar and the like. The smaller the Trd value is, the lower the missing report rate of the squall line is, but the false report rate is increased; conversely, the larger the Trd value is, the higher the missing report rate of squall line identification is, but the false report rate is reduced. During actual application, the parameter value of the Trd can be corrected through historical case analysis on the phenomenon that strong convection weather and squall line weather occur, and then the success rate of squall line identification is improved. In the present application example, Trd is 0.75. Calculated from step 8 for RPFBD _ kn (0, x, y), there is and only MAT [80, 460 ]]The value of the cell satisfies

Fig. 9 shows a straight line drawn on the basis of fig. 8 by substituting r (x ', y', ω) into r (x ', y', ω) 460.

Step 9 if the squall line was identified, then a set of data [ α, Filt _ A, FN ] is recorded into the data set Rlt, proceeding to the next step, where FN is a counter with an initial value of 0.

If no squall line is identified, FN +1 and Filt _ A is Filt _ A-FN × Filt _ B, where Filt _ B is a constant value parameter, typically a positive integer from 1 to 10. At the moment, if Filt _ A is not less than Filt _ N, returning to the step 4; otherwise, if Filt _ A < Filt _ N, go to step 10. Where Filt _ N is a fixed parameter, 0< Filt _ N < Filt _ A, and Filt _ N has the physical meaning that a base reflectivity less than Filt _ N is unlikely to occur for the squall line. In the embodiment of this application, Filt _ N is 50. In the present exemplary embodiment, since the squall line is already identified in this step when FN is 0, FN +1 and subsequent computations in this step need not be performed. Currently, Rlt { [0, 60, 0] }.

And step 10, resetting α to α +1, resetting Filt _ A to an initial value, returning to step 2, continuing to calculate in steps 2 to 9 until α is max _ α, and entering step 11.

Step 11, analyze Rlt a dataset, which may be empty or may consist of 1 or more [ α, Filt _ A, FN ] records, Rlt is practiced to indicate that a squall line is not identified by the analyzed radar data when empty, indicate that a squall line is identified when it is not empty, indicate that a squall line is identified when the number of [ α, Filt _ A, FN ] groups is greater, indicate that the squall line is identified on each of a plurality of radar detection elevation planes, indicating that squall line characteristics are strong, and indicate that the squall line characteristics are more obvious when the FN value is smaller from 0, indicating that squall line characteristics are more obvious on each elevation plane, otherwise, indicating that squall line characteristics are relatively weak.

Reference to the literature

[1] Wangxing, Miaochunsheng, Wangzhong, etc. a method for correcting echo attenuation of a Doppler weather radar [ P ]. Jiangsu: CN105388467A, 2016-03-09.

[2] Liu Cheng Zhi, A manual of conversion of polar coordinates to rectangular coordinates [ M ]. Chinese agricultural machinery Press, 1983.

Claims

1. An intelligent squall line characteristic identification early warning method in a radar echo image is characterized in that: the method comprises the following specific steps:

step 7.2 analysis of each pixel value in RPFBD _ k1(α, x, y) in turn, if soIf the previous pixel value is 0, making the central axis identification value of the pixel be 0, and continuing to analyze the next pixel; if the current pixel value is 1, counting the number of pixels C2 with pixel value of 1 in the 8 pixels immediately adjacent to the current pixel value, if C2 is equal to [3, 4]]If so, the central axis identification value of the pixel is made to be 0, and the next pixel value is continuously analyzed; if it is not

step 7.6, performing multiple rounds of iterative processing on the RPFBD _ k5(α, x, y) according to the method from step 7.1 to step 7.5, taking the RPFBD _ k5(α, x, y) of the previous round as the RPFBD (α, x, y) of step 7.1 in the next round, stopping iteration until no new pixel value is set to be 0, and recording the output result as RPFBD _ kn (α, x, y);

if the RPFBD _ kn (α, x, y) exists in the step 8, recording a data group [ α, Filt _ a, FN ] into the data set Rlt, and entering the step 10, if the squall line does not exist in the step 8, firstly making a counter FN +1 and a high-pass filtering threshold Filt _ a as Filt _ a-FN × Filt _ B, wherein Filt _ B is a constant parameter and takes a positive integer from 1 to 10, then judging the sizes of Filt _ a and the constant parameter Filt _ N, if Filt _ a is not less than or equal to Filt _ N, returning to the step 4, otherwise, if Filt _ a < Filt _ N, entering the step 10, wherein the value range of Filt _ N is 0< Filt _ N < Filt _ a;

2. The intelligent squall line feature identification and early warning method for radar echo images according to claim 1, wherein: the specific steps of the noise reduction treatment in the step 6 are as follows:

step 6.2: if Stat _0> k × Stat _1, setting all pixel values in the square area to 0; otherwise, setting all pixel values in the square area as 1, wherein k is an empirical threshold, and belongs to [0.1,10 ];

and 6.3, sliding the square area to the right or downwards by taking SL as a step size, and repeating the processing procedures of the steps 6.1 and 6.2 until all areas in the RPFB (α, x, y) are processed by the steps 6.1 and 6.2, and recording the processed result as RPFBD (α, x, y).

3. The intelligent squall line feature identification and early warning method for radar echo images according to claim 2, wherein: the specific step of determining whether the squall line exists in step 8 is as follows:

step 8.2, constructing a two-dimensional matrix MAT [ p, q ] by taking omega _ i as a horizontal unit step length and r _ i as a vertical unit step length, wherein p represents a row, p belongs to [ -n · r-i, - (n-1) · r _ i, …, -2 · r _ i, -r _ i,0, r _ i,2r _ i, …, (n-1) · r _ i, n · r _ i ], q represents a column, q belongs to [0, ω _ i,2 ω _ i,3 ω _ i, …,180- ω _ i ], initial values of all units in MAT [ p, q ] are 0, ω _ i is a granularity divided by an included angle between a straight line and a positive direction of an X axis, ω _ i is a positive integer divisible by 180, r _ i is a granularity divided by a shortest distance between the straight line and an origin, n · r _ i and-n · r _ i represent granularities divided by a maximum distance between the straight line and the vertical direction of the longitudinal direction of the X-axis, n · r _ i and a maximum graduation, r _ i is larger than a maximum graduation, rpr is equal to a value of a corresponding integer 365, r _ i, rpr _ i is larger than a diagonal, rpd is equal to a value of a maximum value of a corresponding integer, n is equal to a maximum value, n is equal to;

step 8.3, sequentially traversing each pixel in RPFBD _ kn (α, x, y), skipping and continuing traversing the next pixel if the current pixel value is 0, and if the current pixel value is 1, setting the coordinate of the current pixel as (i, j), namely x ═ i, y ═ j, substituting i into x ', j into y' according to the linear expression in step 8.1, and then substituting q, q ∈ [0, ω _ i,2 ω _ i,3 ω _ i, …,180- ω -i ] in the MAT horizontal direction into ω in the linear expression in step 8.1 one by one, thereby obtaining a series of r (i, j, q);