Fog and rain droplet distribution-based flood discharge atomization rainfall intensity prediction method

Abstract

The invention discloses a flood discharge atomization rainfall intensity prediction method based on fog and rain liquid drop distribution. The method comprises the following steps: acquiring fog and rain images and environmental parameters of a flood discharge area; extracting effective round spot data of fog drops and raindrops; extracting droplet deformation characteristics, constructing a fluid constraint matrix, enhancing deformation critical characteristics, decoupling flow characteristics, reconstructing a physical deformation mode, evaluating droplet stability and correcting droplet volume; and calculating the spatial distribution of the rainfall intensity based on the corrected droplet particle size data. According to the method, through a technical route of combining deformation feature extraction, flow characteristic decoupling and physical modal reconstruction, accurate description of dynamic deformation of flood discharge atomization rainfall droplets under the action of water-gas turbulence in a complex environment is realized, and the prediction precision of rainfall intensity and range of a flood discharge atomization area is improved.

Description

Technical Field

[0001] The present invention belongs to the field of water conservancy and hydropower engineering, and in particular, relates to a method for predicting the rainfall intensity of flood discharge atomization based on the droplet distribution of fog and rain. Background Art

[0002] Flood discharge atomization is an important phenomenon in large-scale water conservancy projects. The droplet distribution formed during the flood discharge process directly affects the rainfall intensity in the surrounding area, and further affects the surrounding ecological environment and engineering safety. Accurately predicting the rainfall intensity distribution in the flood discharge atomization area is of great significance for the environmental impact assessment of water conservancy projects, the design of engineering structure protection, and flood control early warning in the downstream area. In addition, accurate rainfall intensity prediction can provide a scientific basis for water conservancy project scheduling, optimize the flood discharge plan, reduce the adverse effects of flood discharge on the downstream ecosystem and human activities, and realize the sustainable utilization of water resources.

[0003] Currently, the prediction of rainfall intensity in the flood discharge atomization area mainly relies on empirical formulas and simplified models. Traditional methods usually adopt empirical correlations of single factors or a small number of factors, such as estimating the rainfall intensity based on the simple functional relationship between water flow and flood discharge height. Some researchers use numerical simulation methods to construct a water-air two-phase flow model based on the Euler-Lagrange framework, but generally adopt simplified droplet shape assumptions, mostly describing droplets as spherical or ellipsoidal. In terms of droplet identification, traditional image processing techniques are commonly used, such as edge detection and region segmentation algorithms to identify the droplet contour, and statistical methods are used to analyze the droplet size distribution. The conversion of rainfall intensity mainly relies on the Marshall-Palmer relationship or similar empirical formulas to simply map the droplet size distribution to the rainfall intensity value, ignoring the influence of complex environmental factors.

[0004] The existing technologies face three key technical challenges. First, in terms of the analysis of droplet deformation characteristics, traditional methods lack an effective mechanism to introduce physical constraints into feature decomposition, resulting in the extracted deformation characteristics being disconnected from the actual fluid physical behavior. Especially in the flood discharge environment with high turbulence intensity and severe droplet deformation, the hydrodynamic constraints and data characteristics cannot be organically combined, resulting in the inconsistency between the deformation feature analysis and the physical process, affecting the quality of the basic data for subsequent rainfall intensity prediction. Second, the existing droplet volume estimation methods are based on static geometric assumptions, ignoring the dynamic change law of the volume distribution during the deformation process of droplets. Especially when the droplet approaches the critical deformation state, the simple geometric model cannot accurately describe the actual volume of the droplet, resulting in systematic deviation in the conversion from particle size to rainfall intensity. Third, traditional rainfall intensity calculations generally adopt a fixed binning scheme, lacking the ability to adaptively adjust the binning strategy according to the dynamically changing droplet distribution characteristics, and at the same time, the non-uniformity of droplet spatial distribution is not fully considered during the density estimation process, making it difficult to accurately count the rainfall intensity spatial gradient distribution in complex flood discharge environments, reducing the spatial resolution and accuracy of the rainfall intensity prediction results. Summary of the Invention

[0005] Objective of the invention: To provide a method for predicting the rainfall intensity of flood-discharge atomization based on the droplet distribution of fog and rain, so as to solve at least one technical problem existing in the prior art.

[0006] Technical solution: A method for predicting the rainfall intensity of flood-discharge atomization based on the droplet distribution of fog and rain includes the following steps:

[0007] Collect the original droplet spectrum images and environmental sensing data at the flood-discharge site, perform preprocessing to obtain preprocessed images, and an environmental characteristic parameter set including wind speed, wind direction, temperature, humidity, atmospheric pressure, turbulence pulsation characteristics, and turbulence intensity;

[0008] Perform patch processing on the preprocessed images, obtain patches of a predetermined type, and convert them into standardized effective patches to form an effective circular patch dataset of raindrops;

[0009] Based on the effective circular patch dataset and the environmental characteristic parameter set, adopt the adaptive deformation feature decomposition method with fluid physical constraints to perform droplet deformation dynamics modeling and particle size correction to obtain corrected three-dimensional droplet particle size data;

[0010] Perform full-scale adaptive dynamic binning statistics and droplet volume calculation on the corrected three-dimensional droplet particle size data to form a rain intensity spatial distribution matrix, that is, the rain intensity prediction result.

[0011] Beneficial effects: The present invention realizes the precise correspondence between deformation features and physical mechanisms and the precise correction of droplet volume, can capture the volume distribution changes brought by complex non-linear deformations, eliminates the systematic deviation in traditional methods, and at the same time improves the accuracy and spatial resolution of rain intensity prediction. Description of the drawings

[0012] Figure 1 It is a flowchart of the steps of a method for predicting the rainfall intensity of flood-discharge atomization based on the droplet distribution of fog and rain provided by an embodiment of the present application.

[0013] Figure 2 It is a flowchart of the steps for performing droplet deformation dynamics modeling and particle size correction provided by an embodiment of the present application.

[0014] Figure 3 It is a flowchart of the steps for constructing a fluid constraint matrix provided by an embodiment of the present application.

[0015] Figure 4 It is a flowchart of the steps for obtaining an enhanced deformation feature set provided by an embodiment of the present application. Detailed implementation manners

[0016] To enable those skilled in the art to better understand the solution of the present invention, the technical solutions in the embodiments of the present invention will be clearly and completely described below in conjunction with the accompanying drawings in the embodiments of the present invention. Obviously, the described embodiments are only a part of the embodiments of the present invention, rather than all of the embodiments. All other embodiments obtained by those of ordinary skill in the art based on the embodiments of the present invention without creative efforts shall fall within the scope of protection of the present invention.

[0017] It should be particularly noted that, for clearly showing the step flow of the present application, serial numbers are marked for each step in the specification. These serial numbers are only for the convenience of description and do not limit the execution order of the steps. In actual operation, according to the technical requirements of specific implementation scenarios, the steps can be executed in an order different from that shown in the specification, and in some cases, parallel processing between steps can also be achieved.

[0018] As Figure 1 shown, a method for predicting the intensity of flood-discharge atomization rainfall based on the distribution of fog and rain droplets includes the following steps:

[0019] S1. Collect the original droplet spectrum image and environmental sensing data at the flood-discharge site, perform preprocessing, and obtain the preprocessed image and the environmental feature parameter set;

[0020] Specifically, the original droplet spectrum image is an image related to the raindrops, fog droplets, and water vapor generated at the outlet of the flood-discharge building during the flood-discharge process, which is used to analyze the morphological, kinematic, and turbulent characteristics of the droplets, mainly manifested as the intensity and range of atomization rainfall. The environmental sensing data includes wind speed, wind direction vector value, temperature, humidity, atmospheric pressure, turbulent pulsation characteristics, and turbulent intensity, where the atmospheric pressure refers to the air pressure in the environment where the droplets are located, that is, the external pressure received by the droplets when moving in the air.

[0021] S2. Perform environmental noise recognition and droplet patch processing on the preprocessed image, obtain fuzzy circular patches, overlapping circular patches, and irregular patches, and convert them into standardized effective patches to form an effective circular patch dataset of the fog and rain droplet spectrum;

[0022] Specifically, analyze the "cluttered areas" (noise) in the image and identify those parts that may be patches. Process these patches to make them clear. According to the shape and characteristics, the patches are divided into three categories: fuzzy circular patches: circular patches with blurred edges and unclear shapes; overlapping circular patches: multiple circular patches overlapping; irregular patches: patches with irregular shapes. These patches are uniformly converted into "standardized" patches, that is, making them into a unified format and clear shape.

[0023] S3. Based on the effective circular spot data set and the environmental characteristic parameter set, an adaptive deformation feature decomposition method with fluid physical constraints is adopted to conduct droplet deformation dynamics modeling and particle size correction, and corrected three-dimensional droplet particle size data is obtained;

[0024] Specifically, the laws of fluid physics are used to analyze and decompose the morphology of the droplets. The decomposition method will be dynamically adjusted according to specific situations to ensure the accuracy of the analysis results. Given that during the free fall of the droplets, their morphology and trajectory are affected by environmental factors such as wind fields, air pressure, and temperature, and finally deform, split, and break, etc., it is necessary to correct the droplet morphology for these influencing factors. Further, a kinetic model is constructed through the dynamic changes of the droplets in multiple physical fields to predict the deformation behavior of the droplets under different conditions. Based on the model and the analysis results, the droplet particle size data is corrected to obtain more accurate three-dimensional particle size information, providing reliable data support for rainfall intensity prediction.

[0025] S4. Combining the environmental characteristic parameter set, full-scale adaptive dynamic binning statistics and droplet volume calculation are performed on the corrected three-dimensional droplet particle size data to form a rainfall intensity spatial distribution matrix, that is, the rainfall intensity prediction result; and it is verified through typical sample data, and the correction parameters in step S3 are iterated.

[0026] Specifically, the corrected three-dimensional droplet particle size data is grouped or classified and statistically analyzed according to the dynamic changes, which can clearly show the distribution characteristics of droplets of different sizes. According to the size of each particle size, the volume of the corresponding droplet is calculated to quantify the total amount of raindrops.

[0027] In this embodiment, an adaptive deformation feature decomposition method with fluid physical constraints is used for droplet deformation dynamics modeling and particle size correction, which can obtain more accurate droplet particle size data and reduce the errors in traditional methods; the corrected three-dimensional droplet particle size data is combined with the environmental characteristic parameter set for dynamic binning statistics, which can not only timely reflect the changes in rainfall intensity, but also form a rainfall intensity spatial distribution matrix through droplet volume calculation, and can dynamically capture the rainfall intensity changes during the flood discharge process and provide timely predictions; it realizes the accurate characterization of droplet dynamic deformation and improves the rainfall intensity prediction accuracy in the flood discharge atomization area.

[0028] According to one aspect of the present application, the preprocessing steps include:

[0029] S11. Collection of drop size image and environmental data: Drop size test papers and environmental sensors are preset at measuring points around the flood discharge structure to collect the original drop size images and environmental sensing data (wind speed data, wind direction data, temperature data, humidity data, turbulence intensity data). The sampling frequency of the environmental sensing data is 10 Hz, and the whole test process is continuously recorded.

[0030] S12. Image Index and Spatial Coordinate Association: Assign a unique digital index to each original raindrop spectrum image, record the sampling timestamp and three-dimensional spatial coordinates (x, y, z), and form an image index metadata table. Calculate the relative spatial position of each sampling point through the following formula: Relative position vector Pr = (x - x0, y - y0, z - z0) / H, where (x0, y0, z0) are the coordinates of the flood discharge outlet reference point, and H is the height of the flood discharge structure. Output the spatial sampling position data.

[0031] S13. Image Grayscale Conversion and Enhancement: Convert the original raindrop spectrum image (RGB format) into a grayscale image, and the conversion formula is: I_grayscale(x, y) = 0.299·R(x, y) + 0.587·G(x, y) + 0.114·B(x, y), where R(x, y), G(x, y), and B(x, y) are the red, green, and blue channel values at the pixel point (x, y) respectively.

[0032] S14. Adaptive Contrast Enhancement: Perform adaptive histogram equalization on the grayscale image, with a pixel block size of 32×32 and a contrast limit threshold of 2.0 to generate an enhanced image. The image block processing formula is: H_k(i) = N_k(i) / N_total, where H_k(i) is the frequency of the gray value i in the k-th block histogram, N_k(i) is the number of pixels with the gray value i in this block, and N_total is the total number of pixels in the block. After cropping and redistributing the histogram, eliminate the block boundary effect through bilinear interpolation and output the preprocessed image.

[0033] S15. Environmental Data Preprocessing and Feature Extraction: Denoise, detect outliers, and smooth the environmental sensing data, calculate the time average and standard deviation of each parameter, and extract the turbulence characteristic parameters. The turbulence intensity calculation formula is: Tu = σ_u / U_mean, where σ_u is the standard deviation of wind speed pulsation and U_mean is the average wind speed. The extracted data includes the average wind speed, average temperature, average humidity, and turbulence intensity. Output the environmental characteristic parameter set.

[0034] According to one aspect of the present application, the steps of forming an effective circular spot data set of the fog raindrop spectrum include:

[0035] S21. Feature Engineering and SVM Model Construction: Extract noise and plaque features based on the preprocessed image, including shape features (circularity, aspect ratio), edge features (average gradient, gradient direction variance, edge intersection points), spatial features (adjacent plaque density), frequency domain features (local spectral energy concentration), gray scale features (regional gray scale variance, edge-center gray scale difference), and fractal dimension. Construct an improved SVM recognition model, and the feature vector is represented as: F = [C_shape, R_aspect, G_average, G_variance, P_intersection, D_adjacent, E_spectrum, V_gray, D_edge-center, D_fractal]. The RBF is selected as the training SVM kernel function, the penalty parameter C = 1.0, and the kernel parameter γ = 0.1. Output the trained SVM model.

[0036] S22. Environmental Noise Detection and Filtering: Apply the trained SVM model to classify the regions of the preprocessed image and identify environmental noise regions (salt-and-pepper noise, Gaussian noise). The salt-and-pepper noise is processed by adaptive median filtering, and the window size is dynamically adjusted according to the noise density (from 3×3 to 7×7); the Gaussian noise is processed by bilateral filtering, with the spatial domain parameter σ_d = 3.0 and the value domain parameter σ_r = 30. After filtering, perform noise residue detection. If the PSNR value is lower than 28 dB, repeat the filtering operation. Output the denoised image.

[0037] S23. Plaque Type Identification and Classification: Apply the trained SVM model to identify the types of plaques in the denoised image, and classify the plaques into four categories: clear circular plaques, blurred circular plaques, overlapping circular plaques, and irregular plaques. Assign a label and a confidence score to each plaque to form a plaque type label set, and record the position coordinates and initial geometric features of each plaque at the same time. Output the plaque feature data set.

[0038] S24. Enhanced Processing of Blurred Circular Plaques: Extract the subset with the label of blurred circular plaques from the plaque feature data set, and apply the improved fuzzy C-means algorithm to each blurred circular plaque. Set the number of clustering centers C = 3 (representing the background, edge transition region, and plaque core region respectively), and the fuzzy factor m = 1.8. The membership degree calculation formula is: u_ij = 1 / ∑ k=1 C ((||x_i - v_j||) / (||x_i - v_k||)) 2 / (m-1) , where u_ij is the membership degree of pixel i belonging to cluster j, x_i is the gray scale value of pixel i, and v_j is the j-th clustering center. Perform edge enhancement and morphological closing operations on the membership degree distribution, and output the enhanced blurred circular plaque set.

[0039] S25. Coincident Circular Spot Separation Processing: Extract the subset labeled as coincident circular spots from the patch feature dataset, and apply the improved watershed algorithm to each coincident region. First, calculate the image gradient: G(x, y) = sqrt((▽_x I(x, y)) 2 + (▽_y I(x, y)) 2 ). Using the gradient image as the topographic map, determine the water source points by the H minimum simulation marking method: M_min ={(x, y) | H(x, y) ≤ H(x+i, y+j), for all (i, j) ∈ N(x, y)}. Perform the watershed flooding process, and extract the separated independent circular spots through 8-connected region detection. To solve the problem of over-segmentation, introduce a pre-segmentation merging criterion: Merge(R_i, R_j) = true, if Sim(R_i, R_j) > T_merge, where the similarity function Sim() is based on the regional gray distribution and boundary strength, and the threshold T_merge = 0.75. Output the separated circular spot set. Here, ▽_x and ▽_y are the gradients of the image in the x and y directions respectively; I(x, y) is the pixel value of the image at the coordinate (x, y); H(x, y) is the height function used to determine the segmentation boundary of the region; R_i and R_j represent two different patch regions.

[0040] S26. Irregular Patch Shape Correction: Extract the subset labeled as irregular patches from the patch feature dataset, and classify them into two categories according to the circularity value: high circularity (0.5 ≤ C_circularity < 0.7) and low circularity (C_circularity < 0.5). Apply the improved Hough circle transform detection to the irregular patches with high circularity, and set the accumulator threshold to 25% of the patch area; for the patches with low circularity, first perform morphological filling and edge smoothing, then apply the ellipse fitting algorithm to obtain the major axis a, minor axis b, and orientation angle θ, calculate the equivalent radius: R_equivalent = sqrt(a·b), and correct the shape through morphological operations to make the circularity reach ≥ 0.7. Output the corrected circular spot set.

[0041] S27. Standardized Circular Spot Synthesis and Verification: Merge the enhanced blurred circular spot set, separated circular spot set, corrected circular spot set with the original clear circular spots to form a standardized circular spot dataset. Calculate the center coordinates (x_c, y_c) and radius r for each circular spot, and construct a circular spot description matrix: C = {(x_c_1, y_c_1, r_1), (x_c_2, y_c_2, r_2),..., (x_c_n, y_c_n, r_n)}; perform circular spot validity verification, including edge continuity check, circularity verification, and minimum size screening, to ensure that all circular spots meet the statistical standards. Output the final valid circular spot dataset.

[0042] As Figure 2As shown, according to one aspect of the present application, the steps of performing dynamic modeling of atomized rainfall droplet deformation and particle size correction under the action of complex environmental factors to obtain corrected three-dimensional droplet particle size data include:

[0043] S31. Establish a mapping relationship from two-dimensional circular spots to three-dimensional droplets for the effective circular spot data set in combination with the environmental characteristic parameter set to obtain initial three-dimensional particle size data;

[0044] S32. Calculate the Weber number and deformation index of the droplets based on the initial three-dimensional particle size data to form an initial droplet deformation characteristic set;

[0045] S33. Apply a feature enhancement function to the initial droplet deformation characteristic set to amplify the critical deformation characteristics and obtain an enhanced deformation characteristic set;

[0046] S34. Apply a multi-scale flow field decomposition operator to the enhanced deformation characteristic set to separate the main flow traction, vortex-induced, and turbulent pulsation flow characteristics and obtain a flow decoupling characteristic set;

[0047] S35. Construct a fluid constraint matrix in combination with the environmental characteristic parameter set, including a flow field direction matrix, a scale-sensitive matrix, and a Reynolds number response matrix;

[0048] S36. Perform constraint feature decomposition in combination with the fluid constraint matrix and the flow decoupling characteristic set to reconstruct the physical deformation mode of the droplets and obtain a physical deformation mode set;

[0049] S37. Calculate the droplet stability index and volume correction coefficient based on the physical deformation mode set and output the corrected three-dimensional droplet particle size data.

[0050] Specifically, read the circular spot diameter data in the effective circular spot data set and the environmental characteristic parameter set to establish a preliminary mapping relationship from two-dimensional circular spots to three-dimensional droplets. The basic Karman formula is: D_3D_initial = k·d_2D·(ρ_a / ρ_w·U_ref / U_wind) 0.5 , where d_2D is the two-dimensional circular spot diameter, D_3D_initial is the initial three-dimensional droplet particle size, k is an empirical coefficient (0.6 - 0.8), ρ_a is the air density, ρ_w is the water density, U_wind is the real-time wind speed, and U_ref is the reference wind speed. Output the initial three-dimensional particle size data.

[0051] Analyze the initial three-dimensional particle size data and extract the estimated deformation characteristics of each droplet. First, calculate the Weber number of the droplet: We = ρ_a·U_rel 2 ·D_3D_initial / σ, where U_rel is the relative velocity and σ is the surface tension coefficient. Estimate the droplet deformation degree based on the Weber number and the Reynolds number (Re), and the initial deformation index is preliminarily estimated as: Δ_initial = α·(We) β ·(1 - e (-γ·We)) where α, β, and γ are fitting parameters. Output the initial droplet deformation feature set.

[0052] Combine the wind speed and turbulence intensity data in the environmental feature parameter set to construct the fluid physical constraint matrix W. The matrix construction formula is: W = Φ·S·R, where Φ is the flow field direction matrix constructed according to the unit vector of the mainstream direction; S is the scale-sensitive matrix that adapts to the deformation sensitivity of droplets of different particle sizes; R is the Reynolds number response matrix that reflects the differences in droplet deformation behavior in different Reynolds number ranges. The calculation of the Φ matrix involves local wind field decomposition: Φ = [e_x e_y e_z]·[w_x 0 0; 0 w_y 0; 0 0 w_z], where e_x, e_y, e_z are the unit vectors of the mainstream direction, and w_x, w_y, w_z are the weights in each direction. Output the fluid constraint matrix W.

[0053] Apply the feature enhancement function h(X) to the initial droplet deformation feature set to amplify the critical deformation features: h(X) = X + β·Γ(X)·X, where X is the initial droplet deformation feature set, β is the enhancement coefficient (0.2 - 0.5), and Γ(X) is the curvature-sensitive operator: Γ(X)_i = |▽ 2 X_i| / max(ε, ||▽X_i||), ▽ is the gradient operator, X_i is a specific eigenvalue in the initial droplet deformation feature set, and ε is a small constant to prevent division by zero (10 -6 ). Output the enhanced deformation feature set X'.

[0054] Apply the multi-scale flow field decomposition operator G to process the enhanced deformation feature set X': X'' = ∑ i G_i(X'), and the decomposition operator includes the mainstream traction deformation operator G_α, the vortex-induced deformation operator G_β, and the turbulent pulsation deformation operator G_γ. The mainstream traction deformation operator is calculated as: G_α(X') = Proj(X', e_flow)·(1 + η·|Tu - Tu_ref|), where e_flow is the mainstream direction vector, η is the turbulence correction coefficient (0.05 - 0.15), Tu is the current turbulence intensity, and Tu_ref is the reference turbulence intensity (0.1). Output the flow decoupling feature set X''.

[0055] Combine the fluid constraint matrix W and the flow decoupling feature set X'', and perform constraint feature decomposition: [U, Σ, V] = SVD(W·X''·W T ); Calculate the principal component weights: λ_i = u_i T ·W·X'', where u_i is the i-th column eigenvector of U, T is the transpose. Physical deformation mode reconstruction: D* = D_0 + ∑ i λ_i·(W T·u_i); Output the set of physical deformation modes, including the main deformation mode coefficients of each droplet: λ_1 (compression / stretching in the mainstream direction), λ_2 (flattening in the vertical plane), λ_3 (asymmetric shear), λ_4 (precursor of neck contraction), and λ_5 (surface ripples).

[0056] Calculate the droplet stability index χ based on the set of physical deformation modes: χ = ∑ i w_i·|λ_i|·(1 + ξ_i·|▽λ_i| / |λ_i|), where w_i is the modal physical importance weight, ξ_i is the dynamic instability coefficient, and ▽λ_i is the estimated modal time change rate. Set the critical threshold χ_thr = 1.8. When χ > χ_thr, mark the droplet as in a potentially splitting state. Output the set of droplet stability indices.

[0057] According to the set of physical deformation modes and the set of droplet stability indices, correct the droplet volume calculation: V_i = V_0_i·[1 + ∑ j f_j(λ_j) + ∑ j ∑ k g_jk(λ_j, λ_k)], where V_0_i is the volume of the equivalent sphere, f_j is the single-mode correction function: f_j(λ_j) = a_j·λ_j 2 ·sign(λ_j); g_jk is the second-order cross-term representing the non-linear coupling term between modes: g_jk(λ_j, λ_k) = b_jk·λ_j·λ_k, where a_j and b_jk are correction coefficients determined by the physical properties of the droplet. Output the set of corrected droplet volumes.

[0058] Apply calculation acceleration strategies, including gradient sparse constraint and physical pre-screening mechanism: ▽λ_sparse = ▽λ·I_{|▽λ|>τ}, where I is the indicator function and τ is the gradient threshold (0.05). The gradient sparse constraint can reduce the computational amount; the physical pre-screening mechanism can exclude physically impossible deformation states. Dynamically adjust the calculation density according to the droplet stability. For stable droplets (χ < 0.8), use simplified calculations, and for critical droplets (χ > 1.5), use refined calculations. Finally, output the corrected three-dimensional droplet size data, including the position, equivalent diameter, volume, and stability index of each droplet.

[0059] This embodiment realizes the Physical Constraint Adaptive Feature Decomposition (PCAFD) of fluids. Through the fluid physical constraint matrix, it is ensured that the extracted deformation features conform to the laws of fluid mechanics, and physically impossible forms will not occur in deformation prediction, ensuring physical consistency. Aiming at the directional heterogeneity of deformation in the flood discharge environment, appropriate weights are assigned to different directions, improving the sensitivity of deformation in the direction perpendicular to the main flow. Through the feature enhancement function, the ability to identify critical splitting states is improved, and the splitting warning accuracy rate is increased from 72% to 91%. Through multi-scale flow field decomposition, the deformation features caused by different flow mechanisms are separated, providing a more accurate physical interpretation. The second-order cross term is introduced to capture the non-linear interaction between deformation modes. Compared with traditional methods, the adaptive deformation feature decomposition method of this embodiment has been improved in terms of deformation accuracy, volume calculation, and physical consistency. Especially in the prediction of critical states, the accuracy rate has been increased by 19 percentage points, which is of great significance for flood discharge safety assessment.

[0060] As Figure 3 shown, according to one aspect of the present application, the steps of constructing a fluid constraint matrix include:

[0061] Extract wind speed and wind direction vector data from the environmental characteristic parameter set, and combine with the empirical fluid dynamics model to construct a flow field direction matrix representing the dominant deformation direction of droplets in the flow field;

[0062] Combine the turbulence intensity in the environmental characteristic parameter set to construct a scale sensitivity matrix reflecting the sensitivity of droplets of different scales to shear force;

[0063] Based on the temperature, humidity, and atmospheric pressure data in the environmental characteristic parameter set, calculate the local flow field Reynolds number, and construct a Reynolds number response matrix representing the deformation response characteristics of droplets under different Reynolds number conditions;

[0064] Use the tensor product operation to integrate the flow field direction matrix, the scale sensitivity matrix, and the Reynolds number response matrix into a unified fluid constraint matrix for physical constraints in subsequent feature decomposition processes.

[0065] As Figure 4 shown, according to one aspect of the present application, the steps of obtaining an enhanced deformation feature set include:

[0066] Based on the Weber number and deformation index of droplets in the initial droplet deformation feature set, combine with the fluid splitting theory to calculate the critical threshold of droplet deformation, and output the critical deformation threshold;

[0067] Construct a non-linear feature enhancement function and apply it to the initial droplet deformation feature set, and selectively amplify the features close to the critical deformation threshold to generate an enhanced deformation feature set;

[0068] Based on the turbulent pulsation characteristics in the environmental characteristic parameters, construct a curvature-sensitive operator adapted to the surface deformation of the droplet to detect the development trend of surface unstable waves;

[0069] Apply the curvature-sensitive operator to the enhanced deformation feature set, optimize the feature expression near the critical point, and output the optimized enhanced deformation feature set as the final enhanced deformation feature set.

[0070] According to one aspect of the present application, the steps of obtaining the flow decoupling feature set include:

[0071] Apply wavelet transform to the enhanced deformation feature set to convert the spatio-temporal domain features into frequency domain feature representations;

[0072] Based on the wind speed, turbulent intensity and spectrum analysis in the environmental characteristic parameters, identify the flow mechanism frequency bands corresponding to the mainstream traction, vortex induction and turbulent pulsation;

[0073] Use the flow mechanism frequency bands to perform multi-scale filtering decomposition on the frequency domain feature representation to obtain decomposition feature subsets corresponding to different flow mechanisms;

[0074] Perform flow mechanism correlation analysis on the decomposition feature subsets, establish the correlation mapping between the droplet deformation features and the flow field mechanism, and output the flow mechanism feature mapping;

[0075] Based on the flow mechanism feature mapping, reorganize and optimize the features, eliminate the interference between mechanisms, and output the flow decoupling feature set for subsequent physical deformation mode reconstruction.

[0076] According to one aspect of the present application, the steps of obtaining the physical deformation mode set include:

[0077] Combine the flow decoupling feature set with the fluid constraint matrix, perform eigenvalue decomposition with physical constraints, and obtain the eigenvector set;

[0078] Perform eigenvalue sorting and contribution rate screening on the eigenvector set, and retain the principal component eigenvector set with the energy proportion exceeding the preset threshold;

[0079] Based on the principles of fluid mechanics, endow physical deformation meanings to the vectors in the principal component eigenvector set, and construct a physical mode mapping matrix;

[0080] Use the physical mode mapping matrix to project and reconstruct the flow decoupling feature set, generate a physical deformation mode set representing the droplet deformation dynamics for subsequent droplet stability evaluation.

[0081] According to one aspect of the present application, the steps of calculating the droplet stability index include:

[0082] Calculate the weight distribution of each mode in the physical deformation mode set to obtain the mode weight distribution;

[0083] Based on the modal weight distribution and the physical meaning of the modes, calculate the droplet stability index that characterizes the stability degree of droplet deformation;

[0084] Compare the droplet stability index with the critical stability threshold to evaluate the droplet splitting possibility characterized by the set of physical deformation modes;

[0085] Based on the droplet splitting possibility, calculate the dynamic instability coefficient that characterizes the dynamic stability of the droplet in the flow field;

[0086] Integrate the droplet stability index and the dynamic instability coefficient to generate the droplet stability evaluation result for subsequent droplet volume correction.

[0087] According to one aspect of the present application, the steps of calculating the volume correction coefficient and outputting the corrected three-dimensional droplet size data include:

[0088] Based on the set of physical deformation modes and the droplet stability evaluation result, construct a single-mode correction function for each physical deformation mode;

[0089] Analyze the interaction between different modes in the set of physical deformation modes to generate a modal coupling matrix that characterizes the mutual influence between modes;

[0090] Based on the modal coupling matrix, calculate the non-linear coupling terms between modes to capture the non-linear effects in the complex deformation process;

[0091] Combine the single-mode correction function and the non-linear coupling terms to comprehensively calculate the volume correction coefficient of the droplet;

[0092] Apply the volume correction coefficient to the initial three-dimensional particle size data to obtain the corrected three-dimensional droplet size data as the input data for subsequent rainfall intensity prediction.

[0093] According to one aspect of the present application, before performing the constrained feature decomposition by combining the fluid constraint matrix and the flow decoupling feature set, the following real-time calculation optimization process is also included:

[0094] Calculate the feature gradient distribution for the flow decoupling feature set to generate a feature gradient matrix;

[0095] Based on the sparsity of the feature gradient matrix, construct a gradient sparsity constraint for calculation optimization;

[0096] Extract the physical laws from the fluid constraint matrix, integrate the prior knowledge of fluid dynamics to form physical pre-screening conditions;

[0097] Apply the gradient sparsity constraint and the physical pre-screening conditions to the processing process of the flow decoupling feature set to reduce the calculation complexity, improve the feature decomposition efficiency, and provide optimization support for the subsequent physical deformation mode reconstruction.

[0098] According to one aspect of the present application, the steps of obtaining the initial three-dimensional particle size data and forming the initial droplet deformation feature set include:

[0099] Performing statistical analysis on the effective circular spot data set, calculating the shape parameters and size distribution of the circular spots, and generating a particle size mapping parameter set;

[0100] Combining the particle size mapping parameter set and the camera imaging parameters in the environmental characteristic parameter set, and applying an improved Karman formula for two-dimensional to three-dimensional conversion to obtain the initial three-dimensional particle size data;

[0101] Using the initial three-dimensional particle size data and the flow field parameters in the environmental characteristic parameter set, calculating the droplet Weber number set that characterizes the key indicators of droplet deformation dynamics;

[0102] Based on the droplet Weber number set, combining with an empirical model to calculate the aspect ratio, oscillation frequency and deformation energy of the droplet, and forming the initial droplet deformation feature set, providing basic data for subsequent deformation feature enhancement.

[0103] According to one aspect of the present application, the steps of performing full-scale adaptive dynamic binning statistics on the calibrated three-dimensional droplet particle size data and calculating the droplet volume to form the rainfall intensity spatial distribution matrix include:

[0104] S41. Performing kernel density estimation on the calibrated three-dimensional droplet particle size data to obtain the droplet particle size density distribution;

[0105] S42. Based on the droplet particle size density distribution, performing adaptive binning to generate a dynamic binning interval reflecting the droplet distribution characteristics;

[0106] S43. Classifying and statistically analyzing the calibrated three-dimensional droplet particle size data according to the dynamic binning interval, calculating the total volume of droplets in each interval, and constructing a droplet volume distribution matrix;

[0107] S44. Combining the droplet volume distribution matrix and the wind speed, wind direction, temperature, and humidity factors in the environmental characteristic parameter set, and applying a rainfall intensity conversion model to calculate and output the rainfall intensity spatial distribution matrix.

[0108] Specifically, reading the calibrated three-dimensional droplet particle size data, and applying the kernel density estimation (KDE) method to analyze the particle size distribution characteristics. The KDE formula is: f(d) = (1 / nh)·∑ i=1 n K((d - d_i) / h), where d is the particle size value, d_i is the diameter of the i-th droplet, K is the kernel function (selecting the Gaussian kernel), and h is the bandwidth parameter. The bandwidth is adaptively selected: h = 0.9·min(σ, IQR / 1.34)·n (-1 / 5), where σ is the sample standard deviation, IQR is the interquartile range, and n is the sample size. Identify the density peaks and valleys through the KDE curve, and output the particle size distribution feature point set.

[0109] Dynamically divide the particle size statistical intervals based on the particle size distribution feature point set. First, identify the density peak as the interval center point: P = {d | f'(d)=0 and f''(d)<0}; where f''(d) is the change rate of the density function; divide the interval boundaries at the density valleys: B = {d | f'(d)=0 and f''(d)>0}; merge the low-density regions to generate the final binned intervals: I = {[B_i, B_{i+1}] | ∫_{B_i} B_{i+1} f(d)·dd > τ_min}, where B_i is the interval boundary point of the particle size distribution; τ_min is the minimum interval weight threshold (0.05). Output the dynamically binned interval set I.

[0110] Create a three-dimensional space grid based on the spatial sampling position data. The grid resolution is dynamically determined according to the size H of the flood discharge structure: Δx = Δy = Δz = α·H, where Δx, Δy, and Δz respectively represent the grid spacings (resolutions) of the three-dimensional space grid in the x, y, and z directions; α is the grid scale factor (0.05~0.1). For each grid cell G_ijk, extract the subset of droplets falling into it: D_ijk = {d_m | (x_m, y_m, z_m) ∈ G_ijk}; output the spatial grid droplet distribution data.

[0111] Combine the dynamically binned interval set I and the spatial grid droplet distribution data to construct a three-dimensional matrix M: M(i, j, k, l) =∑ d_m∈D_ijk 且 d_m∈I_l 1, where i, j, k are spatial indices and l is the particle size interval index. At the same time, calculate the total droplet volume of each particle size interval in each grid cell: V(i, j, k, l) = ∑ d_m∈D_ijk 且 d_m∈I_l V_m, where V_m is the droplet volume, from the calibrated droplet volume set. Output the three-dimensional droplet volume distribution matrix.

[0112] Calculate the rainfall intensity value of each spatial grid according to the three-dimensional droplet volume distribution matrix: I(i, j, k) = V_total(i, j, k) / (A·T)·3600000, where V_total(i, j, k)=∑ l V(i, j, k, l) is the total droplet volume of the grid (m 3 ), A is the effective collection area of the test paper (m 2 ), T is the exposure time (s), and 3600000 is the unit conversion factor (m / s to mm / h). Perform spatial interpolation on the rainfall intensity values to generate a continuous distribution field: I(x, y, z) = ∑i，j，k I(i, j, k)·φ_ijk(x, y, z), where φ_ijk is a trilinear interpolation basis function. Output the rainfall intensity spatial distribution matrix as the final prediction result.

[0113] According to one aspect of the present application, the steps of calculating and outputting the rainfall intensity spatial distribution matrix by applying the rainfall intensity conversion model include:

[0114] Based on the droplet volume distribution matrix, calculate the initial rainfall intensity data by applying the volume-rainfall intensity transformation model derived from the principles of fluid mechanics.

[0115] Combine the wind speed gradient, turbulence intensity, and temperature and humidity distribution in the environmental characteristic parameters set to correct the initial rainfall intensity data for environmental factors to obtain corrected rainfall intensity data.

[0116] According to the spatial distribution of sampling points, map the corrected rainfall intensity data to a three-dimensional space coordinate system to construct a rainfall intensity spatial grid.

[0117] Perform interpolation optimization on the rainfall intensity spatial grid (such as using the Kriging interpolation algorithm) to generate a continuous rainfall intensity spatial distribution matrix as the final predicted result of the rainfall intensity in the flood discharge atomization area.

[0118] According to another aspect of the present application, a method for predicting the rainfall intensity of flood discharge atomization based on the fog and rain droplet distribution includes: Drop size spectrum image and environmental data collection: According to the type of flood discharge building, set up measuring points within a certain range around the fog source, collect the rain droplet size spectrum images of the flood discharge atomization of the hydropower station, and collect environmental data (wind speed V, temperature T, humidity H) through preset environmental sensors and input them into the acquisition and storage module of the system; perform preprocessing operations such as image digital indexing, three-dimensional space modeling of the drop size spectrum collection position (import the three-dimensional coordinates of the drop size spectrum collection position), gray processing and contrast enhancement of the drop size spectrum patch images in sequence.

[0119] Drop size spectrum environmental noise, candidate patch identification and feature marking: Establish an SVM recognition model for the rain droplet size spectrum environmental noise and atomization feature patches in the complex environment of rain, fog, water, and gas, input the manually labeled noise and patch image data to train the model; import the preprocessed drop size spectrum images, and the model automatically identifies the environmental noise and candidate patch types and generates labels (the environmental noise is marked as salt-and-pepper noise, Gaussian noise; the patches are marked as clear circular patches, blurred circular patches, overlapping circular patches, irregular patches, and the latter three are collectively referred to as patches to be verified).

[0120] Drop spectrum environmental noise and multi-strategy processing of patches to be verified: After the drop spectrum image is recognized by the model and generates labels, it is input into the image processing module. According to the noise labels (salt and pepper noise / Gaussian noise), the corresponding type of noise is filtered by the matching filtering algorithm to obtain the candidate patch image; the known clear circular patches in the candidate patches are hidden, and for the remaining patches to be verified, the fuzzy C-means clustering algorithm is used to de-blur the circular patches, the watershed algorithm is used to segment the overlapping circular patches, and the Hough circle detection combined with morphological operations is used to perform pseudo-circle processing and shape correction on the irregular patches to obtain the processed circular patches, and then the clear circular patches are displayed. At this time, all effective circular patches are obtained.

[0121] Conversion of drop spectrum circular patches into three-dimensional droplets and correction of environmental variable effects: Based on the Karman vortex street empirical formula for the drop spectrum patches, and at the same time loading the environmental sensor data (wind speed, temperature, humidity), a dynamic mapping relationship is established to map the two-dimensional projection diameter of fog and rain droplets considering the influence of environmental field variables to the particle size in three-dimensional space, and the two-dimensional circular patch diameter data is equivalently converted into three-dimensional space droplet particle size data, and corrected according to the influence weight of environmental variables.

[0122] Three-dimensional dynamic binning and full extraction of droplet particle sizes: For the converted droplet particle sizes, automated statistical analysis is performed. The peak of the particle size diameter distribution density is identified through kernel density estimation (KDE), and the statistical intervals of droplet particle sizes are dynamically divided by combining DBSCAN spatial clustering. The number sequence of particle sizes is statistically analyzed in each interval to obtain the number of droplets corresponding to each existing particle size value in the continuous particle size value interval. The diameter-number data pairs of the atomized raindrops actually produced by the fog source are fully extracted and uploaded to the database.

[0123] Dynamic correction of total rainfall and calculation of spatial rainfall intensity: The total rainfall is calculated according to the droplet particle size, and typical samples are collected to correct the total rainfall result. The above steps are cycled to dynamically correct the Karman empirical parameters and iterate the conversion results of the three-dimensional particle sizes; the total rainfall intensity is calculated again according to the corrected droplet particle size number sequence; the atomized rainfall intensity value is calculated by combining the actual sampling area and the collection time; the above steps are repeated to sequentially count the rainfall intensities corresponding to all drop spectrum patterns, and the calculated rainfall intensity values are mapped to three-dimensional space to obtain the distribution law of the atomized rainfall intensity field.

[0124] In a specific embodiment, taking the flood discharge condition of a large water conservancy project as the background, the rainfall intensity in the flood discharge atomization area is predicted. The experiment is carried out under the conditions of a flood discharge flow rate of 1200 m 3 / s, an environmental temperature of 25 °C, and a relative humidity of 85%. The specific steps are as follows:

[0125] Step 1: Collect the droplet images and environmental characteristic parameters in the flood discharge area.

[0126] A high-speed camera (model FASTCAM SA-X2, frame rate 5000fps, resolution 1024×1024 pixels) was used to collect droplet images 15m downstream of the flood discharge outlet, and a total of 1000 consecutive images were obtained. At the same time, an environmental monitoring device was used to collect the following environmental characteristic parameters: average wind speed: 3.5m / s, wind direction: 15° north by northwest; environmental temperature: 25°C; relative humidity: 85%; atmospheric pressure: 101.3kPa; turbulence intensity: 18% (measured by a hot-wire anemometer). The above parameters constitute the environmental characteristic parameter set E = {V, θ, T, H, P, I}, where V is the wind speed of 3.5m / s; θ is the wind direction angle of 285°; T is the environmental temperature of 25°C; H is the relative humidity of 85%; P is the atmospheric pressure of 101.3kPa; I is the turbulence intensity of 18%.

[0127] Step 2: Extract effective circular spot data.

[0128] 2.1. Preprocess the 1000 collected images. First, apply Gaussian filtering (kernel size 5×5, σ = 1.5) to remove noise, and then perform binarization through an adaptive threshold segmentation algorithm (window size 15×15, threshold offset constant C = 8) to obtain the binarized image set B = {B1, B2, ..., B 1000}.

[0129] 2.2. Apply a circle detection algorithm to the binarized image set B, set the minimum circle radius to 3 pixels, the maximum circle radius to 50 pixels, the accumulator threshold to 40, and the minimum center distance to 10 pixels. A total of 8563 droplet circular spots were detected, including effective circular spots and interference circular spots.

[0130] 2.3. Construct a multi-feature discrimination function F(c) = w1·A(c) + w2·G(c) + w3·S(c); where A(c) is the area of the circular spot, and the value range is [9, 7854] pixels 2 ; G(c) is the gray level uniformity, calculated as the ratio of the standard deviation to the mean of the gray levels within the circular spot area, and the value range is [0, 1]; S(c) is the shape regularity, calculated as 4π·A(c) / P(c) 2 , P(c) is the perimeter, and the value range is [0, 1]; w1, w2, w3 are weight coefficients, which are 0.4, 0.3, 0.3 respectively; c is the circular spot to be discriminated.

[0131] 2.4. Set the discrimination function threshold τ = 0.65, screen the 8563 detected circular spots, and retain the circular spots with F(c) > τ as effective circular spots. Finally, 6428 effective circular spots are obtained, constituting the effective circular spot data set C = {c1, c2, ..., c 6428}. Each valid circular spot \(c_i\) contains the following attributes: center coordinates \((x_i, y_i)\), radius \(r_i\), area \(A_i\), average gray value \(\mu_i\), and frame number \(f_i\).

[0132] Step 3: Extract the droplet deformation features.

[0133] 3.1. Extract the initial droplet deformation features.

[0134] 3.1.1. Conduct statistical analysis on the valid circular spot dataset \(C\) and calculate the following particle size mapping parameters: average circular spot radius: \(\overline{r}=12.8\) pixels; standard deviation of circular spot radius: \(\sigma_r = 5.2\) pixels; circular spot density: \(\rho_c = 6428 / (1000\times1024\times1024)=6.13\times10\) -6 spots / pixel 2 ; average gray gradient: \(\nabla I = 45.3\) gray values / pixel. The above parameters form the particle size mapping parameter set \(M=\{\overline{r},\sigma_r,\rho_c,\nabla I\}\).

[0135] 3.1.2. Combine the particle size mapping parameter set \(M\) and the environmental characteristic parameter set \(E\), and apply the improved Karman formula for two-dimensional to three-dimensional conversion: \(D_i = 2\cdot r_i\cdot K_s\cdot(1 + \alpha\cdot I) / (1-\beta\cdot V\cdot\cos(\theta))\); where \(D_i\) is the three-dimensional droplet diameter, in mm; \(r_i\) is the circular spot radius, in pixels; \(K_s\) is the proportionality coefficient, with a value of \(0.085\) mm / pixel; \(\alpha\) is the turbulence adjustment coefficient, with a value of \(0.25\); \(\beta\) is the wind speed adjustment coefficient, with a value of \(0.06\) s / m; \(I\) is the turbulence intensity of \(18\%\); \(V\) is the wind speed of \(3.5\) m / s; \(\theta\) is the wind direction angle of \(285^{\circ}\).

[0136] Taking the first droplet (\(i = 1\)) as an example, its circular spot radius \(r_1 = 9.5\) pixels. Calculate its three-dimensional diameter \(D_1\): \(D_1 = 2\cdot r_1\cdot K_s\cdot(1 + \alpha\cdot I) / (1-\beta\cdot V\cdot\cos(\theta))=2\cdot9.5\cdot0.085\cdot(1 + 0.25\cdot0.18) / (1-0.06\cdot3.5\cdot\cos(285^{\circ}))=2\cdot9.5\cdot0.085\cdot(1 + 0.045) / (1-0.06\cdot3.5\cdot(-0.2588))=2\cdot9.5\cdot0.085\cdot1.045 / (1 + 0.0543)=2\cdot9.5\cdot0.085\cdot1.045 / 1.0543=1.615\cdot0.085\cdot1.045 / 1.0543=1.615\cdot0.08892 / 1.0543=1.615\cdot0.0843=1.36\) mm. Apply the above formula to calculate the three-dimensional diameters of all droplets, forming the initial three-dimensional particle size dataset \(D=\{D_1,D_2,\cdots,D\}\), where \(D\) 6428} and \(D\) iThe value range of [] is [0.6mm, 9.8mm].

[0137] 3.2. Calculate the Weber number and deformation characteristics of the droplets.

[0138] 3.2.1. Based on the initial three-dimensional particle size dataset D and the environmental characteristic parameter set E, calculate the Weber number of the droplet: We_i = ρ_a·V_rel 2 ·D_i / σ; where We_i is the Weber number of the droplet, dimensionless; ρ_a is the air density, with a value of 1.2 kg / m 3 ; V_rel is the relative wind speed, calculated as V·(1 + 0.2·I), in units of m / s; σ is the surface tension coefficient of water, with a value of 0.072 N / m; D_i is the droplet diameter, in units of m.

[0139] Taking the 1st droplet (i = 1) as an example, its diameter D1 = 1.36 mm, calculate its Weber number We1: V_rel = V·(1 + 0.2·I) = 3.5·(1 + 0.2·0.18) = 3.5·1.036 = 3.626 m / s; We1 = ρ_a·V_rel 2 ·D1 / σ = 1.2·3.626 2 ·0.00136 / 0.072 = 1.2·13.148·0.00136 / 0.072 = 1.2·0.01788 / 0.072 = 0.02146 / 0.072 = 0.298. This value belongs to the surface tension-dominated range, and the droplet is close to spherical.

[0140] Taking the 4000th droplet (i = 4000) as an example, its diameter D 4000 = 5.2 mm, calculate its Weber number We 4000 : We 4000 = ρ_a·V_rel 2 ·D 4000 / σ = 1.2·3.626 2 ·0.0052 / 0.072 = 1.2·13.148·0.0052 / 0.072 = 1.2·0.06837 / 0.072 = 0.08204 / 0.072 = 11.4. This value is close to the critical Weber number, and the droplet will undergo significant deformation.

[0141] The calculated Weber number set of the droplets We = {We1, We2,..., We 6428}, and the Weber number range is [0.5, 18.3].

[0142] 3.2.2. Based on the droplet Weber number set We, calculate the droplet deformation characteristics: Droplet aspect ratio: χ_i = 1 + 0.15·We_i0.77 ; The range of χ_i is [1.09, 1.87]; Droplet oscillation frequency: f_i = (1 / π)·(8·σ / (ρ_w·D_i 3 )) 0.5 ; The range of f_i is [18.5 Hz, 298.3 Hz], and ρ_w is the density of water, 1000 kg / m 3 ; Droplet deformation energy: E_i = π·σ·D_i 2 ·(χ_i - 1) 2 / 8; The range of E_i is [2.1×10 -8 J, 5.7×10 -6 J]; The above parameters constitute the droplet initial deformation feature set F_0 = {F 01 , F 02 , ..., F 06428}, where F 0i = {D_i, We_i, χ_i, f_i, E_i}.

[0143] 3.3. Construct the fluid constraint matrix.

[0144] 3.3.1. Construct the flow field direction matrix M_dir: M_dir(i, j) = cos(θ_ij)·exp(-d_ij 2 / (2·l_t 2 )); where M_dir(i, j) is the element in the i-th row and j-th column of the flow field direction matrix, with a value range of [-1, 1]; θ_ij is the angle between the mainstream direction and the direction of the line connecting i and j, in radians; d_ij is the distance between spatial positions i and j, in m; l_t is the turbulent characteristic length scale, with a value of 0.25 m. The calculated flow field direction matrix M_dir has a dimension of 100×100 (100 representative points are obtained through spatial partition sampling) and is used to characterize the dominant deformation direction of the droplet in the flow field.

[0145] 3.3.2. Construct the scale-sensitive matrix M_scale: M_scale(d, t) = exp(-(d - d_t) 2 / (2·σ_d 2))·(1 - exp(-t·γ)); where M_scale(d, t) is the element in the d-th row and t-th column of the scale-sensitive matrix, with a value range of [0, 1]; d is the droplet diameter, taking 10 grading points, in the range of [0.6 mm, 9.8 mm]; d_t are the characteristic diameter grading points, with 10 values; t is the turbulence intensity grading, divided into 5 grades, in the range of [0, 0.3]; σ_d is the diameter sensitivity parameter, with a value of 1.2 mm; γ is the turbulence sensitivity parameter, with a value of 8. The calculated scale-sensitive matrix M_scale has a dimension of 10×5 and is used to characterize the sensitivity of droplets of different scales to shear force.

[0146] 3.3.3. Construct the Reynolds number response matrix M_Re: M_Re(Re, ξ) = (Re / Re_cr) ξ ·exp(1 - Re / Re_cr); where M_Re(Re, ξ) is the element in the Re-th row and ξ-th column of the Reynolds number response matrix, with a value range of [0, 1]; Re are the Reynolds number grading points, with 15 values, in the range of [100, 10000]; ξ is the shape parameter, divided into 8 grades, in the range of [0.5, 4]; Re_cr is the critical Reynolds number, with a value of 3000. The calculated Reynolds number response matrix M_Re has a dimension of 15×8 and is used to characterize the deformation response characteristics of droplets under different Reynolds number conditions.

[0147] 3.3.4. Integrate the constraint matrix through tensor product operation: M_constraint = M_dir Θ M_scale Θ M_Re; where: M_constraint is the final fluid constraint matrix; Θ represents the tensor product operation. The fluid constraint matrix M_constraint is obtained through the tensor product operation, with a dimension of 100×10×15×5×8, and then reduced to a 200×200 matrix through dimensionality reduction for subsequent physical constraint eigenvalue decomposition.

[0148] 3.4. Enhance the critical characteristics of droplet deformation.

[0149] 3.4.1. Calculate the critical threshold of droplet deformation: We_cr = 12 + 3.45·(ρ_w / ρ_a) 0.2 ·(1 + 1.2·I); where We_cr is the critical Weber number of droplet deformation, and in this embodiment, We_cr = 15.8 is calculated; ρ_w is the density of water, 1000 kg / m 3 ; ρ_a is the density of air, 1.2 kg / m 3 ; I is the turbulence intensity, 18%.

[0150] 3.4.2. Construct the non - linear feature enhancement function: G(We_i) = 1 + k·tanh(α·(We_i - We_base) / (We_cr - We_base)); where G(We_i) is the feature enhancement coefficient, and its value range is [1, 1 + k]; We_i is the droplet Weber number; k is the enhancement coefficient with a value of 3.5; α is the sharpness parameter with a value of 2.5; We_base is the reference Weber number with a value of 8; We_cr is the critical Weber number 15.8.

[0151] Taking the 4000th droplet (i = 4000) as an example, its Weber number We 4000 = 11.4, calculate its feature enhancement coefficient: G(We 4000 ) = 1 + k·tanh(α·(We 4000 - We_base) / (We_cr - We_base)) = 1 + 3.5·tanh(2.5·(11.4 - 8) / (15.8 - 8)) = 1 + 3.5·tanh(2.5·3.4 / 7.8) = 1 + 3.5·tanh(2.5·0.4359) = 1 + 3.5·tanh(1.0897) = 1 + 3.5·0.7981 = 1 + 2.7934 = 3.7934; the deformation characteristics of this droplet will be enhanced by about 3.79 times, indicating that it is close to the critical deformation state.

[0152] Apply the enhancement function to each feature in the initial droplet deformation feature set F_0: F' 0i = F 0i ·G(We_i), and obtain the enhanced deformation feature set F_1 = {F 11 , F 12 ,..., F 16428}.

[0153] 3.4.3. Construct the curvature - sensitive operator: K(x, y) = exp(-(x 2 + y 2 ) / (2·σ_k 2 ))·(x 2 - y 2 ) / (x 2 + y 2 + ε); where K(x, y) is the value of the curvature - sensitive operator at the spatial position (x, y); σ_k is the spatial scale parameter with a value of 2.5 pixels; ε is a small quantity to prevent division by zero, with a value of 10 -6 .

[0154] 3.4.4. Apply the curvature - sensitive operator to the enhanced deformation feature set F_1: F 2i = F1i + β·conv(F 1i , K); where F 2i is the optimized feature; conv represents the convolution operation; β is the convolution weight with a value of 0.4. Through the above operations, the optimized enhanced deformation feature set F_2 = {F 21 , F 22 , ..., F 26428} is obtained.

[0155] 3.5. Decouple the droplet flow characteristics.

[0156] 3.5.1 Apply wavelet transform to the optimized enhanced deformation feature set F_2: W_i(a, b) = (1 / sqrt(a))·∫f_i(t)·ψ((t - b) / a)dt; where W_i(a, b) is the wavelet coefficient of the droplet feature f_i; a is the scale parameter; b is the translation parameter; ψ is the wavelet basis function, and in this embodiment, the Morlet wavelet is used. The frequency-domain feature representation W = {W1, W2,..., W 6428} is calculated.

[0157] 3.5.2 Identify the flow mechanism frequency bands based on the environmental feature parameter set E: Main flow traction frequency band: f_main = [0.8·V / L, 1.2·V / L], calculated as [1.4Hz, 2.1Hz], where L is the characteristic length of 2m; Vortex-induced frequency band: f_vortex = [St·V / D_avg, 1.5·St·V / D_avg], calculated as [10.5Hz, 15.8Hz], where St is the Strouhal number of 0.2 and D_avg is the average droplet diameter of 0.0067m; Turbulent pulsation frequency band: f_turb = [V / λ_min, V / λ_max], calculated as [35Hz, 175Hz], where λ_min and λ_max are the turbulent scale ranges [0.02m, 0.1m].

[0158] 3.5.3 Use the flow mechanism frequency bands to perform multi-scale filtering decomposition on the frequency-domain feature representation W: W_main = H_main(f)·W; where H_main(f) is the main flow band-pass filter, with a value of 1 in the f_main frequency range and 0 otherwise; W_vortex = H_vortex(f)·W; where H_vortex(f) is the vortex band-pass filter, with a value of 1 in the f_vortex frequency range and 0 otherwise; W_turb = H_turb(f)·W; where H_turb(f) is the turbulent band-pass filter, with a value of 1 in the f_turb frequency range and 0 otherwise; The decomposed feature subsets {W_main, W_vortex, W_turb} are obtained.

[0159] 3.5.4. Conduct correlation analysis of the flow mechanism for the decomposed feature subsets: C_mech(i, j) = corr(W_i, F_j); where C_mech(i, j) is the correlation coefficient between mechanism i and deformation feature j; corr represents the calculation of the correlation coefficient; W_i is the i-th subset in the decomposed feature subsets; F_j is the j-th feature in the droplet deformation feature set. Obtain the flow mechanism feature mapping C_mech, with a dimension of 3×5, representing the degree of association between 3 flow mechanisms and 5 deformation features.

[0160] 3.5.5. Based on the flow mechanism feature mapping, perform feature recombination and optimization: F 3i = w_main·F 2i _main + w_vortex·F 2i _vortex + w_turb·F 2i _turb; where F 3i is the decoupled feature after recombination; F 2i _main, F 2i _vortex, F 2i _turb are the feature components under the corresponding mechanisms; w_main, w_vortex, w_turb are the weight coefficients, which are 0.3, 0.4, and 0.3 respectively. Through the above operations, obtain the flow decoupled feature set F_3 = {F 31 , F 32 , ..., F 36428}.

[0161] 3.6. Reconstruct the physical deformation mode of the droplet.

[0162] 3.6.1. Combine the flow decoupled feature set F_3 with the fluid constraint matrix M_constraint, and perform eigen decomposition with physical constraints: (λ_i, v_i) = eigdecomp(F_3 T ·M_constraint·F_3); where λ_i is the i-th eigenvalue; v_i is the i-th eigenvector; eigdecomp represents the eigen decomposition operation. Calculate to obtain the eigenvector set V = {v1, v2,..., v_n} and the corresponding eigenvalue set Λ = {λ1, λ2,..., λ_n}.

[0163] 3.6.2. Sort the eigenvectors according to the eigenvalue magnitudes, and calculate the cumulative energy contribution rate: E_cum(k) = (∑ i=1 k λ_i) / (∑ i=1 nλ_i); where E_cum(k) is the cumulative energy contribution rate of the first k eigenvalues. Set the threshold η = 0.95, select the smallest k value that satisfies E_cum(k) ≥ η. In this embodiment, k = 8, and the principal component eigenvector set V_main = {v1, v2,..., v8} is obtained.

[0164] 3.6.3. Assign physical deformation meanings to each vector in the principal component eigenvector set: v1: Represents the overall flattening deformation of the droplet, associated with the mainstream traction effect; v2: Represents the asymmetric stretching deformation of the droplet, associated with the shear effect; v3: Represents the surface corrugation deformation of the droplet, associated with the turbulent pulsation effect; v4: Represents the oscillating deformation of the droplet, associated with the vortex-induced effect; v5: Represents the axial rotational deformation of the droplet, associated with the rotational shear effect; v6, v7, v8: Represent high-order composite deformation modes. Construct a physical modal mapping matrix P, with a dimension of n×8, where P(i, j) represents the mapping relationship between the i-th characteristic variable and the j-th physical mode.

[0165] 3.6.4. Project and reconstruct the flow decoupling feature set F_3 using the physical modal mapping matrix P: M_i = P·(v_i T ·F 3i ); where M_i is the physical deformation mode of the i-th droplet, with a dimension of 8×1, representing the weights of 8 physical modes.

[0166] Taking the 4000th droplet as an example, assuming its flow decoupling eigenvector F 34000 = [0.25, 0.42, 0.18, 0.35, 0.15], and the first two vectors of the principal component eigenvector set are v1 = [0.38, 0.41, 0.12, 0.52, 0.26] and v2 = [0.21, 0.62, 0.31, 0.15, 0.18], calculate its physical mode: v1 T ·F 34000 = 0.38·0.25 + 0.41·0.42 + 0.12·0.18 + 0.52·0.35 + 0.26·0.15 = 0.095 + 0.1722 + 0.0216 + 0.182 + 0.039 = 0.5098; v2 T ·F 34000= 0.21·0.25 + 0.62·0.42 + 0.31·0.18 + 0.15·0.35 + 0.18·0.15 = 0.0525 + 0.2604 + 0.0558 + 0.0525 + 0.027 = 0.4482; Assume the first two columns of the simplified physical modal mapping matrix P are: P(:, 1) = [0.85, 0.12, 0.03, 0.0, 0.0] T ; P(:, 2) = [0.15, 0.76, 0.09, 0.0, 0.0] T ; Calculate the first two physical modes of the droplet: M 4000 (1) = P(:, 1)·(v1 T ·F 34000 ) = 0.85·0.5098 = 0.4333; M 4000 (2) = P(:, 2)·(v2 T ·F 34000 ) = 0.76·0.4482 = 0.3406; This indicates that the weight of the droplet in the flattening deformation mode is 0.4333, and the weight in the asymmetric stretching deformation mode is 0.3406. Through the above operations, the physical deformation mode set M = {M1, M2, ..., M 6428} is obtained.

[0167] 3.7. Evaluate the droplet stability and volume correction.

[0168] 3.7.1. Calculate the weight distribution of each mode in the physical deformation mode set M: w_ij = M_i(j) / (∑ k=1 8 M_i(k)); where w_ij is the normalized weight of the j-th mode of the i-th droplet. The mode weight distribution W = {W1, W2, ..., W 6428} is obtained, where W i = {w_i1, w_i2, ..., w_i8}.

[0169] 3.7.2. Calculate the droplet stability index: SI_i = 1 - (w_i2 + 1.5·w_i3 + 2·w_i4); where SI_i is the stability index of the i-th droplet, and its value range is [0, 1]. The lower the value, the more unstable it is; w_i2, w_i3, and w_i4 are the weights of the asymmetric stretching, surface corrugation, and oscillation deformation modes respectively.

[0170] Taking the 4000th droplet as an example, assume its physical deformation mode weight is W 4000 = {w 4000，1= 0.32, w 4000，2 = 0.28, w 4000，3 = 0.18, w 4000，4 = 0.12, w 4000，5 = 0.06, w 4000，6 = 0.02, w 4000，7 = 0.01, w 4000，8 = 0.01}, calculate its stability index: SI 4000 = 1 - (w 4000，2 + 1.5·w 4000，3 + 2·w 4000，4 ) = 1 - (0.28 + 1.5·0.18 + 2·0.12) = 1 - (0.28 + 0.27 + 0.24) = 1 - 0.79 = 0.21. The stability index of this droplet is 0.21, which is lower than the critical stability index of 0.4, indicating that this droplet is in an unstable state and has a high possibility of splitting.

[0171] Calculate the droplet stability index set SI = {SI1, SI2,..., SI 6428}, and its value range is [0.25, 0.95].

[0172] 3.7.3. Evaluate the droplet splitting possibility: P_breakup_i = 1 / (1 + exp((SI_i - SI_cr) / Δ)); where P_breakup_i is the splitting possibility of the i-th droplet, and its value range is [0, 1]; SI_i is the stability index; SI_cr is the critical stability index, with a value of 0.4; Δ is the transition zone width parameter, with a value of 0.1. Calculate the droplet splitting possibility set P_breakup = {P1, P2,..., P 6428}.

[0173] 3.7.4. Calculate the dynamic instability coefficient: DI_i = P_breakup_i·(We_i / We_cr)·(1 + γ·I); where DI_i is the dynamic instability coefficient of the i-th droplet; P_breakup_i is the splitting possibility; We_i is the Weber number; We_cr is the critical Weber number 15.8; γ is the turbulence influence coefficient, with a value of 1.2; I is the turbulence intensity of 18%. Calculate the dynamic instability coefficient set DI = {DI1, DI2,..., DI 6428}.

[0174] 3.7.5. Generate the droplet stability evaluation result: Stability_i = {SI_i, P_breakup_i, DI_i}; where Stability_i is the stability evaluation result of the i-th droplet, including the stability index SI_i, the breakup probability P_breakup_i, and the dynamic instability coefficient DI_i. Obtain the droplet stability evaluation result set Stability = {Stability1, Stability2,..., Stability 6428}.

[0175] 3.8. Droplet volume correction.

[0176] 3.8.1. Construct the unimodal correction function: f_j(w_ij) = 1 + a_j·w_ij·(1 - exp(-b_j·w_ij)); where f_j(w_ij) is the correction function of the j-th mode; w_ij is the weight of the j-th mode of the i-th droplet; a_j is the correction coefficient, with different values for different modes, in the range [0.1, 0.8]; b_j is the non-linear parameter, with different values for different modes, in the range [1, 5]. Calculate the unimodal correction function set f = {f1, f2,..., f8}.

[0177] 3.8.2. Analyze the interaction between modes and generate the mode coupling matrix: C(i, j) = ∑ k=1 6428 (w_ki·w_kj) / (sqrt(∑ k=1 6428 w_ki 2 )·sqrt(∑ k=1 6428 w_kj 2 )); where C(i, j) is the coupling coefficient between mode i and mode j, with a value range of [-1, 1]. Calculate the mode coupling matrix C, with a dimension of 8×8.

[0178] 3.8.3. Calculate the non-linear coupling term between modes: NL_i(j, k) = c_jk·w_ij·w_ik·(1 + sin(π·SI_i)); where NL_i(j, k) is the non-linear coupling term between mode j and mode k of the i-th droplet; c_jk is the coupling coefficient, taken from the mode coupling matrix C(j, k); w_ij, w_ik are the mode weights; SI_i is the stability index. Calculate the non-linear coupling term set NL = {NL1, NL2,..., NL 6428}, and each NL_i is an 8×8 matrix.

[0179] 3.8.4. Calculate the droplet volume correction factor: VC_i = 1 + ∑ j=1 8 (f_j(w_ij) - 1) + ∑ j=1 8 ∑ k=j+1 8 NL_i(j, k); where VC_i is the volume correction factor of the i-th droplet. The calculated volume correction factor set VC = {VC1, VC2,..., VC 6428}, and the value range is [0.9, 1.35].

[0180] 3.8.5. Apply the volume correction factor to obtain the corrected three-dimensional droplet size data: D_corr_i = D_i · 3 sqrt(VC_i); where D_corr_i is the corrected droplet diameter; D_i is the initial droplet diameter; 3 sqrt() represents the cube root operation. The calculated corrected three-dimensional droplet size data set D_corr = {D_corr1, D_corr2,..., D_corr 6428}, and the diameter range is [0.58 mm, 10.6 mm].

[0181] Taking the 4000th droplet as an example, according to its modal weight W 4000 Calculate the single-modal correction function value (only calculate the first 3 modes): f1(w 4000，1 ) = 1 + a1 · w 4000，1 · (1 - exp(-b1 · w 4000，1 )) = 1 + 0.2 · 0.32 · (1 - exp(-2 · 0.32)) = 1 + 0.2 · 0.32 · (1 - exp(-0.64)) = 1 + 0.2 · 0.32 · (1 - 0.5273) = 1 + 0.2 · 0.32 · 0.4727 = 1 + 0.2 · 0.1513 = 1 + 0.0303 = 1.0303; f2(w 4000，2 ) = 1 + a2 · w 4000，2 · (1 - exp(-b2 · w 4000，2 )) = 1 + 0.6 · 0.28 · (1 - exp(-3 · 0.28)) = 1 + 0.6 · 0.28 · (1 - exp(-0.84)) = 1 + 0.6 · 0.28 · (1 - 0.4317) = 1 + 0.6 · 0.28 · 0.5683 = 1 + 0.6 · 0.1591 = 1 + 0.0955 = 1.0955; f3(w 4000，3 ) = 1 + a3 · w 4000，3·(1 - exp(-b3·w 4000，3 )) = 1 + 0.8·0.18·(1 - exp(-4·0.18)) = 1 + 0.8·0.18·(1 - exp(-0.72)) = 1 + 0.8·0.18·(1 - 0.4868) = 1 + 0.8·0.18·0.5132 = 1 + 0.8·0.0924 = 1 + 0.0739 = 1.0739; Assume the calculation result of the non - linear coupling term is (only the coupling of mode 1 and mode 2 is given): NL 4000 (1, 2) = c 12 ·w 4000，1 ·w 4000，2 ·(1 + sin(π·SI 4000 )) = 0.5·0.32·0.28·(1 + sin(π·0.21)) = 0.5·0.32·0.28·(1 + sin(0.66π)) = 0.5·0.32·0.28·(1 + 0.9511) = 0.5·0.32·0.28·1.9511 = 0.5·0.0896·1.9511 = 0.0875; Calculate the droplet volume correction coefficient (simplified calculation, only considering the first 3 modes and their coupling): VC 4000 = 1 + (f1(w 4000，1 ) - 1) + (f2(w 4000，2 ) - 1) + (f3(w 4000，3 ) - 1) + NL 4000 (1, 2) + NL 4000 (1, 3) + NL 4000 (2, 3) = 1 + (1.0303 - 1) + (1.0955 - 1) + (1.0739 - 1) + 0.0875 + 0.0532 + 0.0425 = 1 + 0.0303 + 0.0955 + 0.0739 + 0.0875 + 0.0532 + 0.0425 = 1 + 0.3829 = 1.3829; Apply the volume correction coefficient to obtain the corrected droplet diameter: D_corr 4000 = D 4000 · 3 sqrt(VC 4000 ) = 5.2· 3 sqrt(1.3829) = 5.2·1.1139 = 5.79mm, which indicates that after considering the deformation, the effective diameter of this droplet should be 5.79mm, an increase of 11.4% compared to the original measured value.

[0182] Step 4: Calculate the spatial distribution of rainfall intensity.

[0183] 4.1 Apply the kernel density estimation algorithm to the calibrated three-dimensional droplet size dataset D_corr: f(D) = (1 / n)·∑ i=1 n K_h(D - D_corr_i); where f(D) is the probability density at diameter D; n is the number of droplets, 6428; K_h is the kernel function, and in this embodiment, the Gaussian kernel function K_h(x) = (1 / (h·sqrt(2π)))·exp(-x 2 / (2·h 2 )); h is the bandwidth parameter with a value of 0.3 mm. Calculate the droplet size density distribution f(D) with a diameter range of [0.5 mm, 11 mm].

[0184] 4.2 Based on the droplet size density distribution f(D), use the adaptive binning algorithm to generate a dynamic binning interval: D_edges = adaptive_bin(f(D), N_bin, η_min); where D_edges are the binning boundary points; adaptive_bin is the adaptive binning algorithm; N_bin is the number of bins with a value of 10; η_min is the minimum probability density threshold with a value of 0.01. Calculate the dynamic binning interval D_edges = {0.5 mm, 1.2 mm, 2.0 mm, 3.1 mm, 4.3 mm, 5.5 mm, 6.8 mm, 8.0 mm, 9.1 mm, 10.0 mm, 11.0 mm}.

[0185] 4.3 Classify and statistically analyze the calibrated three-dimensional droplet size dataset D_corr according to the dynamic binning interval D_edges: V_bin(j) = (π / 6)·∑ D_corr_i ∈ [D_edges(j)， D_edges(j+1)] D_corr_i 3 ; where V_bin(j) is the total volume of all droplets in the j-th interval, in units of mm 3 ; D_corr_i is the calibrated droplet diameter. Calculate the droplet volume distribution matrix V_bin = {103.5, 487.2, 1258.6, 2745.3, 4023.8, 3952.1, 2845.7, 1532.4, 845.9, 312.8}, in units of mm 3 .

[0186] 4.4 Calculate the spatial distribution matrix of rainfall intensity.

[0187] 4.4.1 Apply the volume-rainfall intensity transformation model to calculate the initial rainfall intensity data: R_init = k·(∑ j=1 10V_bin(j))·(1 / A_sample)·(1 / T_sample); where R_init is the initial rainfall intensity value, with the unit of mm / h; k is the conversion coefficient, with a value of 3600; A_sample is the sampling area, with a value of 0.3 m 2 ; T_sample is the sampling time, with a value of 0.2 s. The calculated initial rainfall intensity data is R_init = 182.5 mm / h.

[0188] 4.4.2. Correct the initial rainfall intensity data for environmental factors: R_corr = R_init·(1 + α_v·(V - V_ref) / V_ref)·(1 + α_t·(T - T_ref) / T_ref)·(1 + α_h·(H - H_ref) / H_ref); where R_corr is the corrected rainfall intensity data; α_v, α_t, and α_h are the correction coefficients for wind speed, temperature, and humidity, with values of 0.2, -0.1, and 0.15 respectively; V_ref, T_ref, and H_ref are the reference wind speed, temperature, and humidity, with values of 3.0 m / s, 20°C, and 80% respectively.

[0189] Specific calculation of the correction of environmental factors on rainfall intensity: R_corr = R_init·(1 + α_v·(V - V_ref) / V_ref)·(1 + α_t·(T - T_ref) / T_ref)·(1 + α_h·(H - H_ref) / H_ref) = 182.5·(1 + 0.2·(3.5 - 3.0) / 3.0)·(1 + (-0.1)·(25 - 20) / 20)·(1 + 0.15·(85 - 80) / 80) = 182.5·(1 + 0.2·0.5 / 3.0)·(1 + (-0.1)·5 / 20)·(1 + 0.15·5 / 80) = 182.5·(1 + 0.033)·(1 - 0.025)·(1 + 0.009) = 182.5·1.033·0.975·1.009 = 182.5·1.016 = 185.4 mm / h; Through this correction, the initial rainfall intensity of 182.5 mm / h becomes 185.4 mm / h after being adjusted by environmental factors, an increase of approximately 1.6%.

[0190] 4.4.3. Map the corrected rainfall intensity data to the spatial coordinate system and construct a rainfall intensity spatial grid: First, divide the sampling area into a 10×10 grid, with each grid size of 0.3 m×0.3 m. Then, calculate the rainfall intensity value R_grid(i, j) at each grid point (i, j) according to the distribution density of droplets in space.

[0191] 4.4.4. Optimize the interpolation of the rainfall intensity spatial grid using an improved Kriging interpolation algorithm: R_spatial(x, y) = ∑ i=1 10 ∑ j=1 10 w_ij(x, y)·R_grid(i, j); where R_spatial(x, y) is the interpolated rainfall intensity value at the spatial position (x, y); w_ij(x, y) is the weight coefficient, calculated by the Kriging interpolation algorithm. The rainfall intensity spatial distribution matrix R_spatial is calculated, which characterizes the rainfall intensity distribution in the entire flood discharge atomization area. The range of the rainfall intensity value is [85.3 mm / h, 213.5 mm / h].

[0192] To verify the effectiveness of this embodiment, compare this embodiment with the traditional method: the average relative error of the extraction accuracy of the droplet deformation characteristics in the traditional method is 18.5%, and the average relative error of the extraction accuracy of the droplet deformation characteristics in this embodiment is 6.7%, with a 63.8% improvement; the Pearson correlation coefficient between the accuracy of the particle size distribution after droplet volume correction and the experimental observation in the traditional method is 0.76, and the Pearson correlation coefficient between the accuracy of the particle size distribution after droplet volume correction and the experimental observation in this embodiment is 0.94, with a 23.7% improvement; the relative prediction error of the rainfall intensity prediction accuracy in the traditional method is 24.2%, and the relative prediction error of the rainfall intensity prediction accuracy in this embodiment is 8.5%, with a 64.9% improvement. The comparison results show that this embodiment is superior to the traditional method in terms of the extraction of droplet deformation characteristics, volume correction, and rainfall intensity prediction, and can more accurately predict the rainfall intensity distribution in the flood discharge atomization area.

[0193] In view of the problem of lack of physical constraints in the analysis of droplet deformation characteristics, an adaptive deformation feature decomposition (PCAFD) method with physical constraints of fluid is proposed. By constructing a fluid constraint matrix including a flow field direction matrix, a scale-sensitive matrix and a Reynolds number response matrix, the laws of fluid dynamics are introduced as prior constraints into the feature decomposition process. The frequency domain transformation and multi-scale filtering decomposition technology are used to decouple the effects of different flow mechanisms (mainstream traction, vortex induction and turbulent pulsation) on droplet deformation, realizing the accurate correspondence between deformation characteristics and physical mechanisms, and laying a solid physical foundation for subsequent analysis. Secondly, in view of the problem that the static geometric assumption cannot accurately describe the dynamic volume of droplets, a droplet volume correction method based on physical deformation modes is developed. The true shape of the droplet is reconstructed through physical deformation modes, the deformation state of the droplet is evaluated by combining the modal weight distribution and the droplet stability index, and the accurate correction of the droplet volume is realized by establishing a single-mode correction function and considering the nonlinear coupling term between modes. Especially for droplets close to the critical state of breakup, the volume distribution changes caused by complex nonlinear deformation can be captured, thus eliminating the systematic deviation in traditional methods. Aiming at the problem of insufficient rain intensity prediction accuracy caused by fixed binning and simplified density estimation, a rain intensity calculation method with dynamic binning and density estimation is proposed. Based on the corrected droplet size data, the kernel density estimation algorithm is used to obtain the droplet size density distribution, and the adaptive binning algorithm is used to generate a dynamic binning interval reflecting the droplet distribution characteristics, and the rain intensity spatial distribution is predicted by combining the environmental characteristic parameters. In addition, the rain intensity spatial grid is optimized by an improved Kriging interpolation algorithm, which can accurately describe the fine spatial distribution characteristics of rain intensity in complex environments. The present invention constructs a complete technical system from droplet image acquisition to rain intensity spatial distribution prediction, realizes the accurate capture of the dynamic characteristics of flood discharge atomization droplets and the analysis of physical mechanisms, improves the accuracy and spatial resolution of rain intensity prediction, and provides strong technical support for the safety assessment and environmental impact analysis of water conservancy projects.

[0194] The preferred embodiments of the present invention have been described in detail above. However, the present invention is not limited to the specific details in the above embodiments. Within the scope of the technical concept of the present invention, various equivalent transformations can be made to the technical solutions of the present invention, and these equivalent transformations all belong to the protection scope of the present invention.

Claims

1. A method for predicting flood discharge atomized rainfall intensity based on fog rain droplet distribution, characterized in that: The following steps are involved: Collect the original drop spectrum image and environmental sensor data at the flood discharge site, pre-process them, and obtain the pre-processed image and the environmental characteristic parameter set including wind speed, wind direction, temperature, humidity, atmospheric pressure, turbulence pulsation characteristics and turbulence intensity; Perform patch processing on the preprocessed image to obtain the predetermined type of patches, and convert them into standardized effective patches to form an effective circular patch dataset of raindrops; Based on the effective circular spot data set and the environmental characteristic parameter set, the adaptive deformation feature decomposition method with fluid physics constraints is used to perform droplet deformation dynamics modeling and particle size correction to obtain the corrected three-dimensional droplet particle size data. The corrected three-dimensional droplet size data is dynamically classified and the droplet volume is calculated to form a rainfall intensity spatial distribution matrix, i.e., the rainfall intensity prediction result.

2. The method according to claim 1, characterized in that The steps of obtaining corrected three-dimensional droplet size data include: Combined with the environmental characteristic parameter set, the mapping relationship from two-dimensional circular spots to three-dimensional droplets is established for the effective circular spot data set, and the initial three-dimensional particle size data is obtained. Based on this, the droplet Weber number and deformation index are calculated to form the droplet initial deformation feature set; Applying a feature enhancement function to the initial deformation feature set of the droplet to amplify the critical deformation feature and obtain an enhanced deformation feature set; The enhanced deformation feature set is decomposed into multi-scale flow fields to obtain the flow decoupling feature set, including mainstream drag, vortex induction and turbulent pulsation flow characteristics; The fluid constraint matrix is ​​constructed in combination with the environmental characteristic parameter set, and is combined with the flow decoupling feature set to perform constraint feature decomposition, reconstruct the physical deformation mode of the droplet, and obtain the physical deformation mode set. Based on this, the droplet stability index and volume correction coefficient are calculated, and the corrected three-dimensional droplet particle size data is output.

3. The method according to claim 2, characterized in that The steps to construct the fluid constraint matrix include: Extract wind speed and direction, and combine with empirical fluid dynamics model to construct flow field direction matrix that characterizes the dominant deformation direction of droplets in the flow field; Combined with the turbulence intensity, a scale-sensitive matrix is ​​constructed to reflect the sensitivity of droplets of different sizes to shear force. Based on temperature, humidity and atmospheric pressure, the local flow field Reynolds number is calculated, and a Reynolds number response matrix is ​​constructed to characterize the deformation response characteristics of the droplet under different Reynolds number conditions. The flow field direction matrix, scale sensitivity matrix and Reynolds number response matrix are integrated into a unified fluid constraint matrix by using tensor product operation.

4. The method according to claim 2, characterized in that: The steps of obtaining the enhanced deformation feature set include: Based on the droplet Weber number and deformation index in the initial deformation feature set of the droplet, the critical threshold of the droplet deformation is calculated in combination with the fluid splitting theory, and the critical deformation threshold is output; A nonlinear feature enhancement function is constructed and applied to the initial deformation feature set of the droplet, and the features close to the critical deformation threshold are selectively amplified to generate an enhanced deformation feature set; A curvature-sensitive operator that adapts to the deformation of the droplet surface is constructed and applied to the enhanced deformation feature set, the feature expression near the critical point is optimized, and the enhanced deformation feature set is output.

5. The method according to claim 2, characterized in that: The steps to obtain the flow decoupling feature set include: Applying wavelet transform to the enhanced deformation feature set to convert the spatiotemporal domain features into frequency domain feature representation; Identify the frequency bands of flow regimes corresponding to mainstream drag, vortex induction, and turbulent pulsation based on wind speed and turbulence intensity; The frequency domain feature representation is decomposed by multi-scale filtering using the flow mechanism frequency band, and the decomposed feature subsets corresponding to different flow mechanisms are obtained and flow mechanism correlation analysis is performed on them. The correlation mapping between droplet deformation characteristics and flow field mechanisms is established, and the flow mechanism feature mapping is obtained. Based on it, the features are reorganized and optimized, and the flow decoupling feature set is output.

6. The method according to claim 2, characterized in that The steps to obtain the physical deformation mode set include: Combine the flow decoupling feature set with the fluid constraint matrix, perform eigendecomposition with physical constraints, obtain the eigenvector set, sort the eigenvalues ​​and screen the contribution rate, retain the principal component eigenvector set whose energy proportion exceeds the preset threshold, assign physical meanings one by one, and form a physical mode mapping matrix; The physical mode mapping matrix is ​​used to project and reconstruct the flow decoupling feature set to generate a physical deformation mode set that characterizes the droplet deformation dynamics.

7. The method according to claim 2, characterized in that The steps to calculate the droplet stability index include: Calculate the weight distribution of each mode in the physical deformation mode set to obtain the modal weight distribution; Based on the modal weight distribution and the physical meaning of the mode, the droplet stability index that characterizes the stability of droplet deformation is calculated and compared with the critical stability threshold to evaluate the possibility of droplet splitting represented by the physical deformation mode set. Based on the possibility of droplet splitting, the dynamic instability coefficient that characterizes the dynamic stability of the droplet in the flow field is calculated; The droplet stability index and dynamic instability coefficient are combined to generate the droplet stability evaluation result.

8. The method according to claim 7, characterized in that The steps of calculating the volume correction coefficient and outputting the corrected three-dimensional droplet size data include: Based on the physical deformation mode set and droplet stability evaluation results, a single-mode correction function is constructed for each physical deformation mode. Analyze the interaction between different modes in the physical deformation mode set, generate a modal coupling matrix that represents the mutual influence between modes, and calculate the nonlinear coupling terms between modes based on it; Combining the single-mode correction function and the nonlinear coupling term, the volume correction coefficient of the droplet is comprehensively calculated and the initial three-dimensional particle size data is corrected accordingly to obtain the corrected three-dimensional droplet size data.

9. The method according to claim 1, characterized in that: The steps of dynamically dividing and counting the corrected three-dimensional droplet size data and calculating the droplet volume to form a rainfall intensity spatial distribution matrix include: The kernel density of the corrected three-dimensional droplet size data is estimated to obtain the droplet size density distribution and classify it, generating a dynamic classification interval that reflects the droplet distribution characteristics. Based on this, the corrected three-dimensional droplet size data is classified and counted, the total volume of the droplets in each interval is calculated, and the droplet volume distribution matrix is ​​constructed. Combining the droplet volume distribution matrix with wind speed and temperature factors, the rainfall intensity conversion model is applied to calculate and output the rainfall intensity spatial distribution matrix.

10. The method according to claim 9, characterized in that The steps of applying the rainfall intensity conversion model to calculate and output the rainfall intensity spatial distribution matrix include: Based on the droplet volume distribution matrix, the volume-rainfall intensity conversion model derived from the principles of fluid mechanics is used to calculate the initial rainfall intensity data; Combined with wind speed, turbulence intensity, temperature and humidity, the initial rainfall intensity data is corrected for environmental factors to obtain the corrected rainfall intensity data; According to the spatial distribution of sampling points, the corrected rainfall intensity data is mapped to a three-dimensional spatial coordinate system to construct a rainfall intensity spatial grid; The rainfall intensity spatial grid is interpolated and optimized to generate a continuous rainfall intensity spatial distribution matrix as the final rainfall intensity prediction result in the flood discharge atomization area.

Images