Raindrop size distribution estimation method based on phased array radar and weather unmanned aerial vehicle

Abstract

This invention relates to a raindrop size distribution estimation method based on phased array radar and meteorological UAVs, belonging to the field of radio. It employs a multi-resolution model design strategy, using different models for near- and far-range data to address the resolution differences between near and far targets, providing more accurate droplet size estimation. By constructing a disturbance intensity index and combining it with UAV flight status data, the influence of UAV flight disturbances on droplet size distribution is dynamically corrected, improving the accuracy of droplet size estimation under conditions of high disturbance intensity. Through model fusion and fluctuation loss function optimization, the dynamic changes in precipitation droplet size are effectively captured, enhancing the robustness and accuracy of the estimation model. This invention can estimate droplet size distribution in real time, adapting to heavy precipitation and rapidly changing precipitation systems, exhibiting high real-time response capability and adaptability.

Description

Technical Field

[0001] This invention relates to the field of radio, and more particularly to a method for estimating raindrop size distribution based on phased array radar and meteorological drones. Background Technology

[0002] The estimation of raindrop number density across different droplet size ranges, known as raindrop size distribution (DSD), is a crucial parameter in meteorology, widely applied in precipitation intensity estimation, radar meteorology, and climate prediction. Accurate estimation of raindrop size distribution is of significant value for weather forecasting, flood prediction, and climate change research. However, existing methods for estimating droplet size distribution suffer from insufficient accuracy, high computational complexity, and poor adaptability to dynamic precipitation environments.

[0003] Currently, methods for estimating raindrop size distribution mainly include: 1. Direct observation by ground-based raindrop spectrometers. This method can directly obtain the size distribution and usually has high temporal resolution in actual measurements. However, limited by the density of instrument deployment, it cannot provide dynamic, global precipitation size distribution information like radar, unlike radar's wide-area coverage. 2. The ZR relationship-based method, the most traditional method for size estimation based on radar reflectivity factor, uses the empirical relationship between radar reflectivity factor Z and precipitation intensity R for estimation. This method is simple and computationally efficient, but it has drawbacks such as ignoring the diversity of precipitation size and relatively low accuracy. 3. Methods combining reflectivity factor and differential reflectivity factor. This method improves the estimation of precipitation size distribution and can provide better differentiation between different precipitation phases (such as raindrops and hail). However, this method usually uses fixed mathematical models, which may not be applicable under different weather and meteorological conditions. Therefore, the applicability of this method is limited and cannot flexibly cope with various precipitation situations. 4. Physical model method: This method estimates particle size distribution by simulating the physical processes involved in precipitation particle size distribution (such as raindrop formation, collision, and evaporation). This method typically has high accuracy, especially under specific precipitation scenarios. However, it suffers from high computational complexity, making rapid particle size distribution estimation difficult. Furthermore, this method often relies on fixed physical assumptions, making it challenging to handle diverse precipitation environments and non-uniform precipitation conditions. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of the prior art and provide a method for estimating raindrop size distribution based on phased array radar and meteorological drones, thus solving the deficiencies of the prior art.

[0005] The objective of this invention is achieved through the following technical solution: a raindrop size distribution estimation method based on phased array radar and meteorological UAV, wherein the estimation method includes:

[0006] S1. Collect phased array radar echo data, radar operation status data, meteorological UAV detection data and flight status data. Perform time alignment based on the collected phased array radar echo data and meteorological UAV detection data. Perform spatial alignment based on the radar operation status data and meteorological UAV detection data by converting spatial geographic coordinates to radar polar coordinates model. Obtain raindrop particle size distribution sequence.

[0007] S2. Based on the raindrop size distribution sequence obtained in S1 and the collected flight state data, construct a set of steady-state time segments, and construct a set of disturbance intensity indices based on the collected phased array radar echo data, meteorological UAV detection data and flight state data.

[0008] S3. Based on the raindrop size distribution sequence, steady-state time segment set, and disturbance intensity index set obtained in S1, construct the fluctuation loss function, and optimize the disturbance intensity index set and the raindrop size distribution sequence after disturbance correction.

[0009] S4. Divide the phased array radar echo data into near and far ranges, construct an overall model through a parallel deep convolutional neural network model and a simple convolutional neural network model, and input the phased array radar echo data as the input dataset into the overall model to obtain the output dataset.

[0010] S5. Optimize the loss function of the overall model and train the overall model. The trained overall model is used to obtain the raindrop size distribution estimate based on phased array radar and meteorological UAV.

[0011] The phased array radar echo data includes: radar echo data Data0_par detected at polar coordinate spatial positions (r0_par, phi0_par, seat0_par) under the phased array radar observation time series t0_par, which includes reflectivity factor data Z0_par, differential reflectivity factor data Zdr0_par, and differential propagation phase shift rate KDP0_par, where r0_par is the radial distance, phi0_par is the azimuth, and seat0_par is the antenna elevation angle;

[0012] The radar operating status data includes: the spatial location information of the radar station as P_par = (Lon_par, Lat_par, H_par), where Lon_par, Lat_par, and H_par are the longitude, latitude, and altitude of the antenna feed, respectively; the radial range resolution delt_R_par and beamwidth delt_sita_par of the radar; and the radar volume scanning period delt_T_par.

[0013] The meteorological drone detection data includes: raindrop particle size distribution sequence obtained by the particle imaging probe carried by the drone under the meteorological drone observation time series t0_uav, that is, raindrop number density in different particle size intervals, and the spatial location of the raindrop P_uav=(Lon_uav, Lat_uav, H_uav), where Lon_uav, Lat_uav and H_uav are the longitude, latitude and altitude of the observation data, respectively;

[0014] The flight status data includes: angular velocity vector w_uav=[wx_uav, wy_uav, wz_uav], velocity vector v_uav=[vx_uav, vy_uav, vz_uav], and propeller rotational speed RPM_uav provided by the UAV's onboard inertial navigation and flight control system, where wx_uav, wy_uav, and wz_uav are the components of w_uav in the x, y, and z axes, respectively, and vx_uav, vy_uav, and vz_uav are the components of v_uav in the x, y, and z axes, respectively, with a sampling period of delta_T_uav.

[0015] The time alignment based on collected phased array radar echo data and meteorological UAV detection data, and the spatial alignment based on radar operational status data and meteorological UAV detection data using a spatial geographic coordinate to radar polar coordinate model, include:

[0016] A1. Based on the collected phased array radar observation time series t0_par and meteorological UAV observation time series t0_uav, for a certain moment t0_uav_1 in the meteorological UAV observation time series t0_uav, calculate the time difference between that moment t0_uav_1 and each moment in the phased array radar observation time series t0_par, and form a time difference sequence delta_t.

[0017] A2. Set the maximum permissible error delta_t_max for the alignment of observation time between the phased array radar and the UAV, and find the phased array radar observation time corresponding to the minimum time difference of the elements in delta_t that satisfies the condition of being less than or equal to delta_t_max.

[0018] A3. Repeat steps A1 and A2 to traverse all observation times in the meteorological UAV observation time series t0_uav, and obtain the phased array radar observation time series after time alignment with all observation times in the meteorological UAV observation time series t0_uav, denoted as t1_par. Based on the collected phased array radar echo data of the phased array radar observation time series t0_par, extract the phased array radar echo data Data1_par of the phased array radar observation time series t1_par, including reflectivity factor data Z1_par, differential reflectivity factor data Zdr1_par and differential propagation phase shift rate KDP1_par;

[0019] A4. Input the collected P_par, P_uav, delt_R_par, and delt_sita_par into the spatial geographic coordinate to radar polar coordinate conversion model, run the model, and obtain the raindrop size distribution sequence of the meteorological UAV under the meteorological UAV observation time series t0_uav after conversion to radar polar coordinates, denoted as ND_uav=[N_uav1, N_uav2, …, N_uav B ], where B is the total number of particle size partitions, N_uav b Represents the b-th particle size interval [D b-1 D b The raindrop number density within the range, b=1,2,…,B, and the corresponding spatial position information in radar polar coordinates, are denoted as P0_uav=(r0_uav,phi0_uav,sita0_uav), where r0_uav is the radial distance, phi0_uav is the azimuth, and sata0_uav is the antenna elevation angle.

[0020] A5. Based on the P0_uav information and combined with the spatial location information of Data1_par, extract the phased array radar echo data after spatial alignment with P0_uav according to the principle of shortest spatial distance, and denot it as Data2_par, which includes reflectivity factor data Z2_par, differential reflectivity factor data Zdr2_par and differential propagation phase shift rate KDP2_par.

[0021] The construction of the steady-state time segment set based on the raindrop size distribution sequence obtained from S1 and the collected flight state data includes:

[0022] According to the raindrop size distribution sequence ND_uav obtained by the meteorological UAV at the meteorological UAV observation time series t0_uav, and the sampling period delta_T_uav, set the sliding window length as Wind_move = L×delta_T_uav. Take the start time of the time series t0_uav as the starting position of the sliding window, and perform sliding processing on the time series t0_uav in sequence according to the window width Wind_move. Finally, decompose the time series t0_uav into K sliding time windows, denoted as Wk, k = 1, 2, …, K;

[0023] Extract the raindrop size distribution sequences corresponding to the K sliding time windows from the raindrop size distribution sequence ND_uav respectively, and perform logarithmic transformation, denoted as log_N_uav_k. Then calculate the standard deviation of the log_N_uav_k sequence within each sliding time window, denoted as STDk;

[0024] Set the steady-state threshold as STD0, and the condition for steady-state segment screening is: STDk < STD0. Select the windows that meet the condition from the K sliding time windows, and collect the observation times corresponding to these windows that meet the screening condition to obtain the steady-state segment time set, denoted as STEADY.

[0025] The construction of the disturbance intensity index set according to the collected phased array radar echo data, meteorological UAV detection data, and flight state data includes:

[0026] B1. According to the meteorological UAV at the meteorological UAV observation time series t0_uav, flight state data w_uav, v_uav, and RPM_uav, for a certain moment t0_uav_i in the meteorological UAV observation time series t0_uav, first calculate the angular velocity modulus ||w_uav|| and velocity modulus ||v_uav|| at this moment according to the Euclidean norm calculation formula. Then find the maximum value RPM_uav_max in the propeller speed RPM_uav, and calculate the propeller speed normalization term m_PRM_uav = RPM_uav / RPM_uav_max;

[0027] B2. Construct atmospheric disturbance feature terms u1_uav=||w_uav|| / w_ref for angular velocity, u2_uav=||v_uav|| / v_ref for velocity, and u3_uav=m_PRM_uav for propeller speed. Then, construct the disturbance intensity index dis_i=a1×u1_uav+a2×u2_uav+a3×u3_uav for particles in the sampling area at time t0_uav_i, where w_ref and v_ref are the calibration values ​​of the angular velocity and velocity of the UAV system, respectively, and a1, a2 and a3 are the initial disturbance weight coefficients of feature terms u1_uav, u2_uav and u3_uav, respectively.

[0028] B3. Update the value of t0_uav_i at a certain time step to the value of each time step in the time series t0_uav. Repeat steps B1 and B2 to generate the disturbance intensity index for each time step in the time series t0_uav. Set these indices together to obtain the disturbance intensity index set, denoted as DIS.

[0029] The construction of the fluctuation loss function based on the raindrop size distribution sequence, steady-state time segment set, and disturbance intensity index set obtained from S1 includes:

[0030] C1. Based on the raindrop size distribution sequence ND_uav obtained by the meteorological drone under the observation time series t0_uav, combined with the steady-state time segment set STEADY, firstly extract any observation time t1 and its adjacent time t2 from the time set STEADY, where t2=t1+delta_T_uav, and extract the corresponding raindrop size distribution sequences N_uav_t1 and N_uav_t2.

[0031] C2. Calculate the logarithmic particle size distributions of N_uav_t1 and N_uav_t2, denoted as log_N_uav_t1=log(N_uav_t1+ext_pos) and log_N_uav_t2=log(N_uav_t2+ext_pos), respectively. Then calculate the difference between the two, denoted as delta_log_N_uav= log_N_uav_t1- log_N_uav_t2, where ext_pos is a very small positive value.

[0032] C3. Based on the disturbance intensity index set DIS, extract the disturbance intensity index corresponding to time t1, denoted as dis_t1, and construct the disturbance weight function R_t1=1+alpha0_dis×dis_t1, where alpha0_dis is the disturbance amplification coefficient constant. Then, construct the adjacent time fluctuation loss L_t1= R_t1×sum(delta_log_N_uav) for time t1. 2 ), where sum() is the summation across all intervals of the particle size partition;

[0033] C4. Update any observation time t1 sequentially to all observation times in the steady-state time segment set STEADY and repeat steps C1 to C3 to generate the fluctuation loss of all observation times in STEADY. Sum these fluctuation losses to generate the fluctuation loss sum in the steady-state segment, denoted as LOSS.

[0034] The optimized set of disturbance intensity indices and the raindrop size distribution sequence corrected for disturbance effects include:

[0035] Let the discretization step size delta_a = 1 / N, where N is the discretization parameter. By enumerating non-negative integer triples (i,j,k) satisfying i+j+k=N, we construct perturbation weight coefficients a1=i×delta_a, a2=j×delta_a, a3=k×delta_a, generating a set of G candidate weight coefficients Ag=[a1,a2,a3]. g g = 1, 2, ..., G, where G is the total number of groups;

[0036] Each set of weights in the candidate weight coefficient set Ag is sequentially input into the fluctuation loss and LOSS within the steady-state segment. The fluctuation loss sum of all observation times within STEADY is calculated for each set of weights. The fluctuation loss sums of G different weight coefficient sets are aggregated and denoted as LOSS_g.

[0037] Find the minimum value in the LOSS_g set, and use the combination of weight coefficients corresponding to the minimum value as the optimal weight parameters for calculating the perturbation intensity index, denoted as [a1_opt, a2_opt, a3_opt]. Replace the initial [a1, a2, and a3] with [a1_opt, a2_opt, a3_opt], ​​and re-execute steps B1-B3 to generate an optimized perturbation intensity index set, denoted as DIS_opt;

[0038] Construct a function to determine the effect of perturbation intensity on particle size distribution: F_DIS = 1 + apha1_dis × DIS_opt, where apha1_dis is a correction coefficient, based on [D b-1 D b Constructing particle size D bThe interaction function with the perturbation G_DIS_D b =1 / (1+apha2_dis×D b ), b=1,2,…,B, where apha2_dis are the relation coefficients, and the B interaction functions are set together and denoted as G_DIS;

[0039] Based on the obtained ND_uav, combined with the already obtained F_DIS and G_DIS, the raindrop size distribution sequence after perturbation correction is generated as ND_uav_c = ND_uav × F_DIS × G_DIS.

[0040] S4 includes:

[0041] Set the radial distance distinction threshold to Rang_set. Based on the radial distance parameter in the spatial information of the phased array radar echo data, the echo data that is less than or equal to Rang_set is divided into the near-range echo area, and the radar reflectivity factor data of the corresponding area is denoted as Z2_par_near. The echo data that is greater than Rang_set is divided into the far-range echo area, and the radar reflectivity factor data of the corresponding area is denoted as Z2_par_far.

[0042] For radar echo data in the near range, a deep convolutional neural network model is used, denoted as MODEL1. For radar echo data in the far range, a simple convolutional neural network model is used, denoted as MODEL2. The loss function for both models is the mean square error function, denoted as MSE_MODEL1 and MSE_MODEL2.

[0043] MODEL1 and MODEL2 are connected in parallel to form a whole model MODEL. For the near echo region, the weight is set as weight_near=mean(Z2_par_near) / max(Z2_par_near). For the far echo region, the weight is set as weight_far=min(Z2_par_far) / mean(Z2_par_far). Then, the loss function of the whole model MODEL is set as loss_MODEL=weight_near×MSE_MODEL1+ weight_far×MSE_MODEL2.

[0044] Based on the obtained phased array radar echo data, including reflectivity factor data Z2_par, differential reflectivity factor data Zdr2_par, and differential propagation phase shift rate KDP2_par, the ratio of reflectivity factor to differential reflectivity factor Zdr_Z=Zdr2_par / Z2_par is constructed, and the combination of differential propagation phase shift rate and reflectivity factor KDP_Z=KDP2_par×Z2_par is constructed.

[0045] Z2_par, Zdr2_par, KDP2_par together with Zdr_Z and KDP_Z constitute the input dataset of the model MODEL, while the raindrop size distribution sequence ND_uav_c, which is corrected for the influence of perturbation, is used as the output dataset of the model MODEL.

[0046] S5 includes:

[0047] Based on the corrected raindrop size distribution sequence ND_uav_c and the optimized perturbation intensity index set DIS_opt, for each particle diameter interval [D] in ND_uav_c... b-1 D b Construct the corresponding weights, denoted as weight_D. b =beta1×DIS_opt+ beta2×[D b / D max ], where beta1 and beta2 are both dynamically adjusted parameters, D max This represents the maximum particle diameter in ND_uav_c;

[0048] Based on the obtained loss_MODEL, combined with weight_D b The optimized overall model loss function is loss_MODEL_opt=sum(weight_D b ×loss_MODEL), where sum() is the summation over all particle diameter distribution intervals;

[0049] During the optimization process of beta1 and beta2, the Bayesian optimization method was selected to adjust the parameters and train the model, and finally the optimal parameters of beta1 and beta2 were obtained.

[0050] After the model training is completed, Z2_par, Zdr2_par, KDP2_par and Zdr_Z, KDP_Z are input into the trained model again for calculation, thereby realizing the estimation of raindrop number density in different particle size ranges, that is, completing the estimation of raindrop particle size distribution based on phased array radar and meteorological UAV.

[0051] This invention has the following advantages: The raindrop size distribution estimation method based on phased array radar and meteorological UAVs employs a multi-resolution model design strategy, which solves the problem of resolution differences between near and far targets by using different models for near and far distance data, providing more accurate size estimation. By constructing a disturbance intensity index and combining it with UAV flight status data, the influence of UAV flight disturbances on the size distribution is dynamically corrected, improving the size estimation accuracy under conditions of high disturbance intensity. Through model fusion and fluctuation loss function optimization, the dynamic changes in precipitation size are effectively captured, improving the robustness and accuracy of the estimation model. It can estimate the size distribution in real time, adapting to heavy precipitation and rapidly changing precipitation systems, exhibiting high real-time response capability and adaptability. Attached Figure Description

[0052] Figure 1 This is a schematic diagram of the process of the present invention;

[0053] Figure 2 This is a flowchart illustrating the implementation of the present invention for generating phased array radar echo data after spatiotemporal alignment with meteorological drone detection data;

[0054] Figure 3 This is a flowchart illustrating the implementation of the steady-state time segment set generation method of the present invention.

[0055] Figure 4 This is a flowchart illustrating the implementation process of generating the disturbance intensity index set of the present invention.

[0056] Figure 5 This is a flowchart illustrating the implementation of the steady-state segment fluctuation loss and calculation of the present invention;

[0057] Figure 6 A flowchart illustrating the implementation of the raindrop size distribution sequence after perturbation correction according to the present invention;

[0058] Figure 7 This is a flowchart illustrating the overall model training set construction and loss function setting implementation of the present invention.

[0059] Figure 8 This is a flowchart illustrating the implementation of the raindrop size distribution estimation based on phased array radar and meteorological drones according to the present invention. Detailed Implementation

[0060] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. The components of the embodiments of this application described and shown in the accompanying drawings can generally be arranged and designed in various different configurations. Therefore, the detailed description of the embodiments of this application provided below with reference to the accompanying drawings is not intended to limit the scope of protection of the claimed application, but merely represents selected embodiments of this application. All other embodiments obtained by those skilled in the art based on the embodiments of this application without inventive effort are within the scope of protection of this application. The present invention will be further described below with reference to the accompanying drawings.

[0061] like Figure 1 As shown, this invention specifically relates to a method for estimating raindrop size distribution based on phased array radar and a weather drone, which specifically includes the following:

[0062] Step (1): Collect phased array radar echo data and radar operation status data; collect meteorological UAV detection data and flight status data;

[0063] Among them, the phased array radar echo data refers to the radar echo data Data0_par detected at the polar coordinate spatial position (r0_par, phi0_par, seat0_par) under the phased array radar observation time series t0_par, including reflectivity factor data Z0_par, differential reflectivity factor data Zdr0_par, and differential propagation phase shift rate KDP0_par, where r0 is the radial distance, phi0 is the azimuth, and seat0 is the antenna elevation angle;

[0064] The radar operation status data mainly includes: the spatial location information of the radar site as P_par = (Lon_par, Lat_par, H_par), where Lon_par, Lat_par, and H_par are the longitude, latitude, and altitude of the antenna feed, respectively; the radial range resolution of the radar as delt_R_par = 30 meters and the beamwidth as delt_sita_par = 1.0 degree; and the radar volume scan period as delt_T_par = 30 seconds.

[0065] Meteorological drone detection data refers to the raindrop size distribution sequence obtained by the particle imaging probe carried by the drone under the observation time series t0_uav of the meteorological drone, that is, the raindrop number density in different particle size intervals, and the spatial location of the raindrop P_uav=(Lon_uav, Lat_uav, H_uav), where Lon_uav, Lat_uav and H_uav are the longitude, latitude and altitude of the observation data, respectively;

[0066] The flight status data of meteorological UAVs refers to the angular velocity vector w_uav=[wx_uav, wy_uav, wz_uav], velocity vector v_uav=[vx_uav, vy_uav, vz_uav], and propeller speed RPM_uav provided by the UAV's onboard inertial navigation and flight control system. Here, wx_uav, wy_uav, and wz_uav are the components of w_uav in the x, y, and z axes, respectively, and vx_uav, vy_uav, and vz_uav are the components of v_uav in the x, y, and z axes, respectively. The sampling period is delta_T_uav=5 seconds.

[0067] Step (2): Time alignment and spatial alignment;

[0068] like Figure 2 As shown, firstly, based on the phased array radar observation time series t0_par and the meteorological UAV observation time series t0_uav collected in step (1), for a certain moment t0_uav_1 in the meteorological UAV detection time series t0_uav, the time difference between this moment t0_uav_1 and each moment in the phased array radar observation time series t0_par is calculated to form a time difference sequence delta_t; then, the maximum permissible error for the alignment of the phased array radar and UAV observation time is defined as delta_t_max = delt_T_par, and the phased array radar observation moment corresponding to the minimum time difference of the elements in delta_t that satisfies the condition of being less than or equal to delta_t_max is found. By iterating through all observation times in the meteorological UAV observation time series t0_uav from the above steps, the phased array radar observation time series aligned with all observation times in the meteorological UAV observation time series t0_uav can be obtained, denoted as t1_par; then, based on the phased array radar echo data of the phased array radar observation time series t0_par collected in step (1), the phased array radar echo data Data1_par of the phased array radar observation time series t1_par is extracted, including reflectivity factor data Z1_par, differential reflectivity factor data Zdr1_par and differential propagation phase shift rate KDP1_par; then, P_par, P_uav, delt_R_par, and delt_sita_par collected in step (1) are input together into the spatial geographic coordinate to radar polar coordinate model, the model is run, and the model output is converted to radar polar coordinates and the raindrop size distribution sequence of the meteorological UAV under the meteorological UAV observation time series t0_uav is denoted as ND_uav=[N_uav1, N_uav2, …, N_uav B ], where B is the total number of particle size partitions, N_uav bIndicates the number density of raindrops in the b-th particle size interval [D b-1 , D b , where b = 1, 2, …, B, and the corresponding spatial position information in the radar polar coordinates, denoted as P0_uav = (r0_uav, phi0_uav, sita0_uav), where r0_uav is the radial distance, phi0_uav is the azimuth, and sita0_uav is the antenna elevation angle. Then, according to the P0_uav information, combined with the spatial position information of Data1_par that has been obtained, the phased array radar echo data aligned with P0_uav in space is extracted according to the principle of the shortest spatial distance, denoted as Data2_par, including the reflectivity factor data Z2_par, the differential reflectivity factor data Zdr2_par, and the differential propagation phase shift rate KDP2_par.

[0069] Step (3): Construction of the steady-state segment time set;

[0070] As Figure 3 shown, according to the raindrop size distribution sequence ND_uav obtained by the meteorological UAV at the meteorological UAV observation time series t0_uav output in step (2), and the sampling period delta_T_uav collected in step (1), first set the sliding window length as Wind_move = L × delta_T_uav, L = 3, take the start moment of the time series t0_uav as the starting position of the sliding window, and perform sliding processing on the time series t0_uav in sequence according to the window width Wind_move. Finally, the time series t0_uav is decomposed into K sliding time windows, denoted as Wk, k = 1, 2, …, K; secondly, extract the raindrop size distribution sequences corresponding to the K sliding time windows from the raindrop size distribution sequence ND_uav respectively, and perform logarithmic transformation, denoted as log_N_uav_k. Then, calculate the standard deviation of the log_N_uav_k sequence within each sliding time window, denoted as STDk; after that, set the steady-state threshold as STD0 = 0.1, and set the steady-state segment screening condition as: STDk < STD0. Select the windows that meet the conditions from the K sliding time windows, and collect the observation moments corresponding to these windows that meet the screening conditions to obtain the steady-state segment time set, denoted as STEADY.

[0071] Step (4): Construction of the disturbance intensity index set;

[0072] As Figure 4As shown, based on the meteorological UAV observation time series t0_uav, flight status data w_uav, v_uav and RPM_uav collected in step (1), for a certain moment t0_uav_i in the meteorological UAV observation time series t0_uav, firstly, the angular velocity modulus ||w_uav|| and velocity modulus ||v_uav|| at that moment are calculated according to the Euclidean norm calculation formula. Then, the maximum value RPM_uav_max in the propeller speed RPM_uav is found, and the propeller speed normalization term m_PRM_uav=RPM_uav / RPM_uav_max is calculated. Secondly, the angular velocity atmospheric disturbance characteristic term u1 is constructed. _uav=||w_uav|| / w_ref, velocity atmospheric disturbance feature term u2_uav=||v_uav|| / v_ref, propeller speed atmospheric disturbance feature term u3_uav=m_PRM_uav, where w_ref and v_ref are the calibration values ​​of the UAV system's angular velocity and velocity, respectively. This constructs the disturbance intensity index dis_i=a1×u1_uav+a2×u2_uav+a3×u3_uav for particles in the sampling region at time t0_uav_i, where a1, a2, and a3 are the initial disturbance weight coefficients for feature terms u1_uav, u2_uav, and u3_uav, respectively. By sequentially updating time t0_uav_i in this processing step to the value of each time in the time series t0_uav, the disturbance intensity index for each time in the time series t0_uav can be generated. These indices are then set together to obtain the disturbance intensity index set, denoted as DIS.

[0073] The Euclidean norm is calculated as follows: ||w_uav||=sqrt(wx_uav) 2 +wy_uav 2 +wz_uav 2 ), ||v_uav||=sqrt(vx_uav) 2 +vy_uav 2 +vz_uav 2 ), where sqrt() is the square root function.

[0074] Step (5): Constructing the volatility loss function;

[0075] like Figure 5As shown, based on the raindrop size distribution sequence ND_uav obtained by the meteorological UAV under the observation time series t0_uav of the meteorological UAV output in step (2), combined with the steady-state time segment set STEADY obtained in step (3), firstly, any observation time t1 and its adjacent time t2 are extracted from the time set STEADY, where t2=t1+delta_T_uav, and the corresponding raindrop size distribution sequences N_uav_t1 and N_uav_t2 are extracted; then, the logarithmic size distributions of N_uav_t1 and N_uav_t2 are calculated, denoted as log_N_uav_t1=log(N_uav_t1+ext_pos) and log_N_uav_t2=log(N_uav_t2+ext_pos), respectively, where ext_pos is a very small positive value, and then the difference between the two is calculated, denoted as delta_log_N_uav= log_N_uav_t1- log_N_uav_t2; then, according to step (4), the disturbance intensity index set DIS is obtained, and the disturbance intensity index corresponding to time t1 is extracted and denoted as dis_t1. Thus, the disturbance weight function R_t1=1+alpha0_dis×dis_t1 is constructed, where alpha0_dis=0.2 is the disturbance amplification coefficient constant. Then, the adjacent time fluctuation loss L_t1= R_t1×sum(delta_log_N_uav_t2) is constructed. 2 ), where sum() is the summation across all intervals of the particle size partition. Finally, any observation time t1 in this step is sequentially updated to all observation times in the steady-state time segment set STEADY, thereby generating the fluctuation loss for all observation times within STEADY, and these fluctuation losses are summed to generate the fluctuation loss sum within the steady-state segment, denoted as LOSS.

[0076] Step (6): Optimize the set of disturbance intensity indices and the raindrop size distribution sequence after being corrected for the influence of disturbance;

[0077] like Figure 6 As shown, firstly, let the discretization step size delta_a = 1 / N, where N is the discretization parameter. Then, by enumerating the non-negative integer triples (i,j,k) that satisfy i+j+k=N, and constructing the perturbation weight coefficients set in step (4) accordingly: a1=i×delta_a, a2=j×delta_a, a3=k×delta_a, generating a set of G candidate weight coefficients Ag=[a1,a2,a3] g, g=1,2,…,G, where G is the total number of groups; next, input each weight in the candidate weight coefficient set Ag into the fluctuation loss and LOSS in the steady-state segment obtained in step (5), and calculate the fluctuation loss of all observation times in STEADY for each weight, and set the fluctuation loss of different weight coefficient sets in G, denoted as LOSS_g; then, find the minimum value in the LOSS_g set, and use the weight coefficient combination corresponding to the minimum value as the optimal weight parameter for calculating the disturbance intensity index, denoted as [a1_opt,a2_opt,a3_opt]; then, replace the [a1,a2 anda3] initially set in step (4) with [a1_opt,a2_opt,a3_opt], ​​and re-execute step (4) to generate the optimized disturbance intensity index set, denoted as DIS_opt; then, construct the influence function of disturbance intensity on particle size distribution F_DIS=1+ apha1_dis×DIS_opt, where apha1_dis is the correction coefficient; then, based on [D] obtained in step (2) b-1 D b ], constructing particle size D b Interaction function G_DIS_D with the perturbation b =1 / (1+apha2_dis×D b ), b=1,2,…,B, where apha2_dis is the relation coefficient, and the B interaction functions are set together and denoted as G_DIS; then, based on ND_uav obtained in step (2), combined with the obtained F_DIS and G_DIS, the raindrop size distribution sequence after perturbation correction is generated ND_uav_c=ND_uav×F_DIS×G_DIS.

[0078] Step (7): Connect the two sub-models in parallel and divide them into near and far distances;

[0079] like Figure 7As shown, based on the phased array radar reflectivity factor data Z2_par obtained in step (2), the radial distance distinction threshold is first set to Rang_set = 10 km. According to the radial distance parameter in the spatial information of the phased array radar echo data, echo data less than or equal to Rang_set are divided into near-range echo areas, and the radar reflectivity factor data of the corresponding area is denoted as Z2_par_near. Echo data greater than Rang_set are divided into far-range echo areas, and the radar reflectivity factor data of the corresponding area is denoted as Z2_par_far. Then, for the radar echo data in the near-range area, a deep convolutional neural network model is used, denoted as MODEL1. For the radar echo data in the far-range area, a simple convolutional neural network model is used. Let's denote them as MODEL2. The loss functions for both models are the mean squared error function, denoted as MSE_MODEL1 and MSE_MODEL2. Then, MODEL1 and MODEL2 are concatenated in parallel to form the overall model MODEL. Next, for the near-field echo region, the weights are set as weight_near = mean(Z2_par_near) / max(Z2_par_near), and for the far-field echo region, the weights are set as weight_far = min(Z2_par_far) / mean(Z2_par_far). Finally, the loss function for the overall model MODEL is set as loss_MODEL = weight_near × MSE_MODEL1 + ... weight_far×MSE_MODEL2; then, based on the phased array radar echo data obtained in step (2), including reflectivity factor data Z2_par, differential reflectivity factor data Zdr2_par, and differential propagation phase shift rate KDP2_par, construct the ratio of reflectivity factor to differential reflectivity factor Zdr_Z=Zdr2_par / Z2_par, and construct the combination of differential propagation phase shift rate and reflectivity factor KDP_Z=KDP2_par×Z2_par; finally, Z2_par, Zdr2_par, KDP2_par together with Zdr_Z and KDP_Z constitute the input dataset of the overall model MODEL, and at the same time, the raindrop size distribution sequence ND_uav_c after disturbance correction obtained in step (6) is used as the output dataset of the model MODEL;

[0080] Step (8): Optimize the model loss function;

[0081] like Figure 8 As shown, based on the corrected raindrop size distribution sequence ND_uav_c obtained in step (6) and the optimized perturbation intensity index set DIS_opt, for each particle diameter interval [D] in ND_uav_c... b-1 D bConstruct the corresponding weights, denoted as weight_D. b =beta1×DIS_opt+ beta2×[D b / D max ], where beta1 and beta2 are dynamically adjusted parameters, both initially set to 0.1, D max The maximum value of the particle diameter in ND_uav_c; then, based on the loss_MODEL obtained in step (7), combined with weight_D b The optimized overall model loss function is loss_MODEL_opt=sum(weight_D b ×loss_MODEL), where sum() is the summation across all particle diameter distribution intervals. Then, during the optimization of beta1 and beta2, the Bayesian optimization method is selected to adjust the parameters and train the model, finally obtaining the optimal beta1 and beta2 parameters; finally, after the model training is completed, Z2_par, Zdr2_par, KDP2_par and Zdr_Z, KDP_Z output in step (7) are input into the trained model for calculation, thereby realizing the estimation of raindrop number density in different particle size intervals, that is, completing the estimation of raindrop particle size distribution based on phased array radar and meteorological UAV.

[0082] The above description is merely a preferred embodiment of the present invention. It should be understood that the present invention is not limited to the forms disclosed herein and should not be construed as excluding other embodiments. It can be used in various other combinations, modifications, and improvements, and can be altered within the scope of the concept described herein through the above teachings or related technologies or knowledge. Modifications and variations made by those skilled in the art that do not depart from the spirit and scope of the present invention should be within the protection scope of the appended claims.

Claims

1. A method for estimating raindrop size distribution based on phased array radar and meteorological UAV, characterized in that: The estimation method includes: S1. Collect phased array radar echo data, radar operation status data, meteorological UAV detection data and flight status data. Perform time alignment based on the collected phased array radar echo data and meteorological UAV detection data. Perform spatial alignment based on the radar operation status data and meteorological UAV detection data by converting spatial geographic coordinates to radar polar coordinates model. S2. Based on the data obtained in S1 and the collected radar operating status data, construct a set of steady-state time segments, and based on the collected phased array radar echo data, meteorological UAV detection data and flight status data, construct a set of disturbance intensity indicators. S3. Based on the data obtained in S1, the set of steady-state time segments, and the set of disturbance intensity indices, construct the fluctuation loss function, and optimize the set of disturbance intensity indices and the raindrop size distribution sequence after being corrected for the disturbance. S4. Divide the phased array radar echo data into near and far ranges, construct an overall model through a parallel deep convolutional neural network model and a simple convolutional neural network model, and input the data obtained in S1 into the overall model to obtain the output dataset. S5. Optimize the loss function of the overall model and train the overall model. The raindrop size distribution estimate based on phased array radar and meteorological UAV is obtained through the trained overall model. The phased array radar echo data includes: radar echo data Data0_par detected at polar coordinate spatial positions (r0_par, phi0_par, seat0_par) under the phased array radar observation time sequence t0_par, which includes reflectivity factor data Z0_par, differential reflectivity factor data Zdr0_par, and differential propagation phase shift rate KDP0_par, where r0_par is the radial distance, phi0_par is the azimuth, and seat0_par is the antenna elevation angle; The radar operating status data includes: the spatial location information of the radar station as P_par = (Lon_par, Lat_par, H_par), where Lon_par, Lat_par, and H_par are the longitude, latitude, and altitude of the antenna feed, respectively; the radial range resolution delt_R_par and beamwidth delt_sita_par of the radar; and the radar volume scanning period delt_T_par. The meteorological drone detection data includes: raindrop size distribution sequence obtained by the particle imaging probe carried by the drone under the meteorological drone observation time sequence t0_uav, that is, raindrop number density in different size intervals, and the spatial location of the raindrop P_uav=(Lon_uav, Lat_uav, H_uav), where Lon_uav, Lat_uav and H_uav are the longitude, latitude and altitude of the observation data, respectively; The flight status data includes: angular velocity vector w_uav=[wx_uav, wy_uav, wz_uav], velocity vector v_uav=[vx_uav, vy_uav, vz_uav], and propeller rotational speed RPM_uav provided by the UAV's onboard inertial navigation and flight control system, where wx_uav, wy_uav, and wz_uav are the components of w_uav in the x, y, and z axes, respectively, and vx_uav, vy_uav, and vz_uav are the components of v_uav in the x, y, and z axes, respectively, with a sampling period of delta_T_uav.

2. The raindrop size distribution estimation method based on phased array radar and meteorological UAV according to claim 1, characterized in that: The time alignment based on collected phased array radar echo data and meteorological UAV detection data, and the spatial alignment based on radar operational status data and meteorological UAV detection data using a spatial geographic coordinate to radar polar coordinate model, include: A1. Based on the collected phased array radar observation time series t0_par and meteorological UAV observation time series t0_uav, for a certain moment t0_uav_1 in the meteorological UAV observation time series t0_uav, calculate the time difference between that moment t0_uav_1 and each moment in the phased array radar observation time series t0_par, and form a time difference sequence delta_t. A2. Set the maximum permissible error delta_t_max for the alignment of observation time between the phased array radar and the UAV, and find the phased array radar observation time corresponding to the minimum time difference of the elements in delta_t that satisfies the condition of being less than or equal to delta_t_max. A3. Repeat steps A1 and A2 to traverse all observation times in the meteorological UAV observation time series t0_uav, and obtain the phased array radar observation time series after time alignment with all observation times in the meteorological UAV observation time series t0_uav, denoted as t1_par. Based on the phased array radar echo data of the collected phased array radar observation time series t0_par, extract the phased array radar echo data Data1_par of the phased array radar observation time series t1_par, including reflectivity factor data Z1_par, differential reflectivity factor data Zdr1_par and differential propagation phase shift rate KDP1_par; A4. Input the collected P_par, P_uav, delt_R_par, and delt_sita_par into the spatial geographic coordinate to radar polar coordinate conversion model, run the model, and obtain the raindrop size distribution sequence of the meteorological UAV under the meteorological UAV observation time series t0_uav after conversion to radar polar coordinates, denoted as ND_uav=[N_uav1, N_uav2, …, N_uav B ], where B is the total number of particle size partitions, N_uav b Represents the b-th particle size interval [D b-1 D b The raindrop number density within the area, b=1,2,…,B, and the corresponding spatial position information in radar polar coordinates, are denoted as P0_uav=(r0_uav, phi0_uav, sita0_uav), where r0_uav is the radial distance, phi0_uav is the azimuth, and sita0_uav is the antenna elevation angle. A5. Based on the P0_uav information and combined with the spatial location information of Data1_par, extract the phased array radar echo data after spatial alignment with P0_uav according to the principle of shortest spatial distance, and denot it as Data2_par, which includes reflectivity factor data Z2_par, differential reflectivity factor data Zdr2_par and differential propagation phase shift rate KDP2_par.

3. The raindrop size distribution estimation method based on phased array radar and meteorological UAV according to claim 2, characterized in that: The process of constructing a steady-state time segment set based on the data obtained in S1 and the collected radar operating status data includes: According to the raindrop size distribution sequence ND_uav obtained by the meteorological UAV at the meteorological UAV observation time series t0_uav, and the sampling period delta_T_uav, set the sliding window length as Wind_move = L×delta_T_uav. Take the start time of the time series t0_uav as the starting position of the sliding window, and perform sliding processing on the time series t0_uav in turn according to the window width Wind_move. Finally, decompose the time series t0_uav into K sliding time windows, denoted as Wk, k = 1, 2, …, K; Extract the raindrop size distribution sequences corresponding to the K sliding time windows from the raindrop size distribution sequence ND_uav respectively, and perform logarithmic transformation, denoted as log_N_uav_k. Then calculate the standard deviation of the log_N_uav_k sequence within each sliding time window, denoted as STDk; Set the steady-state threshold as STD0, and the condition for steady-state segment screening is set as: STDk < STD0. Select the windows that meet the condition from the K sliding time windows, and collect the observation times corresponding to these windows that meet the screening condition to obtain the steady-state segment time set, denoted as STEADY.

4. The raindrop size distribution estimation method based on phased array radar and meteorological UAV according to claim 2, characterized in that: The construction of the disturbance intensity index set according to the collected phased array radar echo data, meteorological UAV detection data and flight state data includes: B1. According to the meteorological UAV at the meteorological UAV observation time series t0_uav, flight state data w_uav, v_uav and RPM_uav, for a certain moment t0_uav_i in the meteorological UAV observation time series t0_uav, first calculate the angular velocity modulus ||w_uav|| and velocity modulus ||v_uav|| at this moment according to the Euclidean norm calculation formula. Then find the maximum value RPM_uav_max in the propeller speed RPM_uav, and calculate the propeller speed normalization term m_PRM_uav = RPM_uav / RPM_uav_max; B2. Construct the angular velocity atmospheric disturbance characteristic term u1_uav = ||w_uav|| / w_ref, the velocity atmospheric disturbance characteristic term u2_uav = ||v_uav|| / v_ref, and the propeller speed atmospheric disturbance characteristic term u3_uav = m_PRM_uav. Furthermore, construct the disturbance intensity index of the UAV operating state on the particles in the sampling area at the moment t0_uav_i as dis_i = a1×u1_uav + a2×u2_uav + a3×u3_uav, where w_ref and v_ref are the calibration values of the UAV system angular velocity and velocity respectively, and a1, a2 and a3 are the initially set disturbance weight coefficients of the characteristic terms u1_uav, u2_uav and u3_uav respectively; B3. Update the time t0_uav_i at a certain time to the value of each time in the time series t0_uav. Repeat steps B1 and B2 to generate the disturbance intensity index for each time in the time series t0_uav. Set these indices together to obtain the disturbance intensity index set, denoted as DIS.

5. The raindrop size distribution estimation method based on phased array radar and meteorological UAV according to claim 2, characterized in that: The construction of the fluctuation loss function based on the data obtained from S1, the set of steady-state time segments, and the set of disturbance intensity indices includes: C1. Based on the raindrop size distribution sequence ND_uav obtained by the meteorological drone under the observation time series t0_uav, combined with the steady-state time segment set STEADY, firstly extract any observation time t1 and its adjacent time t2 from the time set STEADY, where t2=t1+delta_T_uav, and extract the corresponding raindrop size distribution sequences N_uav_t1 and N_uav_t2. C2. Calculate the logarithmic particle size distributions of N_uav_t1 and N_uav_t2, denoted as log_N_uav_t1=log(N_uav_t1+ext_pos) and log_N_uav_t2=log(N_uav_t2+ext_pos), respectively. Then calculate the difference between the two, denoted as delta_log_N_uav= log_N_uav_t1- log_N_uav_t2, where ext_pos is a very small positive value. C3. Based on the disturbance intensity index set DIS, extract the disturbance intensity index corresponding to time t1, denoted as dis_t1, and construct the disturbance weight function R_t1=1+alpha0_dis×dis_t1, where alpha0_dis is the disturbance amplification coefficient constant. Then, construct the adjacent time fluctuation loss L_t1= R_t1×sum(delta_log_N_uav) for time t1. 2 ), where sum() is the summation across all intervals of the particle size partition; C4. Update any observation time t1 sequentially to all observation times in the steady-state segment time set STEADY and repeat steps C1 to C3 to generate the fluctuation loss for all observation times in STEADY. Then sum these fluctuation losses to generate the fluctuation loss sum within the steady-state segment, denoted as LOSS.

6. The raindrop size distribution estimation method based on phased array radar and meteorological UAV according to claim 4, characterized in that: The optimized set of disturbance intensity indices and the raindrop size distribution sequence corrected for disturbance effects include: Let the discretization step size delta_a = 1 / N, where N is the discretization parameter. By enumerating non-negative integer triples (i,j,k) satisfying i+j+k=N, we construct perturbation weight coefficients a1=i×delta_a, a2=j×delta_a, a3=k×delta_a, generating a set of G candidate weight coefficients Ag=[a1,a2,a3]. g g = 1, 2, ..., G, where G is the total number of groups; Each set of weights in the candidate weight coefficient set Ag is sequentially input into the fluctuation loss and LOSS within the steady-state segment. The fluctuation loss sum of all observation times within STEADY is calculated for each set of weights. The fluctuation loss sums of G different weight coefficient sets are aggregated and denoted as LOSS_g. Find the minimum value in the LOSS_g set, and use the combination of weight coefficients corresponding to the minimum value as the optimal weight parameters for calculating the perturbation intensity index, denoted as [a1_opt, a2_opt, a3_opt]. Replace the initial [a1, a2, and a3] with [a1_opt, a2_opt, a3_opt], ​​and re-execute steps B1-B3 to generate an optimized perturbation intensity index set, denoted as DIS_opt; Construct a function to determine the effect of perturbation intensity on particle size distribution: F_DIS = 1 + apha1_dis × DIS_opt, where apha1_dis is a correction coefficient, based on [D b-1 D b Constructing particle size D b The interaction function with the perturbation G_DIS_D b =1 / (1+apha2_dis×D b ), b=1,2,…,B, where apha2_dis are the relation coefficients, and the B interaction functions are set together and denoted as G_DIS; Based on the obtained ND_uav, combined with the already obtained F_DIS and G_DIS, the raindrop size distribution sequence after perturbation correction is generated as ND_uav_c = ND_uav × F_DIS × G_DIS.

7. The raindrop size distribution estimation method based on phased array radar and meteorological UAV according to claim 6, characterized in that: S4 includes: Set the radial distance distinction threshold to Rang_set. Based on the radial distance parameter in the spatial information of the phased array radar echo data, the echo data that is less than or equal to Rang_set is divided into the near-range echo area, and the radar reflectivity factor data of the corresponding area is denoted as Z2_par_near. The echo data that is greater than Rang_set is divided into the far-range echo area, and the radar reflectivity factor data of the corresponding area is denoted as Z2_par_far. For radar echo data in the near range, a deep convolutional neural network model is used, denoted as MODEL1. For radar echo data in the far range, a simple convolutional neural network model is used, denoted as MODEL2. The loss function for both models is the mean square error function, denoted as MSE_MODEL1 and MSE_MODEL2. MODEL1 and MODEL2 are connected in parallel to form a whole model MODEL. For the near echo region, the weight is set as weight_near=mean(Z2_par_near) / max(Z2_par_near). For the far echo region, the weight is set as weight_far=min(Z2_par_far) / mean(Z2_par_far). Then, the loss function of the whole model MODEL is set as loss_MODEL=weight_near×MSE_MODEL1+ weight_far×MSE_MODEL2. Based on the obtained phased array radar echo data, including reflectivity factor data Z2_par, differential reflectivity factor data Zdr2_par, and differential propagation phase shift rate KDP2_par, the ratio of reflectivity factor to differential reflectivity factor Zdr_Z=Zdr2_par / Z2_par is constructed, and the combination of differential propagation phase shift rate and reflectivity factor KDP_Z=KDP2_par×Z2_par is constructed. Z2_par, Zdr2_par, KDP2_par together with Zdr_Z and KDP_Z constitute the input dataset of the model MODEL, while the raindrop size distribution sequence ND_uav_c, which is corrected for the influence of perturbation, is used as the output dataset of the model MODEL.

8. The raindrop size distribution estimation method based on phased array radar and meteorological UAV according to claim 7, characterized in that: S5 includes: Based on the corrected raindrop size distribution sequence ND_uav_c and the optimized perturbation intensity index set DIS_opt, for each particle diameter interval [D] in ND_uav_c... b-1 D b Construct the corresponding weights, denoted as weight_D. b =beta1×DIS_opt+beta2×[D b / D max ], where beta1 and beta2 are both dynamically adjusted parameters, D max This represents the maximum particle diameter in ND_uav_c; Based on the obtained loss_MODEL, combined with weight_D b The optimized overall model loss function is loss_MODEL_opt=sum(weight_D b ×loss_MODEL), where sum() is the summation over all particle diameter distribution intervals; During the optimization process of beta1 and beta2, the Bayesian optimization method was selected to adjust the parameters and train the model, and finally the optimal parameters of beta1 and beta2 were obtained. After the model training is completed, Z2_par, Zdr2_par, KDP2_par and Zdr_Z, KDP_Z are input into the trained model again for calculation, thereby realizing the estimation of raindrop number density in different particle size ranges, that is, completing the estimation of raindrop particle size distribution based on phased array radar and meteorological UAV.

Images