Ground precipitation intensity estimation method based on phased array radar echo and geographic factors

Abstract

The present application relates to the method for estimating ground precipitation intensity based on phased array radar echo and geographic factors, belongs to the field of radio, first polar coordinate multivariate radar echo is mapped with three-dimensional wind field and is interpolated to unified three-dimensional grid, forms consistent space-time data foundation, then based on local density clustering, a variety of rainfall type regions are automatically divided, and time lag field is constructed according to precipitation type, adaptive time alignment of radar echo and ground station observation is realized, three-dimensional wind field and raindrop terminal speed are introduced to calculate precipitation falling drift, the position of the grid point matched with the ground station is corrected, and accurate spatial matching of radar echo and ground station observation is completed, terrain lifting index is constructed by fusing into ground layer wind direction, DEM, slope, aspect and slope position type, radar nonlinear precipitation potential vector and terrain characteristics are coupled through physical constraint coupling model, quantitative estimation of ground precipitation intensity is realized, and the reliability and applicability of ground precipitation estimation under complex environment are improved.

Description

Technical Field

[0001] This invention relates to the field of radio, and more particularly to a method for estimating surface precipitation intensity based on phased array radar echoes and geographical factors. Background Technology

[0002] Surface precipitation intensity is a key parameter for monitoring, early warning, and risk assessment of rainstorms, floods, geological disasters, and urban waterlogging. Phased array weather radars, with their fast volume scanning speed, high temporal resolution, and ability to acquire multi-elevation angle and multivariable three-dimensional echoes, provide an important data source for high spatiotemporal resolution quantitative precipitation estimation. However, radar observations reflect the scattering characteristics of precipitation particles within a certain altitude range. Surface precipitation is also affected by the falling process, wind transport, and orographic uplift, resulting in significant spatiotemporal inconsistencies. Therefore, there is an urgent need for high-precision surface precipitation estimation methods tailored to the characteristics of phased array radars.

[0003] Current methods for estimating ground precipitation using phased array radar mainly include: 1. The first type is quantitative precipitation estimation based on empirical relationships. These methods typically utilize empirical relationships between reflectivity factors and precipitation intensity (such as the ZR relationship) or combine dual-polarization parameters such as differential reflectivity factors and differential phase shift rates to construct empirical or semi-empirical formulas. These methods are simple in structure and computationally efficient, but their parameters usually depend on specific regions, specific precipitation types, and statistical samples. Once the microphysical properties of precipitation or environmental conditions change, the estimation error increases significantly, making it difficult to adapt to complex precipitation processes. 2. The second type is physical inversion methods based on dual-polarization variables. These methods enhance the characterization of raindrop particle size distribution and water content by introducing variables such as differential reflectivity factors and differential phase shift rates, improving estimation accuracy to some extent under heavy precipitation conditions. However, these methods still mainly target precipitation characteristics at the radar observation altitude and do not fully consider the spatiotemporal evolution of precipitation particles during their descent and the influence of topography, limiting their applicability to mountainous or complex terrain areas. 3. The third category of methods consists of precipitation estimation methods based on machine learning or deep learning. In recent years, some studies have attempted to use models such as neural networks and random forests to directly map radar observation variables to ground precipitation intensity. These methods have improved the fitting ability to some extent, but most studies only use radar variables as input, ignoring the systematic influence of orographic lifting, windward slope effect, and wind transport on precipitation distribution during the precipitation process. Furthermore, existing methods often perform simple temporal or spatial matching between radar and ground stations, failing to adequately address the time lag and horizontal drift between radar echoes and ground precipitation, resulting in high sample noise and limited model generalization ability. In addition, the high temporal resolution advantage of phased array radar has not been fully utilized in existing research. Most methods still follow the processing approach of traditional volume scanned radar, lacking adaptive time-delay modeling and spatial matching mechanisms for different precipitation intensity regions, making it difficult to accurately characterize the dynamic relationship between radar echoes and ground precipitation responses under different precipitation types. Summary of the Invention

[0004] The purpose of this invention is to overcome the shortcomings of the prior art and to provide a method for estimating ground precipitation intensity based on phased array radar echoes and geographical factors, thus solving the deficiencies of the prior art.

[0005] The objective of this invention is achieved through the following technical solution: a method for estimating surface precipitation intensity based on phased array radar echoes and geographic factors, wherein the estimation method includes:

[0006] S1. Collect phased array radar echo data, radar operation status data, ground automatic station observation data, and geographic factor data. Input the radar operation status data into the radar beam center approximate height estimation model and map projection model. Combine the results with the radar echo data to construct a target 3D mesh. Process the points in the 3D mesh to generate a 3D wind field dataset for the target 3D mesh at different times, as well as a phased array radar reflectivity factor dataset, differential reflectivity factor dataset, and differential phase shift rate dataset. After normalization, generate a normalized reflectivity factor dataset, differential reflectivity factor dataset, and differential phase shift rate dataset. Combine the local clustering model and local density calculation model to divide the region into several rainfall types.

[0007] S2. Based on the collected radar operation status data, ground automatic station observation data, reflectivity factor dataset, and determined rainfall type areas, obtain the regional average reflectivity factor sequence and regional average rainfall sequence for the rainfall type areas. Based on the set time lag sequence, regional average reflectivity factor sequence, and regional average rainfall sequence, obtain the optimal lag time difference between the ground automatic station observation data and the phased array radar observation. Delay the reflectivity factor dataset by the optimal lag time difference to obtain a unified radar reflectivity factor dataset for several rainfall types. Combine the differential reflectivity factor dataset and differential phase shift rate dataset to generate a time-aligned and integrated differential reflectivity factor dataset and differential phase shift rate dataset that are consistent with the ground automatic station observation dataset.

[0008] S3. Based on the collected ground automatic station observation data, three-dimensional wind field dataset, and merged unified radar reflectivity factor dataset, the average level drift distance of raindrops after falling to the ground is calculated by combining the radar reflectivity factor and raindrop terminal velocity relationship model, and then a phased array radar echo dataset is generated after spatiotemporal alignment and spatial matching with the ground automatic station observation dataset.

[0009] S4. Generate a basic topographic lifting index based on the three-dimensional wind field dataset and the collected geographic factor data. Construct a topographic feature vector based on the data output from S3, the collected geographic factor data, and the basic topographic lifting index.

[0010] S5. Construct a radar-terrain coupled physical constraint model based on the ground automatic station observation dataset and the data output from S4, and train it. Finally, use this model to estimate the ground precipitation intensity.

[0011] The phased array radar echo data includes: radar reflectivity factor data Z0_par, differential reflectivity factor data Zdr0_par, differential phase shift rate data KDP0_par, and three-dimensional wind field data (U0, V0, W0) obtained by four-dimensional variational assimilation and inversion at different times in polar coordinate spatial locations (r0, phi0, seat0). Here, r0 is the radial distance, phi0 is the azimuth, seat0 is the antenna elevation angle, and the three-dimensional wind field data refers to (U0, V0, W0), where U0 is the eastward component of the horizontal wind speed, V0 is the northward component of the horizontal wind speed, and W0 is the vertical wind speed.

[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; the radar volume scanning period delt_T; and the start time t0 of the radar scan.

[0013] The ground automatic station observation data includes: the rainfall rate R0_sta observed at different times, whose spatial geographic coordinates are (x_R0, y_R0), where x_R0 and y_R0 are the position information in the east-pointing X-axis and north-pointing Y-axis directions, respectively;

[0014] The geographic factor data includes: Digital Elevation Data (DEM), Slope Data (SLOPE), Aspect Data (ASPECT), and Slope Type Data (TPI).

[0015] The process of inputting radar operating status data into the radar beam center approximate height estimation model and map projection model, and combining the results with radar echo data to construct a target 3D mesh, and processing the points in the 3D mesh to obtain a 3D wind field dataset for generating the target 3D mesh at different times, as well as a phased array radar reflectivity factor dataset, a differential reflectivity factor dataset, and a differential phase shift rate dataset, specifically includes:

[0016] A1. Input the collected radar site spatial location information P_par, (r0, phi0, sita0), delt_R_par, and delt_sita_par into the radar beam center approximate height estimation model, run the model, and output the radar beam center approximate height z0.

[0017] A2. Input the collected spatial location information P_par of the radar station into the map projection model, run the model, and output the planar projection coordinates (X_par, Y_par) of the radar station. Based on the collected (r0, phi0, sita0), delt_R_par, and delt_sita_par, construct a ray with an azimuth angle of phi0 in the planar coordinates. Take the horizontal distance as s0 = r0 × cos(sita0). Then the projection point of the ray centerline on the horizontal plane is x0 = X_par + s0 × sin(phi0) and y0 = Y_par + s0 × cos(phi0). Combine the point (x0, y0) with the acquired z0 to convert the polar coordinates of the radar echo data into three-dimensional spatial coordinates (x0, y0, z0).

[0018] A3. Set the spatial resolution of the horizontal grid in the X and Y directions to delt_x and delt_y respectively, and set the resolution in the vertical height to delt_z, thereby constructing the target 3D grid. For a point M in the target 3D grid, its spatial coordinates are (x_m, y_m, z_m). Calculate the radial distance r_m and azimuth angle phi_m from this point to the radar station's planar projection coordinates (X_par, Y_par). For any point I in (x0, y0, z0), its coordinates are (x_i, y_i, z_i). Based on the collected radar echo data, extract the reflectivity factor data and 3D wind field data corresponding to point I, denoted as Z0_par_i and (U0_i, V0_i, ..., ...) respectively. W0_i), where U0_i is the eastward component of the horizontal wind speed at point I, V0_i is the northward component of the horizontal wind speed at point I, and W0_i is the vertical wind speed at point I. Calculate the coordinate difference between point I and point M, denoted as (delta_x_i, delta_y_i, delta_z_i). Then decompose delta_x_i and delta_y_i into the transverse direction along the radial direction, which is consistent with the beamwidth direction, denoted as (delta_r_i, delta_b_i).

[0019] A4. Set the radial, lateral, and vertical weights to apha_r, apha_b, and apha_z, respectively. Calculate the anisotropic distance at point I as distance_i = apha_r × delta_r_i. 2 +apha_b×delta_b_i 2 +apha_z×delta_z_i 2Then, the inverse distance weighting index is calculated as W_i=1 / distance_i. The reflectivity factor Z_par_m and the three-dimensional wind field data (U0_i_m, V0_i_m, W0_i_m) at point M in the target three-dimensional grid are calculated using the weighted average method. U0_i_m is the eastward component of the horizontal wind speed at point M, V0_i_m is the northward component of the horizontal wind speed at point M, and W0_i_m is the vertical wind speed at point M.

[0020] A5. Apply steps A3 and A4 to all points in the target 3D mesh at different times to generate phased array radar reflectivity factor dataset Z_par and 3D wind field dataset (U, V, W) of the target 3D mesh at different times, where U, V and W are the eastward component of the horizontal wind speed, the northward component of the horizontal wind speed and the vertical wind speed of the target 3D mesh at different times, respectively.

[0021] A6. Replace the reflectivity factor datasets in steps A2-A5 with the differential reflectivity factor dataset Zdr0_par and the differential phase shift rate dataset KDP0_par in sequence to generate phased array radar differential reflectivity factor dataset Zdr_par and differential phase shift rate dataset KDP_par for the target 3D mesh at different times.

[0022] The process of generating normalized reflectivity factor datasets, differential reflectivity factor datasets, and differential phase shift rate datasets, and then dividing the regions into several rainfall types using local clustering models and local density calculation models, specifically includes:

[0023] Find the maximum value of the generated reflectivity factor dataset Z_par, and divide each reflectivity factor data by the maximum reflectivity factor to generate normalized reflectivity factor data. Apply the same processing method to the differential reflectivity factor dataset and the differential phase shift rate dataset to generate normalized differential reflectivity factor dataset and differential phase shift rate dataset.

[0024] The normalized reflectance factor, differential reflectance factor, and differential phase shift rate datasets are input together into a density-based local clustering model. The initial search radius r_orig around the data points and the minimum point value MinPts of the cluster are set. The local clustering model is run to generate regions with different precipitation intensity types that are initially divided. The average reflectance factor of each region is calculated. Based on the relationship between the average reflectance factor and the precipitation intensity classification, each precipitation region is initially labeled as a light rain, moderate rain, heavy rain, or rainstorm region.

[0025] The generated reflectivity factor dataset Z_par is input into the local density calculation model. The local density calculation model is run to generate the local density value of each reflectivity factor data point. Based on the local density value of each reflectivity factor data point, the local density data of the initially marked light rain, moderate rain, heavy rain and rainstorm areas are statistically analyzed to obtain the average local density ρ of each precipitation area.

[0026] Based on the optimization calculation formula of the search radius based on local density, the optimized search radii r_opt of light rain, moderate rain, heavy rain and rainstorm areas are obtained respectively. Then, the initial search radius r_orig around the data points set in the local clustering model is updated to r_opt, and the local clustering model is run again. The output is the area that is accurately divided into light rain, moderate rain, heavy rain and rainstorm types.

[0027] S2 specifically includes the following:

[0028] B1. Based on the collected delt_T and t0, construct the phased array radar observation time series ts=t0+k×delt_T, k=0,1,2,…,K, where K is the total number of radar scanning cycles. Based on the phased array radar reflectivity factor dataset Z_par at the target 3D grid at different output times, extract the phased array radar reflectivity factor dataset corresponding to the time series ts, denoted as Z_par_ts. For areas identified as light rain type, extract the phased array radar reflectivity factor dataset corresponding to the area from Z_par_ts, denoted as Z_par_ts_light. Extract the regional average reflectivity factor sequence of the light rain type area corresponding to the time series ts from Z_par_ts_light using the spatial averaging method, denoted as Z_par_ts_light_m.

[0029] B2. Based on the collected automatic weather station observation data R0_sta at different times, a linear interpolation method is used in time to generate automatic weather station rainfall rate observation data on the time series ts, denoted as R0_sta_ts. For areas identified as light rain, the corresponding automatic weather station observation dataset for that area is extracted from R0_sta_ts, denoted as R0_sta_ts_light. Then, the regional average rainfall rate sequence for the light rain type area corresponding to the time series ts is extracted from R0_sta_ts_light using the spatial averaging method, denoted as R0_sta_ts_ligh_m.

[0030] B3. Set a time lag sequence with step size delt_T: time_delay=m×delt_T, where m=0,1,2,…,M-1, and M is the total number of lags. For a certain lag value τ in the time lag sequence, select the sequence corresponding to the time lag value τ from the R0_sta_ts_ligh_m sequence, and denote it as R0_sta_ts_ligh_m1. Calculate the correlation coefficient between Z_par_ts_light_m and R0_sta_ts_ligh_m1.

[0031] B4. Update the lag value τ sequentially with all M values ​​in time_delay, and repeat step B3 to calculate the correlation coefficient between Z_par_ts_light_m and the rainfall rate data of the ground automatic station after M lag values, forming a correlation coefficient sequence, denoted as CC. Identify the maximum value in the CC sequence and extract the time lag value corresponding to the maximum value, which is the optimal lag time difference between the ground automatic station observation data and the phased array radar observation, denoted as time_delay_light. Delay Z_par_ts_light by time_delay_light to form a phased array radar reflectivity factor dataset that is time-aligned with the ground automatic station observation data R0_sta_ts_light in the light rain area, denoted as Z1_par_ts_light.

[0032] B5. Applying steps B1-B4 to the moderate rain, heavy rain, and storm areas respectively, the optimal time lag differences (time_delay_middle, time_delay_heavy, time_delay_storm) corresponding to these areas can be calculated. This leads to the generation of ground automatic station observation datasets (R0_sta_ts_middle, R0_sta_ts_heavy, R0_sta_ts_storm) and a time-aligned phased array radar reflectivity factor dataset (Z1_par_ts_mi). ddle, Z1_par_ts_heavy, Z1_par_ts_storm), where R0_sta_ts_middle, R0_sta_ts_heavy, and R0_sta_ts_storm are the ground automatic station observation data in the moderate rain, heavy rain, and rainstorm areas, respectively, and Z1_par_ts_middle, Z1_par_ts_heavy, and Z1_par_ts_storm are the phased array radar reflectivity factor datasets that are time-aligned with the ground automatic station observation datasets in the moderate rain, heavy rain, and rainstorm areas, respectively;

[0033] B6. Merge the radar reflectivity factor datasets [Z1_par_ts_light, Z1_par_ts_middle, Z1_par_ts_heavy, Z1_par_ts_storm] under different precipitation intensities into a unified radar reflectivity factor dataset Z1_par_ts, and integrate the corresponding automatic weather station observation data [R0_sta_ts_light, R0_sta_ts_middle, R0_sta_ts_heavy, R0_sta_ts_storm] into an automatic weather station precipitation dataset R1_sta_ts arranged by precipitation intensity;

[0034] B7. Based on the identified areas of light rain, moderate rain, heavy rain, and torrential rain, and combined with the phased array radar differential reflectivity factor dataset Zdr_par and differential phase shift rate dataset KDP_par of the target 3D mesh at different times, extract the differential reflectivity factor and differential phase shift rate of the corresponding rainfall area respectively.

[0035] B8. Using the optimal time lags (time_delay_middle, time_delay_heavy, time_delay_storm) corresponding to the acquired moderate rain, heavy rain, and storm areas, the differential reflectivity factor and differential phase shift rate of each rainfall area are delayed. Then, the processed data are merged and unified to finally generate a differential reflectivity factor dataset Zdr1_par_ts and a differential phase shift rate dataset KDP1_par_ts that are time-aligned with and integrated with the ground automatic station observation dataset.

[0036] S3 specifically includes the following:

[0037] C1. Select a station G from the collected ground automatic station observation data. Its geographical coordinates are (x_R0_G, y_R0_G). Combine the output Z1_par_ts and use the nearest neighbor spatial matching method to find the phased array radar echo grid point P that is directly above station G. The position information on the horizontal plane is recorded as (x_par_P, y_par_P).

[0038] C2. Taking point P as the spatial center point, set the horizontal search window for inclined precipitation to N_h×N_h, where N_h is the width of the window. The number of vertical matching layers for inclined precipitation is N_v. Determine the adjacent spatial domain above station G as Vol_space=N_h×N_h×N_v. Based on the three-dimensional wind field dataset (U, V, W) and Z1_par_ts, combined with the range of the adjacent spatial domain Vol_space, calculate the average values ​​of the U and V components at the i-th height layer, denoted as U_i and V_i, and the average value of the reflectivity factor data at the i-th height layer, denoted as Z_i, i=1, 2, …, N_v;

[0039] C3. Input Z_i into the model relating reflectivity factor and raindrop terminal velocity, run the model, and generate the terminal velocity of the raindrop at the i-th height layer, denoted as Vt_i. Based on the set delt_z, obtain the height of the i-th height layer grid, Height_i = i × delt_z. Combined with the obtained Vt_i, obtain the time difference between the raindrop falling to the ground at the i-th height layer, delta_t_i = Height_i / Vt_i. Based on the obtained U_i and V_i, calculate the distances of raindrop drift in the X and Y directions within the time difference delta_t_i, respectively: delta_X_i = U_i × delta_t_i, delta_Y_i = V_i × delta_t_i. Based on the average reflectivity factor data Z_i of all N_v height layers, calculate the weighting coefficient of reflectivity factor at the i-th height layer, apha_i = Z_i / sum(Z_i), where sum() is the summation function.

[0040] C4. Based on the obtained apha_i, delta_X_i, and delta_Y_i, calculate the weighted average distance of the raindrops drifting in the X and Y axes in the adjacent airspace above the automatic ground station G: delta_X = sum(delta_X_i × apha_i) and delta_Y = sum(delta_Y_i × apha_i). Based on delta_X and delta_Y, and combined with the determined position information (x_par_P, y_par_P) of the nearest phased array radar echo grid point P directly above the automatic ground station G on the horizontal plane, adjust the position of point P to generate a phased array radar echo grid point Q that precisely matches the airspace above the automatic ground station G, with the horizontal position information being (x_par_P + delta_X, y_par_P + delta_Y).

[0041] C5. Apply steps C1-C4 to all automatic ground stations to generate phased array radar reflectivity factor, differential reflectivity factor, and differential phase shift rate datasets that have been time-aligned and spatially matched with the automatic ground station observation dataset R1_sta_ts, denoted as Z1_par_ts_match, Zdr1_par_ts_match, and KDP1_par_ts_match, respectively.

[0042] The generation of the basic topographic lifting index based on the three-dimensional wind field dataset and collected geographic factor data specifically includes:

[0043] Based on the target 3D wind field dataset (U, V, W), the horizontal wind field (U_low, V_low) of the lowest altitude layer is extracted. U_low is the eastward component of the horizontal wind speed at the lowest altitude layer, and V_low is the northward component of the horizontal wind speed at the lowest altitude layer. The near-surface wind direction W_dir=atan2(U_low, V_low) is calculated, where atan2() returns the clockwise angle of the wind direction relative to due north.

[0044] Based on the collected slope aspect data ASPECT, the angle between the slope aspect and the wind direction is calculated as delta_angle = ASPECT - W_dir. Then, the windward slope factor F_wind = max(0, cos(delta_angle × π / 180)) is calculated, where max() is the function to find the maximum value.

[0045] Find the maximum value of the collected slope data SLOPE, denoted as max_SLOPE, and calculate the slope factor F_slope = SLOPE / max_SLOPE. Based on the collected digital elevation data DEM, set the neighborhood search radius R_DEM, and calculate the average digital elevation within each search radius of R_DEM centered on each DEM grid point. Then, calculate the elevation difference of each grid point relative to the average digital elevation within its respective search radius, denoted as delta_DEM. Normalize delta_DEM using the hyperbolic tangent function to generate F_DEM = 0.5 × (1 + tanh(delta_DEM / DEM0)), where tanh() is the hyperbolic tangent function and DEM0 is the scale parameter.

[0046] Based on the collected slope type data TPI, and according to the slope type mapping weight rule, the slope weight coefficient W_TPI is generated. Based on the generated F_wind, F_slope, F_DEM and W_TPI, the basic terrain lift index terrain_lift=W_TPI×F_wind×F_slope×F_DEM is generated.

[0047] The construction of the terrain feature vector based on the data output by S3, the collected geographic factor data, and the basic terrain lifting index specifically includes the following:

[0048] Based on Z1_par_ts_match, Zdr1_par_ts_match, and KDP1_par_ts_match, we first construct the joint reflectivity and phase feature F1 = Z1_par_ts_match × KDP1_par_ts_match, and then construct the particle size distribution index feature F2 = Zdr1_par_ts_match / Z1_par_ts_match, where Z1_par_ts_match, Zdr1_par_ts_match, and KDP1_par_ts_match are the phased array radar reflectivity factor, differential reflectivity factor, and differential phase shift rate datasets after time alignment and spatial matching, respectively.

[0049] A nonlinear radar precipitation potential vector R_radar=[Z1_par_ts_match, Zdr1_par_ts_match, KDP1_par_ts_match, F1, F2] is constructed. Based on the collected DEM, SLOPE, ASPEC, TPI and terrain_lift datasets, the digital elevation, slope, aspect, slope position and basic terrain lifting index datasets that spatially match the R1_sta_ts dataset are extracted according to the nearest neighbor principle and denoted as DEM1, SLOPE1, ASPEC1, TPI1 and terrain_lift1, respectively.

[0050] Construct slope and aspect dynamic features T_d=[SLOPE1, sin(ASPECT1), cos(ASPECT1)], construct terrain morphology and lift features T_m=[TPI1, terrain_lift1], and then construct terrain feature vector T_terrain=[DEM1, T_d, T_m].

[0051] S5 specifically includes the following:

[0052] Based on the output R1_sta_ts, as well as R_radar and T_terrain, an end-to-end radar terrain coupled physical constraint model MODEL consisting of three sub-models is first constructed. Sub-model 1 is composed of a multilayer perceptron model. The input dataset and output dataset of sub-model 1 during training are R_radar and R1_sta_ts, respectively. Sub-model 1 is trained, and the output after training is the radar precipitation potential P_radar ignoring the influence of terrain.

[0053] Sub-model 2 is composed of an XGBoost model. The input and output datasets for training sub-model 2 are T_terrain and R1_sta_ts, respectively. After training sub-model 2, the output is the precipitation-terrain modulation factor M_terrain.

[0054] Sub-model 3 is composed of physical constraint function mapping. The output of sub-model 3 is R_out=P_radar×exp(M_terrain), where exp() is the logarithmic function. The loss function of the radar terrain coupled physical constraint model MODEL is set to the mean square error loss function. The radar terrain coupled physical constraint model MODEL is trained.

[0055] After the radar-terrain coupled physical constraint model MODEL is trained, R1_sta_ts, R_radar, and T_terrain are input into the trained radar-terrain coupled physical constraint model MODEL again, and the model is run to realize the estimation of ground precipitation intensity based on phased array radar echo and geographical factors.

[0056] This invention has the following advantages: The ground precipitation intensity estimation method based on phased array radar echoes and geographic factors first maps and interpolates polar coordinate multivariate echoes with three-dimensional wind fields to a unified three-dimensional grid, forming a consistent spatiotemporal data foundation. Then, it automatically divides areas into light, moderate, heavy, and torrential rain zones based on local density clustering, and constructs time lag fields according to precipitation type to achieve adaptive time alignment between radar and ground station observations. Furthermore, it introduces three-dimensional wind fields and raindrop terminal velocities to calculate precipitation drift, correcting the spatial matching grid point positions with ground station data, thus achieving accurate spatial matching between radar echoes and ground station observations. It integrates DEM, slope, aspect, and slope type to construct a topographic lifting index, and jointly models the radar nonlinear precipitation potential vector and topographic features through a physical constraint coupling model, achieving quantitative estimation of ground precipitation intensity. By employing "zone identification—zone time lag correction—drift correction—radar-topographic coupling," it reduces key error sources and improves the reliability and applicability of ground precipitation estimation in complex environments. Attached Figure Description

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

[0058] Figure 2 A flowchart illustrating the implementation of this invention for generating phased array radar echoes and three-dimensional wind field data of a target three-dimensional mesh at different times;

[0059] Figure 3 A flowchart illustrating the implementation of accurate regional division for light rain, moderate rain, heavy rain, and rainstorm types in this invention;

[0060] Figure 4This is a flowchart illustrating the implementation of the present invention for generating a phased array radar echo dataset that is time-aligned with the ground automatic station observation dataset.

[0061] Figure 5 A flowchart illustrating the implementation of the present invention for generating a phased array radar echo dataset that has been time-aligned and spatially matched with a ground automatic station observation dataset;

[0062] Figure 6 A flowchart illustrating the implementation of the basic terrain lifting index for this invention;

[0063] Figure 7 This is a flowchart illustrating the implementation of the nonlinear radar precipitation potential vector and terrain feature vector generation method of the present invention.

[0064] Figure 8 This is a flowchart illustrating the implementation of the surface precipitation intensity estimation based on phased array radar echo and geographic factors according to the present invention. Detailed Implementation

[0065] 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 a part of the embodiments of this application, and not all of the 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.

[0066] like Figure 1 As shown, this invention specifically relates to a method for estimating surface precipitation intensity based on phased array radar echoes and geographic factors, which specifically includes the following:

[0067] Step (1): Collect phased array radar echo data and radar operation status data; collect encrypted ground automatic station observation data; collect high spatial resolution geographic factor data;

[0068] Among them, the phased array radar echo data refers to the radar reflectivity factor data Z0_par, differential reflectivity factor data Zdr0_par, differential phase shift rate data KDP0_par, and three-dimensional wind field data (U0, V0, W0) obtained by four-dimensional variational assimilation and inversion at different times in the spatial location (r0, phi0, seat0) in polar coordinates. Here, r0 is the radial distance, phi0 is the azimuth, seat0 is the antenna elevation angle, and the three-dimensional wind field data refers to (U0, V0, W0), where U0 is the eastward component of the horizontal wind speed, V0 is the northward component of the horizontal wind speed, and W0 is the vertical wind speed.

[0069] 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 delt_R_par = 30 meters and the beamwidth delt_sita_par = 1.0 degree; the radar volume scanning period delt_T = 1 minute and the start time t0 of the radar scan;

[0070] Encrypted-level ground automatic station observation data refers to the rainfall rate R0_sta observed at different times, with spatial geographic coordinates (x_R0, y_R0), where x_R0 and y_R0 are the location information in the east-pointing X-axis and north-pointing Y-axis directions, respectively.

[0071] High spatial resolution geographic factor data refers to digital elevation data (DEM), slope data (SLOPE), aspect data (ASPECT), and slope position type data (TPI).

[0072] Step (2): As Figure 2As shown, firstly, the spatial location information P_par, (r0, phi0, sita0), delt_R_par, and delt_sita_par of the radar station collected in step (1) are input into the radar beam center approximate height estimation model. The model is run, and the output is the radar beam center approximate height z0. Secondly, the spatial location information P_par of the radar station collected in step (1) is input into the map projection model. The model is run, and the output is the plane projection coordinates (X_par, Y_par) of the radar station. Then, based on (r0, phi0, sita0), delt_R_par, delt_sita_par, and Z0_par collected in step (1), a ray with an azimuth angle of phi0 is constructed in the plane coordinates. The horizontal distance is approximately taken as s0 = r0 × cos(sita0). Then, the projection point of the ray centerline on the horizontal plane is x0 = X_par + s0 × sin(phi0) and y0 = Y_par + s0 × cos(phi0). Then, the point (x0, Combining y0 with the acquired z0, the polar coordinates of the radar echo data are converted into three-dimensional spatial coordinates (x0, y0, z0). The spatial resolution of the horizontal grid in the X and Y directions is set to delt_x = 30 meters and delt_y = 30 meters, respectively, and the resolution in the vertical height is set to delt_z = 30 meters, thus constructing the target's three-dimensional grid. Then, for a point M in the target's three-dimensional grid, whose spatial coordinates are (x_m, y_m, z_m), the radial distance r_m and azimuth angle phi_m from this point to the radar station's projected plane coordinates (X_par, Y_par) are calculated, where r_m = [(x_m - X_par)]. 2 +(y_m-Y_par) 2 ] 1 / 2, phi_m=arctan[(x_m-X_par) / ( y_m-Y_par)], arctan() is the arctangent function; then, for any point I in the coordinates (x0, y0, z0) in the three-dimensional space, its coordinates are (x_i, y_i, z_i), according to the radar echo data collected in step (1), the reflectivity factor data and three-dimensional wind field data corresponding to point I are extracted, and are respectively denoted as Z0_par_i and (U0_i,V0_i, W0_i), where U0_i is the eastward component of the horizontal wind speed at point I, V0_i is the northward component of the horizontal wind speed at point I, and W0_i is the vertical wind speed at point I. The coordinate difference between point I and point M is calculated and denoted as (delta_x_i, delta_y_i, delta_z_i), where delta_x_i= x_i-x_m, delta_y_i= Given y_i-y_m and delta_z_i= z_i-z_m, decompose delta_x_i and delta_y_i into transverse directions along the radial direction that are consistent with the beamwidth direction, denoted as (delta_r_i, delta_b_i), where delta_r_i=delta_x_i×sin(phi_m)+delta_y_i×cos(phi_m) and delta_b_i=delta_x_i×cos(phi_m)-delta_y_i×sin(phi_m); then, set the radial, transverse, and vertical weights to apha_r=1, apha_b=4, and apha_z=9 respectively, and calculate the anisotropic distance at point I: distance_i=apha_r×delta_r_i 2 +apha_b×delta_b_i 2 +apha_z×delta_z_i 2Then, the inverse distance weighting index W_i = 1 / distance_i is calculated. Next, the reflectivity factor Z_par_m and the three-dimensional wind field data (U0_i_m, V0_i_m, W0_i_m) at point M in the target 3D mesh are calculated using a weighted average method. Finally, the above processing flow is applied to all points in the target 3D mesh at different times to generate the phased array radar reflectivity factor dataset Z_par and the three-dimensional wind field dataset (U, V, W) for the target 3D mesh at different times. Replacing the reflectivity factor dataset in this step with the differential reflectivity factor dataset Zdr0_par and the differential phase shift rate dataset KDP0_par in sequence generates the phased array radar differential reflectivity factor dataset Zdr_par and the differential phase shift rate dataset KDP_par for the target 3D mesh at different times. Here, U0_i_m is the eastward component of the horizontal wind speed at point M, V0_i_m is the northward component of the horizontal wind speed at point M, and W0_i_m is the vertical wind speed at point M.

[0073] Furthermore, the weighted average method is as follows: , , , N is the total number of points in the dataset Z0_par;

[0074] Step (3): As Figure 3As shown, firstly, the maximum value of the reflectivity factor dataset Z_par generated in step (2) is found, and each reflectivity factor data is divided by the maximum value of the reflectivity factor to generate normalized reflectivity factor data. The same processing method is applied to the differential reflectivity factor dataset and the differential phase shift rate dataset to generate normalized differential reflectivity factor dataset and differential phase shift rate dataset. Secondly, the normalized reflectivity factor, differential reflectivity factor, and differential phase shift rate datasets are input together into the density-based local clustering model, and the initial search radius around the data points r_orig=0.12 and the minimum point value of the cluster MinPts=16 are set. The model is run to generate regions with different precipitation intensity types that are initially divided, and the average reflectance factor of each region is calculated. Based on the relationship between the average reflectance factor and the precipitation intensity classification, each precipitation region is initially labeled as light rain, moderate rain, heavy rain, and rainstorm. Then, the reflectance factor dataset Z_par generated in step (2) is input into the local density calculation model, and the model is run to generate the local density value of each reflectance factor data point. Then, based on the local density value of each reflectance factor data point, the local density data of the regions initially labeled as light rain, moderate rain, heavy rain, and rainstorm are calculated to obtain the average local density ρ of each precipitation region. Finally, based on the optimization calculation formula of the search radius based on local density, the optimized search radius r_opt of the light rain, moderate rain, heavy rain, and rainstorm regions is obtained respectively. Then, the initial search radius r_orig around the data points set in the local clustering model is updated to r_opt, the model is run again, and the output is the regions that are accurately divided into light rain, moderate rain, heavy rain, and rainstorm types.

[0075] Furthermore, the relationship between the average reflectance factor and precipitation intensity classification is as follows: if the average reflectance factor is less than 30 dBZ, the precipitation category is light rain; if the average reflectance factor is greater than or equal to 30 dBZ and less than 40 dBZ, the precipitation category is moderate rain; if the average reflectance factor is greater than or equal to 40 dBZ and less than 50 dBZ, the precipitation category is heavy rain; and if the average reflectance factor is greater than or equal to 50 dBZ, the precipitation category is torrential rain.

[0076] The optimization formula for the search radius based on local density is: r_opt=(3×MinPts / (4×π×ρ)) 1 / 3 ;

[0077] Step (4): As Figure 4As shown, firstly, based on delt_T and t0 collected in step (1), a phased array radar observation time series ts = t0 + k × delt_T is constructed, k = 0, 1, 2, ..., K, where K is the total number of radar scanning cycles; secondly, based on the phased array radar reflectivity factor data Z_par at different times of the target 3D grid output in step (2), the phased array radar reflectivity factor dataset corresponding to the time series ts is extracted and denoted as Z_par_ts; then, for the area determined to be of the light rain type in step (3), the phased array radar reflectivity factor dataset corresponding to the area is extracted from Z_par_ts and denoted as Z_par_ts_light. Then, the area average reflectivity factor sequence of the light rain type area corresponding to the time series ts is extracted from Z_par_ts_light using the spatial averaging method and denoted as Z_par_ts_l. ight_m; then, based on the ground automatic station observation data R0_sta collected in step (1) at different times, the ground automatic station rainfall rate observation data on the time series ts is generated by using the linear interpolation method in time, and is denoted as R0_sta_ts; then, for the area determined to be light rain type in step (3), the ground automatic station observation dataset corresponding to the area is extracted from R0_sta_ts, and is denoted as R0_sta_ts_light. Then, the regional average rainfall rate sequence of the light rain type area corresponding to the time series ts is extracted by R0_sta_ts_light according to the spatial averaging method, and is denoted as R0_sta_ts_ligh_m; then, a time lag sequence time_delay=m×delt_T is set with a step size of delt_T=1 minute, where m=0,1,2,…,M-1, and M is the total number of lags. For a certain lag value τ in the time lag sequence, select the sequence corresponding to the time lag value τ from the R0_sta_ts_ligh_m sequence, denoted as R0_sta_ts_ligh_m1, and calculate the correlation coefficient between Z_par_ts_light_m and R0_sta_ts_ligh_m1.The lag value τ is updated sequentially with all M values ​​in time_delay, and this process is repeated. The correlation coefficient between Z_par_ts_light_m and the rainfall rate data from the automatic weather station after M lag values ​​can be calculated, forming a correlation coefficient sequence, denoted as CC. Then, the maximum value is identified in the CC sequence, and the time lag value corresponding to the maximum value is extracted. This is the optimal lag time difference between the automatic weather station observation data and the phased array radar observation, denoted as time_delay_light. Finally, Z_par_ts_light is delayed by time_delay_light to form a phased array radar reflectivity factor dataset that is time-aligned with the automatic weather station observation data R0_sta_ts_light in the light rain area, denoted as Z1_par_ts_light. Applying the above processing steps to the areas of moderate rain, heavy rain, and torrential rain respectively, the optimal time lag (time_delay_middle, time_delay_heavy, time_delay_storm) corresponding to the areas of moderate rain, heavy rain, and torrential rain can be calculated. This leads to the generation of ground automatic station observation datasets (R0_sta_ts_middle, R0_sta_ts_heavy, R0_sta_ts_storm) and a phased array radar reflectivity factor dataset (Z1_par_ts_middle, Z1_par_ts_heavy, Z1_par_ts_storm) that are time-aligned with them. Further, the radar reflectivity factor datasets [Z1_par_ts_light, Z1_par_ts_middle, Z1_par_ts_heavy, Z1_par_ts_storm] under different precipitation intensities are merged into a unified radar reflectivity factor dataset Z1_par_ts, and the corresponding automatic ground station observation data [R0_sta_ts_light, R0_sta_ts_middle, R0_sta_ts_heavy, R0_sta_ts_storm] are integrated into an automatic ground station rainfall dataset R1_sta_ts arranged by precipitation intensity. Based on the determined light rain, moderate rain, heavy rain, and rainstorm areas, combined with the phased array radar differential reflectivity factor dataset Zdr_par and differential phase shift rate dataset KDP_par of the target three-dimensional grid at different times output in step (2), the differential reflectivity factor and differential phase shift rate of the corresponding rainfall areas are extracted respectively.Using the optimal time lags (time_delay_middle, time_delay_heavy, time_delay_storm) corresponding to the acquired moderate rain, heavy rain, and storm areas, the differential reflectivity factor and differential phase shift rate of each rainfall area are delayed. Then, the processed data are merged and unified to finally generate a differential reflectivity factor dataset Zdr1_par_ts and a differential phase shift rate dataset KDP1_par_ts that are time-aligned and integrated with the ground automatic station observation dataset. Among them, R0_sta_ts_middle, R0_sta_ts_heavy, and R0_sta_ts_storm are the ground automatic station observation data in the moderate rain, heavy rain, and storm areas, respectively. Z1_par_ts_middle, Z1_par_ts_heavy, and Z1_par_ts_storm are the phased array radar reflectivity factor datasets in the moderate rain, heavy rain, and storm areas that are time-aligned with the ground automatic station observation dataset.

[0078] Step (5): As Figure 5As shown, firstly, a station G is selected from the encrypted ground automatic station observation data collected in step (1), with its geographical coordinates being (x_R0_G, y_R0_G). Then, combined with Z1_par_ts output in step (4), the nearest neighbor spatial matching method is used to find the phased array radar echo grid point P directly above station G, and its position information on the horizontal plane is recorded as (x_par_P, y_par_P). Secondly, with point P as the spatial center point, the horizontal search window for inclined precipitation is set to N_h×N_h=7×7, where N_h is the width of the window size, and the number of vertical matching layers for inclined precipitation is N_v=10. The adjacent airspace above station G is determined to be Vol_space=N_h×N_h×N_v. After that, according to the three-dimensional wind field dataset (U, V, ...) output in step (2), ... The average values ​​of the U and V components at the i-th height layer are calculated using the Z1_par_ts output from W) and step (4), combined with the range of the adjacent spatial domain Vol_space. These average values ​​are denoted as U_i and V_i, respectively. The average value of the reflectivity factor data at the i-th height layer is also denoted as Z_i, where i = 1, 2, …, N_v. Then, Z_i is input into the reflectivity factor and raindrop terminal velocity relationship model. The model is run to generate the terminal velocity of the raindrop at the i-th height layer, denoted as Vt_i. Then, based on delt_z set in step (2), the height of the i-th height layer grid is obtained as Height_i = i × delt_z. Combined with the obtained Vt_i, the time difference between the raindrops falling to the ground at the i-th height layer is delta_t_i = Height_i / Vt_i; then, based on the obtained U_i and V_i, calculate the distances of raindrop drift in the X and Y directions within the time difference delta_t_i, respectively: delta_X_i=U_i×delta_t_i, delta_Y_i=V_i×delta_t_i; then, based on the average reflectance factor data Z_i of all N_v height layers, calculate the weighted coefficient of reflectance factor at the i-th height layer: apha_i=Z_i / sum(Z_i), where sum() is the summation function; then, based on the obtained apha_i, delta_X_i, and delta_Y_i, calculate the weighted average distances of raindrop drift in the X and Y directions of the adjacent airspace above point G of the ground automatic station: delta_X=sum(delta_X_i×apha_i), delta_Y=sum(delta_Y_i×apha_i);Finally, based on delta_X and delta_Y, and combined with the horizontal position information (x_par_P, y_par_P) of the nearest phased array radar echo grid point P directly above the ground automatic station G, the position of point P is adjusted to generate a phased array radar echo grid point Q that precisely matches the sky above the ground automatic station G, with horizontal position information of (x_par_P+delta_X, y_par_P+delta_Y). Applying this processing flow to all ground automatic stations generates a phased array radar echo dataset that is time-aligned and spatially matched with the ground automatic station observation dataset R1_sta_ts, including reflectivity factor, differential reflectivity factor, and differential phase shift rate datasets, denoted as Z1_par_ts_match, Zdr1_par_ts_match, and KDP1_par_ts_match, respectively.

[0079] Step (6): As Figure 6As shown, firstly, based on the target three-dimensional wind field dataset (U, V, W) output in step (2), the horizontal wind field (U_low, V_low) at the lowest altitude layer is extracted, and the near-surface wind direction W_dir=atan2(U_low,V_low) is calculated, where atan2() returns the clockwise angle of the wind direction relative to due north; secondly, based on the slope aspect data ASPECT collected in step (1), the angle delta_angle=ASPECT-W_dir between the slope aspect and the wind direction is calculated, and then the windward slope factor F_wind=max(0, cos(delta_angle×π / 180)), where max() is the maximum value function; then, find the maximum value of the slope data SLOPE collected in step (1), denoted as max_SLOPE, and calculate the slope factor F_slope=SLOPE / max_SLOPE; then, according to the digital elevation data DEM collected in step (1), set the neighborhood search radius R_DEM=9 kilometers, calculate the average digital elevation within the search radius of each DEM grid point with radius R_DEM as the center, and then calculate the elevation difference of each grid point relative to the average digital elevation within its respective search radius, denoted as delta_DEM, and use the hyperbolic tangent function to normalize delta_DEM to generate F_DEM=0.5×(1+tanh(delta_DEM / DEM0)), where tanh() is the hyperbolic tangent function and DEM0 is the scale parameter. Then, based on the slope type data TPI collected in step (1), the slope weight coefficient W_TPI is generated according to the slope type mapping weight rule; finally, based on the generated F_wind, F_slope, F_DEM and W_TPI, the basic terrain lift index terrain_lift=W_TPI×F_wind×F_slope×F_DEM is generated, where U_low is the eastward component of the horizontal wind speed at the lowest height layer and V_low is the northward component of the horizontal wind speed at the lowest height layer.

[0080] Furthermore, the slope type mapping weight rule is as follows: when TPI=0, W_TPI=0.2; when TPI=1, W_TPI=0.7; when TPI=2, W_TPI=1.0;

[0081] Step (7): As Figure 7As shown, based on Z1_par_ts_match, Zdr1_par_ts_match, and KDP1_par_ts_match output in step (5), the joint reflectivity and phase feature F1 = Z1_par_ts_match × KDP1_par_ts_match is first constructed, and the particle size distribution index feature F2 = Zdr1_par_ts_match / Z1_par_ts_match is constructed; then, the nonlinear radar precipitation potential vector R_radar = [Z1_par_ts_match, Zdr1_par_ts_match, KDP1_par_ts_match, F1, F2] is constructed; then, based on the DEM, SLOPE, ASPEC, and TPI collected in step (1), First, based on the nearest neighbor principle, extract the digital elevation, slope, aspect, slope position and basic terrain lift index datasets that spatially match the R1_sta_ts dataset output in step (4), and denote them as DEM1, SLOPE1, ASPEC1, TPI1 and terrain_lift1, respectively; second, construct the slope and aspect dynamic features T_d=[SLOPE1, sin(ASPECT1), cos(ASPECT1)], construct the terrain morphology and lift features T_m=[TPI1, terrain_lift1], and then construct the terrain feature vector T_terrain=[DEM1, T_d, T_m];

[0082] Step (8): As Figure 8As shown, based on R1_sta_ts output in step (4) and R_radar and T_terrain output in step (7), an end-to-end radar terrain coupling physical constraint model MODEL consisting of three sub-models is first constructed. Sub-model 1 is composed of a multilayer perceptron model. The input and output datasets of sub-model 1 during training are R_radar and R1_sta_ts, respectively. Sub-model 1 is trained, and the output after training is the radar precipitation potential P_radar ignoring the influence of terrain. Sub-model 2 is composed of an XGBoost model. The input and output datasets of sub-model 2 during training are T_terrain and R1_sta_ts, respectively. Sub-model 2 is trained... Training is performed, and the output after training is the precipitation terrain modulation factor M_terrain; Sub-model 3 is composed of physical constraint function mapping, and the output of sub-model 3 is R_out=P_radar×exp(M_terrain), where exp() is the logarithmic function, and the loss function of the entire model MODEL is set to the mean square error loss function, and the model MODEL is trained; Finally, after the model MODEL is trained, the R1_sta_ts output in step (4) and the R_radar and T_terrain output in step (7) are input into the trained model MODEL again, and the model is run, so as to realize the estimation of ground precipitation intensity based on phased array radar echo and geographic factors.

[0083] 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 surface precipitation intensity based on phased array radar echoes and geographic factors, characterized in that: The estimation method includes: S1. Collect phased array radar echo data, radar operation status data, ground automatic station observation data, and geographic factor data. Input the radar operation status data into the radar beam center approximate height estimation model and map projection model. Combine the results with the radar echo data to construct a target 3D mesh. Process the points in the 3D mesh to generate a 3D wind field dataset for the target 3D mesh at different times, as well as a phased array radar reflectivity factor dataset, differential reflectivity factor dataset, and differential phase shift rate dataset. After normalization, generate a normalized reflectivity factor dataset, differential reflectivity factor dataset, and differential phase shift rate dataset. Combine the local clustering model and local density calculation model to divide the region into several rainfall types, including light rain, moderate rain, heavy rain, and rainstorm. S2. Based on the collected radar operation status data, ground automatic station observation data, reflectivity factor dataset, and determined rainfall type areas, obtain the regional average reflectivity factor sequence and regional average rainfall sequence for the rainfall type areas. Based on the set time lag sequence, regional average reflectivity factor sequence, and regional average rainfall sequence, obtain the optimal lag time difference between the ground automatic station observation data and the phased array radar observation. Delay the reflectivity factor dataset by the optimal lag time difference to obtain a unified radar reflectivity factor dataset for several rainfall types. Combine the differential reflectivity factor dataset and differential phase shift rate dataset to generate a time-aligned and integrated differential reflectivity factor dataset and differential phase shift rate dataset that are consistent with the ground automatic station observation dataset. S3. Based on the collected ground automatic station observation data, three-dimensional wind field dataset, and merged unified radar reflectivity factor dataset, the average level drift distance of raindrops after falling to the ground is calculated by combining the radar reflectivity factor and raindrop terminal velocity relationship model, and then a phased array radar echo dataset is generated after spatiotemporal alignment and spatial matching with the ground automatic station observation dataset. S4. Generate a basic topographic lifting index based on the three-dimensional wind field dataset and the collected geographic factor data. Construct a topographic feature vector based on the data output from S3, the collected geographic factor data, and the basic topographic lifting index. S5. Construct a radar-terrain coupled physical constraint model based on the ground automatic station observation dataset and the data output from S4, and train it. Finally, use this model to estimate the ground precipitation intensity.

2. The method for estimating surface precipitation intensity based on phased array radar echo and geographic factors according to claim 1, characterized in that: The phased array radar echo data includes: radar reflectivity factor data Z0_par, differential reflectivity factor data Zdr0_par, differential phase shift rate data KDP0_par, and three-dimensional wind field data (U0, V0, W0) obtained by four-dimensional variational assimilation and inversion at different times in polar coordinate spatial locations (r0, phi0, seat0). Here, r0 is the radial distance, phi0 is the azimuth, seat0 is the antenna elevation angle, and the three-dimensional wind field data refers to (U0, V0, W0), where U0 is the eastward component of the horizontal wind speed, V0 is the northward component of the horizontal wind speed, and W0 is the vertical wind speed. 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; the radar volume scanning period delt_T; and the start time t0 of the radar scan. The ground automatic station observation data includes: the rainfall rate R0_sta observed at different times, whose spatial geographic coordinates are (x_R0, y_R0), where x_R0 and y_R0 are the position information in the east-pointing X-axis and north-pointing Y-axis directions, respectively; The geographic factor data includes: Digital Elevation Data (DEM), Slope Data (SLOPE), Aspect Data (ASPECT), and Slope Type Data (TPI).

3. The method for estimating surface precipitation intensity based on phased array radar echo and geographic factors according to claim 2, characterized in that: The process of inputting radar operating status data into the radar beam center approximate height estimation model and map projection model, and combining the results with radar echo data to construct a target 3D mesh, and processing the points in the 3D mesh to obtain a 3D wind field dataset for generating the target 3D mesh at different times, as well as a phased array radar reflectivity factor dataset, a differential reflectivity factor dataset, and a differential phase shift rate dataset, specifically includes: A1. Input the collected radar site spatial location information P_par, (r0, phi0, sita0), delt_R_par, and delt_sita_par into the radar beam center approximate height estimation model, run the model, and output the radar beam center approximate height z0. A2. Input the collected spatial location information P_par of the radar station into the map projection model, run the model, and output the planar projection coordinates (X_par, Y_par) of the radar station. Based on the collected (r0, phi0, sita0), delt_R_par, and delt_sita_par, construct a ray with an azimuth angle of phi0 in the planar coordinates. Take the horizontal distance as s0 = r0 × cos(sita0). Then the projection point of the ray centerline on the horizontal plane is x0 = X_par + s0 × sin(phi0) and y0 = Y_par + s0 × cos(phi0). Combine the point (x0, y0) with the acquired z0 to convert the polar coordinates of the radar echo data into three-dimensional spatial coordinates (x0, y0, z0). A3. Set the spatial resolution of the horizontal grid in the X and Y directions to delt_x and delt_y respectively, and set the resolution in the vertical height to delt_z, thereby constructing the target 3D grid. For a point M in the target 3D grid, its spatial coordinates are (x_m, y_m, z_m). Calculate the radial distance r_m and azimuth angle phi_m from this point to the radar station's planar projection coordinates (X_par, Y_par). For any point I in (x0, y0, z0), its coordinates are (x_i, y_i, z_i). Based on the collected radar echo data, extract the reflectivity factor data and 3D wind field data corresponding to point I, denoted as Z0_par_i and (U0_i, V0_i, ..., ...) respectively. W0_i), where U0_i is the eastward component of the horizontal wind speed at point I, V0_i is the northward component of the horizontal wind speed at point I, and W0_i is the vertical wind speed at point I. Calculate the coordinate difference between point I and point M, denoted as (delta_x_i, delta_y_i, delta_z_i). Then decompose delta_x_i and delta_y_i into the transverse direction along the radial direction, which is consistent with the beam broadening direction, denoted as (delta_r_i, delta_b_i). A4. Set the radial, lateral, and vertical weights to apha_r, apha_b, and apha_z, respectively. Calculate the anisotropic distance at point I as distance_i = apha_r × delta_r_i. 2 +apha_b×delta_b_i 2 +apha_z×delta_z_i 2 Then, the inverse distance weighting index is calculated as W_i=1 / distance_i. The reflectivity factor Z_par_m and the three-dimensional wind field data (U0_i_m, V0_i_m, W0_i_m) at point M in the target three-dimensional grid are calculated using the weighted average method. U0_i_m is the eastward component of the horizontal wind speed at point M, V0_i_m is the northward component of the horizontal wind speed at point M, and W0_i_m is the vertical wind speed at point M. A5. Apply steps A3 and A4 to all points in the target 3D mesh at different times to generate phased array radar reflectivity factor dataset Z_par and 3D wind field dataset (U, V, W) of the target 3D mesh at different times, where U, V and W are the eastward component of the horizontal wind speed, the northward component of the horizontal wind speed and the vertical wind speed of the target 3D mesh at different times, respectively. A6. Replace the reflectivity factor datasets in steps A2-A5 with the differential reflectivity factor dataset Zdr0_par and the differential phase shift rate dataset KDP0_par in sequence to generate phased array radar differential reflectivity factor dataset Zdr_par and differential phase shift rate dataset KDP_par for the target 3D mesh at different times.

4. The method for estimating surface precipitation intensity based on phased array radar echo and geographic factors according to claim 3, characterized in that: The process of generating normalized reflectivity factor datasets, differential reflectivity factor datasets, and differential phase shift rate datasets, and then dividing the regions into several rainfall types using local clustering models and local density calculation models, specifically includes: Find the maximum value of the generated reflectivity factor dataset Z_par, and divide each reflectivity factor data by the maximum reflectivity factor to generate normalized reflectivity factor data. Apply the same processing method to the differential reflectivity factor dataset and the differential phase shift rate dataset to generate normalized differential reflectivity factor dataset and differential phase shift rate dataset. The normalized reflectance factor, differential reflectance factor, and differential phase shift rate datasets are input together into a density-based local clustering model. The initial search radius r_orig around the data points and the minimum point value MinPts of the cluster are set. The local clustering model is run to generate regions with different precipitation intensity types that are initially divided. The average reflectance factor of each region is calculated. Based on the relationship between the average reflectance factor and the precipitation intensity classification, each precipitation region is initially labeled as a light rain, moderate rain, heavy rain, or rainstorm region. The generated reflectivity factor dataset Z_par is input into the local density calculation model. The local density calculation model is run to generate the local density value of each reflectivity factor data point. Based on the local density value of each reflectivity factor data point, the local density data of the initially marked light rain, moderate rain, heavy rain and rainstorm areas are statistically analyzed to obtain the average local density ρ of each precipitation area. Based on the optimization calculation formula of the search radius based on local density, the optimized search radii r_opt of light rain, moderate rain, heavy rain and rainstorm areas are obtained respectively. Then, the initial search radius r_orig around the data points set in the local clustering model is updated to r_opt, and the local clustering model is run again. The output is the area that is accurately divided into light rain, moderate rain, heavy rain and rainstorm types.

5. The method for estimating surface precipitation intensity based on phased array radar echo and geographic factors according to claim 4, characterized in that: S2 specifically includes the following: B1. Based on the collected delt_T and t0, construct the phased array radar observation time series ts=t0+k×delt_T, k=0,1,2,…,K, where K is the total number of radar scanning cycles. Based on the phased array radar reflectivity factor dataset Z_par at the target 3D grid at different output times, extract the phased array radar reflectivity factor dataset corresponding to the time series ts, denoted as Z_par_ts. For areas identified as light rain type, extract the phased array radar reflectivity factor dataset corresponding to the area from Z_par_ts, denoted as Z_par_ts_light. Extract the regional average reflectivity factor sequence of the light rain type area corresponding to the time series ts from Z_par_ts_light using the spatial averaging method, denoted as Z_par_ts_light_m. B2. Based on the collected automatic weather station observation data R0_sta at different times, a linear interpolation method is used in time to generate automatic weather station rainfall rate observation data on the time series ts, denoted as R0_sta_ts. For areas identified as light rain, the corresponding automatic weather station observation dataset for that area is extracted from R0_sta_ts, denoted as R0_sta_ts_light. Then, the regional average rainfall rate sequence for the light rain type area corresponding to the time series ts is extracted from R0_sta_ts_light using the spatial averaging method, denoted as R0_sta_ts_ligh_m. B3. Set a time lag sequence with step size delt_T: time_delay=m×delt_T, where m=0,1,2,…,M-1, and M is the total number of lags. For a certain lag value τ in the time lag sequence, select the sequence corresponding to the time lag value τ from the R0_sta_ts_ligh_m sequence, and denote it as R0_sta_ts_ligh_m1. Calculate the correlation coefficient between Z_par_ts_light_m and R0_sta_ts_ligh_m1. B4. Update the lag value τ sequentially with all M values ​​in time_delay, and repeat step B3 to calculate the correlation coefficient between Z_par_ts_light_m and the rainfall rate data of the ground automatic station after M lag values, forming a correlation coefficient sequence, denoted as CC. Identify the maximum value in the CC sequence and extract the time lag value corresponding to the maximum value, which is the optimal lag time difference between the ground automatic station observation data and the phased array radar observation, denoted as time_delay_light. Delay Z_par_ts_light by time_delay_light to form a phased array radar reflectivity factor dataset that is time-aligned with the ground automatic station observation data R0_sta_ts_light in the light rain area, denoted as Z1_par_ts_light. B5. Applying steps B1-B4 to the moderate rain, heavy rain, and storm areas respectively, the optimal time lag differences (time_delay_middle, time_delay_heavy, time_delay_storm) corresponding to these areas can be calculated. This leads to the generation of ground automatic station observation datasets (R0_sta_ts_middle, R0_sta_ts_heavy, R0_sta_ts_storm) and a time-aligned phased array radar reflectivity factor dataset (Z1_par_ts_mi). ddle, Z1_par_ts_heavy, Z1_par_ts_storm), where R0_sta_ts_middle, R0_sta_ts_heavy, and R0_sta_ts_storm are the ground automatic station observation data in the moderate rain, heavy rain, and rainstorm areas, respectively, and Z1_par_ts_middle, Z1_par_ts_heavy, and Z1_par_ts_storm are the phased array radar reflectivity factor datasets that are time-aligned with the ground automatic station observation datasets in the moderate rain, heavy rain, and rainstorm areas, respectively; B6. Merge the radar reflectivity factor datasets [Z1_par_ts_light, Z1_par_ts_middle, Z1_par_ts_heavy, Z1_par_ts_storm] under different precipitation intensities into a unified radar reflectivity factor dataset Z1_par_ts, and integrate the corresponding automatic weather station observation data [R0_sta_ts_light, R0_sta_ts_middle, R0_sta_ts_heavy, R0_sta_ts_storm] into an automatic weather station precipitation dataset R1_sta_ts arranged by precipitation intensity; B7. Based on the identified areas of light rain, moderate rain, heavy rain, and torrential rain, and combined with the phased array radar differential reflectivity factor dataset Zdr_par and differential phase shift rate dataset KDP_par of the target 3D mesh at different times, extract the differential reflectivity factor and differential phase shift rate of the corresponding rainfall area respectively. B8. Using the optimal time lags (time_delay_middle, time_delay_heavy, time_delay_storm) corresponding to the acquired moderate rain, heavy rain, and storm areas, the differential reflectivity factor and differential phase shift rate of each rainfall area are delayed. Then, the processed data are merged and unified to finally generate a differential reflectivity factor dataset Zdr1_par_ts and a differential phase shift rate dataset KDP1_par_ts that are time-aligned with and integrated with the ground automatic station observation dataset.

6. The method for estimating surface precipitation intensity based on phased array radar echo and geographic factors according to claim 5, characterized in that: S3 specifically includes the following: C1. Select a station G from the collected ground automatic station observation data. Its geographical coordinates are (x_R0_G, y_R0_G). Combine the output Z1_par_ts and use the nearest neighbor spatial matching method to find the phased array radar echo grid point P that is directly above station G. The position information on the horizontal plane is recorded as (x_par_P, y_par_P). C2. Taking point P as the spatial center point, set the horizontal search window for inclined precipitation to N_h×N_h, where N_h is the width of the window. The number of vertical matching layers for inclined precipitation is N_v. Determine the adjacent spatial domain above station G as Vol_space=N_h×N_h×N_v. Based on the three-dimensional wind field dataset (U, V, W) and Z1_par_ts, combined with the range of the adjacent spatial domain Vol_space, calculate the average values ​​of the U and V components at the i-th height layer, denoted as U_i and V_i, and the average value of the reflectivity factor data at the i-th height layer, denoted as Z_i, i=1, 2, …, N_v; C3. Input Z_i into the model relating reflectivity factor and raindrop terminal velocity, run the model, and generate the terminal velocity of the raindrop at the i-th height layer, denoted as Vt_i. Based on the set delt_z, obtain the height of the i-th height layer grid, Height_i = i × delt_z. Combined with the obtained Vt_i, obtain the time difference between the raindrop falling to the ground at the i-th height layer, delta_t_i = Height_i / Vt_i. Based on the obtained U_i and V_i, calculate the distances of raindrop drift in the X and Y directions within the time difference delta_t_i, respectively: delta_X_i = U_i × delta_t_i, delta_Y_i = V_i × delta_t_i. Based on the average reflectivity factor data Z_i of all N_v height layers, calculate the weighting coefficient of reflectivity factor at the i-th height layer, apha_i = Z_i / sum(Z_i), where sum() is the summation function. C4. Based on the obtained apha_i, delta_X_i, and delta_Y_i, calculate the weighted average distance of the raindrops drifting in the X and Y axes in the adjacent airspace above the automatic ground station G: delta_X = sum(delta_X_i × apha_i) and delta_Y = sum(delta_Y_i × apha_i). Based on delta_X and delta_Y, and combined with the determined position information (x_par_P, y_par_P) of the nearest phased array radar echo grid point P directly above the automatic ground station G on the horizontal plane, adjust the position of point P to generate a phased array radar echo grid point Q that precisely matches the airspace above the automatic ground station G, with the horizontal position information being (x_par_P + delta_X, y_par_P + delta_Y). C5. Apply steps C1-C4 to all automatic ground stations to generate phased array radar reflectivity factor, differential reflectivity factor, and differential phase shift rate datasets that have been time-aligned and spatially matched with the automatic ground station observation dataset R1_sta_ts, denoted as Z1_par_ts_match, Zdr1_par_ts_match, and KDP1_par_ts_match, respectively.

7. The method for estimating surface precipitation intensity based on phased array radar echo and geographic factors according to claim 5, characterized in that: The generation of the basic topographic lifting index based on the three-dimensional wind field dataset and collected geographic factor data specifically includes: Based on the target 3D wind field dataset (U, V, W), the horizontal wind field (U_low, V_low) of the lowest altitude layer is extracted. U_low is the eastward component of the horizontal wind speed at the lowest altitude layer, and V_low is the northward component of the horizontal wind speed at the lowest altitude layer. The near-surface wind direction W_dir=atan2(U_low, V_low) is calculated, where atan2() returns the clockwise angle of the wind direction relative to due north. Based on the collected slope aspect data ASPECT, the angle between the slope aspect and the wind direction is calculated as delta_angle = ASPECT - W_dir. Then, the windward slope factor F_wind = max(0, cos(delta_angle × π / 180)) is calculated, where max() is the function to find the maximum value. Find the maximum value of the collected slope data SLOPE, denoted as max_SLOPE, and calculate the slope factor F_slope = SLOPE / max_SLOPE. Based on the collected digital elevation data DEM, set the neighborhood search radius R_DEM, and calculate the average digital elevation within each search radius of R_DEM centered on each DEM grid point. Then, calculate the elevation difference of each grid point relative to the average digital elevation within its respective search radius, denoted as delta_DEM. Normalize delta_DEM using the hyperbolic tangent function to generate F_DEM = 0.5 × (1 + tanh(delta_DEM / DEM0)), where tanh() is the hyperbolic tangent function and DEM0 is the scale parameter. Based on the collected slope type data TPI, and according to the slope type mapping weight rule, the slope weight coefficient W_TPI is generated. Based on the generated F_wind, F_slope, F_DEM and W_TPI, the basic terrain lift index terrain_lift=W_TPI×F_wind×F_slope×F_DEM is generated.

8. The method for estimating surface precipitation intensity based on phased array radar echo and geographic factors according to claim 7, characterized in that: The construction of the terrain feature vector based on the data output by S3, the collected geographic factor data, and the basic terrain lifting index specifically includes the following: Based on Z1_par_ts_match, Zdr1_par_ts_match, and KDP1_par_ts_match, we first construct the joint reflectivity and phase feature F1 = Z1_par_ts_match × KDP1_par_ts_match, and then construct the particle size distribution index feature F2 = Zdr1_par_ts_match / Z1_par_ts_match, where Z1_par_ts_match, Zdr1_par_ts_match, and KDP1_par_ts_match are the phased array radar reflectivity factor, differential reflectivity factor, and differential phase shift rate datasets after time alignment and spatial matching, respectively. A nonlinear radar precipitation potential vector R_radar=[Z1_par_ts_match, Zdr1_par_ts_match,KDP1_par_ts_match, F1, F2] is constructed. Based on the collected DEM, SLOPE, ASPECT, TPI and terrain_lift datasets, the digital elevation, slope, aspect, slope position and basic terrain lifting index datasets that spatially match the R1_sta_ts dataset are extracted according to the nearest neighbor principle and denoted as DEM1, SLOPE1, ASPECT1, TPI1 and terrain_lift1, respectively. Construct slope and aspect dynamic features T_d=[SLOPE1, sin(ASPECT1), cos(ASPECT1)], construct terrain morphology and lift features T_m=[TPI1, terrain_lift1], and then construct terrain feature vector T_terrain=[DEM1, T_d, T_m].

9. The method for estimating surface precipitation intensity based on phased array radar echo and geographic factors according to claim 8, characterized in that: S5 specifically includes the following: Based on the output R1_sta_ts, as well as R_radar and T_terrain, an end-to-end radar terrain coupled physical constraint model MODEL consisting of three sub-models is first constructed. Sub-model 1 is composed of a multilayer perceptron model. The input dataset and output dataset of sub-model 1 during training are R_radar and R1_sta_ts, respectively. Sub-model 1 is trained, and the output after training is the radar precipitation potential P_radar ignoring the influence of terrain. Sub-model 2 is composed of an XGBoost model. The input and output datasets for training sub-model 2 are T_terrain and R1_sta_ts, respectively. After training sub-model 2, the output is the precipitation-terrain modulation factor M_terrain. Sub-model 3 is composed of physical constraint function mapping. The output of sub-model 3 is R_out=P_radar×exp(M_terrain), where exp() is the logarithmic function. The loss function of the radar terrain coupled physical constraint model MODEL is set to the mean square error loss function. The radar terrain coupled physical constraint model MODEL is trained. After the radar-terrain coupled physical constraint model MODEL is trained, R1_sta_ts, R_radar, and T_terrain are input into the trained radar-terrain coupled physical constraint model MODEL again, and the model is run to realize the estimation of ground precipitation intensity based on phased array radar echo and geographical factors.

Images