Artificial precipitation enhancement effect test method based on radar tracking and raindrop spectrum

Abstract

The invention discloses an artificial precipitation enhancement effect test method based on radar tracking and a raindrop spectrum. The method comprises the steps of collecting Doppler weather radar data and weather phenomenon instrument drop spectrum data of a preset area; segmenting the Doppler weather radar data by using a multi-scale morphological cloud body to obtain an influence area and a contrast area, and performing adaptive tracking on the cloud body according to the Doppler weather radar data to obtain a catalytic operation influence area; performing dynamic matching on the catalytic operation influence area and the comparison area by adopting multi-dimensional dynamic similarity matching to obtain a comparison cloud body, constructing a Z-I relation library according to the weather phenomenon instrument drop spectrum data, and obtaining radar estimation rainfall and high temporal-spatial resolution rainfall according to the Z-I relation library; and carrying out space-time alignment on the comparison cloud body and the Z-I relation library to obtain absolute increased rainfall and relative increased rainfall, carrying out physical inspection on the comparison cloud body, carrying out statistical inspection according to the absolute increased rainfall, the relative increased rainfall and high-space-time-resolution rainfall, and outputting an inspection result.

Description

Technical Field

[0001] This invention relates to the field of artificial rain enhancement effect testing technology, and in particular to a method for testing the effect of artificial rain enhancement based on radar tracking and raindrop spectrum. Background Technology

[0002] Verifying the effectiveness of artificial rain enhancement is a core aspect of weather modification operations, and its scientific validity directly determines the optimization of operational plans and the rationality of resource allocation. However, the uneven spatial and temporal distribution and high variability of natural precipitation, coupled with the incomplete understanding of the influence mechanism of catalysts on the physical processes of cloud precipitation, make verifying the effectiveness of artificial rain enhancement a recognized technical challenge both domestically and internationally.

[0003] Existing verification methods are divided into two categories: statistical verification and physical verification. However, they have significant limitations: statistical verification relies on ground rain gauge data, which is constrained by insufficient observation station density and difficulty in data transmission and storage, resulting in limited accuracy in precipitation estimation; physical verification can obtain macro- and micro-changes in clouds through radar tracking, but it is not synchronized with statistical verification in time and space, and the spatial misalignment between cloud areas and rain areas is prominent, making it difficult to form a unified verification.

[0004] Traditional radar precipitation estimation uses a generic ZI relationship without incorporating localized raindrop spectral characteristics, resulting in significant errors. Furthermore, the selection of comparison areas is often static, ignoring dynamic cloud evolution and leading to insufficient similarity among the compared clouds. These issues severely impact the accuracy and reliability of rain enhancement effect verification, hindering the standardized development of weather modification operations. Therefore, a comprehensive verification method integrating dynamic tracking, localized parameters, and multi-dimensional validation is urgently needed. Summary of the Invention

[0005] The purpose of this invention is to provide a method for verifying the effectiveness of artificial rain enhancement based on radar tracking and raindrop spectrum.

[0006] To achieve the above objectives, the present invention is implemented according to the following technical solution: This invention includes the following steps: Collect Doppler weather radar data and weather phenomenon instrument drop spectrum data of a preset area, and preprocess the Doppler weather radar data and the weather phenomenon instrument drop spectrum data. The Doppler weather radar data is segmented by multi-scale morphological cloud body to obtain the affected area and the comparison area. Based on the Doppler weather radar data, the cloud body is adaptively tracked to obtain the catalytic operation's affected area. Multidimensional dynamic similarity matching is used to dynamically match the catalytic operation influence area and the comparison area to obtain the comparison cloud body. A ZI relationship database is constructed based on the drop spectrum data of the weather phenomenon instrument. Radar-estimated precipitation and high spatiotemporal resolution precipitation are obtained based on the ZI relationship database. The absolute and relative rainfall increases are obtained by spatiotemporally aligning the comparison cloud body with the ZI relational database. The comparison cloud body is then subjected to physical verification. Statistical verification is performed based on the absolute and relative rainfall increases, and the verification results are output.

[0007] Furthermore, the method for obtaining the affected region and the contrast region includes: Data standardization was performed on Doppler weather radar data, and cloud edge features were extracted using a morphological gradient operator. The expression is as follows:

[0008] in For expansion operations, For erosion calculation, It is a 1*1 flat plate structural element. This is the raw radar data. The morphological gradient of the raw radar data; Select 3 groups of structural elements of different sizes , , These correspond to cloud features at small, medium, and large scales, respectively. The results of weighted fusion of filtering at each scale are expressed as follows:

[0009] in The result of the opening-closing operation. The weighting coefficients for the first scale. The weighting coefficients for the second scale. The weighting coefficients for the third scale. This is the result of the opening-closing operation at the first scale. This is the result of the opening-closing operation at the second scale. This is the result of the opening-closing operation at the third scale; Multi-scale filtering is applied to the morphological gradient image to eliminate peak features with significant cloud contours while preserving them. A watershed transform is then performed on the filtered gradient image, dividing it into several initial regions. The merging cost of adjacent regions is calculated using the difference in mean grayscale value and area as indicators.

[0010] in Let V be the average reflectance of the v-th region. Let be the average reflectance of the r-th region. Let v be the area of ​​the v-th region. Let r be the area of ​​the r-th region. The cost of merging the v-th region and the r-th region; The merging priority is managed using a region adjacency graph and a min-heap data structure. Each time, the region pair with the lowest global cost is merged until the preset number of regions or cost threshold is met. Regions with smaller areas and similar gray levels are merged first. Multidimensional features were calculated for the merged cloud region to establish a cloud attribute database. The multidimensional features include mean reflectivity, echo top height, volume, and vertical cumulative liquid water content. Based on the location of the operation and the trajectory of the cloud, cloud units that contain the catalytic operation point and meet the preset conditions are marked as the influence area; cloud units with the highest feature similarity in the upwind direction of the influence area are selected as the comparison area; the similarity is calculated by the dynamic time warping algorithm to match the life history characteristics of the cloud in the influence area before the operation; the life history characteristics are the volume change trend and the evolution of vertical cumulative liquid water content.

[0011] Furthermore, the method for obtaining the catalytic operation-affected zone includes: Spatiotemporal registration is performed on Doppler radar reflectivity data at continuous time intervals to unify spatial resolution and time intervals, and range attenuation correction and elevation angle normalization are applied. The affected area and the comparison area are respectively used as the region of interest (ROI) for tracking. Noise points and non-cloud echoes outside the region are removed, and the data within the region of interest is Gaussian smoothed. Set time The radar echo matrix is ,time The matrix is For each grid point within the ROI Define size as × template window ,exist Slide within the search window to calculate the cross-correlation coefficient:

[0012] in For a sub-region within the search window, This is the reflectivity value of the template. The reflectance value of the sub-region. For local pixel indexes within the window, To explore the pixel values ​​of a sub-region of the window, For grid point coordinates, It is a displacement vector. This is the cross-correlation coefficient between the template window and the sub-region at the displacement within the search window; Define the confidence index as follows:

[0013] in The highest cross-correlation coefficient, The second largest correlation coefficient, The confidence index is used to evaluate the motion vector. ; An improved Lucas-Cannard method is used, assuming that the radar echo grayscale value is conserved over a short period of time. Spatiotemporal differentiation is performed on the reflectivity data to construct the optical flow equation:

[0014] in For time gradient, For the lateral spatial gradient, For the longitudinal spatial gradient; The initial optical flow vector is obtained by solving the overdetermined equations in a 3×3 neighborhood using the least squares method. ; Introducing multi-scale pyramid optical flow: Gaussian pyramid downsampling is performed on the echo data, optical flow is calculated starting from the top layer, and the result is used as the initial value for the lower layer. Accumulated error is reduced through iterative optimization, and the optical flow motion vector field of each grid point in the ROI is output. Based on the complementary characteristics of radar echo tracking and optical flow methods, dynamic weighting coefficients are designed to fuse optical flow motion vectors and motion vectors, with the following expression:

[0015] in To fuse motion vectors, For dynamic weighting coefficients, The magnitude of the optical flow gradient; Construct a global optimization objective function that minimizes the following energy term, expressed as:

[0016] in Let be the objective function. For the neighborhood, For smoothing coefficients, Grid points to be optimized Motion vector, The fused motion vector of the grid points. For the neighborhood points of the grid point, The motion vector of the neighboring points, For data fidelity items, For smoothing terms; The optimized motion vector field is subjected to spatiotemporal filtering, which includes spatial smoothing and temporal smoothing. Spatial smoothing uses median filtering to remove isolated outlier vectors. Temporal smoothing is a weighted average of the motion vectors at three consecutive time points. Based on the optimized motion vector field, a region growing method is used to track cloud clusters: the cloud cluster outline of the initial influence area is used as the seed region; according to the motion vector... Predict the location of the cloud grid points at the next moment; perform morphological dilation on the predicted region and compare it with... The radar echoes at any given time are matched to obtain the overlap. When the overlap is greater than 70%, the match is successful. The trajectory of the cloud formation as it evolves over time is generated iteratively, which is the dynamic range of the catalytic operation's influence area.

[0017] Furthermore, the method for obtaining the comparative cloud body includes: The three-dimensional structural features of influencing clouds and potential contrasting clouds are extracted from Doppler radar data, and a multi-dimensional feature vector is constructed, expressed as follows:

[0018] in For a moment The multidimensional feature vectors, For a moment echo peak height, For a moment echo volume, This is the transpose of the vector. For vertically accumulated liquid water content, The density of liquid water, This is the base value in the vertical direction. This represents the upper limit in the vertical direction; Select a time window before the start time of the task. By sampling according to the radar data temporal resolution, we obtained the time series of influencing cloud features and the time series of potential comparative cloud features:

[0019]

[0020] in The start time of the operation. The time window length, To influence cloud feature sequences, For the k-th potential contrast cloud feature sequence, To influence the time-by-time feature values ​​of the cloud feature sequence, Let be the time-series feature value of the k-th candidate cloud feature sequence. The radar data sampling interval; the length of the cloud feature sequence is affected by... The potential contrast cloud feature sequence length is ; Zero-mean standardization of features is used to screen candidate clouds that meet spatial conditions from the segmented cloud bodies. Spatial conditions include location constraints and range constraints. Location constraint: located upwind of the influencing cloud, with a horizontal distance of ≥20km. Range constraint: located in the same radar scanning area as the influencing cloud, and with an elevation difference of ≤2km between the cloud body centers. Calculate the static feature similarity between the candidate cloud and the influencing cloud at the start of the time window:

[0021] in The static feature similarity between candidate clouds and influencing clouds. Let be the number of dimensions of the static features. For at any time The m-th static eigenvalue of the cloud is affected. For at any time The m-th static feature value of the h-th candidate cloud; Candidate clouds with static feature similarity greater than or equal to a similarity threshold are retained. A multidimensional cost matrix is ​​defined based on the standardized impact cloud feature sequence and the candidate cloud feature sequence, expressed as:

[0022] in Let m be the dynamic feature value of the m-th dimension at time t in the standardized cloud sequence. Let tu be the dynamic feature value at time m in the standardized candidate cloud sequence. The optimal time alignment path is found using dynamic programming to minimize the cumulative distance between the two sequences.

[0023] in For from (1,1) to The total cumulative distance of the optimal path, For path selection operators, for The cost matrix; Boundary conditions are , Introduce the Sakoe-Chiba window constraint, allowing windows to be opened only on both sides of the diagonal. Search path within the range; where ; Dynamic time warping distance Define dynamic similarity: For candidate cloud computing dynamic similarity, the cloud with the highest dynamic similarity is selected as the initial comparison cloud. If the dynamic similarity is less than the comparison similarity threshold, it is determined that there is no qualified comparison cloud, and the cloud search range is expanded or the time window is adjusted. The matching quality is verified by trend consistency and morphological similarity. Trend consistency: Calculate the Pearson correlation coefficient of each characteristic time series of the influencing cloud and the comparison cloud. If the echo top height, volume, and vertical cumulative liquid water content are all greater than or equal to 0.7, the result is positive. Morphological similarity: Compare the three-dimensional morphological changes of the influencing cloud and the comparison cloud within the time window. Based on the optimal path of dynamic time warping, the feature sequence of the comparison cloud is aligned with the sequence of the influencing cloud, and the comparison cloud body is output.

[0024] Furthermore, the method for constructing the ZI relational database includes: Given the radar reflectivity factor and rainfall intensity, the expression is:

[0025]

[0026] in Let be the diameter of the precipitation particles at level i. The inner diameter per cubic meter The number of particles, For reflectivity factor, Let be the mass of the i-th order raindrop. Let be the final velocity of the i-th order raindrop. Rainfall intensity; The median volume diameter was obtained, and the droplet spectral samples were classified by rainfall intensity and median volume diameter. The samples were divided into rainfall intensity categories and median diameter categories. The rainfall intensity categories were divided into 16 categories with unequal intervals, covering the main precipitation intensity range of 0~5 mm / h. The median diameter categories were divided into small, medium and large droplets. Drop spectrum samples falling within the same grid are averaged to generate 48 distinct average drop spectra based on averaging rules. Averaging rule: Take the arithmetic mean of the drop spectrum data of all samples within the same grid. For each type of average drop spectrum, the least squares method is used for fitting. coefficients in Sum of Indices Based on the average droplet spectrum, the radar reflectivity factor and rainfall intensity are calculated, and an objective function is constructed. ,in The rainfall intensity for the k-th sample is calculated from the drop spectrum. The radar reflectivity factor of the k-th sample is calculated from the drop spectrum; ZI relationships are classified and labeled according to precipitation cloud system types to form localized ZI relationships; precipitation cloud system types include stratiform clouds, cumulus clouds, and mixed stratocumulus clouds.

[0027] Furthermore, the expressions for the absolute rainfall increase and the relative rainfall increase are as follows:

[0028]

[0029] in This refers to the absolute increase in rainfall. This refers to the relative increase in rainfall. To account for the total precipitation in the affected area, The total precipitation in the comparison area.

[0030] Furthermore, the method for physically examining the comparison cloud includes: Based on radar tracking, time-series curves of echo top height, echo volume, maximum reflectivity, and vertical cumulative liquid water content of clouds in the affected and control areas were plotted. The effectiveness of the catalytic operation was physically verified based on five physical quantities: echo top height, echo volume, maximum reflectivity, vertical cumulative liquid water content, and precipitation flux. Echo top height refers to the echo top height of a given dBZ threshold cell, used to identify radar echoes. The echo top height of a cell refers to the maximum height of radar echoes greater than or equal to the dBZ threshold. Echo volume refers to the echo volume of a given dBZ threshold cell. Maximum reflectivity refers to the maximum dBZ value within a given dBZ threshold cell. Vertical cumulative liquid water content refers to the vertical cumulative liquid water content calculated from the maximum reflectivity factor of each layer from the zero-degree layer to the top of the cell within a given dBZ threshold cell. Precipitation flux refers to the sum of precipitation flux calculated from the vertical maximum reflectivity factor of each cell in the horizontal projection within a given dBZ threshold cell. The curves of the changes of five physical quantities over time were plotted for the operational cloud unit and the comparison cloud unit, respectively. The operational information was marked, and the time series changes of the physical quantities of the operational cloud unit and the comparison cloud unit were compared and analyzed.

[0031] Furthermore, the statistical test method includes: The area affected by the catalytic operation was designated as the treatment group, and the dynamically matched control area was designated as the control group. The time window covered 1 hour before the operation, 0-2 hours during the operation, and 1 hour after the operation. High spatiotemporal resolution precipitation data were aggregated into a spatiotemporal scale consistent with radar-estimated precipitation. Calculate the absolute and relative rainfall increases during the time period of the catalytic influence of cloud bodies in the affected area and the control area based on radar tracking; Calculate the Pearson correlation coefficient between the precipitation sequences of the affected area and the control area before the operation. If the Pearson correlation coefficient is greater than 0.6, the correlation requirement is met. Use the two independent samples t test to verify whether there is no significant difference between the mean precipitation of the two areas before the operation. If the p value of the t test is greater than 0.05, the mean difference test is passed. When precipitation data from the affected area and the control area are independent, a two-sample t-test is used to verify whether the rainfall increase is significant: Null hypothesis : Alternative hypothesis : ; Calculate statistics : ; in The average precipitation in the affected area. To compare the average precipitation in the region, The sample size of the affected area. The sample size for the comparison region. To represent the variance of precipitation samples in the affected area, The variance of precipitation samples in the comparison area; When there is a spatiotemporal correlation between the affected area and the control area, a paired t-test is used: null hypothesis : Alternative Hypothesis : ; Calculate the statistics: ,in This represents the average absolute rainfall increase. The standard deviation of absolute rainfall increase This represents the number of paired samples; Monte Carlo simulation test: Precipitation data from the catalytic operation's impact area and the control area were merged into a total sample set. Subsamples with the same sample size as the original impact area and control area were randomly selected from the total sample set to serve as virtual impact area and virtual control area, respectively. Virtual rainfall increase was calculated, and the probability distribution of virtual rainfall increase was generated by repeated sampling. The 95% confidence interval of the virtual rainfall increase was calculated using the percentile method. If the actual absolute rainfall increase is greater than the upper quantile of the 95% confidence interval, the rainfall increase effect is significant. Statistical test of relative rainfall increase rate: For relative rainfall increase, take the logarithm of the precipitation data to transform the relative difference into an absolute difference, and construct a statistical measure. ,in To determine the standard deviation of logarithmic precipitation, a standard normal distribution test is performed on the rainfall increase. A significance level of 0.05 is set, and the one-sided critical value is obtained from the standard normal distribution table. If the calculated statistic is greater than 1.645 and... If the value is greater than 0, the rain enhancement effect is considered significant at a 95% confidence level. After aligning the cloud bodies with the ZI relation database in time and space, the time series and spatial range of the cloud body physical quantities from the physical test are simultaneously matched to the spatiotemporal scale of the precipitation data from the statistical test.

[0032] The beneficial effects of this invention are: This invention is a method for verifying the effectiveness of artificial rain enhancement based on radar tracking and raindrop spectrum. Compared with existing technologies, this invention has the following technical advantages: This invention integrates radar dynamic tracking and raindrop spectrum analysis to achieve precise spatiotemporal matching of physical and statistical verification, solving the problems of spatiotemporal disconnect and cloud-rain area misalignment in traditional methods, thus improving the scientific rigor of verification. It constructs a localized cloud type ZI relation database and optimizes raindrop spectrum processing using the SIFT method, reducing radar precipitation estimation errors and overcoming the limitations of insufficient ground observation station density. It employs multi-dimensional dynamic similarity matching to select comparison clouds, combined with adaptive cloud tracking technology, improving the similarity and tracking accuracy of comparison clouds and reducing interference from natural variability. By integrating multiple statistical and physical verification indicators, it forms a comprehensive verification system, resulting in more reliable verification results. This provides precise support for optimizing artificial rain enhancement operations and is suitable for operational deployment needs. Attached Figure Description

[0033] Figure 1 This is a flowchart illustrating the steps of the artificial rain enhancement effect verification method based on radar tracking and raindrop spectrum of the present invention. Figure 2 This is a schematic diagram illustrating the selection of the affected area and the comparison area in the embodiments of this specification. Detailed Implementation

[0034] The present invention will be further described below through specific embodiments. The illustrative embodiments and descriptions herein are used to explain the present invention, but are not intended to limit the present invention.

[0035] The present invention provides a method for verifying the effectiveness of artificial rain enhancement based on radar tracking and raindrop spectrum, comprising the following steps: like Figure 1 As shown, this embodiment includes the following steps: Collect Doppler weather radar data and weather phenomenon instrument drop spectrum data of a preset area, and preprocess the Doppler weather radar data and the weather phenomenon instrument drop spectrum data. In the actual assessment, taking Province A1 as the research object, affected by the upper-level trough and the low-level shear line, the low-level jet stream in front of the shear line was established and guided the northward transport of water vapor. On April 25, moderate rain and local heavy rain occurred in the central and northern parts of Province A1. Among them, sleet or light snow occurred in the eastern part of A2, the eastern part of A3, most of A4, most of A5, the northwestern part of A6 and Dengta, and moderate snow occurred in A6 County. Using silver iodide as a de-icing agent as a catalyst, this sortie operated at an altitude of 4500m, with the main operating period between 13:00 and 14:00, near the -16°C level. The actual operating altitude and temperature were reasonable and within radar echo range. The aircraft's operational information is shown in Table 1, and the ground operation information is shown in Table 2. Table 1. Aircraft Rain Enhancement Operations on April 25, 2023

[0036] Table 2. Ground-based rain enhancement operations on April 25, 2023

[0037] The Doppler weather radar data is segmented by multi-scale morphological cloud body to obtain the affected area and the comparison area. Based on the Doppler weather radar data, the cloud body is adaptively tracked to obtain the catalytic operation's affected area. Multidimensional dynamic similarity matching is used to dynamically match the catalytic operation influence area and the comparison area to obtain the comparison cloud body. A ZI relationship database is constructed based on the drop spectrum data of the weather phenomenon instrument. Radar-estimated precipitation and high spatiotemporal resolution precipitation are obtained based on the ZI relationship database. The absolute and relative rainfall increases are obtained by spatiotemporally aligning the comparison cloud body with the ZI relational database. The comparison cloud body is then subjected to physical verification. Statistical verification is performed based on the absolute and relative rainfall increases, and the verification results are output.

[0038] In this embodiment, the method for obtaining the affected region and the contrast region includes: Data standardization was performed on Doppler weather radar data, and cloud edge features were extracted using a morphological gradient operator. The expression is as follows:

[0039] in For expansion operations, For erosion calculation, It is a 1*1 flat plate structural element. This is the raw radar data. The morphological gradient of the raw radar data; Select 3 groups of structural elements of different sizes , , These correspond to cloud features at small, medium, and large scales, respectively. The results of weighted fusion of filtering at each scale are expressed as follows:

[0040] in The result of the opening-closing operation. The weighting coefficients for the first scale. The weighting coefficients for the second scale. The weighting coefficients for the third scale. This is the result of the opening-closing operation at the first scale. This is the result of the opening-closing operation at the second scale. This is the result of the opening-closing operation at the third scale; Multi-scale filtering is applied to the morphological gradient image to eliminate peak features with significant cloud contours while preserving them. A watershed transform is then performed on the filtered gradient image, dividing it into several initial regions. The merging cost of adjacent regions is calculated using the difference in mean grayscale value and area as indicators.

[0041] in Let V be the average reflectance of the v-th region. Let be the average reflectance of the r-th region. Let v be the area of ​​the v-th region. Let r be the area of ​​the r-th region. The cost of merging the v-th region and the r-th region; The merging priority is managed using a region adjacency graph and a min-heap data structure. Each time, the region pair with the lowest global cost is merged until the preset number of regions or cost threshold is met. Regions with smaller areas and similar gray levels are merged first. Multidimensional features were calculated for the merged cloud region to establish a cloud attribute database. The multidimensional features include mean reflectivity, echo top height, volume, and vertical cumulative liquid water content. Based on the operation location and cloud movement trajectory, cloud units containing catalytic operation points and meeting preset conditions are marked as the influence area; cloud units with the highest feature similarity in the upwind direction of the influence area are selected as the comparison area; the similarity is calculated by dynamic time warping algorithm to match the life history characteristics of the cloud in the influence area before the operation; the life history characteristics are the volume change trend and the evolution of vertical cumulative liquid water content. In the actual assessment, the preset conditions were reflectivity >25dBZ and top height >5km.

[0042] In this embodiment, the method for obtaining the catalytic operation-affected zone includes: Spatiotemporal registration is performed on Doppler radar reflectivity data at continuous time intervals to unify spatial resolution and time intervals, and range attenuation correction and elevation angle normalization are applied. The affected area and the comparison area are respectively used as the region of interest (ROI) for tracking. Noise points and non-cloud echoes outside the region are removed, and the data within the region of interest is Gaussian smoothed. Set time The radar echo matrix is ,time The matrix is For each grid point within the ROI Define size as × template window ,exist Slide within the search window to calculate the cross-correlation coefficient:

[0043] in For a sub-region within the search window, This is the reflectivity value of the template. The reflectance value of the sub-region. For local pixel indexes within the window, To explore the pixel values ​​of a sub-region of the window, For grid point coordinates, It is a displacement vector. This is the cross-correlation coefficient between the template window and the sub-region at the displacement within the search window; Define the confidence index as follows:

[0044] in The highest cross-correlation coefficient, The second largest correlation coefficient, The confidence index is used to evaluate the motion vector. ; An improved Lucas-Cannard method is used, assuming that the radar echo grayscale value is conserved over a short period of time. Spatiotemporal differentiation is performed on the reflectivity data to construct the optical flow equation:

[0045] in For time gradient, For the lateral spatial gradient, For the longitudinal spatial gradient; The initial optical flow vector is obtained by solving the overdetermined equations in a 3×3 neighborhood using the least squares method. ; Introducing multi-scale pyramid optical flow: Gaussian pyramid downsampling is performed on the echo data, optical flow is calculated starting from the top layer, and the result is used as the initial value for the lower layer. Accumulated error is reduced through iterative optimization, and the optical flow motion vector field of each grid point in the ROI is output. Based on the complementary characteristics of radar echo tracking and optical flow methods, dynamic weighting coefficients are designed to fuse optical flow motion vectors and motion vectors, with the following expression:

[0046] in To fuse motion vectors, For dynamic weighting coefficients, The magnitude of the optical flow gradient; Construct a global optimization objective function that minimizes the following energy term, expressed as:

[0047] in Let be the objective function. For the neighborhood, For smoothing coefficients, Grid points to be optimized Motion vector, The fused motion vector of the grid points. For the neighborhood points of the grid point, The motion vector of the neighboring points, For data fidelity items, For smoothing terms; The optimized motion vector field is subjected to spatiotemporal filtering, which includes spatial smoothing and temporal smoothing. Spatial smoothing uses median filtering to remove isolated outlier vectors. Temporal smoothing is a weighted average of the motion vectors at three consecutive time points. Based on the optimized motion vector field, a region growing method is used to track cloud clusters: the cloud cluster outline of the initial influence area is used as the seed region; according to the motion vector... Predict the location of the cloud grid points at the next moment; perform morphological dilation on the predicted region and compare it with... The radar echoes at any given time are matched to obtain the overlap. When the overlap is greater than 70%, the match is successful. The trajectory of the cloud formation as it evolves over time is generated iteratively, which is the dynamic range of the catalytic operation's influence area. In actual assessments, noise points outside the region and non-cloud echoes (i.e., reflectivity) are less than 15 dBZ. The search window is 2-3 times the template window. The Gaussian pyramid has 4 layers, with the top layer being low resolution and the bottom layer being high resolution. The median filter uses a 3×3 window. Multi-scale morphology was used to segment Doppler radar data into cloud units, and a cloud attribute database was established. Adaptive spatiotemporal tracking was performed on the segmented cloud units based on Doppler radar data to obtain the cloud's motion trajectory and life history characteristics. Combining the operation location and cloud motion trajectory, cloud units containing the catalytic operation point and meeting preset conditions were marked as the initial influence area. The spatiotemporal evolution range of the initial influence area, i.e., the catalytic operation influence area, was obtained through cloud tracking.

[0048] In this embodiment, the method for obtaining the comparison cloud body includes: The three-dimensional structural features of influencing clouds and potential contrasting clouds are extracted from Doppler radar data, and a multi-dimensional feature vector is constructed, expressed as follows:

[0049] in For a moment The multidimensional feature vectors, For a moment echo peak height, For a moment echo volume, This is the transpose of the vector. For vertically accumulated liquid water content, The density of liquid water, This is the base value in the vertical direction. This represents the upper limit in the vertical direction; Select a time window before the start time of the task. By sampling according to the radar data temporal resolution, we obtained the time series of influencing cloud features and the time series of potential comparative cloud features:

[0050]

[0051] in The start time of the operation. The time window length, To influence cloud feature sequences, For the k-th potential contrast cloud feature sequence, To influence the time-by-time feature values ​​of the cloud feature sequence, Let be the time-series feature value of the k-th candidate cloud feature sequence. The radar data sampling interval; the length of the cloud feature sequence is affected by... The potential contrast cloud feature sequence length is ; Zero-mean standardization of features is used to screen candidate clouds that meet spatial conditions from the segmented cloud bodies. Spatial conditions include location constraints and range constraints. Location constraint: located upwind of the influencing cloud, with a horizontal distance of ≥20km. Range constraint: located in the same radar scanning area as the influencing cloud, and with an elevation difference of ≤2km between the cloud body centers. Calculate the static feature similarity between the candidate cloud and the influencing cloud at the start of the time window:

[0052] in The static feature similarity between candidate clouds and influencing clouds. Let be the number of dimensions of the static features. For at any time The m-th static eigenvalue of the cloud is affected. For at any time The m-th static feature value of the h-th candidate cloud; Candidate clouds with static feature similarity greater than or equal to a similarity threshold are retained. A multidimensional cost matrix is ​​defined based on the standardized impact cloud feature sequence and the candidate cloud feature sequence, expressed as:

[0053] in Let m be the dynamic feature value of the m-th dimension at time t in the standardized cloud sequence. Let tu be the dynamic feature value at time m in the standardized candidate cloud sequence. The optimal time alignment path is found using dynamic programming to minimize the cumulative distance between the two sequences.

[0054] in For from (1,1) to The total cumulative distance of the optimal path, For path selection operators, for The cost matrix; Boundary conditions are , Introduce the Sakoe-Chiba window constraint, allowing windows to be opened only on both sides of the diagonal. Search path within the range; where ; Dynamic time warping distance Define dynamic similarity: For candidate cloud computing dynamic similarity, the cloud with the highest dynamic similarity is selected as the initial comparison cloud. If the dynamic similarity is less than the comparison similarity threshold, it is determined that there is no qualified comparison cloud, and the cloud search range is expanded or the time window is adjusted. The matching quality is verified by trend consistency and morphological similarity. Trend consistency: Calculate the Pearson correlation coefficient of each characteristic time series of the influencing cloud and the comparison cloud. If the echo top height, volume, and vertical cumulative liquid water content are all greater than or equal to 0.7, the result is positive. Morphological similarity: Compare the three-dimensional morphological changes of the influencing cloud and the comparison cloud within the time window. Based on the dynamic time warping optimal path, the comparison cloud feature sequence is aligned with the influencing cloud sequence to output the comparison cloud body; In actual assessments, echo volume ,in , , For radar data spatial resolution; The time is 30 minutes, adjusted according to the cloud's life cycle. Based on the spatiotemporal resolution of radar echoes, 70% was selected as the similarity threshold for cloud cluster matching. If the cloud cluster is a strong convective cloud, the similarity threshold was adjusted to 80%. The dynamic similarity value ranges from 0 to 1. The closer the value is to 1, the higher the similarity of cloud features. In combination with the cloud feature matching requirements of artificial rain enhancement operations, 0.5 is selected as the comparison similarity threshold. If there are few cloud resources in the area, the threshold is lowered to 0.4. The cloud unit with the highest feature similarity in the upwind direction of the catalytic operation influence area was selected as the comparison area; Based on surface precipitation data, the rain enhancement effect of this operation was statistically tested. Combined with radar echo images, it was determined that the target clouds were primarily stratocumulus-mixed clouds suitable for cloud seeding and rain enhancement operations. Aircraft rain enhancement operations were conducted in Panjin on April 25th. Daily rainfall was used as a statistical variable to test the rain enhancement effect of this aircraft operation. Based on aircraft operation information, the operational sample was determined to be the regional average daily rainfall on April 25th, and the historical sample was the regional average daily rainfall of rainy days from April 20th to April 30th, 1961-1990. Combined with radar echo images, it was determined that the target clouds were primarily stratocumulus-mixed clouds, and the echoes showed multiple strong centers. Considering the diffusion of the aircraft catalyst itself, this analysis selected the area within the dashed frame as the operational area and the area within the solid frame as the control area. Figure 2 As shown.

[0055] In this embodiment, the method for constructing the ZI relational database includes: Quality control was performed on the raindrop spectrometer observation data, outlier particle diameter values ​​were removed, missing spatiotemporal data were supplemented, and the spatiotemporal scale of the raindrop spectrum samples was matched to the observation scale of the Doppler radar to obtain standardized raindrop spectrum samples. Given the radar reflectivity factor and rainfall intensity, the expression is:

[0056]

[0057] in Let be the diameter of the precipitation particles at level i. The inner diameter per cubic meter The number of particles, For reflectivity factor, Let be the mass of the i-th order raindrop. Let be the final velocity of the i-th order raindrop. Rainfall intensity; The median volume diameter was obtained, and the droplet spectral samples were classified by rainfall intensity and median volume diameter. The samples were divided into rainfall intensity categories and median diameter categories. The rainfall intensity categories were divided into 16 categories with unequal intervals, covering the main precipitation intensity range of 0~5 mm / h. The median diameter categories were divided into small, medium and large droplets. Drop spectrum samples falling within the same grid are averaged to generate 48 distinct average drop spectra based on averaging rules. Averaging rule: Take the arithmetic mean of the drop spectrum data of all samples within the same grid. For each type of average drop spectrum, the least squares method is used for fitting. coefficients in Sum of Indices Based on the average droplet spectrum, the radar reflectivity factor and rainfall intensity are calculated, and an objective function is constructed. ,in The rainfall intensity for the k-th sample is calculated from the drop spectrum. The radar reflectivity factor of the k-th sample is calculated from the drop spectrum; ZI relationships are classified and labeled according to precipitation cloud system types to form localized ZI relationships; precipitation cloud system types include stratiform clouds, cumulus clouds, and mixed stratocumulus clouds; In actual assessment, the median diameter is divided into small, medium, and large droplets, with small droplets being greater than 0.812 mm, medium droplets being greater than or equal to 0.812 mm and less than 0.937 mm, and large droplets being greater than or equal to 0.937 mm. The particle diameter anomaly is: Less than 0.1mm or Greater than 8mm.

[0058] In this embodiment, the expressions for the absolute rainfall increase and the relative rainfall increase are:

[0059]

[0060] in This refers to the absolute increase in rainfall. This refers to the relative increase in rainfall. To account for the total precipitation in the affected area, The total precipitation in the comparison area; In actual assessment, high spatiotemporal resolution precipitation data is used as the basic data for statistical verification. After being aggregated to a spatiotemporal scale consistent with radar-estimated precipitation, the total precipitation in the affected area and the comparison area is calculated separately, thereby obtaining the absolute and relative rainfall increases, which serve as the core indicators for statistical verification.

[0061] In this embodiment, the method for physical examination of the comparison cloud includes: Based on radar tracking, time-series curves of echo top height, echo volume, maximum reflectivity, and vertical cumulative liquid water content of clouds in the affected and control areas were plotted. The effectiveness of the catalytic operation was physically verified based on five physical quantities: echo top height, echo volume, maximum reflectivity, vertical cumulative liquid water content, and precipitation flux. Echo top height refers to the echo top height of a given dBZ threshold cell, used to identify radar echoes. The echo top height of a cell refers to the maximum height of radar echoes greater than or equal to the dBZ threshold. Echo volume refers to the echo volume of a given dBZ threshold cell. Maximum reflectivity refers to the maximum dBZ value within a given dBZ threshold cell. Vertical cumulative liquid water content refers to the vertical cumulative liquid water content calculated from the maximum reflectivity factor of each layer from the zero-degree layer to the top of the cell within a given dBZ threshold cell. Precipitation flux refers to the sum of precipitation flux calculated from the vertical maximum reflectivity factor of each cell in the horizontal projection within a given dBZ threshold cell. Plot the curves of the changes of five physical quantities over time for the operational cloud unit and the comparison cloud unit respectively, identify the operational information, and compare and analyze the time series changes of physical quantities in the operational cloud unit and the comparison cloud unit. In practical assessments, the vertical cumulative liquid water content is determined as follows: after identifying a given dBZ threshold, for each identified unit cell, the vertical direction is measured from the zero-degree layer. Initially, find the grid point with the largest dBZ value in this unit cell. Based on the relationship between water density M and reflectivity factor Z (M-Z), convert it to the liquid water content per unit grid. Then, find the unit cell layer upwards using the same method. To obtain the vertically cumulative liquid water content, the liquid water content of each layer from the zero-degree layer height to the top layer height is accumulated. The specific calculation steps are as follows: a. At any layer height within the identified unit cell, find the grid point with the largest dBZ value for that layer.

[0062] b. Convert the maximum dBZ value into a reflectivity factor: ; c. Calculate the density of water using the reflectivity factor Z: ; d. Integrate the density M of each water layer within the unit cell in the vertical direction from the zero-degree layer height to the top of the layer: ; e. Convert units to (kg / m2): In practical applications, all reflectance values ​​greater than 55dBZ are taken as 55dBZ; Precipitation flux: a. For each cell identified by a given dBZ threshold, find the maximum dBZ value in the vertical direction of the cell point by point on the horizontal projection plane grid, and convert the maximum dBZ value into a reflectivity factor: ; b. Calculate the precipitation rate based on the empirical relationship (Z—R) between the reflectance factor Z and the precipitation rate R: ; c. Convert the units to (m / s); d. Multiply the precipitation rate R by the grid area. The precipitation flux over the grid area is obtained: Unit: m 3 / s; e. Add up the precipitation fluxes at all grid points on the horizontal projection surface of the unit cell to obtain the total precipitation flux of the unit cell; Aircraft physical inspection: The intuitive comparative physical inspection method was adopted to compare the radar parameters of the affected area and the comparison area before and after the operation. After the operation, the radar echo intensity of the affected area showed a continuous upward trend, and the operation was effective. Using the radar combined reflectivity regional dynamic comparative analysis method (K-value method), where K = the ratio of the average radar combined reflectivity of the operational area to the average radar combined reflectivity of the comparison area, the radar combined reflectivity of the operational area and the comparison area were compared and calculated. Before the operation, k = 0.93; after the operation, k = 1.02; 60 minutes after the operation, k = 1.30; 120 minutes after the operation, k = 1.14; 180 minutes after the operation, k = 1.30; and 240 minutes after the operation, k = 1.19. After the operation, the radar echo intensity increased in both the affected area and the comparison area, with the increase being slightly higher in the affected area than in the comparison area.

[0063] In this embodiment, the statistical test method includes: The area affected by the catalytic operation was designated as the treatment group, and the dynamically matched control area was designated as the control group. The time window covered 1 hour before the operation, 0-2 hours during the operation, and 1 hour after the operation. High spatiotemporal resolution precipitation data were aggregated into a spatiotemporal scale consistent with radar-estimated precipitation. Calculate the absolute and relative rainfall increases during the time period of the catalytic influence of cloud bodies in the affected area and the control area based on radar tracking; Calculate the Pearson correlation coefficient between the precipitation sequences of the affected area and the control area before the operation. If the Pearson correlation coefficient is greater than 0.6, the correlation requirement is met. Use the two independent samples t test to verify whether there is no significant difference between the mean precipitation of the two areas before the operation. If the p value of the t test is greater than 0.05, the mean difference test is passed. When precipitation data from the affected area and the control area are independent, a two-sample t-test is used to verify whether the rainfall increase is significant: Null hypothesis : Alternative hypothesis : ; Calculate statistics : ; in The average precipitation in the affected area. To compare the average precipitation in the region, The sample size of the affected area. The sample size for the comparison region. To represent the variance of precipitation samples in the affected area, The variance of precipitation samples in the comparison area; When there is a spatiotemporal correlation between the affected area and the control area, a paired t-test is used: null hypothesis : Alternative Hypothesis : ; Calculate the statistics: ,in This represents the average absolute rainfall increase. The standard deviation of absolute rainfall increase This represents the number of paired samples; Monte Carlo simulation test: Precipitation data from the catalytic operation's impact area and the control area were merged to form a total sample set. Subsamples with the same sample size as the original impact area and control area were randomly selected from the total sample set to serve as virtual impact area and virtual control area, respectively. Virtual rainfall increase was calculated, and the probability distribution of virtual rainfall increase was generated by repeated resampling. The 95% confidence interval was calculated using the percentile method. If the actual absolute rainfall increase is greater than the upper quantile of the 95% confidence interval, the rainfall increase effect is significant. Statistical test of relative rainfall increase rate: For relative rainfall increase, take the logarithm of the precipitation data to transform the relative difference into an absolute difference, and construct a statistical measure. ,in To determine the standard deviation of logarithmic precipitation, a standard normal distribution test is performed on the rainfall increase. A significance level of 0.05 is set, and the one-sided critical value is obtained from the standard normal distribution table. If the calculated statistic is greater than 1.645 and... If the value is greater than 0, the rain enhancement effect is considered significant at a 95% confidence level. After aligning the cloud bodies with the ZI relation database in time and space, the time series and spatial range of the cloud body physical quantities of the physical test are simultaneously matched to the time and space scale of the precipitation data of the statistical test. In actual assessments, the lower 2.5% quantile and the upper 97.5% quantile are within the 95% confidence interval. The sampling was repeated 1000 times. The daily rainfall data of the provincial A1 national-level ground meteorological station were used to statistically test and analyze the effect of the ground rain enhancement operation. The time series of the daily rainfall dataset used is from 1961 to 2015, and the data format is standard and general. After determining the stations in the operation's impact area, region B11 was selected as the comparison area based on the selection principles. The reasons are: a. The comparison area is located on the crosswind side of the operation's impact area, more than 20km away, and is unaffected by the catalytic operation; b. The topography and area of ​​the comparison area are roughly similar to those of the operation's impact area; c. The comparison area and the operation's impact area are affected by the same or similar weather systems; d. The precipitation in the comparison area and the operation's impact area shows a good correlation. Aircraft statistical test: The rain enhancement effect of this operation was analyzed based on daily rainfall data using the regional historical regression statistical test method. The absolute rainfall increase of Xinzhou 60-B3435 was 0.13 mm, the relative rainfall increase rate was 11%, and the precipitation increased by 0.116 billion cubic meters. The significance level was less than 0.05. Ground statistical test: The rain enhancement effect of the operation was analyzed based on daily rainfall data using the regional historical regression statistical test method. The absolute rainfall increase was 0.05 mm and the relative rainfall increase rate was 8%, with a significance level of less than 0.05.

[0064] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for verifying the effect of artificial rain enhancement based on radar tracking and raindrop spectrum, characterized in that, Includes the following steps: Collect Doppler weather radar data and weather phenomenon instrument drop spectrum data of a preset area, and preprocess the Doppler weather radar data and the weather phenomenon instrument drop spectrum data. The Doppler weather radar data is segmented by multi-scale morphological cloud body to obtain the affected area and the comparison area. Based on the Doppler weather radar data, the cloud body is adaptively tracked to obtain the catalytic operation's affected area. Multidimensional dynamic similarity matching is used to dynamically match the catalytic operation influence area and the comparison area to obtain the comparison cloud body. A ZI relationship database is constructed based on the drop spectrum data of the weather phenomenon instrument. Radar-estimated precipitation and high spatiotemporal resolution precipitation are obtained based on the ZI relationship database. The absolute and relative rainfall increases are obtained by spatiotemporally aligning the comparison cloud body with the ZI relation database. The comparison cloud body is then subjected to physical verification. Statistical verification is performed based on the absolute rainfall increase, the relative rainfall increase, and high spatiotemporal resolution precipitation, and the verification results are output.

2. The method for verifying the effect of artificial rain enhancement based on radar tracking and raindrop spectrum according to claim 1, characterized in that, The method for obtaining the affected region and the contrast region includes: Doppler weather radar data is standardized, and cloud edge features are extracted using morphological gradient operators. Select 3 groups of structural elements of different sizes , , These correspond to cloud features at small, medium, and large scales, respectively, and the filtering results at each scale are weighted and fused together. Multi-scale filtering is performed on the morphological gradient image to eliminate the peak features of the significant cloud contours while preserving them. Watershed transformation is then performed on the filtered gradient image to divide the image into several initial regions. The merging cost of adjacent regions is calculated using the difference in the mean gray level and the area of ​​the regions as indicators. The merging priority is managed using a region adjacency graph and a min-heap data structure. Each time, the region pair with the lowest global cost is merged until the preset number of regions or cost threshold is met. Regions with smaller areas and similar gray levels are merged first. Multidimensional features were calculated for the merged cloud region to establish a cloud attribute database. The multidimensional features include mean reflectivity, echo top height, volume, and vertical cumulative liquid water content. Based on the location of the operation and the trajectory of the cloud, cloud units that contain the catalytic operation point and meet the preset conditions are marked as the influence area; cloud units with the highest feature similarity in the upwind direction of the influence area are selected as the comparison area; the similarity is calculated by the dynamic time warping algorithm to match the life history characteristics of the cloud in the influence area before the operation; the life history characteristics are the volume change trend and the evolution of vertical cumulative liquid water content.

3. The method for verifying the effect of artificial rain enhancement based on radar tracking and raindrop spectrum according to claim 1, characterized in that, The method for obtaining the catalytic operation-affected zone includes: Spatiotemporal registration is performed on Doppler radar reflectivity data at continuous time intervals to unify spatial resolution and time intervals, and range attenuation correction and elevation angle normalization are applied. The affected area and the comparison area are respectively used as the region of interest (ROI) for tracking. Noise points and non-cloud echoes outside the region are removed, and the data within the region of interest is Gaussian smoothed. Set time The radar echo matrix is ,time The matrix is For each grid point within the ROI Define size as × template window ,exist Slide within the search window to calculate the cross-correlation coefficient; A confidence index is defined, and a modified Lucas-Cornard method is used, assuming that the radar echo grayscale value is conserved over a short period of time. Spatiotemporal differentiation is performed on reflectivity data to construct optical flow equations; Introducing multi-scale pyramid optical flow: Gaussian pyramid downsampling is performed on the echo data, optical flow is calculated starting from the top layer, and the result is used as the initial value for the lower layer. Accumulated error is reduced through iterative optimization, and the optical flow motion vector field of each grid point in the ROI is output. Based on the complementary characteristics of radar echo tracking and optical flow methods, dynamic weighting coefficients are designed, optical flow motion vectors and motion vectors are fused, and a global optimization objective function is constructed to minimize the following energy terms; The optimized motion vector field is subjected to spatiotemporal filtering, which includes spatial smoothing and temporal smoothing. Spatial smoothing uses median filtering to remove isolated outlier vectors. Temporal smoothing is a weighted average of the motion vectors at three consecutive time points. Based on the optimized motion vector field, a region growing method is used to track cloud clusters: the cloud cluster outline of the initial influence area is used as the seed region; according to the motion vector... Predict the location of the cloud grid points at the next moment; perform morphological dilation on the predicted region and compare it with... The radar echoes at any given time are matched to obtain the overlap. When the overlap is greater than 70%, the match is successful. The trajectory of the cloud formation as it evolves over time is generated iteratively, which is the dynamic range of the catalytic operation's influence area.

4. The method for verifying the effect of artificial rain enhancement based on radar tracking and raindrop spectrum according to claim 1, characterized in that, The method for obtaining the comparison cloud body includes: Three-dimensional structural features of influencing clouds and potential contrasting clouds were extracted from Doppler radar data to construct multi-dimensional feature vectors. A time window prior to the operation was selected, using the start time of the operation as a benchmark. By sampling radar data at the time resolution, time series of influencing cloud features and time series of potential comparative cloud features are obtained. Zero-mean standardization of features is used to screen candidate clouds that meet spatial conditions from the segmented cloud bodies. Spatial conditions include location constraints and range constraints. Location constraint: located upwind of the influencing cloud, with a horizontal distance of ≥20km. Range constraint: located in the same radar scanning area as the influencing cloud, and with an elevation difference of ≤2km between the cloud body centers. Calculate the static feature similarity between candidate clouds and impact clouds at the start of the time window, retain candidate clouds with static feature similarity greater than or equal to the similarity threshold, and define a multidimensional cost matrix based on the standardized impact cloud feature sequence and candidate cloud feature sequence. The optimal time alignment path is found through dynamic programming to minimize the cumulative distance between the two sequences; Boundary conditions are , Introduce the Sakoe-Chiba window constraint, allowing windows to be opened only on both sides of the diagonal. Search path within the range; where ; Dynamic time warping distance Define dynamic similarity: For candidate cloud computing dynamic similarity, the cloud with the highest dynamic similarity is selected as the initial comparison cloud. If the dynamic similarity is less than the comparison similarity threshold, it is determined that there is no qualified comparison cloud, and the cloud search range is expanded or the time window is adjusted. The matching quality is verified by trend consistency and morphological similarity. Trend consistency: Calculate the Pearson correlation coefficient of each characteristic time series of the influencing cloud and the comparison cloud. If the echo top height, volume, and vertical cumulative liquid water content are all greater than or equal to 0.7, the result is positive. Morphological similarity: Compare the three-dimensional morphological changes of the influencing cloud and the comparison cloud within the time window. Based on the optimal path of dynamic time warping, the feature sequence of the comparison cloud is aligned with the sequence of the influencing cloud, and the comparison cloud body is output.

5. The method for verifying the effect of artificial rain enhancement based on radar tracking and raindrop spectrum according to claim 1, characterized in that, The method for constructing the ZI relational database includes: Given the radar reflectivity factor and rainfall intensity, the expression is: ； ； in Let be the diameter of the precipitation particles at level i. The inner diameter per cubic meter The number of particles, For reflectivity factor, Let be the mass of the i-th order raindrop. Let be the final velocity of the i-th order raindrop. Rainfall intensity; The median volume diameter was obtained, and the droplet spectral samples were classified by rainfall intensity and median volume diameter. The samples were divided into rainfall intensity categories and median diameter categories. The rainfall intensity categories were divided into 16 categories with unequal intervals, covering the main precipitation intensity range of 0~5 mm / h. The median diameter categories were divided into small, medium and large droplets. Drop spectrum samples falling within the same grid are averaged to generate 48 distinct average drop spectra based on averaging rules. Averaging rule: Take the arithmetic mean of the drop spectrum data of all samples within the same grid. For each type of average drop spectrum, the least squares method is used for fitting. coefficients in Sum of Indices Based on the average droplet spectrum, the radar reflectivity factor and rainfall intensity are calculated, and an objective function is constructed. ,in The rainfall intensity for the k-th sample is calculated from the drop spectrum. The radar reflectivity factor of the k-th sample is calculated from the drop spectrum; ZI relationships are classified and labeled according to precipitation cloud system types to form localized ZI relationships; precipitation cloud system types include stratiform clouds, cumulus clouds, and mixed stratocumulus clouds.

6. The method for verifying the effect of artificial rain enhancement based on radar tracking and raindrop spectrum according to claim 1, characterized in that, The expressions for the absolute rainfall increase and the relative rainfall increase are: ； ； in This refers to the absolute increase in rainfall. This refers to the relative increase in rainfall. To account for the total precipitation in the affected area, The total precipitation in the comparison area.

7. The method for verifying the effect of artificial rain enhancement based on radar tracking and raindrop spectrum according to claim 1, characterized in that, A method for physically examining the comparison cloud body includes: Based on radar tracking, time-series curves of echo top height, echo volume, maximum reflectivity, and vertical cumulative liquid water content of clouds in the affected and control areas were plotted. The effectiveness of the catalytic operation was physically verified based on five physical quantities: echo top height, echo volume, maximum reflectivity, vertical cumulative liquid water content, and precipitation flux. Echo top height refers to the echo top height of a given dBZ threshold cell, used to identify radar echoes. The echo top height of a cell refers to the maximum height of radar echoes greater than or equal to the dBZ threshold. Echo volume refers to the echo volume of a given dBZ threshold cell. Maximum reflectivity refers to the maximum dBZ value within a given dBZ threshold cell. Vertical cumulative liquid water content refers to the vertical cumulative liquid water content calculated from the maximum reflectivity factor of each layer from the zero-degree layer to the top of the cell within a given dBZ threshold cell. Precipitation flux refers to the sum of precipitation flux calculated from the vertical maximum reflectivity factor of each cell in the horizontal projection within a given dBZ threshold cell. The curves of the changes of five physical quantities over time were plotted for the operational cloud unit and the comparison cloud unit, respectively. The operational information was marked, and the time series changes of the physical quantities of the operational cloud unit and the comparison cloud unit were compared and analyzed.

8. The method for verifying the effect of artificial rain enhancement based on radar tracking and raindrop spectrum according to claim 1, characterized in that, The statistical test methods include: The area affected by the catalytic operation was designated as the treatment group, and the dynamically matched control area was designated as the control group. The time window covered 1 hour before the operation, 0-2 hours during the operation, and 1 hour after the operation. High spatiotemporal resolution precipitation data were aggregated into a spatiotemporal scale consistent with radar-estimated precipitation. Calculate the absolute and relative rainfall increases during the time period of the catalytic influence of cloud bodies in the affected area and the control area based on radar tracking; Calculate the Pearson correlation coefficient between the precipitation sequences of the affected area and the control area before the operation. If the Pearson correlation coefficient is greater than 0.6, the correlation requirement is met. Use the two independent samples t test to verify whether there is no significant difference between the mean precipitation of the two areas before the operation. If the p value of the t test is greater than 0.05, the mean difference test is passed. When precipitation data from the affected area and the control area are independent, a two-sample t-test is used to verify whether the rainfall increase is significant: Null hypothesis : Alternative hypothesis : ; Calculate statistics : ； in The average precipitation in the affected area. To compare the average precipitation in the region, The sample size of the affected area. The sample size for the comparison region. To represent the variance of precipitation samples in the affected area, The variance of precipitation samples in the comparison area; When there is a spatiotemporal correlation between the affected area and the control area, a paired t-test is used: null hypothesis : Alternative Hypothesis : ; Calculate the statistics: ,in This represents the average absolute rainfall increase. The standard deviation of absolute rainfall increase This represents the number of paired samples; Monte Carlo simulation test: Precipitation data from the catalytic operation's impact area and the control area were merged into a total sample set. Subsamples with the same sample size as the original impact area and control area were randomly selected from the total sample set to serve as virtual impact area and virtual control area, respectively. Virtual rainfall increase was calculated, and the probability distribution of virtual rainfall increase was generated by repeated sampling. The 95% confidence interval of the virtual rainfall increase was calculated using the percentile method. If the actual absolute rainfall increase is greater than the upper quantile of the 95% confidence interval, the rainfall increase effect is significant. Statistical test of relative rainfall increase rate: For relative rainfall increase, take the logarithm of the precipitation data to transform the relative difference into an absolute difference, and construct a statistical measure. ,in To determine the standard deviation of logarithmic precipitation, a standard normal distribution test is performed on the rainfall increase. A significance level of 0.05 is set, and the one-sided critical value is obtained from the standard normal distribution table. If the calculated statistic is greater than 1.645 and... If the value is greater than 0, the rain enhancement effect is considered significant at a 95% confidence level. After aligning the cloud bodies with the ZI relation database in time and space, the time series and spatial range of the cloud body physical quantities from the physical test are simultaneously matched to the spatiotemporal scale of the precipitation data from the statistical test.

Images