Three-dimensional storm cell identification and boundary refinement method and system based on multi-source data

Abstract

The application provides a three-dimensional storm cell identification and boundary refinement method and system based on multi-source data, which comprises the following steps: obtaining multi-source radar data and fusing the data into three-dimensional radar grid data; selecting candidate seed points in the fused data, performing spatial clustering analysis and aggregation into initial growth cores, and performing three-dimensional competitive growth according to multi-constraint growth criteria with the initial growth cores as starting points; when the distance between the growth cores is less than a threshold, triggering a competitive segmentation mechanism to realize adjacent cell decoupling and separation, and outputting an initial three-dimensional region label graph; taking the growth boundary in the initial region label graph as an initial contour, calculating a multi-scale three-dimensional morphological gradient field, defining an energy functional of a three-dimensional dynamic contour model, and performing three-dimensional curve evolution by using a level set method to output a three-dimensional contour of the storm cell. Through integrated processing of fusion, growth and refinement, the application improves the accuracy of storm cell identification and the degree of boundary refinement.

Description

Technical Field

[0001] This invention relates to the field of meteorological radar data processing and severe convective weather monitoring technology, and in particular to a method and system for three-dimensional storm cell identification and boundary refinement based on multi-source data. Background Technology

[0002] In radar detection technology, meteorological radars are mainly classified into S-band (wavelength approximately 10 cm), C-band (wavelength approximately 5 cm), and X-band (wavelength approximately 3 cm) based on their operating wavebands. S-band radars have low electromagnetic wave attenuation and a large detection range, but their resolution is relatively low. X-band radars have high resolution and can detect the fine structure of storms, but their electromagnetic wave attenuation is severe. Dual-polarization radar technology can simultaneously transmit and receive electromagnetic waves in both horizontal and vertical polarization directions, acquiring multi-parameter data such as reflectivity factor (Z), differential reflectivity (ZDR), and correlation coefficient (CC), providing richer information for the microphysical characteristic analysis of precipitation particles. Multi-radar network observation has become a trend in meteorological operations. How to effectively integrate radar observation data from different wavebands, balancing large-area coverage with local high-resolution capabilities, is one of the core technical challenges in the application of multi-source radar data.

[0003] In terms of storm identification algorithms, representative technical solutions include the TITAN algorithm (Thunderstorm Identification, Tracking, Analysis, and Nowcasting algorithm) and the SCIT algorithm (Storm Cell Identification and Tracking algorithm). The TITAN algorithm sets a reflectivity threshold to filter voxels, performs three-dimensional connected component analysis to aggregate them into storm clusters, and then performs secondary segmentation along the local minimum reflectivity locations to separate cohesive cells. The SCIT algorithm searches for reflectivity segments exceeding a threshold along the radar radial direction, merges them into two-dimensional components in the azimuth direction, and then vertically links them between different elevation layers to form three-dimensional storm cells. However, the above-mentioned existing technical solutions still have the following shortcomings.

[0004] First, the accuracy of multi-source radar data fusion is insufficient. Existing technologies often employ simple maximum value mosaicking or single range weighting methods when fusing data from multiple radar sources. These methods treat multi-source radar data fusion merely as a spatial interpolation problem, failing to consider the fundamental differences in attenuation characteristics and resolution between the S-band and X-band. This results in fusion results that struggle to simultaneously meet the requirements of wide-area coverage and local high resolution.

[0005] Secondly, the physical constraints of three-dimensional spatial growth and segmentation are weak. Existing methods mainly rely on a single reflectivity threshold for determining storm boundaries and cell segmentation, treating storm boundaries as isosurfaces determined by reflectivity differences. However, from a storm dynamics perspective, the natural boundary of a storm cell is essentially a convergence-divergence conversion surface of the airflow field, manifesting as a velocity shear extremum region in the radar radial velocity field. Existing technologies do not utilize this physical fact, making it difficult to effectively separate adjacent cells in dense storm areas, and also making it difficult to distinguish between meteorological echoes and non-meteorological echoes, resulting in false targets in the identification results.

[0006] Third, the storm boundaries are coarse and lack three-dimensional morphological refinement methods. Existing methods directly output the voxel-level labeling results as the storm boundaries after completing connected domain analysis or region growing, lacking a post-processing step to transform the coarse boundaries into smooth, continuous three-dimensional surfaces. This results in jagged boundaries, which directly affects the accuracy of subsequent quantitative analyses such as volume calculations and vertical integration liquid water content estimation.

[0007] In summary, existing technologies have significant shortcomings in terms of the physical adaptability of data fusion, the dynamic constraints of growth segmentation, and the geometric smoothness of boundary output. Summary of the Invention

[0008] This invention addresses the technical problems existing in the prior art by providing a method and system for three-dimensional storm cell identification and boundary refinement based on multi-source data. It aims to achieve high-precision, fully automated identification and three-dimensional contour refinement of storm cells by constructing an integrated processing chain of "adaptive fusion of multi-source data - three-dimensional competitive growth - level set boundary refinement", thereby solving the shortcomings of the prior art in terms of data fusion, physical segmentation and boundary smoothness.

[0009] To solve the above-mentioned technical problems, the technical solution proposed by this invention is as follows: A method for identifying and refining the boundaries of 3D storm cells based on multi-source data includes: In the multi-source radar data 3D mesh fusion processing stage, multi-source radar data is acquired and fused into 3D radar mesh data; In the three-dimensional competitive growth stage guided by convergence boundary constraints and physical phase state, local maxima with reflectivity factors exceeding a threshold are searched in the fused three-dimensional radar grid data as candidate seed points. Spatial clustering analysis is performed on these candidate seed points, and those satisfying the spatial distance threshold and reflectivity gradient direction consistency are aggregated into an initial growth kernel. Starting from the initial growth kernel, candidate voxels are expanded outward layer by layer according to multi-constraint growth criteria for three-dimensional competitive growth. When the distance between growth kernels is less than a set threshold, a competitive segmentation mechanism is triggered to achieve automatic decoupling and separation of adjacent storm cells, and the initial three-dimensional region label map of each storm cell is output. The multi-constraint growth criteria include reflectivity continuity constraints, dual polarization phase state constraints, and radial velocity convergence boundary termination constraints. In the multi-scale morphological gradient and dynamic contour model 3D boundary refinement stage, the growth boundary in the initial 3D region label map is used as the initial contour. The multi-scale 3D morphological gradient field is calculated in the boundary neighborhood. The energy functional of the 3D dynamic contour model is defined. The energy functional includes external energy terms and internal energy terms. The level set method is used to perform 3D surface evolution on the energy functional. The iteration continues until the energy functional converges, and the 3D contour of the storm cell is output.

[0010] As a further improvement to the method of the present invention: the multi-source radar data three-dimensional mesh fusion processing stage includes: Acquire volume scan observation data from at least one S-band weather radar and at least one X-band weather radar, select an S-band radar that meets the preset requirements as the reference radar, and establish a three-dimensional polar coordinate reference grid with the reference radar site as the origin. The volume scan observation data of other S-band and X-band weather radars are mapped to the three-dimensional polar coordinate reference grid to obtain multi-source radar observation data with the three-dimensional polar coordinate reference grid as a unified spatial framework. An adaptive weighting function is constructed to perform weighted fusion of volume scan observation data from multiple weather radars at each grid point, resulting in the fusion results of volume scan observation data from S-band weather radar and X-band weather radar, respectively. The volume scan observation data fusion results of the S-band weather radar and the volume scan observation data fusion results of the X-band weather radar are fused in the frequency domain using a frequency-band synthesis strategy, and the fused three-dimensional fused radar grid data is output.

[0011] As a further improvement to the method of the present invention: the adaptive weight function is jointly determined by a distance dependence factor, a physical confidence factor, and a beam obstruction perception factor; the distance dependence factor is used to adjust the weight according to the distance between the weather radar and the observed target, the physical confidence factor is used to adjust the weight according to the measurement reliability of different physical quantities on radars of different bands, and the beam obstruction perception factor is used to adjust the weight according to the terrain obstruction.

[0012] As a further improvement to the method of the present invention: the frequency-band synthesis strategy includes: performing three-dimensional Gaussian pyramid decomposition on the fusion results of volume scan observations from the S-band weather radar and the X-band weather radar to obtain Laplace components at each level; performing selective weighted synthesis according to the spatial frequency level, with the high-frequency layer dominated by the X-band Laplace component and the low-frequency layer dominated by the S-band Laplace component; reconstructing the Laplace components of each layer after weighted synthesis layer by layer through the pyramid to output the fused three-dimensional mesh data after frequency-band synthesis.

[0013] As a further improvement to the method of the present invention: in the three-dimensional competitive growth stage guided by convergence boundary constraints and physical phase state, the method for judging the consistency of the reflectivity gradient direction is as follows: calculate two candidate seed points. and The angle between the three-dimensional gradient directions of reflectivity When the included angle When the reflectivity is less than the threshold, it is determined that the reflectivity gradient direction is consistent; the included angle The calculation formula is:

[0014] in, and Representing candidate seed points and The reflectance three-dimensional gradient vector at a given location, where || represents the magnitude of the vector; the reflectance three-dimensional gradient vector It is composed of the partial derivatives of reflectivity along the x, y, and z directions, and is calculated using the central difference method. Its three components are as follows:

[0015]

[0016]

[0017] in, voxels The reflectivity factor value at that location, Seed point spatial coordinates, , , The three-dimensional radar mesh is respectively in Spatial resolution in the direction.

[0018] As a further improvement to the method of the present invention: in the three-dimensional competitive growth stage guided by convergence boundary constraints and physical phase states, the dual-polarization phase state constraints include: extracting the multi-polarization feature set corresponding to the candidate growth voxels encountered when the growth nucleus expands outward; inputting the multi-polarization feature set into a fuzzy logic classifier, which calculates the comprehensive membership score of each phase state type based on the membership function corresponding to each phase state type and the preset feature weights, and taking the phase state with the highest score as the classification result; when the classification result corresponding to the candidate growth voxel is the preset phase state, the growth of the candidate growth voxel is prohibited; the comprehensive membership score The calculation formula is:

[0019] in, For the first Each input feature For feature weights, phase state Features The membership function.

[0020] As a further improvement to the method of the present invention: in the three-dimensional competitive growth stage guided by the convergence boundary constraint and physical phase state, the radial velocity convergence boundary termination constraint includes: for each candidate voxel extending one layer beyond the boundary of the growth body. In the candidate voxels Calculate the candidate voxel in the azimuth direction Radial velocity shear variable: When the absolute value of the radial velocity shear variable is greater than a set threshold, the candidate voxel is determined to be a convergent boundary voxel, the convergent boundary voxel is prohibited from being included in the growth region, and the growth front terminates its expansion in the direction of the convergent boundary voxel; the radial velocity shear variable... The calculation formula is:

[0021] in, Radial distance, Angle of elevation, It is the azimuth angle. voxels The fusion radial velocity, The radial distance where the voxel is located. This is the azimuth resolution.

[0022] As a further improvement to the method of the present invention: in the three-dimensional competitive growth stage guided by the convergence boundary constraint and physical phase state, the competitive segmentation mechanism includes: when the distance between growth nuclei is less than a set threshold, a conflict region is determined; in the conflict region, the weighted distance from voxels to each growth nucleus is calculated; the weighted distance takes the minimum value at the location with a large reflectivity gradient and the maximum value at the location with a small reflectivity gradient; a three-dimensional distance difference field is constructed; and the isosurface with zero difference field is used as the watershed segmentation surface to divide the voxels in the conflict region into the corresponding growth regions.

[0023] As a further improvement to the method of the present invention: in the stage of refining the three-dimensional boundary of the multi-scale morphological gradient and dynamic contour model, the calculation method of the multi-scale three-dimensional morphological gradient field includes: Centered on the boundary surface of each storm cell, calculate the Euclidean distance from each voxel in the three-dimensional radar grid to the boundary surface, and add voxels with a distance of less than 5 voxel spaces to the narrowband region. Within the narrow band region, spherical structure elements with radius scales of 3, 5, and 7 voxel distances were used to calculate the difference between the three-dimensional expansion and three-dimensional erosion of each voxel at each scale, thus obtaining the morphological gradient at each scale. :

[0024] in, A spherical structural element with scale s. Represents the reflectivity factor, and ⊕ represents the three-dimensional dilation operation. This represents a three-dimensional erosion operation, where the three-dimensional dilation is the maximum value of the reflectivity factor within a spherical neighborhood, and the three-dimensional erosion is the minimum value of the reflectivity factor within a spherical neighborhood. The morphological gradients at each scale are weighted and averaged to fuse, with the small-scale gradient having the largest weight and the large-scale gradient having the smallest weight. The fusion result is normalized to the [0,1] interval to obtain the final gradient field.

[0025] The present invention also provides a three-dimensional storm cell identification and boundary refinement system based on multi-source data, including a microprocessor and a memory interconnected thereto, wherein the microprocessor is programmed or configured to execute the three-dimensional storm cell identification and boundary refinement method based on multi-source data.

[0026] Compared with the prior art, the beneficial effects of the present invention are as follows: This invention automates the entire process of storm cell identification by constructing an integrated processing chain of "fusion-growth-refinement". At the data fusion level, it combines the advantages of wide coverage of S-band radar and high spatial resolution of X-band radar through three-dimensional mesh fusion processing of multi-source radar data, providing a high-quality data foundation for subsequent identification. At the physical segmentation level, by introducing radial velocity convergence boundary termination constraints and dual polarization phase constraints, it endows the growth process with physical termination conditions and guidance mechanisms from the perspectives of fluid dynamics and microphysics, respectively, effectively solving the problems of easy adhesion in dense storm areas and inability to distinguish non-meteorological echoes in traditional reflectivity threshold methods. At the boundary refinement level, through the construction of multi-scale morphological gradient fields and the solution of dynamic contour models of level sets, the initial voxel-level rough boundary is refined into a smooth and continuous three-dimensional surface, so that the final boundary converges to the physical boundary position where the reflectivity space changes most drastically, filling the technical gap between the voxel level and the smooth physical boundary in the existing technology. The synergistic effect of these three levels ultimately enabled fully automated, high-precision identification of storm cells, from data fusion and three-dimensional competitive growth to the refinement of their three-dimensional boundaries. This significantly improved the accuracy, physical rationality, and boundary refinement of the identification, providing reliable technical support for the monitoring and early warning of severe convective weather. Attached Figure Description

[0027] Figure 1 This is a flowchart of the method for identifying and refining the boundaries of a three-dimensional storm cell based on multi-source data in an embodiment of the present invention.

[0028] Figure 2 This is a detailed flowchart of the algorithm for multi-source radar three-dimensional mesh data fusion in an embodiment of the present invention.

[0029] Figure 3 This is a detailed flowchart of the algorithm for the convergence boundary termination mechanism in three-dimensional competitive growth in this embodiment of the invention.

[0030] Figure 4 This is the refinement process of the three-dimensional boundary refinement algorithm for multi-scale morphological gradient and dynamic contour model in the embodiments of the present invention. Detailed Implementation

[0031] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present invention, and not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those of ordinary skill in the art without creative effort are within the scope of protection of the present invention.

[0032] like Figure 1 As shown, this embodiment provides a method for three-dimensional storm cell identification and boundary refinement based on multi-source data, including: Step 1: Multi-source radar data 3D mesh fusion processing stage, acquire multi-source radar data, and fuse the multi-source radar data into 3D radar mesh data.

[0033] In this embodiment, the three-dimensional mesh fusion processing stage of multi-source radar data includes: Step 101: Acquire volume scan observation data from at least one S-band weather radar and at least one X-band weather radar, select an S-band radar that meets the preset requirements as the reference radar, and establish a three-dimensional polar coordinate reference grid with the reference radar site as the origin.

[0034] like Figure 2 As shown, this embodiment acquires volume scan observation data from N S-band dual-polarization weather radars and M X-band dual-polarization weather radars, where N and M are both greater than or equal to one. Radar data includes four physical quantities: reflectivity factor (denoted by the symbol Z, proportional to the sixth power of precipitation particle number concentration and particle size, measured in dBZ; higher values ​​indicate stronger precipitation. Generally, above 35 dBZ is considered convective precipitation), radial velocity (denoted by the symbol V, measured in m / s, reflecting the speed of meteorological targets along the radar radial direction), differential reflectivity (denoted by ZDR, the difference between horizontal and vertical polarization reflectivity, measured in dB; positive values ​​indicate that the particle size is larger in the horizontal direction than in the vertical direction (e.g., large raindrops are oblate spheroids), while hail, due to its tumbling, is nearly spherical and has a ZDR close to 0 dB), and correlation coefficient (denoted by CC, the degree of correlation between horizontal and vertical polarization echo signals, dimensionless, ranging from 0 to 1. Pure precipitation areas are usually greater than 0.95, while ground clutter and bio-echoes are usually less than 0.85).

[0035] When multiple S-band radars exist, one is selected as the reference radar to establish a three-dimensional polar coordinate grid. The selection of the reference radar comprehensively considers the effective detection coverage volume, the geographical center location relative to the observation network, data quality, and calibration accuracy. Typically, the S-band radar located at the geometric center of the observation network, possessing the largest effective detection coverage volume, and having dual-polarization observation capabilities is preferred. After selecting the reference S-band radar, a three-dimensional polar coordinate reference grid is established with the radar site as the origin. This grid consists of three dimensions: radial distance r, azimuth angle θ, and elevation angle φ. The observation data from the reference S-band radar directly correspond to the grid coordinates.

[0036] Step 102: Map the volume scan observation data of other S-band and X-band weather radars to a three-dimensional polar coordinate reference grid to obtain multi-source radar observation data with the three-dimensional polar coordinate reference grid as a unified spatial framework.

[0037] Specifically, for all non-reference radars (including non-reference S-band radars and all X-band radars), their observation data are mapped to a three-dimensional polar coordinate reference grid. That is, for any grid point P in the three-dimensional polar coordinate reference grid, the corresponding sampling position of each non-reference radar at that grid point is calculated. The observation values ​​of each non-reference radar at grid point P are obtained through spatial coordinate transformation and bilinear interpolation, thus ensuring that both non-reference radars and reference radars have corresponding observation data on the same set of three-dimensional polar coordinate grids. To achieve the above mapping, the spatial rectangular coordinate system of grid point P in the atmosphere is first calculated. The specific calculation formula is as follows:

[0038]

[0039]

[0040] in, The radial distance is the location of grid point P in the reference grid. Let P be the elevation angle at grid point P. Let P be the azimuth angle of the grid point P. For the equivalent Earth radius, This item is for Earth curvature correction.

[0041] Secondly, calculate the first Non-reference radar observation of this grid point The formulas for slant range, azimuth, and elevation are:

[0042]

[0043]

[0044] in,( , , ) is the first Spatial position coordinates of non-reference radar The coordinates are in the east-west direction. The coordinates are in the north-south direction. This refers to altitude. The reflectivity factor needs to be converted to linear units first. ,in To convert the reflectivity factor value to linear units (where z is the reflectivity factor value), interpolation is performed, and the interpolation result needs to be converted back to dBz units. Radial velocity, differential reflectivity, and correlation coefficient can be directly interpolated using bilinear interpolation. Through the above mapping, both non-reference radar and reference radar have corresponding observation data on the same set of three-dimensional polar coordinate grids, resulting in multi-source radar observation data with a unified spatial framework of a three-dimensional polar coordinate reference grid.

[0045] To achieve effective fusion of observation data from multiple radars, it is necessary to determine the fusion weight of each radar at each grid point. Taking the reflectivity factor as an example, if a direct weighted summation method is used, the fused reflectivity factor is theoretically equal to the sum of the weighted contributions of all S-band radars and all X-band radars at that grid point, i.e.: (7) in, The fusion weights for the nth S-band data are... Let p be the observed reflectivity factor of the nth S-band radar. The fusion weights for the m-th X-band radar data are: Let be the observed reflectivity factor of the m-th X-band radar at point p. And it satisfies the normalization constraint: .

[0046] However, in order to further utilize the frequency-band synthesis strategy to achieve complementary advantages between the S-band and X-band in subsequent steps, this embodiment does not directly perform one-time fusion according to the above formula. Instead, it first calculates the fusion results of the S-band radar and the X-band radar separately, and then performs frequency domain fusion through the frequency-band synthesis strategy. For this purpose, it is necessary to determine the comprehensive weight of each radar at each grid point.

[0047] Step 103: Construct an adaptive weighting function to perform weighted fusion of volume scan observation data from multiple weather radars at each grid point, and obtain the fusion results of volume scan observation data from S-band weather radar and X-band weather radar respectively.

[0048] Taking reflectivity factor as an example, the fused reflectivity factor Theoretically, it equals the sum of the weighted contributions of all S-band radars at that grid point plus the sum of the weighted contributions of all X-band radars at that grid point. To achieve optimal fusion, it is necessary to calculate the combined weight of each S-band radar and each X-band radar at grid point P.

[0049] In this embodiment, the adaptive weighting function is jointly determined by three parts: distance dependence factor, physical confidence factor, and beam occlusion perception factor.

[0050] (1) Distance-dependent factor

[0051] For the m-th X-band radar, a range-dependent weighting factor is defined. This is a function of the distance *d* between the radar and grid point P. The principle behind this function is that X-band radar has high data quality at close range and should be assigned a higher weight; as the distance increases, X-band electromagnetic waves attenuate significantly, leading to a decrease in data reliability, so the weight should gradually shift towards S-band data. The function's form is: ={ ,

[0052] ,

[0053] ,

[0054] } in, The value is the attenuation coefficient of X-band data due to distance, then 1- The confidence level to be retained for X-band data. The function achieves a linear and smooth transition of confidence from X-band to S-band between 15 and 40 kilometers.

[0055] For the nth S-band radar, its electromagnetic wave attenuation is small, and its effective detection range is long. A range-dependent weighting factor is defined. The function representing the distance d between the radar and grid point P is as follows: ={ ,

[0056] ,

[0057] ,

[0058] } In this context, the distance d of the reference radar is directly equal to the radial distance of the grid point, while the distance d of the non-reference S-band radar is equal to the Euclidean spatial distance of the radar to point p. Although the distances to the same target grid point p are different, the attenuation law is the same.

[0059] (2) Physical confidence factor

[0060] Different physical quantities exhibit varying degrees of dependence on different wavelengths. The reflectivity factor is more accurately measured in the S-band due to less attenuation, while differential reflectivity and correlation coefficient provide a more refined characterization of particle physics in the X-band due to higher spatial resolution. A confidence adjustment coefficient for the physical quantity is introduced. As shown in Table 1.

[0061] Table 1 Confidence Adjustment Coefficients for Physical Quantities

[0062] (3) Beam blocking sensing factor γ Using digital elevation model data in advance, the beam obstruction rate of each radar in each azimuth and elevation direction is calculated. This involves using publicly available terrain height grid data to compare the beam center height with the terrain height point-by-point along the beam propagation path, determining the proportion of beam cross-section truncated by the terrain. For radars including all radars in the azimuth direction... Angle of elevation Beam blocking rate in the direction Define the occlusion perception factor: ={ , (Take 60%) (Take 60%) } When the occlusion rate Exceeding the threshold hour, The fusion weight of all grid points in this direction of the radar is reset to 0.

[0063] Based on the above three factors, the original weight of the q-th radar (non-reference radar or X-band radar) at grid point P is calculated. for: (8) in, This represents the attenuation coefficient due to distance. The confidence level of the physical factor. Where P is located i Obstruction rate at azimuth k-elevation angle.

[0064] The original weights of the benchmark S-band radar are : (9) The distance from the reference radar to point P is the radial distance. Therefore Its distance attenuation coefficient, For the confidence level of the physical factor, Where the occlusion rate P is located i Obstruction rate at azimuth k-elevation angle.

[0065] The original integrated weights of all radars (reference S-band radar, N-1 non-reference S-band radars, and M X-band radars) at grid point P for:

[0066] in, This represents the original weight of the nth S-band non-reference radar at point P. Let represent the original weight of the m-th X-band radar at point P. The original weight is obtained by formula (8).

[0067] After normalization, the overall weight of each S-band radar is obtained as follows: The overall weight of each X-band radar is :

[0068]

[0069] Step 104: Using the above-mentioned normalized comprehensive weights, the volume scan observation data fusion results of the S-band weather radar and the volume scan observation data fusion results of the X-band weather radar are fused in the frequency domain using a frequency-band synthesis strategy, and the fused three-dimensional fused radar grid data is output.

[0070] In this embodiment, the frequency-band synthesis strategy includes: performing three-dimensional Gaussian pyramid decomposition on the fusion results of volume scan observations from S-band weather radar and X-band weather radar to obtain Laplace components at each level; selectively weighting synthesis according to the spatial frequency level, with the high-frequency layer weighted by the X-band Laplace component and the low-frequency layer weighted by the S-band Laplace component; and reconstructing the Laplace components of each layer after weighted synthesis layer by layer through the pyramid to output the fused three-dimensional mesh data after frequency-band synthesis.

[0071] Taking the reflectivity factor as an example, after obtaining the fusion weight of each radar at each grid point through formulas (11) and (12), the fused reflectivity observations of N S-band radars at point P are calculated respectively. Combined reflectivity observations at point P from the M-band X-band radar :

[0072]

[0073] At this point, adaptive weighted fusion was completed, yielding S-band and X-band fusion results. The next step is frequency-band synthesis. and A 3D Gaussian pyramid decomposition is performed separately. This embodiment sets up a 4-level pyramid, with level 0 representing the original resolution.

[0074]

[0075] Next, downsampling is performed step by step. Operation, which includes performing The process involves two steps: 3D Gaussian kernel convolution and 3D downsampling with a stride of 2. Each step of downsampling halves the 3D resolution.

[0076]

[0077] in, This indicates that the S-band fusion result is in the [number]th [number]. Downsampling results of the pyramid (i.e., after...) (resolution hierarchical data after subsampling) This indicates that the X-band fusion result is in the [number]th [band]. Downsampling results of the pyramid.

[0078] Then, progressively upgraded sampling was performed. That is, the inverse process of downsampling mentioned above. Upsampling involves inserting zero values ​​between the existing voxels (the smallest volume unit in three-dimensional space, the corresponding concept of two-dimensional pixels in three-dimensional space) before performing upsampling. 3D Gaussian kernel convolution smoothing aims to lose high-frequency details (including fine textures, edges, and local abrupt changes), that is:

[0079]

[0080] in, This indicates that the S-band fusion result is in the [number]th [number]. Upsampling results of the pyramid level This indicates that the X-band fusion result is in the [number]th [band]. Upsampling results of the pyramid.

[0081] After obtaining the downsampling and upsampling results, calculate the Laplace component for each level, that is, take the difference between the downsampling and upsampling results to obtain the high-frequency residual, which is the lost high-frequency details:

[0082]

[0083] To more intuitively illustrate the calculation process of each layer of Laplace components, the overall formula is simplified as follows:

[0084] After the pyramid decomposition was completed, the Laplace components of the S-band and X-band radars at each level were obtained. and Subsequently, selective frequency band synthesis is performed, and frequency band allocation thresholds are set. =2. The weighted Plas components at each level are determined by the following formula: ={

[0085] ,

[0086] } in, The value ranges from 0.75 to 0.92, with an empirically optimized value of 0.85, indicating that the X-band contributes 85% in the high-frequency layer. The value ranges from 0.05 to 0.30, with an empirically preferred value of 0.20, indicating that the X-band accounts for 20% in the low-frequency layer. The above parameter values ​​can be adjusted according to actual data conditions and application requirements.

[0087] After obtaining the weighted Laplacian components for each level, pyramid reconstruction is performed. Reconstruction is carried out layer by layer, starting with the lowest resolution layer:

[0088]

[0089]

[0090]

[0091] This is the final fused three-dimensional mesh data after frequency band synthesis. This data retains the large-scale reliability of the S-band at low spatial frequencies and absorbs the fine structural information of the X-band at high spatial frequencies, achieving the complementary advantages of the two bands.

[0092] The above fusion process is performed on reflectivity factor, radial velocity, differential reflectivity, and correlation coefficient respectively, and three-dimensional fused radar grid data containing these four physical quantities is output for subsequent storm cell identification.

[0093] Alternatively, in addition to the frequency-band synthesis strategy based on Gaussian pyramid decomposition mentioned above, a multi-scale fusion method based on wavelet transform can be used as an alternative. Specifically, three-dimensional wavelet decomposition is performed on the S-band and X-band fused data respectively, and the fusion weights are adaptively selected according to the spatial frequency at each wavelet coefficient layer. After reconstruction by inverse wavelet transform, the complementary advantages of S-band and X-band in the frequency domain can also be achieved.

[0094] Step 2: The 3D competitive growth stage guided by convergence boundary constraints and physical phase states. Local maxima points with reflectivity factors exceeding a threshold are searched in the fused 3D radar grid data as candidate seed points. Spatial clustering analysis is performed on these candidate seed points, aggregating those satisfying the spatial distance threshold and reflectivity gradient direction consistency into initial growth kernels. Starting from these initial growth kernels, candidate voxels are expanded layer by layer outward according to multi-constraint growth criteria, performing 3D competitive growth. When the distance between growth kernels is less than a set threshold, a competitive segmentation mechanism is triggered, achieving automatic decoupling and separation of adjacent storm cells, and outputting the initial 3D region label map for each storm cell. The multi-constraint growth criteria include reflectivity continuity constraints, dual-polarization phase state constraints, and radial velocity convergence boundary termination constraints. The physical basis of the radial velocity convergence boundary termination constraint is that in the radar radial velocity field, the region where positive velocity (indicating the target is moving away from the radar) and negative velocity (indicating the target is moving towards the radar) appear adjacently is called the velocity convergence zone. This region corresponds to the convergence or upward motion of airflow in the atmosphere. The outflow boundary leading edge of severe convective storms typically exhibits a distinct velocity convergence zone, with a horizontal shear rate reaching 0.01 s. - ¹Above. This convergence zone constitutes the natural dynamic boundary between the storm cell and its surrounding ambient air, and in this embodiment, it is used as the physical boundary constraint for storm cell identification.

[0095] like Figure 3 As shown, the method for identifying 3D seeds in this embodiment includes: traversing all voxels of the fused 3D mesh and marking voxels that simultaneously meet the following two conditions as candidate seed points: Condition 1: Reflectivity factor exceeds 35 dBz; Condition 2: The reflectivity factor is greater than its The 26 adjacent grid points within the three-dimensional neighborhood are the three-dimensional local maxima.

[0096] After obtaining the candidate seed point set, spatial clustering analysis is performed on the candidate seed point set, and two seed points are defined. and Distance metric between and the measure of the difference in reflectivity gradient direction :

[0097] in, Seed point The coordinates in a spatial rectangular coordinate system with the reference radar site as the origin represent the east-west coordinate values ​​(east is positive). Seed point The coordinates in a spatial rectangular coordinate system with the reference radar site as the origin represent the north-south coordinate values ​​(north is positive) and the altitude values.

[0098] In this embodiment, the method for determining the consistency of reflectivity gradient direction is as follows: calculate two candidate seed points. and The angle between the three-dimensional gradient directions of reflectivity When the included angle When the reflectivity is less than the threshold, it is determined that the reflectivity gradient direction is consistent; the included angle The calculation formula is: (29) in, and Representing candidate seed points and The reflectance three-dimensional gradient vector at a given location, where || represents the magnitude of the vector; the reflectance three-dimensional gradient vector... It is composed of the partial derivatives of reflectivity along the x, y, and z directions, and is calculated using the central difference method. Its three components are as follows: (30) (31) (32) in, voxels The reflectivity factor value at that location, Seed point spatial coordinates, , , The three-dimensional radar mesh is respectively in Spatial resolution in the direction.

[0099] When two seed points simultaneously satisfy <2km and When the angle is less than 45°, the seed points are aggregated into the same initial growth nucleus. The above dual conditions ensure that the seed points aggregated into the same initial growth nucleus are not only spatially adjacent to each other, but also belong to the same continuous storm entity in terms of physical structure.

[0100] In this embodiment, dual polarization phase state constraint is used to automatically exclude non-meteorological targets such as ground clutter, biological echoes, and clear-sky echoes during the growth process. The specific steps are as follows: Step 201: Extract the multi-polarization feature set corresponding to the candidate growth voxels encountered when the growth nucleus expands outward, wherein the multi-polarization feature set includes reflectivity factor, differential reflectivity, correlation coefficient, local standard deviation of reflectivity factor and local standard deviation of differential reflectivity; Step 202: Input the multi-polarization feature set into the fuzzy logic classifier. The fuzzy logic classifier calculates the comprehensive membership score for each phase type based on the membership function corresponding to each phase type and the preset feature weights. The phase type c is defined as follows: c=1 for ground clutter / bioecho, c=2 for rain, c=3 for wet snow, c=4 for dry snow, c=5 for graupel, and c=6 for hail. The comprehensive membership score for each phase type is shown below. Calculated by weighted summation: (33) in, For the first Each input feature (reflectance factor, differential reflectance, correlation coefficient, local standard deviation of reflectance factor, local standard deviation of differential reflectance) phase state Features The membership function.

[0101] Step 203: Obtain the phase with the highest score as the classification result for the candidate voxel. When the classification result is ground clutter or bioecho (c=1), the growth of the candidate voxel is prohibited; when the classification result is other phase types (rain, wet snow, dry snow, graupel, or hail), the candidate voxel is allowed to participate in growth, and subsequent extended judgments are performed.

[0102] In this embodiment, the radial velocity convergence boundary termination constraint includes: for each candidate voxel of the outermost layer of the growth body boundary. In candidate voxels The radial velocity shear variable is calculated along the azimuth direction. When the absolute value of the radial velocity shear variable is greater than a set threshold, the candidate voxel is determined to be a convergent boundary voxel, and the convergent boundary voxel is prohibited from being included in the growth region. Furthermore, the growth front terminates its expansion in the direction of the convergent boundary voxel. The calculation formula is: (34) in, Radial distance, Angle of elevation, It is the azimuth angle. voxels The fusion radial velocity, The radial distance where the voxel is located. The azimuth resolution is expressed in radians.

[0103] The criteria for determining the convergence boundary are:

[0104] in, This is an empirical threshold, with a value of 0.008s. - ¹, corresponding to velocity convergence of moderate to high intensity.

[0105] The physical basis for radial velocity convergence boundary detection lies in the fact that the outflow boundary front of a strong convective storm exhibits significant velocity convergence in the radar radial velocity field, meaning that positive and negative velocities appear adjacent to each other, and the horizontal shear variable increases significantly. This convergence zone constitutes the natural dynamic boundary between the storm cell and the surrounding ambient air, and growth should automatically stop when it reaches this boundary. Based on this principle, when a candidate voxel is determined to be a convergence boundary voxel, it is prohibited from being included in the growth region, and the growth front terminates its expansion in the direction of that voxel.

[0106] Alternatively, in addition to using the horizontal shear criterion, the radial velocity convergence boundary detection described above can also use the velocity direction difference criterion as an alternative. Specifically, the directional consistency of the radial velocity in the neighborhood of the growth front is calculated, and when the standard deviation of the direction exceeds a set threshold, it is determined to be a convergence boundary, which can also achieve effective detection of the convergence boundary.

[0107] The aforementioned convergence boundary termination mechanism enables the directional stopping of a single growth nucleus upon reaching a physical boundary. Building upon this, when multiple growth nuclei extend outwards in parallel in three-dimensional space, it is also necessary to address the ownership conflict that arises when the boundaries of different growth nuclei are close together.

[0108] In this embodiment, the competitive segmentation mechanism includes: determining a conflict region when the distance between growth nuclei is less than a set threshold; calculating the weighted distance from voxels to each growth nucleus in the conflict region; taking a smaller value for the weighted distance at locations with a large reflectivity gradient and a larger value at locations with a small reflectivity gradient; constructing a three-dimensional distance difference field; using an isosurface with zero difference field as a watershed segmentation surface; and dividing voxels in the conflict region into corresponding growth regions.

[0109] Specifically, the implementation methods for competitive growth and three-dimensional watershed segmentation are as follows: When the distance between the boundaries of two or more growth nuclei is less than 3 voxels, a competitive segmentation mechanism is triggered.

[0110] First, the conflict region is determined. The conflict region is defined as a set of spatial voxels that satisfy the condition that "the nearest distance to at least two different growth nuclei is less than 3 voxels". If three or more growth nuclei are competing at the same time, a strategy of pairwise competition and pairwise partitioning is adopted: all growth nucleus pairs that meet the starting conditions are detected, and partitioning is performed in each pair of growth nuclei. If a voxel is located in the conflict region of multiple pairs of competing nuclei, it is allocated according to the principle of nearest distance.

[0111] After identifying the conflict region, calculate the weighted distance from each voxel p in the conflict region to each relevant growth nucleus k. The weighted distance takes a smaller value at locations with a larger reflectance gradient (clear boundary) and a larger value at locations with a smaller gradient (overly gentle), so that the watershed boundary naturally tends to unfold along the location where the reflectance change is the gentlest. The weighted distance formula is:

[0112] in, Let be the weighted distance from voxel p to growth nucleus k, and q be the voxel distance in the region where growth nucleus k has grown. This is the gradient weight coefficient (an empirical value of 0.5 is chosen). This represents the reflectance gradient mode at voxel q.

[0113] After obtaining the weighted distances from each voxel to each growth nucleus, a three-dimensional distance difference field is constructed. For any two competing growth nuclei A and B, a three-dimensional distance difference field from p to the two growth nuclei is constructed. The definition of is:

[0114] The isosurface where the difference field is 0 is the watershed dividing surface between the two growth nuclei. Specifically, when When p < 0, the weighted distance from p to growth nucleus A is closer, and p belongs to growth nucleus A; when When = 0, the weighted distances from p to the two growth nuclei are equal, and p lies on the watershed dividing surface, marked as the boundary voxel; when When p > 0, the weighted cohesion of p is closer to the growth nucleus B, and p belongs to the growth nucleus B.

[0115] Using the aforementioned three-dimensional distance difference field criterion, each voxel within the conflict region is assigned to a corresponding growth region, thereby achieving decoupling and separation of adjacent storm cells. Each growth nucleus continues to grow within its designated region until convergence, ultimately resulting in mutually independent storm cells with clearly defined boundaries.

[0116] After the above-described 3D competitive growth and watershed segmentation processes, this stage outputs initial 3D region label maps for K storm cells. Each voxel is assigned a unique cell number or marked as background (i.e., not belonging to any storm). The initial boundary of each cell is the outer envelope of the growth region.

[0117] Step 3: Multi-scale morphological gradient and dynamic contour model 3D boundary refinement stage. Using the growth boundary in the initial 3D region label map as the initial contour, calculate the multi-scale 3D morphological gradient field in the boundary neighborhood, define the energy functional of the 3D dynamic contour model. The energy functional includes external energy terms and internal energy terms. Use the level set method to perform 3D surface evolution on the energy functional. Iterate until the energy functional converges and output the 3D contour of the storm cell.

[0118] like Figure 4 As shown, in this embodiment, the multi-scale morphological gradient and dynamic contour model stereo boundary refinement stage is used to refine the initial storm boundary obtained in step 2, so that the boundary is continuous and converges to the physical boundary position where the reflectance space changes most drastically.

[0119] First, a narrowband is constructed. For each storm cell output in step 2, a narrowband processing region is constructed centered on its boundary surface. The purpose is to reduce the computational scope from the entire area to a local region near the boundary surface, thereby significantly reducing the computational load. Specifically, the Euclidean distance from each voxel in the 3D radar mesh to the boundary surface is calculated, and voxels with a distance less than 5 voxel spaces are added to the narrowband region. Subsequent calculations are performed only within this narrowband region.

[0120] After completing the narrowband construction, multi-scale morphological gradient calculations are performed within the narrowband region. This embodiment uses spherical structural elements with radius scales of 3, 5, and 7 voxels as the basic operational units for morphological operations. The physical basis for choosing spherical structures is that the spherical neighborhood is consistent with the isotropic characteristics of electromagnetic wave propagation and precipitation particle distribution. The small-scale radius (3 voxels) is used to capture sharp boundaries and the edges of strong echo cores; the medium-scale radius (5 voxels) is used to capture transition zones with moderate ambiguity; and the large-scale radius (7 voxels) is used to capture blurred, gradually changing boundaries and the periphery of weak echoes.

[0121] For each scale *s*, the morphological gradient at that scale is calculated using the difference between 3D dilation and 3D erosion. The 3D dilation operation takes the maximum reflectivity factor *Z* of all voxels within a spherical neighborhood of radius *s* centered on voxel *p*, causing high-reflectivity regions to expand outwards and low-reflectivity regions to compress. At the boundary, dilation expands strong echoes. The 3D erosion operation takes the minimum reflectivity of all voxels within a spherical neighborhood of radius *s* centered on voxel *p*, causing low-reflectivity regions to expand outwards and high-reflectivity regions to be eroded. At the boundary, erosion causes weak echoes to contract inwards. The dilation results at scale *s* are then analyzed. Corrosion results The morphological gradient at scale s is obtained by subtracting the difference. for: (38) in, A spherical structural element with scale s. Represents the reflectivity factor, and ⊕ represents the three-dimensional dilation operation. This represents a three-dimensional erosion operation. Three-dimensional dilation takes the maximum value of the reflectivity factor within a spherical neighborhood, while three-dimensional erosion takes the minimum value of the reflectivity factor within a spherical neighborhood.

[0122] After obtaining the morphological gradients at three scales, a weighted average method is used for fusion, assigning the largest weight to the small-scale gradient to preserve fine boundaries, and assigning the smallest weight to the large-scale gradient to help capture gradient boundaries:

[0123] in, , , These represent the weights of spherical structural units with radius scales of 3, 5, and 7 voxels, respectively. , , This represents the corresponding morphological gradient value.

[0124] The fusion result is normalized to the [0,1] interval to obtain the final gradient field. The gradient field has a larger value (close to 1) in locations where the reflectivity changes drastically in space, and a smaller value (close to 0) in uniform regions, thus accurately indicating the location of the physical boundary.

[0125] After obtaining the multi-scale morphological gradient field, an energy functional for the 3D dynamic contour model is established using this gradient field. This energy functional is a weighted combination of external and internal energy terms: the external energy term drives the contour to move towards the point of maximum multi-scale gradient, while the internal energy terms include surface tension terms (constraining the smoothness of the surface) and bending stiffness terms (constraining the continuity of the surface). The total energy functional is defined as:

[0126] in, Let the level set function be the function corresponding to the surface. For external energy function, For surface tension, For bending stiffness, and These are the weighting coefficients for controlling the strength of the surface tension term and the bending stiffness term, respectively.

[0127] According to the variational principle (Euler-Lagrange equations), the energy functional... A necessary condition for taking a minimum value is that its first variation is 0. The evolution equation of the level set is derived using the gradient descent method:

[0128]

[0129]

[0130]

[0131]

[0132] in, The average curvature of a surface is used to describe the degree of curvature of the surface at a point in space. This is the regularized form of the Dirac function, used to restrict the evolution to occur only near narrowband boundaries. For the final multi-scale morphological gradient field, Gradient field The spatial gradient vector, level set function The spatial gradient vector, yes Gradient in three-dimensional space, For surface gradient, It is the unit normal vector of the surface. It is the component of the gradient in the normal direction.

[0133] Alternatively, in addition to using the level set method, an explicit surface evolution method based on triangular meshes can be used as an alternative to solve the above-mentioned three-dimensional dynamic contour model. Specifically, the initial growth body boundary is transformed into a triangular mesh surface, and the gradient force and surface internal force derived from multi-scale morphological gradients are calculated at each triangular patch node. The surface evolution converges by updating the node positions through explicit time stepping, which can also achieve the refinement of the boundary surface.

[0134] By iteratively solving the above evolution equations, the initial profile gradually evolves under the combined drive of external gradient forces and internal smoothing forces until the energy functional converges. At this point, the profile stops at the physical boundary position where the reflectance space changes most drastically. After the energy functional converges, the refined three-dimensional boundary surface of each storm is determined, and the relevant parameters of the storm cells are output, including: unique cell number, identification timestamp, coordinates of the three-dimensional outer envelope vertex, volume, three-dimensional coordinates of the centroid, maximum reflectance factor and its location, vertical integral liquid water content, etc.

[0135] This embodiment achieves full-process automation of storm cell identification by constructing an integrated physical model-driven processing chain of "fusion → growth → refinement → identification".

[0136] At the data fusion level, the multi-source radar 3D mesh fusion method in this embodiment incorporates range dependence factors, physical confidence factors, and beam occlusion perception factors into the fusion decision through an adaptive weighting function and a frequency-band synthesis strategy. The range dependence factor enables a smooth transition from near-range high-confidence X-band data to attenuated data at long range in the S-band. The physical confidence factor selects the most suitable band contribution for different physical quantities such as reflectivity factor, differential reflectivity, and correlation coefficient. The beam occlusion perception factor effectively eliminates distorted data caused by terrain occlusion. Furthermore, the frequency-band synthesis strategy further achieves complementary advantages between the S-band and X-band in the frequency domain, ensuring that the fusion result maintains the wide-area coverage reliability of the S-band while fully absorbing the high-resolution characterization capability of the X-band for fine storm boundaries and internal details, thus solving the problem of insufficient utilization of multi-source radar data in existing methods.

[0137] At the physical segmentation level, the 3D competitive growth method in this embodiment fundamentally changes the traditional segmentation mode that relies solely on a single reflectivity threshold by introducing radial velocity convergence boundary termination constraints and dual-polarization phase state classification guidance. Radial velocity convergence boundary detection uses the natural dynamic boundary (i.e., convergence zone) between the storm cell and its surrounding environment as the growth termination condition, ensuring that the segmentation results conform to the hydrodynamic characteristics of the storm and effectively solving the problem of adjacent cells in dense storm areas becoming stuck and unable to separate due to interconnected reflectivity echoes. Dual-polarization phase state classification constraints exclude non-meteorological targets such as ground clutter and biological echoes in real time during the growth stage, ensuring that the segmentation results conform to the microphysical characteristics of the storm and preventing false targets from entering the identification results.

[0138] At the boundary refinement level, the stereo boundary refinement method in this embodiment refines the initial voxel-level rough boundary into a smooth and continuous three-dimensional surface by constructing a multi-scale morphological gradient field and solving a level set dynamic contour model. The multi-scale morphological gradient uses spherical structuring elements to perform multi-scale expansion and erosion operations on the reflectivity field, fully considering the isotropic characteristics of electromagnetic wave propagation and precipitation particle distribution. The level set method, by defining an energy functional containing external and internal energy terms, drives the contour to naturally converge to the physical boundary position where the reflectivity space changes most drastically while maintaining the topological structure, filling the technical gap in the prior art from the voxel level to a smooth physical boundary.

[0139] In summary, this embodiment achieves fully automated storm entity identification through the synergistic effect of the above three levels, from multi-source data fusion and three-dimensional competitive growth to three-dimensional boundary refinement, significantly improving the accuracy, physical rationality, and boundary refinement of the identification.

[0140] This embodiment also provides a three-dimensional storm cell identification and boundary refinement system based on multi-source data, including a microprocessor and a memory interconnected, wherein the microprocessor is programmed or configured to execute a three-dimensional storm cell identification and boundary refinement method based on multi-source data.

[0141] The above description is merely a preferred embodiment of the present invention and is not intended to limit the invention in any way. Although the present invention has been disclosed above with reference to preferred embodiments, it is not intended to limit the invention. Therefore, any simple modifications, equivalent changes, and alterations made to the above embodiments based on the technical essence of the present invention without departing from the scope of the present invention should fall within the protection scope of the present invention.

Claims

1. A method for identifying and refining the boundaries of 3D storm cells based on multi-source data, characterized in that, include: In the multi-source radar data 3D mesh fusion processing stage, multi-source radar data is acquired and fused into 3D radar mesh data; In the three-dimensional competitive growth stage guided by convergence boundary constraints and physical phase state, local maxima with reflectivity factors exceeding a threshold are searched in the fused three-dimensional radar grid data as candidate seed points. Spatial clustering analysis is performed on these candidate seed points, and those satisfying the spatial distance threshold and reflectivity gradient direction consistency are aggregated into an initial growth kernel. Starting from the initial growth kernel, candidate voxels are expanded outward layer by layer according to multi-constraint growth criteria for three-dimensional competitive growth. When the distance between growth kernels is less than a set threshold, a competitive segmentation mechanism is triggered to achieve automatic decoupling and separation of adjacent storm cells, and the initial three-dimensional region label map of each storm cell is output. The multi-constraint growth criteria include reflectivity continuity constraints, dual polarization phase state constraints, and radial velocity convergence boundary termination constraints. In the multi-scale morphological gradient and dynamic contour model 3D boundary refinement stage, the growth boundary in the initial 3D region label map is used as the initial contour. The multi-scale 3D morphological gradient field is calculated in the boundary neighborhood. The energy functional of the 3D dynamic contour model is defined. The energy functional includes external energy terms and internal energy terms. The level set method is used to perform 3D surface evolution on the energy functional. The iteration continues until the energy functional converges, and the 3D contour of the storm cell is output.

2. The method for three-dimensional storm cell identification and boundary refinement based on multi-source data according to claim 1, characterized in that, The multi-source radar data 3D mesh fusion processing stage includes: Acquire volume scan observation data from at least one S-band weather radar and at least one X-band weather radar, select an S-band radar that meets the preset requirements as the reference radar, and establish a three-dimensional polar coordinate reference grid with the reference radar site as the origin. The volume scan observation data of other S-band and X-band weather radars are mapped to the three-dimensional polar coordinate reference grid to obtain multi-source radar observation data with the three-dimensional polar coordinate reference grid as a unified spatial framework. An adaptive weighting function is constructed to perform weighted fusion of volume scan observation data from multiple weather radars at each grid point, resulting in the fusion results of volume scan observation data from S-band weather radar and X-band weather radar, respectively. The volume scan observation data fusion results of the S-band weather radar and the volume scan observation data fusion results of the X-band weather radar are fused in the frequency domain using a frequency-band synthesis strategy or a multi-scale fusion strategy based on wavelet transform, and the fused three-dimensional fused radar grid data is output.

3. The method for identifying and refining the boundaries of three-dimensional storm cells based on multi-source data according to claim 2, characterized in that, The adaptive weighting function is jointly determined by a distance dependence factor, a physical confidence factor, and a beam obstruction perception factor. The distance dependence factor is used to adjust the weights based on the distance between the weather radar and the observed target. The physical confidence factor is used to adjust the weights based on the measurement reliability of different physical quantities on radars in different bands. The beam obstruction perception factor is used to adjust the weights based on terrain obstruction.

4. The method for identifying and refining the boundaries of three-dimensional storm cells based on multi-source data according to claim 2, characterized in that, The frequency-band synthesis strategy includes: performing three-dimensional Gaussian pyramid decomposition on the fusion results of volume scan observations from the S-band weather radar and the X-band weather radar to obtain Laplace components at each level; selectively weighting synthesis according to the spatial frequency level, with the high-frequency layer weighted by the X-band Laplace component and the low-frequency layer weighted by the S-band Laplace component; reconstructing the Laplace components of each layer after weighted synthesis layer by layer through the pyramid to output the fused three-dimensional mesh data after frequency-band synthesis.

5. The method for three-dimensional storm cell identification and boundary refinement based on multi-source data according to claim 1, characterized in that, In the three-dimensional competitive growth stage guided by convergence boundary constraints and physical phase state, the method for determining the consistency of the reflectivity gradient direction is as follows: calculate two candidate seed points. and The angle between the three-dimensional gradient directions of reflectivity When the included angle When the reflectivity is less than the threshold, it is determined that the reflectivity gradient direction is consistent; the included angle The calculation formula is: in, and Representing candidate seed points and The reflectance three-dimensional gradient vector at a given location, where || represents the magnitude of the vector; the reflectance three-dimensional gradient vector It is composed of the partial derivatives of reflectivity along the x, y, and z directions, and is calculated using the central difference method. Its three components are as follows: in, voxels The reflectivity factor value at that location, Seed point spatial coordinates, , , The three-dimensional radar mesh is respectively in Spatial resolution in the direction.

6. The method for three-dimensional storm cell identification and boundary refinement based on multi-source data according to claim 1, characterized in that, In the three-dimensional competitive growth stage guided by convergence boundary constraints and physical phase states, the dual-polarization phase state constraints include: extracting the multi-polarization feature set corresponding to the candidate growth voxels encountered when the growth nucleus expands outward; inputting the multi-polarization feature set into a fuzzy logic classifier, which calculates the comprehensive membership score of each phase state type based on the membership function and preset feature weights, and takes the phase with the highest score as the classification result; when the classification result corresponding to the candidate growth voxel is a preset phase state, the growth of the candidate growth voxel is prohibited; the comprehensive membership score... The calculation formula is: in, For the first Each input feature For feature weights, phase state Features The membership function.

7. The method for three-dimensional storm cell identification and boundary refinement based on multi-source data according to claim 1, characterized in that, In the three-dimensional competitive growth stage guided by convergence boundary constraints and physical phase states, the radial velocity convergence boundary termination constraint includes: for each candidate voxel extending one layer beyond the growth body boundary. In the candidate voxels Calculate the candidate voxel in the azimuth direction Radial velocity shear variable: When the absolute value of the radial velocity shear variable is greater than a set threshold, the candidate voxel is determined to be a convergent boundary voxel, the convergent boundary voxel is prohibited from being included in the growth region, and the growth front terminates its expansion in the direction of the convergent boundary voxel; the radial velocity shear variable... The calculation formula is: in, Radial distance, Angle of elevation, It is the azimuth angle. voxels The fusion radial velocity, The radial distance where the voxel is located. This is the azimuth resolution.

8. The method for three-dimensional storm cell identification and boundary refinement based on multi-source data according to claim 1, characterized in that, In the three-dimensional competitive growth stage guided by convergence boundary constraints and physical phase state, the competitive segmentation mechanism includes: determining a conflict region when the distance between growth nuclei is less than a set threshold; calculating the weighted distance from voxels to each growth nucleus in the conflict region; taking the minimum value at the location with a large reflectivity gradient and the maximum value at the location with a small reflectivity gradient; constructing a three-dimensional distance difference field; using the isosurface with zero difference field as the watershed segmentation surface; and dividing the voxels in the conflict region into the corresponding growth regions.

9. The method for three-dimensional storm cell identification and boundary refinement based on multi-source data according to claim 1, characterized in that, The method for calculating the multi-scale three-dimensional morphological gradient field in the stereo boundary refinement stage of the multi-scale morphological gradient and dynamic contour model includes: Centered on the boundary surface of each storm cell, calculate the Euclidean distance from each voxel in the three-dimensional radar grid to the boundary surface, and add voxels with a distance of less than 5 voxel spaces to the narrowband region. Within the narrow band region, spherical structure elements with radius scales of 3, 5, and 7 voxel distances were used to calculate the difference between the three-dimensional expansion and three-dimensional erosion of each voxel at each scale, thus obtaining the morphological gradient at each scale. : in, A spherical structural element with scale s. Represents the reflectivity factor, and ⊕ represents the three-dimensional dilation operation. This represents a three-dimensional erosion operation, where the three-dimensional dilation is the maximum value of the reflectivity factor within a spherical neighborhood, and the three-dimensional erosion is the minimum value of the reflectivity factor within a spherical neighborhood. The morphological gradients at each scale are weighted and averaged to fuse, with the small-scale gradient having the largest weight and the large-scale gradient having the smallest weight. The fusion result is normalized to the [0,1] interval to obtain the final gradient field.

10. A three-dimensional storm cell identification and boundary refinement system based on multi-source data, comprising a microprocessor and a memory interconnected, characterized in that, The microprocessor is programmed or configured to execute the method for three-dimensional storm cell identification and boundary refinement based on multi-source data as described in any one of claims 1 to 9.

Images