Rainfall analysis method and system based on micro rain radar data

Abstract

The invention discloses a rainfall analysis method and system based on micro-rain radar data, and belongs to the field of meteorological radar signal processing, and the method comprises the steps: carrying out the vertical detection of a target airspace through a micro-rain radar, and obtaining a radar echo signal; performing noise filtering on the original Doppler power spectrum density matrix to obtain a cleaned observation Doppler power spectrum; performing feature extraction on the observation Doppler power spectrum by using a pre-trained neural network model, and outputting prior probability distribution about rainfall types; constructing and solving a precipitation state unified inversion model to obtain a globally optimal unified precipitation state vector; and based on the unified rainfall state vector, calculating rainfall particle number density distribution of each height layer, and carrying out integration to obtain an instantaneous rainfall rate of each height layer. According to the method, the problem of attenuation correction error accumulation in a traditional method is solved, and accurate early warning of rainfall dynamic risks is realized.

Description

Technical Field

[0001] This invention relates to the field of meteorological radar signal processing, and in particular to a method and system for rainfall analysis based on micro-rain radar data. Background Technology

[0002] Rain radar is a vertically pointing Doppler radar that can provide high spatiotemporal resolution information on the vertical profiles of precipitation particles and is commonly used for meteorological observation.

[0003] Currently, traditional analysis methods are easily affected by data quality. The raw spectral data acquired by micro-rain radar is highly susceptible to contamination from non-meteorological echoes (e.g., ground clutter, echoes from flying objects like insects, and electromagnetic noise from both the system's internal and external environments). Furthermore, during heavy rainfall, radar signals experience path integral attenuation (PIA) along their propagation path, with the attenuation becoming more severe with increasing rainfall intensity. Current methods for addressing data contamination and attenuation typically involve simplification or complete neglect; however, this approach often leads to significant distortion of the calculated radar reflectivity Z-value, reducing the reliability and stability of the final rainfall estimate. Moreover, current analyses often assume the detected target is liquid raindrops. However, in reality, a large amount of solid precipitation (e.g., snow, ice crystals) exists at high altitudes during weather events, transforming into liquid raindrops as it descends through the 0°C layer (melting layer). Particularly in the melting layer region, melting ice crystals form a 'bright band' encased in a water film, and the radar scattering characteristics of this bright band differ significantly from those of purely liquid or solid particles. However, current analysis generally applies the ZR relationship established for liquid precipitation to the entire probe profile indiscriminately. This method can lead to serious estimation errors near the height of the melt layer, resulting in an inaccurate reflection of the actual precipitation physical processes.

[0004] Therefore, overcoming the above-mentioned shortcomings and improving the accuracy, stability and environmental adaptability of rainfall analysis based on micro-rain radar data is a technical problem that urgently needs to be solved in this field. Summary of the Invention

[0005] One of the objectives of this invention is to provide a rainfall analysis method based on micro-rain radar data, in order to solve the problems of low accuracy, poor robustness and weak adaptability of rainfall analysis in the prior art, which rely on fixed ZR relationships, ignore data quality issues and cannot handle complex precipitation phases.

[0006] This invention is achieved through the following technical solution: a rainfall analysis method based on micro-rain radar data, comprising the following steps: S100, using micro-rain radar to perform vertical detection of the target airspace and acquire radar echo signals; S200, filtering noise from the original Doppler power spectral density matrix to obtain the cleaned observed Doppler power spectrum; S300, using a pre-trained neural network model to extract features from the observed Doppler power spectrum and outputting a prior probability distribution regarding precipitation type; S400, constructing and solving a unified precipitation state inversion model, wherein the unified precipitation state inversion model includes: constructing a unified precipitation state vector, the state vector including parameters describing the microphysical structure of precipitation and mixing coefficients describing the phase weights of precipitation, and constructing a variational inversion-based model. A global inversion model is constructed, comprising a forward physics operator configured to calculate the theoretical Doppler power spectrum considering path integral attenuation based on a hypothetical unified precipitation state vector, and constructing an objective functional. The objective functional includes at least a data fidelity term for constraining the consistency between the theoretical and observed spectra, and a structured prior term for constraining the inversion solution to approximate the prior probability distribution. The extreme values ​​of the objective functional are solved using a numerical optimization algorithm to obtain the globally optimal unified precipitation state vector. S500: Based on the unified precipitation state vector, the precipitation particle number density distribution at each altitude layer is calculated, and the instantaneous rainfall rate at each altitude layer is obtained by integrating the precipitation particle number density distribution, particle volume, and particle falling velocity.

[0007] Furthermore, the objective functional can be expressed as follows: in, To obtain the globally optimal state vector field. For the overall objective functional, that is, to enable The state vector field that has reached its minimum value. This is a data fidelity term used to measure the Euclidean distance between the predicted spectrum and the observed spectrum; To smooth the regularization term, the smoothness of the vertical gradient is constrained by calculating the squared magnitude of the derivative of the state vector with respect to the height. For structured priors, The regularization hyperparameter for smoothing the regularization term; These are the regularization hyperparameters for the structured priors; these two regularization hyperparameters are used to balance the degree of data fit, the smoothness of the solution, and the weights between prior constraints, and can be set by the L-Curve method or empirically.

[0008] Furthermore, the instantaneous rainfall rate at each altitude level can be calculated using the following formula: in, In order to be at a high altitude Instantaneous rainfall rate at the location; The formula for the volume of a sphere Geometric coefficients in; To obtain the optimal raindrop spectrum, substitute the optimal parameters obtained from the inversion. The subsequent distribution function; This represents the falling speed of liquid raindrops when calculating rainfall rate.

[0009] Furthermore, the rainfall analysis method also includes: S600, constructing a dynamic rainfall risk early warning process, the process including: calculating the time derivative of the instantaneous rainfall rate on a continuous time series to obtain a burst acceleration index characterizing the rate of change of precipitation intensity; and / or, calculating the spatial gradient of the instantaneous rainfall rate on a vertical height layer to obtain a microphysical evolution factor characterizing the vertical evolution characteristics of precipitation particles.

[0010] Furthermore, the rainfall analysis method also includes: S700, calculating a comprehensive risk index, wherein the comprehensive risk index is a dimensionless value constructed based on the real-time ground rainfall rate, the burst acceleration index, and the microphysical evolution factor; wherein the burst acceleration index acts on the comprehensive risk index in an exponential form to nonlinearly amplify the influence of the precipitation enhancement trend.

[0011] Furthermore, the comprehensive precipitation dynamics risk index can be calculated using the following formula: in, For the ground layer ( Real-time rainfall rate; For reference heavy rainfall thresholds; Nonlinear amplification factor (usually) ); The sensitivity coefficient to the outbreak trend. This is to ensure that the risk is amplified only when the rainfall intensifies; The function is used to normalize the effect of the vertical gradient when When it is negative (due to collision and growth), this item is greater than 1, further amplifying the risk value.

[0012] Furthermore, the burst acceleration index is calculated using a weighted difference method within a sliding time window, with greater weight given to values ​​closer to the current moment. When the burst acceleration index exceeds a preset positive threshold, the precipitation is determined to be in a burst enhancement phase, and an early warning is triggered.

[0013] Furthermore, the precipitation burst acceleration index can be calculated using the following formula: in, The result obtained by integration in Example 2 Rainfall rate at any given time; The sampling time interval; The size of the sliding window; The weight is determined by time decay, with a greater weight closer to the current time.

[0014] Furthermore, the microphysical evolution factor is used to correct the precipitation forecast, specifically including: when the microphysical evolution factor shows that the upper-layer rainfall rate is significantly greater than the lower-layer rainfall rate, it is determined that there is evaporation loss, and the surface precipitation forecast is lowered; when the microphysical evolution factor shows that the upper-layer rainfall rate is significantly less than the lower-layer rainfall rate, it is determined that there is collision and growth, and the surface precipitation forecast is maintained or raised.

[0015] Furthermore, the vertical microphysical evolution factor can be calculated using the following formula: in, For height Vertical microphysical evolution factor; For height increments.

[0016] Furthermore, the radar echo signal is the original Doppler power spectral density matrix that dynamically changes over time by setting the detection period to the second level and dividing the vertical detection path into multiple height layer range libraries. Each data point in the matrix represents the echo power density at a specific height and a specific Doppler velocity.

[0017] Furthermore, in S200, the path attenuation characteristics of the signal when penetrating the precipitation area are preserved, and no path attenuation correction is performed; the noise filtering includes: transforming the original Doppler power spectral density matrix to the time-frequency domain; based on the difference between the correlation characteristics of the precipitation signal in the time-frequency domain and the randomness characteristics of the noise, an adaptive threshold is constructed to filter the transformed coefficients; wherein, the denoising process does not include layer-by-layer attenuation compensation operation for the radar beam propagation path.

[0018] Furthermore, the observed Doppler power spectrum can be expressed by the following formula: in, This represents the actual observed Doppler power spectrum with a high signal-to-noise ratio after reconstruction; and These represent the discrete wavelet transform and its inverse transform, respectively. The original power spectrum; This is a thresholding operator that operates based on a preset noise floor level. High-frequency wavelet coefficients below this energy level are removed, thereby outputting a high signal-to-noise ratio observed Doppler power spectrum. Specifically, in this embodiment, the threshold processing operator can be constructed by the following steps: Furthermore, the threshold processing operator can be constructed by the following steps: First, in this embodiment, to avoid interference from heavy precipitation signals in noise estimation, the median absolute deviation estimation method can be used to obtain a threshold for the background noise level. The system first makes an accurate estimate of the noise standard deviation; then, based on the estimated noise standard deviation, it calculates the optimal threshold for the current height layer using a general thresholding criterion; finally, it employs a soft thresholding operator that not only sets coefficients below the preset noise floor level to zero, but also shrinks coefficients above the preset noise floor level toward zero, thereby ensuring the smoothness of the reconstructed signal.

[0019] Further, S300 includes the following sub-steps: S310, extracting the macroscopic feature profile of the preprocessed observation spectrum; S320, inputting the macroscopic feature profile into a neural network with the output layer configured as a probability output; S330, outputting a probability vector containing multiple dimensions, the value of each dimension representing the confidence level of different precipitation phases, wherein the precipitation phases include at least stratiform cloud precipitation, convective precipitation, and solid precipitation.

[0020] Furthermore, the prior probability distribution can be calculated using the following formula: , in, For height The input feature vector at the location; This represents the forward inference function of a convolutional neural network; The function ensures that the sum of the probabilities of the output is 1; This represents the prior probability of precipitation from stratiform clouds at that altitude. This represents the prior probability of convective precipitation at that altitude; This represents the prior probability of solid precipitation at that altitude.

[0021] Furthermore, the mixing coefficient in the unified precipitation state vector in S400 satisfies the normalization condition and is used to characterize the coexistence probability of different precipitation phases at the current altitude layer; when calculating the theoretical Doppler power spectrum, the forward physical operator uses the mixing coefficient to perform a weighted summation of the scattering cross section and falling velocity of different precipitation phases.

[0022] Furthermore, when calculating the theoretical Doppler power spectrum, the forward physics operator is also configured to perform the following operations: calculate the specific attenuation coefficient of the radar wave when penetrating the precipitation area based on the unified precipitation state vector; integrate the specific attenuation coefficient along the radar beam propagation path to obtain the two-way path integral attenuation value for each height layer; and use the two-way path integral attenuation value to simulate the attenuation of the calculated theoretical echo signal.

[0023] Furthermore, the theoretical Doppler power spectrum can be calculated using the following formula: in, Given state Below, the model predicts at altitude The predicted Doppler velocity is The power spectral density; It is the path integral decay operator; It is an index of precipitation types; for type At height Mixed weights; The equivalent sphere diameter of the precipitation particles; These are the lower and upper limits of integration, respectively, defining the physical range of particle size; It is a diameter of Particles in type The radar backscattering cross section is a physical constant function that depends on the particle size. and type ; It is the raindrop spectral distribution function, which is derived from the uniform precipitation state vector. A gamma distribution with defined parameters represents a distribution with a diameter of [missing information] per cubic meter. The number of particles; It is a Gaussian convolution kernel function used to simulate the spectral broadening effect caused by air turbulence; It is the Doppler velocity variable obtained from radar observation; It is the terminal velocity of a particle falling under specific conditions, representing a diameter of... Type is Particles at high altitude The theoretical falling speed.

[0024] Furthermore, the path integral decay operator can be expressed as follows: in, It is a natural exponential function used to describe the exponential decay of electromagnetic waves in lossy media; For vertical integration variables, representing Any height above; For height The specific attenuation coefficient at that location.

[0025] Another aspect of the present invention provides a rainfall analysis system for light rain radar data, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the rainfall analysis method for light rain radar data as described above.

[0026] Compared with the prior art, the present invention has the following advantages and beneficial effects: 1. This invention abandons the traditional attenuation correction in the preprocessing stage, and instead incorporates path integral attenuation into the forward physical model for unified solution. This closed-loop feedback mechanism ensures the physical consistency between attenuation processing and microphysical parameter inversion, and avoids irreversible truncation errors caused by open-loop correction.

[0027] 2. This invention introduces the prior probability distribution predicted by deep learning as a soft constraint into the objective functional of physical inversion through a variational inference framework. This not only utilizes the fast reasoning ability of deep learning models to prevent physical solutions from getting trapped in local extrema, but also preserves the rigor of physical models, achieving a balance between interpretability and accuracy.

[0028] 3. This invention is not limited to the analysis and monitoring of static rainfall. Instead, by introducing burst acceleration in the time dimension and microphysical evolution factors in the spatial dimension, it can identify the trend of rainstorm outbreaks and the authenticity of aerial rainfall earlier, providing more forward-looking decision support for urban disaster prevention and mitigation. Attached Figure Description

[0029] The accompanying drawings, which are included to provide a further understanding of embodiments of the invention and form part of this application, do not constitute a limitation thereof. In the drawings: Figure 1 The above is a flowchart of the overall method provided in Embodiment 1 of the present invention.

[0030] Figure 2 This is a timing diagram of the overall method provided in Embodiment 1 of the present invention.

[0031] Figure 3 This is a vertical cross-sectional view of the micro-rain radar provided in Embodiment 1 of the present invention.

[0032] Figure 4 The vertical profile of the microphysical feature and the prior probability distribution diagram provided in Embodiment 1 of the present invention.

[0033] Figure 5 This is a schematic diagram of the variational inversion fitting principle provided in Embodiment 1 of the present invention.

[0034] Figure 6 This is a schematic diagram illustrating the early warning system provided in Embodiment 2 of the present invention. Detailed Implementation

[0035] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. The components of the embodiments of the present invention described and shown in the accompanying drawings can generally be arranged and designed in various different configurations.

[0036] Currently, rainfall analysis using micro-rain radar data typically involves the following steps: First, acquiring raw Doppler power spectrum data using micro-rain radar equipment; then, calculating basic parameters such as radar reflectivity factor (Z) and average particle falling velocity based on this data; finally, using a fixed empirical ZR relationship (i.e., the power-law relationship between radar reflectivity Z and rainfall rate R, such as Z = aR)... b The ZR relationship (where a and b are empirical coefficients) is used to estimate the surface rainfall rate. However, this method suffers from poor model adaptability, leading to insufficient estimation accuracy. In particular, the empirical coefficients (a, b) of the ZR relationship are usually long-term statistical results for specific regions and rainfall types. In actual rainfall observations, the microphysical properties of precipitation, especially the droplet spectral distribution (DSD), change significantly with rainfall type (e.g., stratiform cloud precipitation versus convective precipitation), geographical location, or seasonal variations. This method, using a single fixed ZR relationship, cannot dynamically adapt to these changes, resulting in a large error between the final estimated rainfall and the actual value, making it difficult to meet the overall accuracy requirements of meteorological services.

[0037] To address the problems encountered in existing technologies that analyze rainfall using the ZR relationship, this application provides a rainfall analysis method based on micro-rain radar data. The specific solutions of this invention will be described in detail below with reference to specific embodiments and accompanying drawings.

[0038] Example 1 Figure 1 The overall method flowchart in this embodiment is shown. Figure 2 The overall method timing diagram of this embodiment is shown. As can be seen from the diagram, this embodiment includes the following steps: Step 1: First, a micro-rain radar is used to vertically probe the target airspace, acquiring high spatiotemporal resolution data and receiving echo signals scattered by atmospheric precipitation particles. Simultaneously, to capture transient changes in the precipitation process, the system is configured to output high-precision raw spectral data. For example, the observation period can be set to the second level (e.g., 10 to 60 seconds), and the vertical probe path can be divided into several altitude layers (range library) to obtain the raw Doppler power spectral density matrix that dynamically changes over time.

[0039] Specifically, in this embodiment, the system can adopt the following configuration: The temporal resolution is adjustable, with an observation period ranging from 10 seconds (severe convective weather) to 60 seconds (stratified cloud precipitation), enabling 6-60 continuous observation samples per minute. Spatial stratification involves dividing the 0-6 km detection altitude into 128 range bins, each 47 meters thick, achieving high-density sampling in the vertical direction. Data output directly outputs the raw Doppler power spectral density matrix, containing 128 (range bins) × 512 (velocity channels) data points, achieving a velocity resolution of 0.1 m / s. Spectral processing involves acquiring 256 complex sampling points from each range bin and obtaining the Doppler spectrum through a 1024-point FFT transform.

[0040] The raw data matrix obtained through this configuration has a three-dimensional structure: time dimension (second level), height dimension (128 layers), and velocity dimension (512 channels), thus providing complete spectral distribution information for subsequent inversion of precipitation microphysical parameters.

[0041] Step 2: Next, the raw data is cleaned. In this step, only noise filtering is performed, without path attenuation correction. This is because the subsequent inversion model contains a complete physical forward operator, which can handle the attenuation problem more accurately in a closed loop. If open-loop iterative correction is performed during the data cleaning stage, it may introduce irreversible truncation errors, leading to large errors in the final inversion process.

[0042] Specifically, in this embodiment, wavelet thresholding denoising technology can be used to clean the original data. The Doppler spectrum is processed using discrete wavelet transform, decomposing the Doppler spectrum signal into the time-frequency domain. Utilizing the strong correlation of precipitation signals in the time-frequency domain, while random thermal noise exhibits no correlation, the irregular noise coefficients are set to zero. Specifically, an adaptive thresholding operator can be constructed that preserves and enhances high-amplitude wavelet coefficients while eliminating low-amplitude noise coefficients. On the other hand, since the path attenuation characteristics of the signal must be preserved, traditional layer-by-layer attenuation correction must be strictly prohibited at this stage. This is because the subsequent inversion model in this embodiment contains a fully physical closed loop that can automatically solve for attenuation; introducing empirical hard correction in the preprocessing stage would disrupt physical consistency and introduce irreversible truncation errors.

[0043] For example, in this embodiment, the denoised observed Doppler power spectrum can be expressed by the following formula: in, This represents the actual observed Doppler power spectrum with a high signal-to-noise ratio after reconstruction; and These represent the discrete wavelet transform and its inverse transform, respectively. The original power spectrum; This is a thresholding operator that operates based on a preset noise floor level. High-frequency wavelet coefficients below this energy level are removed, thereby outputting a high signal-to-noise ratio observed Doppler power spectrum. Specifically, in this embodiment, the threshold processing operator can be constructed by the following steps: First, in this embodiment, to avoid interference from heavy precipitation signals in noise estimation, the median absolute deviation (MAD) estimation method can be used to obtain a value for the background noise level. A precise estimate.

[0044] An accurate estimate of the background noise level can be obtained by the following formula: in, These are the detail coefficients of the first layer (i.e., the highest frequency layer) of wavelet decomposition, which mainly contains noise; It is the adjustment factor for the normal distribution; This is for median calculation.

[0045] Then, based on the estimated noise standard deviation, the optimal threshold for the current height layer is calculated using a general thresholding criterion. Specifically, it can be shown in the following formula: in, This formula ensures that Gaussian white noise is filtered out with a very high probability, given the total number of sampling points in the current Doppler spectrum data.

[0046] Finally, to avoid pseudo-Gibbs oscillations caused by hard thresholding at coefficient cutoff points, this embodiment employs a soft thresholding operator. This soft thresholding operator not only reduces the threshold value below a certain threshold, but also... Setting the coefficient to zero will still result in a value higher than [the coefficient]. The coefficients shrink towards zero, thus ensuring the smoothness of the reconstructed signal. Specifically, this can be calculated using the following formula: in, It is a symbolic function; The process is applied to each detail coefficient after wavelet decomposition, and then an inverse transform is performed. Recovering the time domain .

[0047] It should be noted that the above three sub-steps are used to construct... The operator achieves mathematical unification of noise adaptive estimation and signal smoothing preservation. Compared with simple frequency truncation, this operator based on statistical principles can effectively distinguish between weak rain signals (which have low amplitude but wavelet domain correlation) and system noise (which has low amplitude and no correlation), thereby improving the signal-to-noise ratio and physical accuracy of the input data of the subsequent inversion model. Figure 3 The image shows a vertical cross-sectional view of the rain radar after noise reduction in this embodiment.

[0048] Step 3: This step uses a deep neural network to initially extract the macroscopic features of the observed Doppler power spectrum, generating a prior probability distribution on precipitation type to guide the subsequent physical optimization process and provide prior data for subsequent physical inversion.

[0049] Specifically, the purpose of this step is to leverage the rapid reasoning capabilities of the data-driven model to provide a high-quality initial guide for the subsequent complex physical inversion process, preventing the physical solution from getting trapped in local optima. This is achieved by extracting vertical profile features of radar reflectivity factor, Doppler velocity, and spectral width, and then inputting the extracted data into a lightweight one-dimensional convolutional neural network (1D-CNN). This network is configured not to directly output hard yes / no classification labels, but rather a probability vector describing the likelihood of precipitation phase. This soft constraint approach allows the physical model to fine-tune or even correct this prior information in subsequent steps based on actual observations.

[0050] Specifically, in this embodiment, the preliminary extraction to generate a prior probability distribution about precipitation type may include the following steps: First, key moment features of the observed Doppler power spectrum (e.g., reflectivity factor Z, Doppler velocity v, spectral width σ) are extracted to construct a feature profile.

[0051] It is then fed into a pre-trained lightweight one-dimensional convolutional neural network (1D-CNN), whose output layer is designed as a softmax probability layer, outputting a multi-dimensional probability vector corresponding to the confidence levels of stratiform clouds, convective clouds, and solid precipitation, respectively.

[0052] For example, in this embodiment, the prior probability distribution can be calculated using the following formula: in, For height The input feature vector at that location (e.g., may contain...) ); This represents the forward inference function of a convolutional neural network; The function ensures that the sum of the probabilities of the output is 1. These represent the prior probabilities of precipitation at that altitude being stratiform cloud precipitation, convective precipitation, and solid precipitation, respectively. Figure 4 The diagram shows the vertical profile of the microphysical features and the prior probability distribution in this embodiment. In the diagram, the three small plots on the left show the reflectivity (Z), velocity (V), and spectral width (W) as they change with altitude; the chart on the right shows the probability distribution of the 1D-CNN output, which shows that the probability of rain and snow in the melting layer alternates at around 2250m.

[0053] Step 4: Then construct and solve the unified inversion model of precipitation state based on variational inference. This model can find the precipitation state that best matches the observation data under the dual constraints of physical laws and AI priors.

[0054] Specifically, in this embodiment, the unified inversion model of precipitation state based on variational inference can be constructed through the following steps: 1) First, it is necessary to transform the physical state of complex precipitation at any altitude along the vertical profile of the atmosphere into a standardized state vector that includes mixed phase information and microphysical properties, i.e., a unified precipitation state vector.

[0055] In this embodiment, a state space capable of simultaneously describing the raindrop spectral distribution characteristics and precipitation phase probability weights is constructed by introducing a 'continuous mixing' approach. This eliminates the need for discrete, rigid precipitation type labels in traditional schemes.

[0056] Specifically, a probability mixing coefficient vector can be used to represent complex precipitation phases (e.g., stratiform precipitation, convective precipitation, and solid precipitation). This is because precipitation processes in the real atmosphere are often continuously changing. By defining a probability vector that satisfies the normalization condition, mixed precipitation and smooth transitions in phases along the vertical direction can be naturally expressed (e.g., a smooth transition from snow at high altitudes to rain at low altitudes), thereby expanding the model's ability to express complex weather scenarios.

[0057] For precipitation microphysical structure, we can assume that the particle size distribution follows a gamma distribution at any altitude. We can use a set of parameters to capture the concentration, average size and dispersion of precipitation particles, thereby defining the core element of precipitation microphysical structure.

[0058] Exemplarily, in this embodiment at high The uniform precipitation state vector at a given location can be defined by the following formula: in, For height Precipitation state vector at location; The set of gamma distribution parameters describing the raindrop spectrum (DSD) at this altitude directly determines the microstructure of precipitation; among them, The intercept parameter is used to control the particle concentration level. Shape parameters for controlling the curvature of the spectrum; The slope parameter is used to control the cutoff velocity of large particles.

[0059] This is a vector of precipitation type mixing coefficients.

[0060] It should be noted that in the above formula: Each component in Represents at a high level In this area, precipitation exhibits different types. (For example, Stratoclouds For convection, The vector represents the weights or probabilities of characteristics of upper-level solid precipitation; this vector must satisfy... The physical constraints allow the model to move beyond simple either / or classifications and enable more precise quantification of the solid and liquid components.

[0061] 2) After defining a unified precipitation state vector, it is necessary to establish a mapping relationship between physical state and radar observation data.

[0062] In this embodiment, a forward physics prediction model based on integral equations can be used to map the physical state to radar observation data. The forward physics prediction model based on integral equations in this embodiment acts as a digital twin engine. Given a hypothetical physical state, it accurately calculates the Doppler power spectrum that the micro-rain radar should theoretically observe. To achieve accurate calculations, the microscopic particle distribution needs to be integrally transformed into a macroscopic radar echo signal through a physical scattering mechanism.

[0063] Specifically, in this embodiment, the mapping process can be described using Fredholm's first kind of integral equation. By considering the significant differences in the scattering characteristics (i.e., scattering cross-section) and dynamic characteristics (i.e., falling velocity) of radar waves from different phases (e.g., liquid raindrops and solid snowflakes), this model must utilize the aforementioned precipitation type mixing coefficient. Weighted summation is performed on different types of physical contributions.

[0064] At the same time, the velocity spectrum broadening effect caused by air turbulence and the path integral attenuation (PIA) of electromagnetic waves when penetrating precipitation areas must also be considered to ensure the physical accuracy of the prediction results.

[0065] For example, in this embodiment, the predicted Doppler power spectrum can be calculated using the following formula: in, Given state Below, the model predicts at altitude The predicted Doppler velocity is The power spectral density; It is the path integral decay (PIA) operator; It is an index of precipitation type (e.g., liquid, solid); for type At height The mixed weights (from the uniform precipitation state vector); The equivalent sphere diameter of the precipitation particles (i.e., the integral variable); These are the lower and upper limits of integration, which define the physical range of particle size (e.g., it can be set from 0.1 mm to 8 mm). It is a diameter of Particles in type The radar backscattering cross section is a physical constant function that depends on the particle size. and type (This determines the dielectric constant and shape of the particle), which can be calculated in advance using Mie scattering theory or the T-matrix method; It is the raindrop spectral distribution function, which is derived from the uniform precipitation state vector. A gamma distribution with defined parameters represents a distribution with a diameter of [missing information] per cubic meter. The number of particles; It is a Gaussian convolution kernel function used to simulate the spectral broadening effect caused by air turbulence; It is the Doppler velocity variable obtained from radar observation; It is the terminal velocity of a particle falling under specific conditions (such as air density), representing a diameter of... Type is Particles at high altitude The theoretical falling speed (affected by air density).

[0066] It should be noted that in the above formula... The path integral attenuation operator can be expressed by the following equation: in, It is a natural exponential function used to describe the exponential decay of electromagnetic waves in lossy media; For vertical integration variables, representing Any height above; For height The specific attenuation coefficient at a given location is determined by the precipitation conditions at that location. The physical quantity that determines this is calculated in a similar way to reflectivity, and is also an integral of DSD.

[0067] Since radar waves propagate in two paths (transmission + return), the coefficient is -2, where altitude is represented in the above formula. The signal attenuation at that point depends on the distance from that height to the top of the radar beam. The cumulative attenuation capacity of all media; by introducing this operator, the model can automatically compensate for signal attenuation, which is crucial for the accurate inversion of micro-rain radar signals in the next step.

[0068] 3) Next, in order to invert the highly complex unknown state from the limited observation data, simple data fitting is insufficient; physical constraints and prior knowledge must also be introduced. In this embodiment, the inversion problem is transformed into a functional extremum problem: finding an optimal state field that perfectly explains the observation data, conforms to the continuity of atmospheric physics, and integrates the intelligent recognition results of deep learning. Minimizing the difference between the predicted spectrum and the observed spectrum is taken as the core objective. By constructing a variational objective functional that includes multiple physical constraints and intelligent priors, the globally optimal precipitation state field is solved.

[0069] Specifically, the target functional in this embodiment should include three parts: First, a data fidelity term to ensure that the predicted power spectrum is highly consistent with the actual observed power spectrum, which is the source of high accuracy; second, a spatial smoothing regularization term to reduce noise and ensure image smoothness, which can force the solution to have physical continuity by penalizing sharp jumps in the vertical direction of the state vector, thereby suppressing artifacts caused by noise; and finally, a structured prior regularization term to ensure deep coupling between the physical model and the intelligence, which can use a data-driven intelligence interface to quickly predict the input data using a pre-trained lightweight neural network (e.g., 1D-CNN) to obtain a prior type distribution, and guide the solution of the physical model to approach the intelligent prior through KL divergence (KL divergence) to achieve deep coupling between the physical model and intelligent recognition.

[0070] For example, in this embodiment, the target functional can be expressed as follows: in, To obtain the globally optimal state vector field. For the overall objective functional, that is, to enable The state vector field that has reached its minimum value. This is a data fidelity term used to measure the Euclidean distance between the predicted spectrum and the observed spectrum; To smooth the regularization term, the smoothness of the vertical gradient is constrained by calculating the squared magnitude of the derivative of the state vector with respect to the height. For structured priors, The regularization hyperparameter for smoothing the regularization term; These are the regularization hyperparameters for the structured priors; these two regularization hyperparameters are used to balance the degree of data fit, the smoothness of the solution, and the weights between prior constraints, and can be set by the L-Curve method or empirically.

[0071] Then, using numerical optimization with gradient descent or L-BFGS algorithm, the objective functional is iteratively solved. After several iterations and convergence, the globally optimal state vector field can be obtained. .

[0072] Additionally, it should be noted that: In this embodiment, the data fidelity item can be represented by the following formula: in, These are the lower and upper limits of the vertical integration, representing the starting height (usually the top of the lowest observation blind zone of the radar) and the ending height (usually the top of the radar echo) of the inversion region, respectively, used to limit the spatial range of the physical inversion. The lower and upper limits of the Doppler velocity integral define the effective observation window of the Doppler spectrum (e.g., it can be set to -12 m / s to +12 m / s), ensuring coverage of all possible particle fall velocities and turbulence broadening. To predict the Doppler power spectrum, which is a theoretical observation calculated by the forward physics model based on the current assumed state vector, it includes simulation results of physical processes such as scattering, attenuation, and velocity broadening. To actually observe the Doppler power spectrum, this was done by a micro-rain radar at a high altitude. and speed The actual power spectral density data collected at the site serves as the true reference for the inversion process. It represents the square of the L2 norm (i.e., the square of the Euclidean distance); using the square form not only eliminates the canceling effect of positive and negative errors, but also imposes a stronger penalty on larger fitting errors, forcing the model to prioritize correcting regions with larger deviations. It is the unified precipitation state vector field to be inverted. This state vector field is the independent variable in the formula, and the entire integral calculation result depends on it. The value of . This represents error accumulation across the entire vertical height layer and the entire velocity spectrum. The above equation is a standard least-squares term used to ensure the Doppler power spectrum predicted by the physical model. Approximate as close as possible to the Doppler power spectrum actually observed by radar. .

[0073] It should be noted that the above formula represents the accuracy of the inversion result. If the value of this term is large, it means that the physical mapping model in sub-step 2) cannot explain the current observation data (for example, the model parameters are set incorrectly, or physical phenomena not covered by the model have occurred). By minimizing this term, the state vector can be forced to fit the original radar observation.

[0074] In this embodiment, the smoothing regularization term can be expressed as follows: in, Represents the vertical height variable; Indicates height The state vector representation at that location; The gradient represents the first derivative of the state vector with respect to altitude, which describes the rate of change of precipitation conditions (such as raindrop concentration and mixing coefficient) in the vertical direction.

[0075] It should be noted that the above equation is used to impose physical continuity constraints, preventing non-physical turbulence caused by observation noise and enhancing the stability of the solution; it acts as a low-pass filter, filtering out high-frequency noise; and it ensures the continuity of the vertical profile. If the state vector undergoes a sharp jump between adjacent heights (i.e., a large gradient), this term will increase dramatically, resulting in a large penalty cost in the overall objective function. By minimizing this term, the algorithm can be made to seek solutions with smooth changes, so that the inverted precipitation structure conforms to the basic characteristics of hydrodynamics, thus obtaining solutions that conform to atmospheric hydrodynamics.

[0076] In this embodiment, the structured prior terms can be represented by the following formula: in, The KL divergence, also known as relative entropy, is a measure of probability distribution. Compared with the reference distribution An asymmetric measure of the difference between them, when When the KL divergence is zero, it is a statistical metric used to measure the difference between two probability distributions. In this formula, it measures the information loss of the physical inversion distribution relative to the prior distribution of the deep learning. This represents the mixing coefficient vector of precipitation types currently retrieved; it represents the probability distribution of precipitation types currently calculated in the physical retrieval cycle (e.g., [stratified cloud probability, convective cloud probability, solid cloud probability]). This represents the intelligent prior probability distribution; it is the probability distribution of precipitation type directly predicted by a pre-trained convolutional neural network (CNN) model based on the input radar echo profile; this distribution does not change with the inversion process and is a fixed reference target.

[0077] It should be noted that the structured prior regularization term is used to fuse the intelligent prediction results of the deep learning model, providing global structured guidance for the inversion process and resolving the multiple solutions in the inversion problem. The value calculated by the above formula reflects the degree to which the physical inversion result deviates from the prior judgment of the deep learning. When the physical observation data ( When the signal-to-noise ratio is low or the features are ambiguous, the constraint of the data fidelity term weakens. In this case, the term dominates the optimization direction and guides the solution. Tend to This allows us to use patterns learned from historical data to fill the information gaps in real-time observations. Figure 5 The diagram illustrates the variational inversion fitting principle in this embodiment. It shows how, at a specific altitude, the model in this embodiment uses a forward operator to fit the theoretical spectrum to the observed spectrum (including noise), thereby achieving data fidelity. The gray scatter dots represent the noisy actual observation data, and the red solid line represents the theoretical spectrum derived from the physical model. The objective functional in this embodiment adjusts the physical parameters (raindrop spectrum, phase coefficient) to make the curve fit the gray dots as closely as possible (data fidelity) while maintaining smoothness (regularization).

[0078] Step 5: After obtaining the globally optimal state vector field through the above steps, since the state vector fully describes the microphysical structure (DSD parameters) and phase composition at each altitude, all macroscopic precipitation parameters can be calculated through explicit physical definitions.

[0079] Specifically, since the globally optimal state vector field contains the optimal raindrop spectral parameters and phase mixing coefficient This eliminates the need to rely on any empirical ZR relationship. Instead, by directly integrating the optimal raindrop spectrum, multiplying and integrating the particle volume, falling velocity, and number density distribution, a highly accurate and physically meaningful estimate of the rainfall rate can be obtained. This method not only distinguishes the contributions of liquid water and solid ice to radar echoes but also accurately calculates precipitation flux, especially in mixed-phase regions (such as the zero-degree bright band), where its calculation accuracy is significantly better than traditional methods.

[0080] For example, in this embodiment, height The rainfall rate at a given location can be calculated using the following formula: in, In order to be at a high altitude Instantaneous rainfall rate at a location (usually measured in mm / hr). The formula for the volume of a sphere Geometric coefficients in; To obtain the optimal raindrop spectrum, substitute the optimal parameters obtained from the inversion. The subsequent distribution function; This indicates that when calculating rainfall rate, only the contribution of liquid water is usually considered; therefore, this specifically refers to the falling velocity of liquid raindrops. The above formula allows for the precise derivation of macroscopic physical quantities from the inverted microscopic distribution, thus achieving a complete physical inversion from observational data to high-precision meteorological products.

[0081] Example 2 This embodiment is used to illustrate in detail how, based on the high-precision rainfall rate integral results obtained through the steps in Embodiment 1, a dynamic rainfall risk early warning system capable of capturing the trend of rainfall outbreaks and vertical evolution characteristics is constructed by further introducing differential operators of time and space dimensions.

[0082] Specifically, in this embodiment, the rainfall dynamic risk early warning system may include: 1) By analyzing rainfall rates over a continuous time series The rate of change of the rainfall rate is obtained by taking the time derivative. This rate of change is used as an indicator of the acceleration of precipitation bursts, thereby realizing the transformation from static rainfall monitoring and analysis to dynamic trend early warning. In meteorology, taking the time derivative of the rainfall rate is similar to the concept of acceleration; it focuses not on the current amount of rain, but on how quickly the rain is increasing in intensity.

[0083] Specifically, for urban drainage scheduling or flood warning, a rapidly intensifying moderate precipitation (with a very large positive derivative) often poses a higher potential risk than a weakening heavy precipitation (with a negative derivative).

[0084] Therefore, a 'burst acceleration' can be defined by calculating the first or even second derivative of the rainfall rate. In order to overcome the differential noise caused by the small fluctuations in the observed data, a difference fitting within a sliding time window can be used to robustly estimate this derivative.

[0085] For example, in this embodiment, the precipitation burst acceleration can be calculated using the following formula: in, The result obtained by integration in Example 2 Rainfall rate at any given time; The sampling time interval; The size of the sliding window; The weight is determined by time decay, with a greater weight closer to the current time.

[0086] It should be noted that the above formula is essentially a weighted differential operator, when Furthermore, when the value exceeds a certain threshold, the system can determine that a rainstorm is in its peak phase, even if the current absolute rainfall... Even if the warning line has not been reached, the system should still trigger an early warning and notify the drainage system to start the pumping station in advance.

[0087] 2) By analyzing the rainfall rate along the vertical profile By performing height differentiation, we can reveal whether precipitation particles undergo evaporation loss or collision-coalescence growth as they pass through the atmosphere, thereby quantifying the changes in mass flux of raindrops during their descent and enabling in-depth diagnosis of precipitation landing efficiency and intra-cloud processes.

[0088] Specifically, in arid or semi-arid regions, radar echoes at high altitudes are often very strong, but raindrops evaporate violently as they fall through dry air layers, resulting in very little actual ground rainfall. Conversely, under conditions of sufficient moisture, raindrops absorb cloud water as they fall, causing rainfall to increase as they descend.

[0089] Simple integration cannot distinguish between these two cases, but spatial differentiation can accurately capture this process. By calculating the gradient of rainfall rate with respect to altitude, an evolution factor can be defined to correct the surface precipitation estimate.

[0090] For example, in this embodiment, the vertical microphysical evolution factor can be calculated using the following formula: in, For height Vertical microphysical evolution factor; For height increments.

[0091] It should be noted that, if (That is, the rainfall in the upper layer is greater than that in the lower layer), indicating the presence of evaporation or wind shear causing particles to move out of the beam; if (That is, the rainfall in the upper layer is less than that in the lower layer), indicating the presence of collisional growth or cloud-water conversion processes. This can be monitored in real time. Based on the numerical distribution, the system can automatically determine the real and virtual nature of the radar echo. When the value is significantly positive (severe evaporation), the estimated value of surface precipitation is automatically lowered to avoid false alarms in weather forecasts.

[0092] Example 3 In this embodiment, a comprehensive precipitation dynamics risk index model is constructed based on the integral obtained in Embodiment 1 and the differential component obtained in Embodiment 2 to achieve the final intelligent decision-making.

[0093] Specifically, in this embodiment, the static rainfall intensity obtained through integration, the dynamic rainfall trend obtained through temporal differentiation, and the structural characteristics of rainfall obtained through spatial differentiation are integrated into a unified dimensionless risk index. This risk index can respond most strongly to high-risk scenarios.

[0094] For example, in this embodiment, the comprehensive precipitation dynamics risk index can be calculated using the following formula: in, For the ground layer ( Real-time rainfall rate; For reference heavy rainfall thresholds; Nonlinear amplification factor (usually) ); The sensitivity coefficient to the outbreak trend. This is to ensure that the risk is amplified only when the rainfall intensifies; The function is used to normalize the effect of the vertical gradient when When it is negative (due to collision and growth), this item is greater than 1, further amplifying the risk value.

[0095] It should be noted that this formula will integrate ( ) and differential ( A perfect blend, in which, This ensures that the baseline rainfall amount determines the risk level; this measure can assess the ground level at the current moment ( How heavy is the rain? (In the formula,) Used for normalization processing, for example, it can be set mm / h (rainstorm threshold); if the current rainfall is 25, this value is 0.5; if it is 100, it is 2.0. Nonlinear amplification factor. (for example, take) ) is a nonlinear amplifier; if it is a small rain (0.5), then it is (The overall risk is reduced and can be ignored); if it is heavy rain (2.0), (At this point, the risk is dramatically amplified). This feature allows the system to be highly sensitive to heavy rain, while remaining insensitive to light rain.

[0096] By utilizing the properties of the exponential function, when a positive precipitation acceleration is detected, the risk index increases exponentially, thus achieving a differential-based early warning effect. It is the derivative of the rainfall rate with respect to time (that is, the rate at which the rain intensifies). This indicates that we only care about the situation where the rain intensifies; if the rain intensifies (the derivative is negative), this term becomes 0, and it does not increase the risk; then through... Exponential function, as rainfall increases slowly ( (very small) The overall risk multiplier remains unchanged; however, if the rainfall increases dramatically ( (very large) It will become a huge multiplier (for example, 3 times or 5 times); through this formula, the effect of differential first can be achieved. Even if the rain is not heavy now (A is small), as long as the system detects that the rate of increase in rainfall is increasing, the total risk will also soar instantly, so that the warning can be issued a few minutes earlier than the actual rain peak.

[0097] This ensures that a risk is only identified when the rain actually falls to the ground (rather than evaporating in the air), giving the assessment system a dynamic insight into the precipitation process. It refers to the change in rainfall rate with altitude. When encountering collisions and intensification, the rainfall will be heavier as it goes down (less in the upper layers, more in the lower layers). It will be a negative number; there's a negative sign in the formula, and two negatives make a positive. If it becomes a positive number (e.g., +0.8), the value of that term becomes... The risk increases significantly at this point because the rain not only doesn't decrease but also absorbs moisture from the clouds, resulting in direct impact on the ground. Furthermore, when encountering evaporation / rain banners, the rain decreases as it descends (this is reflected in the spectrum as strong upper-level echoes and dry conditions below). It's a positive number; in the formula it becomes a negative number. If the number becomes negative (e.g., -0.8), the correction term will become... The risk is reduced because the rain evaporates in mid-air. This means that although the radar spectrum shows strong high-altitude echoes, there is actually little rainfall on the ground. This method helps to distinguish between true and false alarms, preventing the system from triggering false alarms due to high-altitude echoes in the radar spectrum, and ensuring that the system only reacts to precipitation that actually falls to the ground. Figure 6 A schematic diagram illustrating the early warning demonstration in this embodiment is shown. The diagram illustrates the temporal relationship between rainfall rate and burst acceleration. It can be seen from the diagram that the scheme in this embodiment can issue early warnings earlier than the traditional threshold method.

[0098] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above description is only a specific embodiment of the present invention and is not intended to limit the scope of protection of the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A rainfall analysis method based on micro-rain radar data, characterized in that, The rainfall analysis method includes: S100. Use micro-rain radar to vertically detect the target airspace and acquire radar echo signals; S200. Noise is filtered out from the original Doppler power spectral density matrix to obtain the cleaned observed Doppler power spectrum. S300. Use a pre-trained neural network model to extract features from the observed Doppler power spectrum and output a prior probability distribution of precipitation type. S400. Construct and solve the unified inversion model for precipitation state. The unified inversion model for precipitation status includes: A unified precipitation state vector is constructed, which includes parameters describing the microphysical structure of precipitation and mixing coefficients describing the phase weights of precipitation. A global inversion model based on variational inference is constructed. This model includes a forward physics operator configured based on a uniform precipitation state vector with assumptions. The theoretical Doppler power spectrum, considering the path integral attenuation effect, is then calculated. Construct a target functional, wherein the target functional includes at least: Data fidelity terms used to constrain the consistency between the theoretical spectrum and the observed spectrum, and structured prior terms used to constrain the approximation of the prior probability distribution by the inversion solution; The extreme values ​​of the objective functional are solved by numerical optimization algorithms to obtain the globally optimal unified precipitation state vector. S500: Based on a unified precipitation state vector, calculate the precipitation particle number density distribution at each altitude level. The instantaneous rainfall rate at each altitude level is obtained by integrating the number density distribution, particle volume, and particle falling velocity of the precipitation particles.

2. The rainfall analysis method based on micro-rain radar data according to claim 1, characterized in that, The rainfall analysis method also includes: S600. Construct a dynamic rainfall risk early warning process, the process including: The time derivative of the instantaneous rainfall rate on a continuous time series is calculated to obtain the burst acceleration index characterizing the rate of change of precipitation intensity; And / or, The spatial gradient of instantaneous rainfall rate at vertical height is calculated to obtain the microphysical evolution factor characterizing the vertical evolution of precipitation particles.

3. The rainfall analysis method based on micro-rain radar data according to claim 2, characterized in that, The rainfall analysis method also includes: S700, Calculate the comprehensive risk index. The comprehensive risk index is a dimensionless value constructed based on real-time ground rainfall rate, the burst acceleration index, and the microphysical evolution factor. The burst acceleration index acts on the comprehensive risk index in an exponential manner, amplifying the influence of the precipitation enhancement trend in a nonlinear way.

4. The rainfall analysis method based on micro-rain radar data according to claim 2, characterized in that, The burst acceleration index is calculated using a weighted difference method within a sliding time window, with the weight increasing as it approaches the current moment. When the burst acceleration index exceeds a preset positive threshold, the precipitation is determined to be in the burst enhancement stage and an early warning is triggered.

5. The rainfall analysis method based on micro-rain radar data according to claim 2, characterized in that, The microphysical evolution factors are used to correct precipitation forecasts, specifically including: When the microphysical evolution factor shows that the upper layer rainfall rate is significantly greater than the lower layer rainfall rate, it is determined that there is evaporation loss, and the estimated surface precipitation is lowered. When the microphysical evolution factor shows that the upper-layer rainfall rate is significantly lower than the lower-layer rainfall rate, it is determined that there is collisional growth, and the estimated surface precipitation is maintained or increased.

6. The rainfall analysis method based on micro-rain radar data according to claim 1, characterized in that, The radar echo signal in S100 is the original Doppler power spectral density matrix that dynamically changes over time by setting the detection period to the second level and dividing the vertical detection path into multiple height layer distance libraries. Each data point in the matrix represents the echo power density at a specific height and a specific Doppler velocity.

7. The rainfall analysis method based on micro-rain radar data according to claim 1, characterized in that, In S200, the path attenuation characteristics of the signal when penetrating the precipitation area are retained, and no path attenuation correction is performed. The noise filtering includes: Transform the original Doppler power spectral density matrix to the time-frequency domain; Based on the difference between the correlation characteristics of precipitation signals in the time-frequency domain and the randomness characteristics of noise, an adaptive threshold is constructed to filter the transformed coefficients. The denoising process does not include layer-by-layer attenuation compensation for the radar beam propagation path.

8. The rainfall analysis method based on micro-rain radar data according to claim 1, characterized in that, S300 includes the following sub-steps: S310. Extract the macroscopic feature profile of the preprocessed observation spectrum; S320. Input the macroscopic feature profile into a neural network whose output layer is configured as a probability output form; S330. Output a probability vector containing multiple dimensions, where the value of each dimension represents the confidence level of different precipitation phases, and the precipitation phases include at least stratiform cloud precipitation, convective precipitation, and solid precipitation.

9. The rainfall analysis method based on micro-rain radar data according to claim 1, characterized in that, The mixing coefficients in the unified precipitation state vector in S400 satisfy the normalization condition and are used to characterize the coexistence probability of different precipitation phases at the current altitude layer. When calculating the theoretical Doppler power spectrum, the forward physics operator uses the mixing coefficient to perform a weighted summation of the scattering cross section and falling velocity for different precipitation phases.

10. A rainfall analysis system based on micro-rain radar data, characterized in that, The rainfall analysis system includes: processor; The memory stores a computer program that, when executed by a processor, implements the rainfall analysis method based on micro-rain radar data as described in any one of claims 1 to 9.

Images