Morphology-based Doppler meteorological radar multi-speed fuzzy correction method

Abstract

The invention discloses a Doppler meteorological radar multi-speed fuzzy correction method based on morphology, and the method specifically comprises the steps: firstly, extracting an authenticity 0-speed line in a speed PPI image through the morphological operation; secondly, according to the distribution rule of the 0-speed lines, constructing an asymmetry score matrix, and grouping and pairing the 0-speed lines; counting the mean value of the azimuth range of each library-by-library sampling point of each 0-speed line, and eliminating an abnormal path through an IQR (Interquartile Range) method; and finally, real 0-speed line detection and extraction are realized. And taking the obtained real 0-speed line as a reliable initial reference value, and continuously de-blurring according to the speed, so that the multi-speed blurring correction method is completed.

Description

Technical Field

[0001] The present invention relates to the technical field of meteorological radar, and in particular, to a method for correcting multiple velocity ambiguities of a Doppler meteorological radar based on morphology. Background Art

[0002] The accuracy of the radial velocity data observed by a Doppler meteorological radar directly affects the accuracy of wind field inversion, nowcasting, and data assimilation. Short-wavelength radars such as X / Ku band have high detection resolution, but the maximum unambiguous velocity V N is too small, and multiple velocity ambiguity phenomena often occur during the observation of severe weather processes, restricting the further development of the above applications. Therefore, the correction of velocity ambiguity, especially multiple velocity ambiguity, is crucial.

[0003] At present, velocity de-aliasing methods are mainly divided into a hardware de-aliasing method that alternately uses two different pulse repetition frequencies and a software de-aliasing method based on the principle of velocity continuity. The former significantly increases the complexity of data processing and may lead to a decline in data quality. Due to the continuity of the wind field, the core idea of the latter de-aliasing is to find an unambiguous reference point in the adjacent area of the currently processed velocity as the reference velocity, and compare the absolute value of the difference between the current velocity and the reference velocity with a preset threshold. If it is less than the threshold, it is determined that there is no ambiguity at this velocity point; otherwise, it is determined that this point is ambiguous, and de-aliasing is achieved through an integer multiple of ±2V N Therefore, a reliable starting reference velocity is the key to successfully achieving de-aliasing. There are roughly three ways to select the starting reference velocity for traditional software de-aliasing algorithms: (1) Assume that the radial velocity data in several range bins closest to the radar in the starting scan radial direction do not have velocity ambiguity, and use the observed velocities in the first several range bins close to the radar in this radial direction as the starting reference velocity to de-alias along the radial direction; (2) Search for a radial direction with a small and gently changing velocity value in the weakest wind area as the starting reference velocity; (3) When a velocity ambiguity occurs, extract the zero-velocity line based on the observed radial velocity change characteristics, and use the extracted zero-velocity line as the starting reference velocity. The zero-velocity line is a line connected by points with a radial velocity of 0, and the possibility of velocity ambiguity is relatively low, but the extraction process is relatively complex.

[0004] However, when multiple velocity ambiguities occur in the radial velocity data, multiple sets of false zero-velocity lines and true zero-velocity lines caused by ambiguity will be generated. Therefore, when facing the problem of multiple ambiguities, multiple sets of initial radial zero-velocity lines meeting the above conditions will appear, and it is difficult for traditional algorithms to distinguish between false zero-velocity lines and true zero-velocity lines. The incorrect selection of the starting reference velocity will directly affect the incorrect de-aliasing result of the entire elevation data.

[0005] The Doppler radial velocity observed by radar is the projection component of the three-dimensional wind field in the radial direction. The wind component away from the radar direction is the positive Doppler velocity, and the wind component towards the radar direction is the negative Doppler velocity. The radar is located at the center of the image displayed on the Plan Position Indicator (PPI) of the radar radial velocity plane. In the radar surveillance area, when large-scale weather generates vertical wind shear but the wind is basically uniform horizontally, the true radar radial velocity PPI image (without velocity ambiguity) will have the following characteristics. Regarding the true zero-velocity line: The true zero-velocity line passes through the center of the image and is centrosymmetric about the radar position. The positive and negative velocity regions are separated by the true zero-velocity line, clearly distinguishable. During the actual radar scanning process, due to the limitation of the minimum detection distance of the radar, there is a blind area in radar observation. Therefore, a complete zero-velocity line is divided into a group of two zero-velocity lines, and these two zero-velocity lines are centrosymmetric. Regarding the velocity isopleth: In the positive and negative value regions, the isopleths with equal numerical values and opposite directions are centrosymmetric to each other. In the positive (negative) value region, ideally, the isopleth passes through both sides of the boundary of the positive (negative) value region, that is, within the same fixed-distance bin, when sampling a single isopleth, two sampling points with significantly different azimuth angles will be obtained. In actual scanning, if the isopleth does not pass through the radar blind area (meteorological targets in the radar blind area do not have this velocity value), that is, in the PPI image, the radar blind area will not destroy the connectivity of the isopleth. The isopleth will be completely retained, otherwise, the isopleth will be divided into two. At this time, when sampling each part of the broken isopleth in the fixed-distance bin, only a single sampling point can be obtained, that is, the azimuth range of the sampling points approaches zero. When multiple velocity ambiguities occur, the isopleths of ±2nV N will be displayed as false zero-velocity lines. From the above description, it can be known that the false zero-velocity lines mainly have the following characteristics: The false zero-velocity lines have the characteristic of central symmetry; The azimuth range of the sampling points of the first type of false zero-velocity line on the fixed-distance bin is close to 0; The azimuth range of the sampling points of the second type of false zero-velocity line on the fixed-distance bin is extremely large. The true zero-velocity line and the first type of false zero-velocity line have similar shape characteristics. This is because the first type of false zero-velocity line evolves from the isopleth divided into two parts by the radar blind area and is grouped according to the degree of symmetry, thus obtaining two groups of false zero-velocity lines. Therefore, the number of groups of the first type of false zero-velocity line is even and is distributed on both sides of the true zero-velocity line. Therefore, the idea of the algorithm of the present invention for extracting the true zero-velocity line is: First, group the zero-velocity lines according to central symmetry; Second, statistically calculate the azimuth range of the sampling points of each group of zero-velocity lines on the fixed-distance bin, and use the outlier detection method to eliminate the second type of false zero-velocity lines; Finally, according to the position distribution of each group of zero-velocity lines, select the zero-velocity line in the middle position as the true zero-velocity line.

[0006] The key part of the algorithm of the present invention involves the statistical analysis of the spatial shape characteristics of the zero-velocity line, while morphological operations focus on the overall shape and structural characteristics of the object, which is very suitable for extracting the zero-velocity line. Therefore, in view of the multiple blurring problems, the present invention proposes an algorithm for realizing the extraction of the true zero-velocity line by means of morphological operations based on the spatial shape characteristics of the true zero-velocity line and velocity contour lines, and on this basis, realizing successful deblurring. Summary of the Invention

[0007] A method for correcting multiple velocity ambiguities of a Doppler weather radar based on morphology provided by the present invention. The main technical problem to be solved is: aiming at the problem of multiple velocity ambiguities, realizing the identification and extraction of the true zero-velocity line, so as to realize velocity correction.

[0008] The technical solution of the present invention is as follows:

[0009] A method for correcting multiple velocity ambiguities of a Doppler weather radar based on morphology, comprising:

[0010] Step 1, preprocessing the radar radial velocity PPI image;

[0011] Step 2, grouping and pairing the zero-velocity lines according to the central symmetry characteristics of the zero-velocity line;

[0012] Step 3, statistically calculating the mean azimuth range of the sampling points of each group of zero-velocity lines for each bin, and using the interquartile range (IQR) outlier detection method to eliminate the corresponding irrelevant paths;

[0013] Step 4, identifying and extracting the true zero-velocity line;

[0014] Step 5, using the extracted true zero-velocity line as the starting reference velocity to perform deblurring operation.

[0015] Furthermore, the Step 1 includes:

[0016] 1) Image binarization: setting a threshold v x =βV N , where V N is the maximum unambiguous velocity, β∈(0.1, 0.4) is the threshold factor, and at the same time ensuring that v x ≥1 m / s. The threshold is to allow for the position deviation of the zero-velocity line caused by Doppler observation errors, and 1 m / s is the accuracy that can be achieved by current meteorological radar observations of Doppler velocity, which is used as the lower limit of the threshold v x . When the absolute value of the radial velocity is greater than the threshold v x , it is set to False; otherwise, it is set to True, thereby extracting the region close to the zero-velocity and converting it into a binary image.

[0017] 2) Image closing operation processing: at vx = βV N Under the threshold setting of, due to the influence of radar observation accuracy and the degree of wind field nonlinearity, there may be a speed magnitude slightly greater than v in the connected "block". x As a result, holes may appear. The closing operation can fill the small cracks in the connected area, smooth the object boundary and remove isolated noise points at the same time.

[0018] 3) Skeletonization operation: After image binarization and closing operation, what is obtained is the connected "block" area. In order to convert the "block" structure into a "line" structure, skeleton extraction, bridging and pruning operations are performed in sequence.

[0019] In 1), the area where the speed magnitude is less than v x = βV N is set to True. Therefore, the connected area of the "block" structure close to 0 speed is extracted. The skeleton extraction operation is used to remove the edge pixels of the "block" object and extract the skeleton of the connected "block". If the closing operation in 2) does not fill the holes thoroughly, the skeleton fracture problem will occur after the skeleton extraction operation. Therefore, the bridging operation is used to connect the adjacent end points and repair the broken pixels. At different altitudes, the magnitude of the wind speed often varies greatly. Generally, near the ground surface, the wind speed is small. Therefore, in the original observed radial velocity PPI image, the area closer to the radar, that is, the area with lower altitude, has a smaller wind speed. Therefore, the area where the speed magnitude is less than v x is larger, and the width of the connected "block" extracted here is wider. As the altitude increases, the wind speed increases, and the connected "block" area extracted here is narrower. Due to the inconsistent width of the connected "block" before and after and the irregularity of the edge, redundant branches other than the main path will be generated after skeleton extraction. Therefore, the pruning operation is adopted to eliminate the small protrusions on the main path and make the contour object smoother.

[0020] Furthermore, the step 2 includes:

[0021] According to the distribution characteristic that the 0-speed lines are mostly centrosymmetric. Group and pair the main paths. Ideally, whether it is the false 0-speed line caused by velocity ambiguity or the real 0-speed line, it has the characteristic of centrosymmetry about the center of the image (radar position). The main idea of this step is: First, search for the main path through the Depth-First Search (DFS) algorithm. Secondly, by selecting a certain main path, obtain the azimuth distribution of the centrosymmetric path of this main path in each range bin, select the main path closest to its azimuth distribution, and construct an asymmetry score matrix by statistically calculating the mean value of the azimuth angle deviation between the symmetric paths of each main path and other main paths in the sampled range bin. Note that the azimuth angle starts from the due north direction and rotates clockwise.

[0022] 1) The DFS algorithm searches for the backbone path. After skeletonization, there will still be a small number of branches at the end of the main path. To avoid the influence of the end branches on the spatial shape features of the subsequent recognition of the zero-speed line, the DFS algorithm is used to further search for the backbone path. The DFS algorithm is a technique widely used for traversing and searching tree or graph structures. The point closest to the center point (the position of the radar) in each connected component is selected as the root node to search for the longest path, and the longest path found is the backbone path. Thus, the backbone paths of each connected component can be obtained. Connected component combinations with too few pixels (1 / 10 of the number of range bin samples) are filtered out according to the number of valid pixel points. These areas are often isolated noise points or are generated by partial breaks in a certain backbone path. Even if there is valid information, due to the too small number of pixel points, the amount of information provided is still limited.

[0023] 2) Grouping and pairing of zero-speed lines (backbone paths)

[0024] First, arbitrarily select a connected path L m (m = 1, 2,..., l, where l represents the number of paths). The radar position is at the center of the image. According to central symmetry, the symmetric path L m ' of L can be obtained. Compare L m ' with the average angular deviation of other paths L m (k = 1, 2,..., l, k ≠ m) in turn, and record it as the asymmetry score. The specific calculation process is as follows:

[0025] a) Calculate the range bin sampling range where L k ' is located, and calculate the azimuth angle of each range bin sampling point m The subscript m n represents the nth (n = 1, 2,...., N) range bin of the connected path L m . At the same time, record the azimuth angles of each range bin of the path L k . The subscript k n represents the nth range bin of the connected path L k . If there are multiple sampling points in the nth range bin, the average value of the azimuth angles of the sampling points is selected.

[0026] b) Calculate the average angular deviation

[0027] At the same range bin n, the azimuth angle deviation between L m ' and L k is calculated as follows:

[0028]

[0029] ​​Calculate the average of the angular deviations of these N distance libraries, which is the average angular deviation, denoted as for L m , L k 's asymmetry score.

[0030] c) Record for L m The path combination L with the lowest asymmetry score x . Only when for L x , the combination with the lowest asymmetry score is L m , is it considered a successful pairing.

[0031] Furthermore, the said step 3 includes:

[0032] According to the difference in the azimuth range of the sampling points of the first type of false zero-speed line and the second type of false zero-speed line on the fixed distance library, statistically calculate the mean of the azimuth range of the sampling points of each group of zero-speed lines per library, and use the IQR method to detect outliers and eliminate the corresponding irrelevant paths. The specific characteristic differences are as follows: The mean of the azimuth range of the sampling points per library of the first type of false zero-speed line and the true zero-speed line is close to 0, and the mean of the azimuth range of the sampling points per library of the second type of false zero-speed line is relatively large. The IQR method is a method for measuring variability and is used to detect outliers. By calculating the upper and lower boundaries, potential outliers in the data can be identified, and the corresponding backbone paths with extremely large mean azimuth range of sampling points per library can be eliminated.

[0033] The calculation process of the mean of the azimuth range of the sampling points per library of each connected path is as follows:

[0034] 1) First, arbitrarily select a connected path L m (m = 1, 2,..., l, where l represents the number of paths).

[0035] 2) Calculate the sampling range of the distance gate where L m is located, and calculate the azimuth angle of each sampling point of the distance gate n represents the nth distance library, r (r = 1, 2,...., q) represents the sampling point of this path in this distance library. If q ≥ 2, count once, and calculate the difference between the azimuth angles of any two sampling points. The formula is the same as (1), and select the largest azimuth angle difference and denote it as the azimuth angle deviation Δθ of the distance library n n .

[0036] 3) Accumulate the Δθ of all distance libraries n , and divide by the count, then the mean of the azimuth range per library of the connected path L m is obtained The mean of the azimuth range per library of each backbone path is in the range (0, π), which is used as the input data for the IQR method to obtain the backbone path corresponding to the mean of the azimuth range per library with abnormal values .

[0037] The calculation steps of the IQR method are as follows: First, calculate the first quartile Q1: the value at the 25% position after the data is arranged in ascending order; calculate the third quartile Q3: the value at the 75% position after the data is arranged in ascending order. Then IQR = Q3 - Q1. The outlier boundaries are: upper limit = Q3 + 1.5 × IQR, lower limit = Q1 - 1.5 × IQR. 1.5 times IQR is a commonly used threshold in statistics, which is sufficient to identify the zero-speed lines with abnormal per-library azimuth range means of the second type. After screening out the outliers, the main path corresponding to the outlier and its paired path are both removed.

[0038] Furthermore, step 4 includes:

[0039] This step will statistically calculate the overall azimuth angle of the remaining main path groups, and the main path group in the middle position is determined as the true zero-speed line group. Assume that the azimuth angle of a certain sampling point of a certain main path group C (C = 1, 2, 3,..., L, where L represents the total number of remaining main path groups after removing the main paths corresponding to the zero-speed lines of the second type. Note that there are two paired and connected main paths in each main path group.) ( i = 1, 2, 3,..., s c , s c being the total number of sampling points on the main path group C), then the overall azimuth angle φ c of the main path group C is:

[0040]

[0041] The selection principle of the main path in the middle position, which is the true zero-speed line group, is as follows:

[0042] 1) Arbitrarily select a certain main path group D (D = 1, 2, 3,..., L, where L represents the total number of remaining main path groups after removing the main paths corresponding to the zero-speed lines of the second type.), calculate the overall azimuth angle φ D of the main path group D, and the overall azimuth angle φ C of other main path groups C (C = 1, 2, 3,..., L, C ≠ D). The azimuth deviation Δφ DC between them is:

[0043]

[0044] 2) Calculate the overall azimuth deviation Δφ D of the main path group D from other groups as:

[0045]

[0046] 3) Select the overall azimuth deviation Δφ I (I = 1, 2, 3, ..., L) The main path group I with the smallest

[0047] Furthermore, the step 5 includes:

[0048] Through steps 1 - 4, the true 0 - speed line of the velocity - blurred image can be detected, identified, and extracted. Taking the true 0 - speed line as the starting reference velocity for de - blurring operation can effectively solve the problem that the traditional algorithm cannot correct multiple velocity blurs. The core formula for de - blurring is as follows:

[0049] |V ref - V i | ≤ V a (5)

[0050] where V i is the current velocity gate being processed, V ref is the reference velocity within the neighborhood of V i The reference velocity can be the velocity value or relevant statistical value without ambiguity in the radial direction, or in the azimuthal direction along the same range bin, or in the vicinity of the elevation angle, or in the nearest - neighbor region of the adjacent time. V a represents the discontinuous difference between the current processed velocity and the reference velocity. In the de - blurring operation of the UNRAVEL algorithm, V a = aV N , V N is the maximum unambiguous velocity, a is the threshold factor, generally taking values from 0.6 to 0.8. When equation (5) holds, it is determined that there is no velocity blur in this velocity gate; otherwise, it is considered that there is velocity blur in this velocity gate, and the following operations are performed:

[0051] V t = V i + 2nV N (6)

[0052] where n = ±1, ±2, ..., n is an integer, V t is the velocity after de - blurring. By choosing an appropriate n, V t can satisfy equation (5). The above is the core operation of de - blurring. The core of this patent lies in how to extract the true 0 - speed line for multiple - blur problems to provide a reliable starting reference velocity for traditional de - blurring algorithms. The specific de - blurring process can refer to the traditional method of the UNRAVEL algorithm, and of course, it is not limited to the UNRAVEL algorithm.

[0053] Beneficial effects:

[0054] For short-wavelength radars such as X / Ku bands, the Nyquist velocity interval of the observed velocity is relatively small, making it easy to generate multiple ambiguity problems. The key to software de-ambiguity algorithms lies in the reliability of the initial reference velocity. However, in the case of multiple ambiguity problems, there will be multiple false zero-velocity lines and real zero-velocity lines generated by the ambiguity problem. Existing software de-ambiguity algorithms are difficult to determine the real zero-velocity line, greatly reducing the reliability of the initial reference velocity, which will lead to incorrect de-ambiguity results for the entire elevation data. Therefore, traditional software de-ambiguity methods have obvious limitations in multiple velocity ambiguity correction. Based on this, the present invention aims to explore the determination of the real zero-velocity line under multiple ambiguity problems and use the real zero-velocity line as the initial reference value to achieve accurate correction of multiple velocity ambiguity data. The pioneering effects and advantages of the present invention are as follows:

[0055] (1) The present invention applies morphological operations to extract the zero-velocity line from the radial velocity PPI image, simplifying the extraction process and reducing the dependence on the artificially set velocity change characteristics under different weather scenarios;

[0056] (2) The present invention analyzes the differences in the spatial morphological characteristics of true and false zero-velocity lines in the ideal scenario and uses them to accurately identify the real zero-velocity line from among numerous false zero-velocity lines;

[0057] (3) Based on the extraction of the real zero-velocity line, the present invention de-ambiguates, and compared with traditional mainstream de-ambiguity algorithms, it achieves accurate correction of multiple velocity ambiguity data, providing reliable data support for applications such as wind field inversion and data assimilation. BRIEF DESCRIPTION OF THE DRAWINGS

[0058] Figure 1 is the algorithm flow chart of the present invention;

[0059] Figure 2 is a typical illustration of the contour distribution of Doppler velocity ±2nV N under different vertical profiles of horizontal uniform winds;

[0060] Figure 3 is, due to the limitation of the minimum detection distance of the radar, in the actual scan, ±2nV N classification schematic diagram of the contour;

[0061] Figure 4 is the process image of the preprocessing part of the radar radial velocity PPI image in the first embodiment of the present invention;

[0062] Figure 5 is the process image of extracting the real zero-velocity line in the first embodiment of the present invention;

[0063] Figure 6 is the comparison diagram of the de-ambiguity results between the first embodiment of the present invention and traditional mainstream algorithms. DETAILED DESCRIPTION OF THE INVENTION

[0064] In order to make the objectives, technical solutions, etc. of the present invention clearer and more understandable, the present invention will be further described in detail below through specific embodiments in combination with the accompanying drawings. It should be understood that the specific embodiments described herein are only used to explain the present invention and are not used to limit the present invention.

[0065] Embodiment 1:

[0066] In this embodiment, Ku-band meteorological radar experimental observation data is used to verify the invention. The observation experiment was carried out in Ningbo, Zhejiang on November 1, 2024. The main task was to observe Typhoon "Kong-rey". During the experiment, data such as radial velocity and reflectivity obtained from the radar's 9.9° elevation angle observation at 3:02 (UTC) on that day were collected. The original observed radial velocity PPI image can be seen in Figure 4 (a). The range resolution of the radial velocity data is 37.5 m, and 1334 range gates are sampled in each azimuth direction. The maximum unambiguous velocity is approximately 5.17 m / s. Table 1 shows the system parameters of the Ku-band meteorological radar, where "AZ" represents the azimuth direction and "EL" represents the elevation direction.

[0067] Table 1 System Parameters of Ku-band Doppler Meteorological Radar

[0068] Parameter Value Parameter Value Operating wavelength λ (m) 0.0222 Polarization mode Alternating polarization Signal bandwidth (MHz) 4 Pulse repetition frequency (Hz) 1860 Beam width (°, AZ / EL) 1.6 / 2.4 Pulse width (μs) 64 Maximum unambiguous velocity (m / s) 5.17 Elevation angle (°) 9.9

[0069] This Embodiment 1 provides a morphological-based multiple velocity ambiguity correction method for Doppler meteorological radar. The algorithm flowchart can be seen in Figure 1 . Figure 2 Several Doppler velocity diagrams under horizontal uniform winds with different vertical profiles are listed, showing the distribution of ±2nV N isoline and the true 0-velocity line. Among them, Figure 2 The function images in the leftmost column and the bottom row show the vertical profiles of the horizontal uniform wind speed and direction. It is assumed that the radial velocity away from the radar is positive Doppler velocity, represented by a solid line in the figure; the radial velocity towards the radar is negative Doppler velocity, represented by a dashed line; the thick long black dashed line is the true 0-velocity line. The blue solid (dashed) line represents the isoline with a radial velocity value of 2V N , and the orange solid (dashed) line represents the isoline with a radial velocity value of 4V N . The radar is located at the center of the PPI image, and the black dots in the figure represent the radar. From Figure 2 it can be seen that the true 0-velocity line and the isolines with equal and opposite velocity magnitudes are each centrosymmetrically distributed. Figure 3 Shows a schematic diagram of the classification of isolines in an actual scan. Figure 3 (a) is an example of the velocity isoline distribution diagram in an ideal state, corresponding to weather conditions where the wind speed and direction do not change significantly with height. InFigure 3 (b) Under the influence of the radar blind area, a complete isopleth is divided into two parts. Taking the orange 4V N isopleth as an example, it is divided into two parts, green and red, by the radar blind area. If the isopleths that are centrosymmetric to each other are grouped into one group, then the originally regarded group of ±2nV N isopleths are divided into two groups by the radar blind area. Such isopleths whose connectivity is damaged by the radar blind area are the first type of isopleths. Grouping the first type of isopleths according to central symmetry, the number of groups of the first type of isopleths (excluding the true 0-speed line) is even, and they are distributed on both sides of the true 0-speed line. The typical distribution of the second type of isopleths whose connectivity is not damaged by the radar blind area can be seen in Figure 3 (c) and (d). If there is a peak in the wind speed profile within the height range on the display, then a pair of 180° opposite closed isopleths will appear. For this type of isopleth, the connectivity remains unchanged, so there is no situation where an isopleth is divided into two parts. Based on the analysis of the number of intersection points between the above isopleths and the scanning ring at a certain height, it is easy to know that the second type of isopleths often intersect the scanning ring at two points, that is, at a fixed range bin, the second type of isopleths can often sample two positions with a large difference in azimuth angles, which is also the difference from the true 0-speed line and the first type of isopleths. The first type and the second type of ±2nV N isopleths correspond to the first type and the second type of false 0-speed lines. Summarizing the distribution characteristics of the true 0-speed line and the first type and the second type of false 0-speed lines, the following prior information can be obtained: (i) Whether it is the true 0-speed line or the false 0-speed line, it has the characteristics of central symmetry; (ii) At a fixed range bin, on the scanning ring with a radial distance of R, the second type of false 0-speed lines can often sample two positions with different azimuth angles, with the characteristic of a large azimuth deviation between the sampling points; (iii) After removing the second type of false 0-speed lines, the first type of false 0-speed lines are distributed on both sides of the true 0-speed line, and the number of groups is even, and the true 0-speed line is in the middle of the first type of false 0-speed lines. According to the above prior information, the specific steps of the present invention are as follows:

[0070] Step 1 specifically includes:

[0071] 1) Binarize the original observed radial velocity PPI image. Select the threshold v x ≈1.04 m / s, which is 0.2 times the maximum unambiguous velocity of 5.17 m / s. When the absolute value is less than the threshold v x the absolute value is set to True, and thus it is converted into a binary image, and the area close to or showing a 0-speed in the original radial PPI image is extracted. The result of the binarization process can be seen in Figure 4 (b)

[0072] 2) Closed operation is adopted to effectively fill the holes in the zero-speed area. In the closed operation process of the first embodiment, the morphological structure element selected is a circular structure with a radius of 3. The principle for selecting the morphological structure element is to fill the holes in each zero-speed area as much as possible while ensuring that the zero-speed areas corresponding to different zero-speed lines are not connected. For the result of the closed operation, see Figure 4 (c)

[0073] 3) Skeletonization processing: Successively perform skeleton extraction, bridging, and pruning operations. Among them, the skeleton extraction and bridging operations are iterated infinitely to ensure complete refinement to a single-pixel-width skeleton while ensuring that possible breaks are bridged. For the result after the skeleton extraction and bridging processing operations, see Figure 4 (d), and for the result after the pruning operation, see Figure 4 (e).

[0074] For the process diagram of Step 1, see Figure 4 .

[0075] Step 2 specifically includes:

[0076] 1) Use the DFS algorithm to search for the main path, find the longest path, and eliminate the influence of the bifurcation at the end of the main path on the subsequent identification of the spatial position characteristics of the zero-speed line. At the same time, remove the main paths with less than 134 pixels to obtain 8 main paths, numbered ①, ②, ③, ④, ⑤, ⑥, ⑦, ⑧.

[0077] 2) Group and pair the zero-speed lines. Successively select the main path L m (m ∈ 1, 2,..., 8). The radar position is at the center of the image. According to the central symmetry, the symmetric path L m of L m ' can be obtained. Assume that the coordinates of the pixel point i on L m are (x i , y i ), and the position coordinates of the center of the image are (x c , y c ), then the corresponding pixel point on the symmetric path L m ' is (2x c - x i , 2y c - y i ). First, determine the range of the range gate sampling where L m ' is located, and record the azimuth angle subscript m n representing the nth (n = 1, 2,...., N) range bin of the connected path L m . At the same time, record the azimuth angles k of each range gate of the path L According to Equation (1), calculate that at the same range gate n, Lm ' and L k Azimuth angle deviation Obtain the azimuth angle deviations of these N distance bins The average value, i.e., the average angle deviation, is denoted as for L m , L k Asymmetry score. Record for L m The path combination L with the lowest asymmetry score x . Only when for L x , the combination with the lowest asymmetry score is L m , is it considered a successful pairing. Traverse the backbone paths and group and pair all backbone paths. The pairing results can be seen in Figure 5 (a), where the numbers with the same color represent a group, and there are four groups in total, namely ①⑦, ②⑧, ③⑥, and ④⑤.

[0078] Step 3 specifically includes:

[0079] Calculate the mean of the azimuth range for each backbone path In this embodiment, the mean of the azimuth range for each backbone path is: 0.0253, 0.0137, 0.6226, 0.0239, 0.0382, 0.0992, 0.0157, 0.0427. Note that the above data is in radians. The mean of the azimuth range for each backbone path As the input data for the IQR method, the upper and lower bound limits are calculated as: -0.0569, 0.1477. Therefore, backbone path ⑥ and its paired backbone path ③ can be excluded. After excluding backbone paths ③ and ⑥, only backbone paths ①, ②, ④, ⑤, ⑦, and ⑧ remain; among them, backbone paths ① and ⑦ are in group 1, backbone paths ② and ⑧ are in group 2, and backbone paths ④ and ⑤ are in group 3. See Figure 5 (b)

[0080] Step 4 specifically includes:

[0081] According to Equation (2), calculate the overall azimuth angle φ for groups 1, 2, and 3 respectively D (D = 1, 2, 3). Calculate the overall azimuth deviation Δφ of each backbone path group D in turn according to Equation (3) and Equation (4) D , select the group with the smallest overall azimuth deviation, i.e., group 2, and determine it as the true 0-speed line. Thus, the identification and extraction of the true 0-speed line are achieved. See Figure 5 (c), Figure 5 (d) is the extracted true 0-speed line.

[0082] Step 5 specifically includes:

[0083] Based on the determination of the true zero Doppler velocity line, the de-aliasing operation is carried out with this as the initial reference velocity. The de-aliasing result of this example will be compared with the de-aliasing result of the traditional algorithm, the UNRAVEL algorithm. Therefore, the de-aliasing module can refer to the UNRAVEL algorithm and perform initialization (preprocessing, azimuth continuity processing), de-aliasing (range de-aliasing, azimuth de-aliasing, and "Box" de-aliasing), and error checking (range and "Box" error checking) in sequence.

[0084] For the multiple-aliasing problem, based on the extraction of the true zero Doppler velocity line, de-aliasing is performed. The de-aliasing results compared with the traditional method, the UNRAVEL algorithm, and the operational method of the current new generation weather radar, the NEXRAD algorithm, can be seen in Figure 6 . The maximum unambiguous velocity in this scenario is approximately 5.17 m / s. The too small Nyquist velocity interval causes the true radial velocity to fold multiple times, forming a typical multiple-aliasing phenomenon. The original observed velocity PPI image can be seen in Figure 6 (a). The core key to such problems lies in the existence of multiple zero Doppler velocity lines. Accurately identifying the most reliable true zero Doppler velocity line becomes the key to solving the problem. The traditional operational de-aliasing method is based on the following basic assumption: the radial velocities detected in the first several range bins close to the radar do not have aliasing. However, in the multiple-aliasing scenario, this assumption often fails. In addition, this method lacks an effective noise and clutter suppression mechanism, resulting in abnormal phenomena such as local negative (positive) velocities still appearing in large positive (negative) velocity regions as shown in Figure 6 (c). Figure 6 (d) shows the result obtained after de-aliasing by the UNRAVEL algorithm. The UNRAVEL algorithm determines the initial reference range by optimizing the initial reference velocity selection strategy, requiring conditions such as the absolute value of the radial average velocity being the smallest and the single velocity value in the range not being close to the maximum unambiguous velocity. However, in the multiple-aliasing scenario, this method faces two main challenges: First, the existence of multiple zero Doppler velocity lines increases the difficulty of determining the initial reference range; second, the relatively small maximum unambiguous velocity easily leads to the situation where there are velocity aliasing points in a single range. Due to the incorrect selection of the initial reference range, an incorrect velocity discontinuity de-aliasing result appears in a single sector. As shown in the lower right corner of the image in Figure 6 (d). The de-aliasing result after performing the de-aliasing operation based on the present invention can be seen in Figure 6 (b), successfully achieving de-aliasing for the multiple-aliasing problem. The present invention provides a reliable initial reference velocity for the de-aliasing operation, demonstrating excellent de-aliasing capabilities.

[0085] The above content is a further detailed description of the present invention in combination with specific embodiments, and it cannot be determined that the specific implementation of the present invention is only limited to these descriptions. For those of ordinary skill in the technical field to which the present invention pertains, without departing from the concept of the present invention, several simple deductions or substitutions can still be made, and all should be regarded as belonging to the protection scope of the present invention.

Claims

1. A morphological-based multiple velocity ambiguity correction method for Doppler weather radar, characterized in that, Including: Step 1, preprocessing the image displayed on the radar radial velocity Plan Position Indicator (PPI); Step 2, grouping and pairing the zero-velocity lines according to the central symmetry feature of the zero-velocity lines; Step 3, statistically calculating the mean azimuth range of sampling points for each group of zero-velocity lines in each range bin, and using the Interquartile Range (IQR) outlier detection method to eliminate the corresponding irrelevant paths; Step 4, identifying and extracting the true zero-velocity lines; Step 5, using the extracted true zero-velocity lines as the starting reference velocity to perform deblurring operations.

2. The method for multi-velocity ambiguity correction of a Doppler weather radar based on morphology according to claim 1, wherein In the above Step 1, morphological operations such as binarization, closing operation, skeleton extraction, bridging, and pruning are used for image preprocessing.

3. The method for correcting multiple velocity ambiguities of a Doppler weather radar based on morphology according to claim 1, wherein The above Step 2 includes two steps:

1. Searching for the main path by the DFS algorithm: Searching for the main path through the Depth-First Search (DFS) algorithm; 2. Grouping and pairing of zero-velocity lines (main paths): By selecting a certain main path, obtaining the azimuth angle distribution of the central symmetric path of this main path in each range bin, selecting the main path with the closest azimuth angle distribution to it, and constructing an asymmetry score matrix by statistically calculating the mean azimuth angle deviation between the symmetric paths of each main path and other main paths in the sampling range bins.

4. The method for multi-velocity ambiguity correction of a Doppler weather radar based on morphology according to claim 3, characterized in that, In the above Step 2, the method for grouping and pairing zero-velocity lines (main paths) is as follows: Arbitrarily select a connected path L m (m = 1, 2, ..., l, where l represents the number of paths); the radar position is at the center of the image. According to central symmetry, L m 's symmetric path L m '; successively compare the average angular deviation between L m ' and other paths L k (k = 1, 2, ..., l, k ≠ m), and record it as the asymmetry score; the specific calculation process is as follows: 1) Calculate L m The distance gate sampling range where' is located, and calculate the azimuth angle of each distance gate sampling point Subscript m n Represents the nth (n = 1, 2,...., N) range bin of the connected path L m ; At the same time, record the azimuth angles of each distance gate of path L k Subscript k n Represents the nth range bin of the connected path L k ​​ If there are multiple sampling points in the nth range bin, then select the mean of the azimuth angles of the sampling points; 2) Calculate the average angle deviation At the same distance from door n, L m ' and L k The azimuth angle deviation The calculation formula is as follows: Calculate the average of the angular deviations of these N distance libraries, which is the average angular deviation, denoted as for L m , L k 's asymmetry score; 3) Record for L m The path combination L with the lowest asymmetry score x ; Only when for L x the combination with the lowest asymmetry score is L m is it considered a successful pairing.

5. The method for multi-velocity ambiguity correction of a Doppler weather radar based on morphology according to claim 1, wherein In the above Step 3, the calculation process of the mean azimuth range of sampling points for each connected path in each range bin is as follows: 1) First, arbitrarily select a connected path L m (m = 1, 2,..., l, where l represents the number of paths); 2) Calculate L m Determine the distance gate sampling range where it is located, and calculate the azimuth angle of each distance gate sampling point n represents the nth range bin, and r (r = 1, 2,...., q) represents the sampling point of this path in this range bin. If q ≥ 2, count once, and calculate the difference between the azimuth angles of any two sampling points. The formula is the same as (1). Select the largest azimuth angle difference and denote it as the azimuth angle deviation Δθ of this range bin n n ; 3) Accumulate the Δθ for all distance libraries, and divide by the count to obtain the connected path L n The average azimuth range difference per library m of the connected path L The average azimuth range difference per library of each main path is in the range (0, π), which is used as the input data for the IQR method to obtain the average azimuth range difference per library with abnormal values and the corresponding main path 6. The method for multi-velocity ambiguity correction of a Doppler weather radar based on morphology according to claim 1, wherein In the said step 4, the overall azimuth angle φ of the main path group C c is as follows: Among them is a certain pixel point c of the main path group C i azimuth angle i ∈ 1, 2, 3,..., s c , s c is the total number of sampling points on the main path C of the group 7. The method for correcting multiple velocity ambiguities of a Doppler weather radar based on morphology according to claim 1, wherein In the above Step 4, the selection principle for the main path at the middle position, which is the true zero-velocity line group, is as follows: 1) Arbitrarily select a certain main path group D (D = 1, 2, 3,..., L, where L represents the total number of remaining main path groups after excluding the main paths corresponding to the second - type 0 - speed lines), and calculate the overall azimuth angle φ of the main path group D D , and the overall azimuth angle φ of other main path groups C (C = 1, 2, 3,..., L, C ≠ D) C . The azimuth deviation Δφ DC is as follows: 2) Calculate the overall azimuth deviation Δφ of the main path group D from other groups D It is: 3) Select the overall azimuth deviation Δφ I (I = 1, 2, 3, ..., L) The main path group I with the smallest

Images