Mountainous terrain radar quantitative precipitation estimation method, system, storage medium, and apparatus

Abstract

The present application belongs to the technical field of quantitative precipitation estimation, and discloses a mountainous terrain radar quantitative precipitation estimation method, system, storage medium and device. The beam blocking rate is calculated by using high-resolution terrain data, and the adaptive elevation angle is selected, a multi-channel input data block containing polarization parameters, beam blocking rate, beam height from the ground and other terrain perception information is constructed, a physical constraint deep learning network introducing a blocking perception gating mechanism and a blocking filling sub-network is constructed, the beam blocking geometric optics principle and the reflectivity factor vertical variation law are embedded into the deep network training process through terrain blocking consistency constraint and vertical profile consistency constraint. The problems of low precipitation estimation accuracy and ineffective reconstruction of the blocking blind area under complex mountainous blocking conditions are solved, the precision and physical reasonableness of quantitative precipitation estimation in complex mountainous terrain areas are improved, high-precision, physically interpretable and strong generalization ability quantitative precipitation estimation is realized, and the present application has engineering practical value.

Description

Technical Field

[0001] This invention belongs to the field of quantitative precipitation estimation technology, specifically relating to a method, system, storage medium, and device for quantitative precipitation estimation using mountain topographic radar. Background Technology

[0002] Quantitative precipitation estimation is a key technology for inverting surface precipitation intensity using weather radar observation data. It plays an irreplaceable role in fields such as urban flooding early warning, mountain torrent and geological disaster monitoring, watershed flood control scheduling, and numerical weather prediction data assimilation.

[0003] X-band dual-polarization weather radar, with its narrow beam and high spatiotemporal resolution (up to 3 minutes for time and 30 meters for space), can finely characterize the internal structure and rapid evolution of severe convective systems, making it a core piece of equipment for refined precipitation monitoring and short-term early warning. In recent years, cities such as Beijing, Shenzhen, Wuhan, Chengdu, and Jiaxing have deployed dozens to hundreds of X-band dual-polarization phased array radars in network configurations.

[0004] However, in complex mountainous terrain areas such as the Longmenshan Mountains in northern Sichuan, X-band radar quantitative precipitation estimation faces a series of unique and severe challenges. Taking the Longmenshan mountain area in Mianyang as a typical example, this region has the following significant topographic features: The Longmenshan Mountains run northeast to southwest, about 500 kilometers long and 30-70 kilometers wide, forming a natural boundary between the Sichuan Basin and the western Sichuan Plateau; the area within Mianyang belongs to the northern section of the Longmenshan Mountains, diagonally traversing most of Beichuan, southeast of Pingwu, Jiangyou, and northwest of Anzhou, with an average elevation of 1000-3000 meters; the Beichuan-Yingxiu fault divides the Longmenshan Mountains into two parts: the front mountain and the back mountain. The front Longmenshan Mountains include hills, low mountains, and medium mountains, with mountain elevations generally below 2500 meters; the back Longmenshan Mountains are mainly medium and high mountains, with elevations mostly above 3500 meters, and the main peak, the Lion King Peak of Jiudingshan, reaches 4989 meters. Within a horizontal distance of only about 60 kilometers, the elevation rises sharply from 450 meters at the edge of the basin to 4989 meters at the main peak of the back mountain, a difference of more than 4500 meters, forming an extreme topographic gradient zone rarely seen in China. Meanwhile, the Longmenshan fault zone consists of three nearly parallel thrust faults (Anxian-Guanxian fault, Beichuan-Yingxiu fault, and Wenchuan-Maoxian fault). The dense fault zones result in fragmented terrain and numerous gullies. The eastern slope of this area is the famous Lutoushan rainstorm area, with an annual precipitation of over 1700 mm, making it a high-risk area for flash floods and geological disasters, thus creating an urgent engineering need for refined QPE (Quality, Physical, and Environmental Protection).

[0005] Under the aforementioned extreme terrain conditions, X-band radar quantitative precipitation estimation faces the following unique and severe challenges: First, the mountainous terrain causes severe radar beam obstruction. The Longmenshan area has a dramatically undulating terrain, with an elevation range of 1000-4800 meters or more. The high mountains and deep valleys create extensive obstruction of X-band radar deployed on the edge of the basin due to the mountain ridges at the front. Particularly towards the rear of the Longmenshan Mountains, the ridges at elevations of 1500-2500 meters almost completely block low-elevation beams, creating a continuous, large-scale blind spot in the rear mountain area. Beam obstruction leads to a significant attenuation of the radar's received echo power. This attenuation, combined with attenuation caused by precipitation, further exacerbates the existing attenuation problem of the X-band radar. In the obstructed areas, the reflectivity factor Z_H can be underestimated by 10-25 dBZ, directly resulting in a severely underestimation or even complete absence of precipitation data.

[0006] Second, the lack of near-surface observation information in the obstructed area is difficult to compensate for using traditional methods. Existing methods for addressing beam obstruction mainly fall into three categories: First, calculating the beam obstruction rate based on a digital elevation model and using linear or kriging interpolation with adjacent unobstructed azimuth data. This method relies solely on geometric optics principles and cannot reflect actual atmospheric refraction conditions and the true spatial distribution characteristics of precipitation echoes. Furthermore, in the Longmenshan Mountains' front-back mountain binary structure, the large-scale continuous obstructed area behind the front ridge cannot be effectively reconstructed using simple azimuth interpolation. Second, using the vertical profile of the reflectivity factor to correct the high-altitude unobstructed elevation angle data to near-surface height. This method assumes the same vertical... The distribution of reflectivity factor along the height within a straight air column satisfies the statistical average profile. However, the vertical structure of the precipitation system in the Longmenshan area exhibits high heterogeneity and spatiotemporal variability due to factors such as topographic forced uplift and valley wind circulation. Correcting the vertical profile of the fixed reflectivity factor will introduce significant bias. Thirdly, beam occlusion correction methods based on deep learning, such as using GAN, DenMerD, and CNN-M to fill echoes in occluded areas, have achieved certain results. However, these methods are purely data-driven and suffer from poor physical interpretability and insufficient cross-regional generalization ability.

[0007] Third, the mountainous terrain exacerbates the spatiotemporal heterogeneity of raindrop spectrum distribution. The formation mechanism of precipitation in the Longmenshan Mountains is complex, including stratiform cloud precipitation from large-scale weather systems, convective precipitation triggered by orographic forcing, and mixed precipitation caused by local circulation of valley winds. Due to the significant elevation differences between the front and rear Longmenshan Mountains, the raindrop spectrum distribution parameters of the same precipitation process vary considerably at different altitudes. The annual precipitation in the Lutoushan rainstorm area on the eastern slope of the Longmenshan Mountains exceeds 1700 mm, with frequent heavy rainfall, but the raindrop spectrum distribution parameters differ significantly between different precipitation types and at different altitudes. Traditional precipitation estimation methods based on fixed Zr relationships are completely inadequate for this highly heterogeneous raindrop spectrum distribution environment.

[0008] Fourth, X-band electromagnetic wave attenuation is compounded by the mountainous terrain. Besides attenuation caused by precipitation, the mountainous terrain of the Longmenshan Mountains further attenuates the radar signal. Especially in the area behind the ridgeline of the front mountain, the beam must pass through the obstructed area before reaching the rear mountain area, making the combined effect of obstruction and precipitation attenuation more pronounced. At long distances or at low elevation angles with severe obstruction, the radar receiving power may drop below the noise level, resulting in a significant reduction in effective detection range. Traditional attenuation correction methods suffer severe performance degradation in complex scenarios where obstruction and precipitation attenuation coexist.

[0009] Fifth, the physical consistency problem of existing deep learning methods is more pronounced in mountainous scenarios. Purely data-driven deep learning methods perform reasonably well in mountainous precipitation processes covered by training data, but their performance drops sharply in new areas, with new terrain conditions, and new precipitation types not covered by the training data. The Longmenshan area has complex and varied terrain. The dual structure of the front and back mountains means that within the same radar coverage area, there are multiple terrain conditions such as flat basins, gentle hills, and steep mountains. Radar stations are sparse, and training samples are limited. The generalization bottleneck of purely data-driven methods is particularly prominent, and the output results cannot guarantee that they meet the basic physical constraints of precipitation, which may produce physically unreasonable inversion results.

[0010] Currently, there is an urgent need to develop a method, system, storage medium, and device for quantitative precipitation estimation using mountain topographic radar. Summary of the Invention

[0011] The first technical problem to be solved by this invention is to provide a quantitative precipitation estimation method for mountain topographic radar. The second technical problem to be solved by this invention is to provide a quantitative precipitation estimation system for mountain topographic radar. The third technical problem to be solved by this invention is to provide a storage medium for the quantitative precipitation estimation system for mountain topographic radar. The fourth technical problem to be solved by this invention is to provide a quantitative precipitation estimation device for mountain topographic radar. These solutions address the technical problems of low precipitation estimation accuracy in complex mountainous terrain, exemplified by the Longmenshan Mountains, caused by severe beam obstruction due to the dual structure of the front and back mountains, strong spatiotemporal heterogeneity of raindrop spectrum distribution, superposition of precipitation and obstruction attenuation, poor physical interpretability of pure data-driven models, and insufficient generalization ability.

[0012] The quantitative precipitation estimation method for mountain topographic radar of the present invention includes the following steps: Step S1: Data collection and preprocessing; Acquire X-band dual-polarization radar data, ground rain gauge precipitation data, and high-resolution digital elevation model (DEM) data for the target mountainous terrain region. Extract four polarization parameters from the X-band dual-polarization radar data: horizontal reflectivity factor Z_H, differential reflectivity factor Z_DR, differential propagation phase shift rate K_DP, and correlation coefficient CC. Perform non-meteorological echo identification and removal, differential propagation phase shift Φ_DP defolding and quality control, and system bias correction. Simultaneously acquire minute-level precipitation data from synchronously observed ground rain gauges as tags, raindrop spectrometer observation data for raindrop spectral distribution DSD parameter verification, and high-resolution DEM data for the target region (spatial resolution better than 30 meters, preferably 12.5 meters or higher accuracy). Step S2: Terrain occlusion pre-analysis and adaptive hybrid elevation angle construction; Using high-resolution digital elevation model data, based on the radar site's latitude and longitude, antenna altitude, and parameters of each elevation layer, the ray tracing method is used to calculate the radar's beam obstruction rate at each azimuth-radial range library-elevation layer. Based on the beam obstruction rate, adaptive elevation angle selection is performed for each azimuth-radial range library to construct hybrid elevation polarization parameter data. Step S3: Create the training dataset; To ensure the diversity of training data, the dataset needs to cover samples from different seasons (flood season / non-flood season), different precipitation types (stratified cloud precipitation, convective precipitation, mixed precipitation), and different degrees of occlusion (no occlusion, partial occlusion, severe occlusion, complete occlusion). The training dataset is divided into training set, validation set, and test set. Centered on the location of the ground rain gauge, mixed elevation polarization parameter data within the horizontal sampling window are extracted. In addition to the polarization parameter channels, beam occlusion rate channel, beam height above ground channel, and terrain elevation gradient channel are added to construct a multi-channel input data block containing X-band dual-polarization radar data polarization parameters and terrain perception information. Step S4: Construct a neural network model; A physically constrained deep learning network is constructed, which includes a terrain occlusion awareness coding module, a spatial feature extraction module, a temporal context aggregation module, a dual-head output module, an occlusion filling sub-network, and a physical consistency verification module. The spatial feature extraction module introduces an occlusion awareness gating mechanism to adaptively adjust the feature contribution weights of regions with different occlusion levels based on the beam occlusion rate. The dual-head output module includes a first output head and a second output head in parallel. The first output head outputs the spatial distribution of minute-level precipitation intensity, and the second output head outputs the spatial distribution of raindrop spectrum DSD parameters. Step S5: Construct the total loss function; Total loss function L total as follows: ; in, For data-driven loss terms, For the physical consistency constraint term of raindrop spectral distribution, These are the constraint terms of the attenuation equation. These are non-negative constraint terms. For spatial smoothing constraints, For terrain occlusion consistency constraints, This is a vertical profile consistency constraint term. ~ These are the weighting coefficients for each constraint term; Step S6: Train the neural network model; The physical constraint deep learning network is trained using the training dataset, and the network parameters are optimized using gradient descent with the total loss function as the objective until the total loss function converges. Step S7: Deploy the neural network model; The trained physical constraint deep learning network is deployed in the business environment to receive radar volume scan data and digital elevation model (DEM) in real time. After preprocessing, terrain occlusion pre-analysis, and multi-channel data block construction, the data is input into the network for forward inference. The network outputs minute-level precipitation intensity R value for each radar grid point and performs time accumulation to generate a quantitative precipitation estimation grid point product with the corresponding time resolution.

[0013] Furthermore, the specific strategy for adaptive elevation angle selection described in step S2 is as follows: When the beam obstruction rate of the lowest available elevation angle in the azimuth-radial distance database is less than or equal to the first threshold T1, the observation data is used directly; where the value of the first threshold T1 ranges from 0.3 to 0.5. When the beam obstruction rate of the azimuth-radial distance coulomb is between the first threshold T1 and 0.1, it is in a partial obstruction state. After correction of the horizontal reflectivity factor Z_H, it is used, and the revised formula is: ,in, For beam blocking rate, The corrected horizontal reflectance. The measured value of the horizontal reflectivity factor before correction. This is the beam obstruction power correction factor; When the beam obstruction rate of the azimuth-radial distance library is greater than the first threshold T1, the search is performed at the next higher elevation angle. The highest available elevation angle with a beam obstruction rate of less than or equal to 0.1 is selected, and the highest available elevation angle is corrected using the vertical profile in subsequent steps and converted to the reference height. When the beam obstruction rate of the azimuth-radial distance library is greater than the first threshold T1 in all elevation layers, it is marked as a completely obstructed blind zone and reconstructed through the obstruction filling mechanism of a deep learning network. The beam obstruction rate is calculated using a Gaussian beam model. Specifically, under standard atmospheric or actual sounding refraction conditions, the beam's azimuth angle θ and elevation angle are calculated. The center height at a distance r is obtained by reading the terrain elevation profile along the beam propagation path from the digital elevation model (DEM). The power ratio of the beam blocked by the terrain is calculated based on the Gaussian distribution model of the antenna pattern. The beam blocking rate ranges from 0 to 1, where 0 represents no blocking and 1 represents complete blocking. The calculation of beam obstruction rate incorporates an echo climate statistical correction factor. Specifically, for each azimuth-radial distance database, the frequency of clear-sky echoes during non-precipitation periods is statistically analyzed to construct the observation obstruction probability (OBP). When the difference between the beam obstruction rate (BBR) calculated by the digital elevation model (DEM) and the observation obstruction probability (OBP) exceeds a set threshold, a correction formula is used. Correction is performed, among which... λ For the credibility weight of the digital elevation model (DEM), The corrected beam obstruction rate, The beam obstruction rate calculated for the Digital Elevation Model (DEM) is used to compensate for non-topographic obstruction factors such as tall buildings and seasonal changes in vegetation that may not be reflected in the DEM data. When the accuracy of the DEM data does not meet the pre-set accuracy requirements or is difficult to obtain (such as earthquakes, landslides, or topographic changes caused by vegetation changes), the obstruction azimuth is identified by statistically analyzing the frequency of radar echoes over a long time series (more than 6 months). The effective minimum available elevation angle for each azimuth-radial distance database is determined using the echo-based climatological frequency method, and the obstruction information is dynamically updated.

[0014] Furthermore, the construction method of the multi-channel input data block mentioned in step S3 is as follows: A rectangular sampling window of M×N azimuth-radial distance libraries is selected in the horizontal dimension to form a three-dimensional input tensor, wherein the values ​​of M and N range from 16 to 64, preferably 32×32; The input channels include four polarization parameter channels and three terrain perception auxiliary channels, totaling seven channels. The information from these seven channels enables the deep learning network to simultaneously perceive the polarization characteristics of precipitation echoes and the impact of terrain on radar beam obstruction. The four polarization parameter channels are Z_H, Z_DR, K_DP, and CC, and the three terrain perception auxiliary channels are beam obstruction rate (BBR), beam height above ground (BH), and terrain elevation gradient (TG). The beam obstruction rate (BBR) records the beam obstruction rate of each spatial location after calculation and correction in step S2. The beam height above ground (BH) records the height of the radar beam center above ground at each spatial location, in km, reflecting the spatial relationship between the beam and the terrain. The terrain elevation gradient (TG) records the local gradient and aspect information of the terrain elevation at each spatial location, reflecting terrain features including valley bottoms, slopes, ridges, and mountain peaks. The 7 channels of three-dimensional data blocks from T consecutive time intervals are stacked according to the time dimension to form a four-dimensional spatiotemporal input sample with dimensions [T,7,M,N], where the value of T ranges from 3 to 12, preferably 6; each channel is standardized or normalized separately.

[0015] Furthermore, the specific structure of the terrain occlusion perception coding module described in step S4 is as follows: The information of beam obstruction rate (BBR) and beam height above ground (BH) is encoded into an obstruction-aware feature vector through two fully connected layers. This vector is then converted into an obstruction-aware feature map (O) with the same resolution as the spatial feature map by a spatial broadcast layer. The dimension of the obstruction-aware feature map O is M×N×D_o, where D_o is the encoding dimension and its value ranges from 16 to 32.

[0016] Furthermore, the specific operation of the occlusion perception gating mechanism described in step S4 is as follows: In each layer of the encoder of the spatial feature extraction module, the occlusion-aware feature map O generated by the terrain occlusion-aware coding module is sampled to the same spatial resolution as the current layer feature map F_l, and a gated weight map G_l is generated through the Sigmoid activation function. The current layer feature map F_l is then gated and weighted as follows: F_l_gated = F_l⊙G_l + F_l⊙(1-G_l)⊙β_l; where ⊙ denotes element-wise multiplication, β_l is the learnable baseline compensation factor for the occluded region, and F_l_gated is the gated and weighted feature map. The learnable weight matrix for gated convolution. This is an occlusion-aware feature map downsampled to the current layer resolution. This is a gated bias term; In occluded regions, the gating weights reduce the contribution of the original echo features and enhance the context compensation learned from the surrounding unoccluded regions; in unoccluded regions, the gating weights are close to 1, preserving the original echo features. The generation of the gating weight map G_l uses at least one of the following information sources: univariate gating based on beam obscuration ratio BBR, bivariate gating based on BBR and BH, trivariate gating based on beam obscuration ratio BBR, beam height above ground BH and terrain elevation gradient TG, or learnable gating learned from the digital elevation model DEM.

[0017] The spatial feature extraction module adopts an encoder-decoder structure. The encoder extracts multi-scale spatial features through multi-layer convolution, with each layer followed by a CBAM attention module and an occlusion-aware gating mechanism. The encoder-decoder structure uses at least one of the following: 3D ResNet, a CBAM attention-enhanced U-Net architecture, or a Swing Transformer module. The temporal context aggregation module performs temporal modeling on the spatial feature maps over T consecutive time intervals, outputting spatiotemporal joint features that incorporate historical information. It uses at least one of the following structures: a temporal convolutional network TCN, a bidirectional long short-term memory network BiLSTM, or a spatiotemporal attention module. The occlusion filling subnetwork is based on partial convolution or spatial attention mechanisms. For grid points marked as completely occluded blind areas, it takes the spatiotemporal joint features of the surrounding unoccluded areas as input, performs convolution operations on the effective unoccluded pixels through partial convolution, and propagates the effective information to the center of the blind area layer by layer through a mask update mechanism. The output is the filled precipitation intensity field and raindrop spectrum distribution DSD parameter field. The occlusion filling subnetwork adopts a U-Net structure decoder, which receives the feature maps after occlusion perception gating in each layer of the encoder through skip connections to ensure the continuity of information. The raindrop spectral distribution DSD parameters output by the second output head are described using a normalized Gamma distribution, which is N(D) = N_w·f(μ)·(D / D_m)^μ·exp[-(4+μ)·D / D_m], where N_w is the normalized intercept parameter, D is the raindrop diameter, D_m is the mass-weighted average diameter, μ is the shape parameter, and f(μ) is the normalization function; the second output head outputs respectively The spatial distribution of the three parameters, D_m and μ.

[0018] Furthermore, the terrain occlusion consistency constraint L_terrain mentioned in step S5 includes the following three sub-items: Sub-item 1 is the occlusion gradient consistency constraint: near the boundary of the occlusion region, the gradient of precipitation intensity R is statistically consistent with the topographic elevation gradient; the cosine similarity between the spatial gradient magnitude of precipitation intensity R and the normalized topographic elevation gradient magnitude is calculated, and a penalty term L_terrain_grad is applied to regions where the two are inconsistent. Where w(BBR) is a weighting function that increases with the beam obstruction rate, R is the estimated precipitation intensity output by the network, and H is the terrain elevation; to ensure that sub-item 1 plays a greater role near the obstruction boundary; the physical basis of sub-item 1 is: mountain precipitation is affected by topographic forced uplift, and there is a systematic spatial distribution pattern on the windward side of the mountain slope and the valley bottom, and there is a learnable physical correlation between the spatial gradient of precipitation intensity and the topographic gradient; Sub-item 2 is the joint consistency constraint for occlusion and attenuation: in the partially occluded area, the polarization parameters corrected by the digital elevation model DEM for occlusion correction and precipitation attenuation correction are constrained to be consistent with the equivalent polarization parameters inversely calculated from the precipitation intensity R output by the physically constrained deep learning network through the Z-R relationship; specifically, in the partially occluded area (0.1<BBR<T1), the estimated value of the precipitation intensity R shall satisfy that the equivalent value of the measured Z_H after DEM occlusion correction and precipitation attenuation correction is consistent with the Z_H value inversely calculated from the precipitation intensity R output by the physically constrained deep learning network through the Z-R relationship, that is: shall be consistent with Z_H(R) inversely calculated from R through the Z-R relationship, wherein A_H_path(R) is the path-integrated precipitation attenuation estimated from R; Sub-item 3 is the spatial continuity constraint for occlusion blind areas: the weighted spatial second-order difference between the precipitation intensity of the occlusion filling area and the adjacent non-occluded area is constrained to be minimized, and the minimized value of the weighted spatial second-order difference L_terrain_continuity=Σw_i·|R_i-2·R_j+R_k|², where i is the grid point inside the blind area, j and k are two adjacent grid points along the terrain gradient direction, w_i is a weight that decreases as the distance between the grid point i and the edge of the blind area increases, R_i is the estimated precipitation intensity value of the grid point i inside the blind area, R_j is the estimated precipitation intensity value of the adjacent grid point j, and R_k is the estimated precipitation intensity value of the adjacent grid point k; the weight w_i in sub-item 3 adopts a distance attenuation function, and the distance attenuation function is at least one of the following: exponential attenuation exp(-d_i / σ), Gaussian attenuation, and attenuation function based on terrain connectivity, where d_i is the distance from the grid point i to the nearest non-occluded grid point at the edge of the blind area, and σ is the attenuation scale; The vertical profile consistency constraint term L_VPR is a physical constraint term introduced for mountain terrain scenarios; the physical basis of this constraint term is: radar beam occlusion forces the network to use data at higher elevation angles for estimation, while the radar sampling volume at higher elevation angles is farther from the ground, and the vertical variation of the reflectivity factor Z_H will lead to precipitation estimation deviation (i.e., the VPR effect). Therefore, this constraint term ensures that the network learns the correct vertical variation rule of Z_H, so that the precipitation estimations retrieved at different elevation angles are consistent at the near-ground reference height; The calculation method of the vertical profile consistency constraint term L_VPR is: For the same horizontal position, forward inference is performed using low-elevation and high-elevation data obtained from the adaptive elevation angle selection results. The difference in precipitation intensity between the two inferences, L_VPR, is calculated as a constraint: L_VPR=(1 / N_pair)·Σ(R_low-R_high)²; where R_low is the precipitation intensity inferred using low-elevation data, R_high is the precipitation intensity inferred using high-elevation data, and N_pair is the effective logarithm of elevation angles. In unobstructed areas where the beam obstruction rate of both elevation angle data is less than or equal to 0.1, the differential L_VPR is applied to avoid obstruction noise affecting the accuracy of the constraint. The physical consistency constraint L_DSD is calculated as follows: Based on the raindrop spectral distribution DSD parameters output by the network, including the normalized intercept parameter N_w, the mass-weighted average diameter D_m, and the shape parameter μ, the equivalent polarization parameter prediction values ​​Z_H^pred, Z_DR^pred, and K_DP^pred are calculated forward using the T-matrix scattering model. These prediction values ​​are then compared with the radar-measured polarization parameters Z_H^obs, Z_DR^obs, and K_DP^obs. The differences are used as constraints. Constraints are applied to grid points in unobstructed regions where the beam obstruction rate is less than or equal to 0.1 to avoid polarization parameter distortion caused by obstruction affecting the reliability of the constraints. The weight coefficients λ1 to λ6 are selected using an adaptive scheduling strategy: in the early stage of training, λ5 and λ6 are selected with smaller values, ranging from 0.01 to 0.05, while λ1 to λ4 are selected with medium values, ranging from 0.1 to 0.5, so that data-driven terms and basic physical constraints dominate the fast convergence; in the middle and later stages of training, λ5 and λ6 are gradually increased to 0.5-1.5, so that terrain-related physical constraints play a dominant role.

[0019] Furthermore, the target mountainous terrain region has a dual structure of front mountain and back mountain. With the regional fault zone as the boundary, the front mountain region consists of hills, low mountains, and medium mountains with an altitude below 2,500 meters, while the back mountain region consists of medium mountains and high mountains with an altitude above 3,500 meters. The ridge of the front mountain partially obstructs the radar deployed on the edge of the basin, while the back mountain region forms a severe obstruction or a complete blind spot. Adaptive elevation angle selection and terrain occlusion consistency constraints employ differentiated threshold parameters and constraint weights for partially occluded areas in the front mountain and severely occluded areas in the back mountain. When S-band or C-band operational weather radars are deployed simultaneously in the study area, the polarization parameters of the S-band or C-band radars are added as an additional input channel to the multi-channel input data block; a cross-band feature alignment module is added to the network to learn the feature correspondence between S-band, C-band and X-band data at corresponding spatial locations through a spatial attention mechanism; a cross-band consistency constraint term L_crossband is added to the total loss function to constrain the consistency of precipitation intensity retrieved from different bands at the same horizontal location.

[0020] The mountain terrain radar quantitative precipitation estimation system of the present invention includes the following modules: The data acquisition and preprocessing module is used to acquire X-band dual-polarization radar base data, ground rain gauge data, and digital elevation model (DEM), extract four polarization parameters: horizontal reflectivity factor Z_H, differential reflectivity factor Z_DR, differential propagation phase shift rate K_DP, and correlation coefficient CC, and perform non-meteorological echo identification and removal, differential propagation phase shift Φ_DP defolding and quality control, and system deviation correction. The terrain obstruction pre-analysis and adaptive elevation angle selection module is used for the digital elevation model (DEM). It uses the ray tracing method to calculate the beam obstruction rate (BBR) of the radar at each azimuth-radial range-elevation layer, and performs adaptive elevation angle selection based on the beam obstruction rate to construct hybrid elevation polarization parameter data. The multi-channel input data construction module is used to construct standardized multi-channel input data blocks from the polarization parameter channel, beam obstruction rate channel, beam height above ground channel, and terrain elevation gradient channel, and stack them in the time dimension to form a four-dimensional spatiotemporal input sample. The physically constrained deep learning mountain QPE network module includes a terrain occlusion-aware encoding submodule, a spatial feature extraction submodule, a temporal context aggregation submodule, a dual-head output submodule, an occlusion filling subnetwork, and a physical consistency verification submodule. The spatial feature extraction submodule introduces an occlusion-aware gating mechanism to adaptively adjust the feature contribution weights of regions with different occlusion levels based on the beam occlusion rate. The dual-head output submodule includes a parallel precipitation intensity output head and a DSD parameter output head. The loss function calculation module is used to construct the total loss function, which includes the data-driven loss term L_data, the physical consistency constraint term L_DSD, the attenuation equation constraint term L_attenuation, the nonnegation constraint term L_nonneg, the spatial smoothness constraint term L_smooth, the terrain occlusion consistency constraint term L_terrain, and the vertical profile consistency constraint term L_VPR. Among them, L_terrain includes the occlusion gradient consistency constraint sub-term, the occlusion attenuation joint consistency constraint sub-term, and the occlusion blind zone spatial continuity constraint sub-term. The model training module is used to optimize network parameters based on training data with the total loss function as the objective. The QPE product generation module is used to receive radar data in real time and output quantitative precipitation estimation grid products. The quantitative precipitation estimation grid products include occlusion probability products, which indicate the beam occlusion rate (BBR) value of each grid point, and occlusion filling confidence products, which indicate the filling confidence score of each grid point for forecasters to refer to and judge.

[0021] The storage medium of the mountain terrain radar quantitative precipitation estimation system of the present invention is used to store the calculation program of the mountain terrain radar quantitative precipitation estimation method. The storage medium of the mountain terrain radar quantitative precipitation estimation system is installed on a computer and the calculation program of the mountain terrain radar quantitative precipitation estimation method is run by the computer.

[0022] The present invention relates to a mountain terrain radar quantitative precipitation estimation device, which is used to install a mountain terrain radar quantitative precipitation estimation system. The mountain terrain radar quantitative precipitation estimation device includes a memory and a processor. The memory stores a calculation program for the mountain terrain radar quantitative precipitation estimation method, and the processor executes the calculation program for the mountain terrain radar quantitative precipitation estimation method for field testing.

[0023] The mountain terrain radar quantitative precipitation estimation method, system, storage medium, and device of the present invention have the following characteristics: (1) Achieve organic integration of terrain occlusion physical constraints and deep learning nonlinear fitting. Through occlusion perception gating mechanism, terrain occlusion consistency constraints and vertical profile consistency constraints, the deep network, when performing precipitation inversion in occluded areas, not only relies on the statistical laws of training data, but also follows the geometric optical principle of beam occlusion and the vertical variation law of reflectivity factor, which fundamentally improves the reliability of precipitation estimation in occluded areas; (2) Simultaneous processing of precipitation attenuation and terrain occlusion attenuation. Existing attenuation correction methods such as ZPHI only address attenuation caused by precipitation, while digital elevation model (DEM) occlusion correction methods only address geometric occlusion caused by terrain. The two are used independently. This invention organically combines the two parts into a unified physical constraint deep learning framework, solving the problem of decoupling and coordinating attenuation sources. (3) Retain effective near-surface observation information. An adaptive hybrid elevation angle construction strategy oriented towards mountainous terrain is adopted. Instead of using a fixed elevation angle or simply taking the maximum value / median, the optimal available elevation angle is dynamically selected based on the BBR of each azimuth angle-radial distance library, and quantitative correction is performed when there is partial obstruction, so as to maximize the retention of effective near-surface observation information; (4) It has the function of intelligent filling of blind spots. Through the shading filling sub-network and spatial continuity constraints, it can intelligently reconstruct the precipitation distribution of completely shading areas from the surrounding effective data, thus filling the detection blind spots; (5) Improve the accuracy of precipitation estimation in mountainous terrain. In tests using actual data from the Longmenshan area in northern Sichuan, the method of this invention reduced the root mean square error (RMSE) of precipitation estimation in occluded areas by 35%-50% compared to pure data-driven deep learning methods; and reduced the RMSE by 25%-40% compared to the traditional digital elevation model (DEM+VPR) method. In severely occluded / blind areas in the back mountains (behind the ridge of the front mountains), the RMSE was reduced by 50%-60% compared to the traditional DEM+VPR method. The RMSE of the entire region (including occluded and unoccluded areas) was reduced by 20%-30% compared to the pure physical constraint deep learning method without terrain processing. (6) It has strong generalization ability and physical interpretability. Physical constraints make the mapping function learned by the network conform to the physical laws of beam propagation and the microphysical laws of precipitation. Even when the radar stations in the mountains are sparse and the training samples are limited, it can still maintain stable estimation performance under new terrain conditions and precipitation processes. The terrain occlusion perception feature map can be used as an intermediate interpretable representation, making the precipitation inversion process in the occluded area transparent and traceable.

[0024] This invention addresses the challenges of severe X-band radar beam obstruction, prominent hierarchical obstruction characteristics due to the front-mountain-back-mountain binary structure, and the superimposed attenuation of precipitation and obstruction in complex mountainous terrain, exemplified by the Longmenshan Mountains. It utilizes high-resolution topographic data to calculate the beam obstruction rate and adaptively selects the elevation angle, constructing a multi-channel input data block containing polarization parameters and topographic perception information such as beam obstruction rate and beam height above ground. A physically constrained deep learning network is constructed, incorporating an obstruction perception gating mechanism and an obstruction filling sub-network. Through topographic obstruction consistency constraints and vertical profile consistency constraints, the geometric optics principles of beam obstruction and the vertical variation law of reflectivity factor are embedded into the deep network training process. This invention solves the problems of low precipitation estimation accuracy and ineffective reconstruction of obstruction blind areas under complex mountainous terrain conditions, improving the accuracy and physical rationality of quantitative precipitation estimation in complex mountainous areas.

[0025] In summary, the quantitative precipitation estimation method, system, storage medium, and device for mountain terrain radar of the present invention can simultaneously handle precipitation attenuation and terrain occlusion attenuation. By integrating the geometric priors of digital elevation models, the physical constraints of radar polarization parameters, and the advantages of deep learning nonlinear fitting, the X-band radar quantitative precipitation estimation method achieves high-precision, physically interpretable, and highly generalizable quantitative precipitation estimation under complex mountain terrain conditions such as Longmenshan, and has practical engineering value. Attached Figure Description

[0026] The accompanying drawings are for illustrative purposes only and are not intended to limit the scope of this disclosure. It is obvious that the drawings described below are merely some embodiments of this disclosure, and those skilled in the art can obtain other drawings based on these drawings without any inventive effort. Furthermore, the same reference numerals denote the same parts throughout the drawings.

[0027] Figure 1 The overall flowchart of the mountain terrain radar quantitative precipitation estimation method is shown in the example. Figure 2 This is a schematic diagram illustrating the beamout rate calculation based on the digital elevation model (DEM) in the embodiment. Figure 3 This is a schematic diagram of the adaptive hybrid elevation angle construction strategy in the embodiment; Figure 4 A schematic diagram of the physically constrained deep learning mountain QPE network architecture provided for this embodiment; Figure 5 This is a detailed structural diagram of the occlusion perception gating mechanism in the embodiment; Figure 6 This is a schematic diagram of the structure of the multi-physical constraint loss function (including the mountain terrain constraint term) in the embodiment. Detailed Implementation

[0028] The present invention will now be described in detail with reference to the accompanying drawings and embodiments.

[0029] Example 1: This example is based on observation data from an X-band dual-polarization radar network deployed in the Longmenshan area of ​​Sichuan, data from 52 ground automatic weather stations, and high-resolution DEM data, to explain in detail the implementation process of the present invention.

[0030] The Longmenshan region in northern Sichuan is located in the northern segment of the Longmenshan fault zone, where the eastern edge of the Qinghai-Tibet Plateau transitions into the Sichuan Basin. Its latitude and longitude range is approximately 31.5°N–33.0°N and 103.5°E–105.5°E, encompassing most of Beichuan Qiang Autonomous County in Mianyang City, the northwestern part of Anzhou District, the northern part of Jiangyou City, and the southeastern part of Pingwu County. This region has the following typical topographic features: The dual structure of the front and back mountains: with the Beichuan-Yingxiu fault (central fault) as the boundary, the front mountain is mainly composed of hills and low mountains (mostly below 2500 meters above sea level), while the back mountain is composed of medium and high mountains (mostly above 3500 meters above sea level), forming a hierarchical shielding characteristic of "partial shielding by the front mountain and severe shielding or complete blind spots by the back mountain"; extreme topographic gradient: within a horizontal distance of 60 kilometers, the elevation rises sharply from about 470 meters to 4989 meters, with a height difference of more than 4500 meters, and three nearly parallel thrust faults result in fragmented terrain and numerous V-shaped canyons, creating an extremely complex shielding environment for radar beam propagation; frequent heavy rainfall: the annual precipitation on the eastern slope reaches more than 1700 mm, with a clear characteristic of increased precipitation on the windward slope of the front mountain and decreased precipitation on the leeward slope of the back mountain.

[0031] like Figure 1 As shown, the quantitative precipitation estimation method for mountain terrain radar in this embodiment includes the following steps: Step S1: Data collection and preprocessing; Radar baseline data for key precipitation events during two flood seasons (May-September) and non-flood seasons from May 2024 to September 2025 were acquired. Four polarization parameters, Z_H, Z_DR, Φ_DP, and CC, were extracted from the baseline data. K_DP was calculated after processing Φ_DP through sliding filtering and defolding.

[0032] Non-meteorological echo identification: Due to the complexity of mountainous terrain, ground object echoes (especially permanent mountain echoes) are exceptionally prominent. This embodiment employs a two-stage identification strategy: the first stage uses a fuzzy logic algorithm based on CC texture, Z_DR texture, and Φ_DP texture for initial screening; the second stage uses the U-Net semantic segmentation network to refine the initial screening results. The U-Net semantic segmentation network takes a four-channel polarization parametric image as input and outputs pixel-by-pixel binary classification labels (meteorological / non-meteorological echoes). The training data comes from manually labeled typical clear sky and precipitation processes.

[0033] Digital elevation model (DEM) processing: Acquire ALOS PALSAR 12.5-meter resolution DEM data, resample to 30-meter resolution to match radar data grid, and generate terrain elevation values ​​corresponding to each radar azimuth-radial distance database through bilinear interpolation.

[0034] Step S2: Terrain occlusion pre-analysis and adaptive hybrid elevation angle construction; For each radar, the beam obstruction rate (BBR) is calculated for each azimuth angle (0°-360°, step 1.0°), each azimuth angle-radial distance library (0-105km, step 30m, total 3500 libraries), and each elevation layer (11 layers in total).

[0035] like Figure 2As shown, the beam obstruction rate (BBR) is calculated using a modified Gaussian beam model. Under standard atmospheric refraction conditions, the beam is calculated at azimuth angle θ and elevation angle θ. The center height h_c(θ, at a distance r) For this sampling volume, the terrain elevation profile along the beam propagation path is read from the DEM data, and the percentage of beam power blocked by the terrain is calculated: ; Where α(θ, ,r) is the difference in elevation angle between the beam center and the highest point of the obstructing terrain, and f(α) is the antenna pattern function (Gaussian form).

[0036] An echo climate statistical correction factor is introduced. For clear-sky periods from May 2024 to September 2025 (selecting periods with more than 30 consecutive days without precipitation records), the observation obstruction probability (OBP) for each azimuth-radial distance database is calculated. When |BBR_DEM-OBP|>0.15, the correction formula BBR=0.7·BBR_DEM+0.3·OBP is used.

[0037] like Figure 3 As shown, adaptive hybrid elevation angle selection is performed: the threshold T1 is set to 0.4. For each azimuth-radial distance library: If the lowest available elevation angle BBR ≤ 0.1, use that elevation angle; If the BBR is between 0.1 and 0.4, use this elevation angle but apply... Correction; If BBR ≥ 0.4, proceed to the next higher elevation angle and repeat the judgment. If the BBR at all 11 elevation angles is ≥0.4, it is marked as a completely obscured blind spot.

[0038] Taking Radar A (2850 meters above sea level, on the ridge of the Houlongmenshan Mountains) as an example, at an elevation angle of 0.5°, the BBR of its western and northwestern azimuth sectors (approximately 120°) is between 0.6 and 1.0, forming a large area of ​​blind spots; the BBR of its eastern and southern sides is close to 0, representing effective detection areas; after adaptive elevation angle selection, the mixed elevation angle of the blind spots is automatically increased to 1.5°-6.0°. Radar C (980 meters above sea level, on the edge of the basin) shows the opposite: the BBR of its northwestern side (approximately 150° azimuth sector towards the Houlongmenshan Mountains) is between 0.4 and 1.0, forming a continuous blind spot in the rear mountain area. The four radars, through complementary observations from stations at different elevations, form a three-dimensional observation network covering both the front and rear mountains.

[0039] Step S3: Create the training dataset; The training dataset is divided into a training set, a validation set and a test set; taking the position of each ground rain gauge as the center, a horizontal sampling window of 32×32 bins is selected. The input data block has 7 channels: Z_H, Z_DR, K_DP, CC, BBR, BH, and TG. Each data item is processed by z-score standardization (subtracting the mean and dividing by the standard deviation) to construct training samples. The data of the BBR and BH channels retain their original physical dimensions and value ranges, and BH is in km.

[0040] Time stacking is performed on the data of 6 consecutive time instants to form a 4-dimensional input sample with a dimension of [6,7,32,32]. The corresponding label is the minute precipitation intensity R (mm / h) of the rain gauge at the current moment.

[0041] A total of 86 effective precipitation processes (hourly cumulative rainfall ≥ 0.5 mm and rain gauge availability ≥ 80%) are screened out, and divided into a training set (60 processes), a validation set (13 processes) and a test set (13 processes) according to a ratio of 7:1.5:1.5. It is ensured that the training set contains balanced samples of each occlusion degree (no occlusion BBR ≤ 0.1, partial occlusion 0.1 < BBR ≤ 0.4, severe occlusion 0.4 < BBR ≤ 0.7, complete occlusion BBR > 0.7) and each precipitation type (stratiform, convective, mixed).

[0042] Step S4: As Figure 4 shown, a neural network model is constructed; Terrain occlusion-aware encoding module: BBR and BH of each azimuth-radial range bin are converted through two fully connected layers (2→32→16) and spatial broadcasting into a 16-channel occlusion-aware feature map O (32×32×16).

[0043] Spatial feature extraction module: A U-Net architecture enhanced by CBAM attention is adopted. The number of input channels of the encoder is 7 (4 polarization parameters + 3 terrain auxiliary channels). The encoder includes 4 downsampling layers, each layer consists of two 3×3 convolution-BatchNorm-ReLU units, and the number of channels is 32, 64, 128, 256 in sequence. Each layer is followed by a CBAM attention module, which is followed by an occlusion-aware gating OAG.

[0044] OAG mechanism: As Figure 5 shown, the 16-channel occlusion-aware feature map O is downsampled through adaptive average pooling to the same spatial size as the feature map F_l of the current layer (for example, 16×16 for the second layer), the dimension is reduced to the same number of channels as F_l through a 1×1 convolution, and the gating map G_l is generated via Sigmoid. The gated feature is F_l_gated=F_l⊙G_l+F_l⊙(1-G_l)⊙β_l, where β_l is a learnable baseline compensation factor for each layer.

[0045] The decoder gradually restores the spatial resolution to 32×32 through bilinear interpolation upsampling and skip connections. The output spatial feature map size is 32×32×64.

[0046] Temporal context aggregation module: A single-layer BiLSTM (hidden layer dimension 256) is used as input. The spatial feature vectors flattened over 6 consecutive time intervals (32×32×64=65536 dimensions) are input, and the forward and backward hidden states at the last time interval are concatenated to form a 512-dimensional spatiotemporal feature. This feature is then mapped back to a 32×32×16 spatiotemporal encoded feature map through a fully connected layer.

[0047] Dual-head output module: The first output head is a precipitation intensity head, which, after two 3×3 convolutions (16→8→1 channels) and Softplus activation (ensuring non-negativity), outputs a 32×32×1 R field. The second output head is a DSD parameter head, which, after two 3×3 convolutions (16→16→3 channels), outputs log values ​​respectively. 10 (N_w), D_m, and μ.

[0048] Occlusion Filling Subnetwork: The region in the spatiotemporal encoded feature map corresponding to the grid points of the completely occluded blind area is taken as the region to be filled. Partial convolution is used to propagate and fill the region layer by layer towards the center of the blind area. Partial convolution only performs convolution operations on valid (non-occluded) pixels, and the occluded region is "swallowed" layer by layer through a mask update mechanism. After 3 layers of partial convolution and 2 layers of upsampling, the filled 32×32×1 precipitation intensity field is output (merged with the output of the first output head).

[0049] Step S5: As Figure 6 As shown, construct the total loss function; The total loss function contains 7 terms. The parameters for the newly added mountain terrain-related constraint term are set as follows: The L_terrain_grad weight λ5 increases from 0.02 to 0.8 during training; The L_terrain_continuity weight increases from 0.01 to 0.5 during training; The L_VPR weight λ6 increases from 0.02 to 1.0; Adaptive weight scheduling: Total training epochs E=250. For the first E / 5=50 epochs, λ5=0.02, λ6=0.02, with data-driven terms and basic physical constraints dominating rapid convergence; from epochs 50-125, λ5 and λ6 linearly increase to 0.8 and 1.0 respectively; from epochs 125-200, they remain high; from epochs 200-250, they gradually decay to stable values. The scheduling of other weights λ1-λ4 is the same as the original scheme.

[0050] Training hyperparameters: AdamW optimizer, initial learning rate 1e-3, cosine annealing to 1e-6, batch size 24 (limited by GPU memory). Training is performed on an NVIDIA A100 GPU, with a total time consumption of approximately 18 hours.

[0051] Step S6: Train the neural network model, deploy the neural network model, and perform effect evaluation; Quantitative evaluation is performed on the method of the present invention and 6 comparison methods for 13 precipitation processes in the test set. The evaluation is carried out in three sub-regions: the overall region, the blocked region (BBR>0.1), and the non-blocked region (BBR≤0.1). The blocked region is further subdivided into the partial pre-mountain blocked area (0.1<BBR≤0.4) and the severe post-mountain blocked / blind area (BBR>0.4).

[0052] Comparison methods: M1: Traditional R(K_DP, Z_DR) empirical formula method; M2: Traditional DEM+VPR method (DEM block correction first + VPR correction then + empirical formula estimation finally); M3: Pure data-driven CNN (7-channel input, including terrain auxiliary channels, without physical constraints and block gating); M4: GAN block inpainting method (Tan&Chen, 2025) + traditional R(K_DP, Z_DR) method; M5: Physically constrained deep learning QPE method (the previous version of the present invention, without the mountain terrain module); M6: CBAM U-Net + TCN method (the method of Chengdu Meteorological Bureau patent CN120468851B); Evaluation metrics: RMSE, MAE, CC, NSE, relative bias (RB), and physical consistency compliance rate.

[0053] The evaluation results are shown in the following table: Table 1 Evaluation results of the overall region ; Table 2 Evaluation results of the blocked region (BBR>0.1) ; Table 3a Evaluation results of the partial pre-mountain blocked area (0.1<BBR≤0.4, corresponding to the area near the ridge of the Front Longmen Mountain) ; Table 3b Evaluation results of the severe post-mountain blocked / blind area (BBR>0.4, corresponding to the area from the rear of the pre-mountain ridge to the Back Longmen Mountain) ; Table 4. Filling effect of completely occluded blind spots (BBR≥0.7)

[0054] Evaluation results show that the method of this invention outperforms all comparative methods across all evaluation metrics. Particularly in the severely shaded / blind areas of the back mountain (BBR>0.4), the RMSE is reduced by 58.4% compared to the traditional DEM+VPR method, and the relative bias decreases from -52.3% to -12.8%, significantly improving the systematic underestimation problem caused by large-area shading behind the front mountain ridge. In completely shaded blind areas, the RMSE of the shading-filling subnetwork is reduced by 37.8% compared to the GAN-filling method.

[0055] Taking a heavy precipitation event caused by a mesoscale convective system (MCS) on August 15, 2023 (maximum hourly rainfall of 68 mm / h) as an example: the precipitation system developed along the windward slope of the eastern slope of Longmenshan Mountain, forming a heavy precipitation belt along the Beichuan-Anxian line. Within an approximately 150° azimuth sector facing the back mountain from radar C, the 0.5°-1.0° elevation beam was severely obstructed by the front ridge of Longmenshan Mountain. Traditional methods and pure deep learning methods underestimated this area by more than 50%-80%. The method of this invention increases the effective elevation angle of the obstructed area to 2.4°-6.0°, and combined with an obstruction perception gating mechanism and an obstruction filling sub-network, successfully restored the reasonable spatial distribution of the precipitation field within the obstructed area (approximately 180 square kilometers), which highly matches the measured data from the rain gauge station on the back mountain (CC=0.86).

[0056] Although the embodiments of the present invention have been disclosed above, they are not limited to the applications listed in the specification and embodiments. For those skilled in the art, all features disclosed in the present invention, or all steps in all methods or processes disclosed, except for mutually exclusive features and / or steps, can be combined in any way without departing from the principles of the present invention. The present invention is not limited to the specific details and illustrations shown and described herein.

Claims

1. A method for quantitative precipitation estimation using mountain topographic radar, characterized in that, Includes the following steps: Step S1: Data collection and preprocessing; Acquire X-band dual-polarization radar data, ground rain gauge precipitation data, and high-resolution digital elevation model data for the target mountainous terrain area. Extract four polarization parameters from the X-band dual-polarization radar data: horizontal reflectivity factor Z_H, differential reflectivity factor Z_DR, differential propagation phase shift rate K_DP, and correlation coefficient CC. Perform non-meteorological echo identification and removal, differential propagation phase shift Φ_DP defolding and quality control, and system bias correction. Step S2: Terrain occlusion pre-analysis and adaptive hybrid elevation angle construction; Using high-resolution digital elevation model data, based on the radar site's latitude and longitude, antenna altitude, and parameters of each elevation layer, the ray tracing method is used to calculate the radar's beam obstruction rate at each azimuth-radial range library-elevation layer. Based on the beam obstruction rate, adaptive elevation angle selection is performed for each azimuth-radial range library to construct hybrid elevation polarization parameter data. Step S3: Create the training dataset; The training dataset is divided into a training set, a validation set, and a test set. Centered on the location of the ground rain gauge, mixed elevation polarization parameter data within the horizontal sampling window are extracted. Based on the polarization parameter channel, a beam obstruction rate channel, a beam height above ground channel, and a terrain elevation gradient channel are added to construct a multi-channel input data block containing polarization parameters of X-band dual-polarization radar data and terrain perception information. Step S4: Construct a neural network model; A physically constrained deep learning network is constructed, which includes a terrain occlusion perception coding module, a spatial feature extraction module, a temporal context aggregation module, a dual-head output module, an occlusion filling sub-network, and a physical consistency verification module. An occlusion perception gating mechanism is introduced in the spatial feature extraction module to adaptively adjust the feature contribution weights of regions with different occlusion levels according to the beam occlusion rate. The dual-head output module includes a first output head and a second output head in parallel. The first output head outputs the spatial distribution of precipitation intensity in minutes, and the second output head outputs the spatial distribution of raindrop spectrum distribution DSD parameters. Step S5: Construct the total loss function; Total loss function L total as follows: ； in, For data-driven loss terms, For the physical consistency constraint term of raindrop spectral distribution, These are constraint terms in the attenuation equation. These are non-negative constraint terms. For spatial smoothing constraints, For terrain occlusion consistency constraints, This is a vertical profile consistency constraint term. ~ These are the weighting coefficients for each constraint term; Step S6: Train the neural network model; The physical constraint deep learning network is trained using the training dataset, and the network parameters are optimized using gradient descent with the total loss function as the objective until the total loss function converges. Step S7: Deploy the neural network model; The trained physical constraint deep learning network is deployed in the business environment to receive radar volume scan data and digital elevation model (DEM) in real time. After preprocessing, terrain occlusion pre-analysis, and multi-channel data block construction, the data is input into the network for forward inference. The network outputs minute-level precipitation intensity R value for each radar grid point and performs time accumulation to generate a quantitative precipitation estimation grid point product with the corresponding time resolution.

2. The method for quantitative precipitation estimation using mountain topographic radar according to claim 1, characterized in that, The specific strategy for adaptive elevation angle selection in step S2 is as follows: When the beam obstruction rate of the lowest available elevation angle in the azimuth-radial distance database is less than or equal to the first threshold T1, the observation data is used directly; where the value of the first threshold T1 ranges from 0.3 to 0.

5. When the beam obstruction rate of the azimuth-radial distance coulomb is between the first threshold T1 and 0.1, it is in a partial obstruction state. After correction of the horizontal reflectivity factor Z_H, it is used, and the revised formula is: ,in, For beam blocking rate, The corrected horizontal reflectance. The measured value of the horizontal reflectivity factor before correction. This is the beam obstruction power correction factor; When the beam obstruction rate of the azimuth-radial distance library is greater than the first threshold T1, the search is performed at the next higher elevation angle. The highest available elevation angle with a beam obstruction rate of less than or equal to 0.1 is selected, and the highest available elevation angle is corrected using the vertical profile in subsequent steps and converted to the reference height. When the beam obstruction rate of the azimuth-radial distance library is greater than the first threshold T1 in all elevation layers, it is marked as a completely obstructed blind zone and reconstructed through the obstruction filling mechanism of a deep learning network. The beam obstruction rate is calculated using a Gaussian beam model. Specifically, under standard atmospheric or actual sounding refraction conditions, the beam's azimuth angle θ and elevation angle are calculated. The center height at a distance r is obtained by reading the terrain elevation profile along the beam propagation path from the digital elevation model (DEM). The power ratio of the beam blocked by the terrain is calculated based on the Gaussian distribution model of the antenna pattern. The beam blocking rate ranges from 0 to 1, where 0 represents no blocking and 1 represents complete blocking. The calculation of beam obstruction rate incorporates an echo climate statistical correction factor. Specifically, for each azimuth-radial distance database, the frequency of clear-sky echoes during non-precipitation periods is statistically analyzed to construct the observation obstruction probability (OBP). When the difference between the beam obstruction rate (BBR) calculated by the digital elevation model (DEM) and the observation obstruction probability (OBP) exceeds a set threshold, a correction formula is used. The correction is performed, among which... λ For the credibility weight of the digital elevation model (DEM), The corrected beam obstruction rate, The beam obstruction rate is calculated for the digital elevation model (DEM). When the accuracy of the DEM data does not meet the preset accuracy requirements or is difficult to obtain, the effective minimum available elevation angle of each azimuth-radial distance library is determined using the echo climate statistical frequency method.

3. The method for quantitative precipitation estimation using mountain topographic radar according to claim 1, characterized in that, The method for constructing the multi-channel input data block in step S3 is as follows: A rectangular sampling window of M×N azimuth-radial distance libraries is selected in the horizontal dimension to form a three-dimensional input tensor, where the values ​​of M and N range from 16 to 64. The input channels include 4 polarization parameter channels and 3 terrain sensing auxiliary channels, for a total of 7 channels; the 4 polarization parameter channels are Z_H, Z_DR, K_DP and CC, and the 3 terrain sensing auxiliary channels are beam obstruction rate BBR, beam height above ground BH and terrain elevation gradient TG. The 7 channels of 3D data blocks from T consecutive time intervals are stacked along the time dimension to form a 4D spatiotemporal input sample with dimensions [T,7,M,N], where T ranges from 3 to 12; each channel is standardized or normalized separately.

4. The method for quantitative precipitation estimation using mountain topographic radar according to claim 1, characterized in that, The specific structure of the terrain occlusion perception coding module mentioned in step S4 is as follows: The information of beam obstruction rate (BBR) and beam height above ground (BH) is encoded into an obstruction-aware feature vector through two fully connected layers. This vector is then converted into an obstruction-aware feature map (O) with the same resolution as the spatial feature map by a spatial broadcast layer. The dimension of the obstruction-aware feature map O is M×N×D_o, where D_o is the encoding dimension and its value ranges from 16 to 32.

5. The method for quantitative precipitation estimation using mountain topographic radar according to claim 1, characterized in that, The specific operation of the occlusion sensing gating mechanism described in step S4 is as follows: In each layer of the encoder of the spatial feature extraction module, the occlusion-aware feature map O generated by the terrain occlusion-aware coding module is sampled to the same spatial resolution as the current layer feature map F_l, and a gated weight map G_l is generated through the Sigmoid activation function. The current layer feature map F_l is then gated and weighted as follows: F_l_gated = F_l⊙G_l + F_l⊙(1-G_l)⊙β_l; where ⊙ denotes element-wise multiplication, β_l is the learnable baseline compensation factor for the occluded region, and F_l_gated is the gated and weighted feature map. The learnable weight matrix for gated convolution. This is an occlusion-aware feature map downsampled to the current layer resolution. This is a gated bias term; In occluded regions, the gating weights reduce the contribution of the original echo features and enhance the context compensation learned from the surrounding unoccluded regions; in unoccluded regions, the gating weights are close to 1, preserving the original echo features. The generation of the gating weight map G_l adopts at least one of the following information sources: univariate gating based on beam obscuration ratio BBR, bivariate gating based on BBR and BH, trivariate gating based on beam obscuration ratio BBR, beam height above ground BH and terrain elevation gradient TG, or learnable gating learned from the digital elevation model DEM. The spatial feature extraction module adopts an encoder-decoder structure. The encoder extracts multi-scale spatial features through multi-layer convolution, and each layer is followed by a CBAM attention module and an occlusion-aware gating mechanism. The encoder-decoder structure adopts at least one of the following: 3D ResNet, a CBAM attention-enhanced U-Net architecture, or a Swing Transformer module. The temporal context aggregation module adopts at least one of the following structures: Temporal Convolutional Network (TCN), Bidirectional Long Short-Term Memory Network (BiLSTM), or a spatiotemporal attention module. The occlusion infilling subnetwork is based on partial convolution or spatial attention mechanisms. For grid points marked as completely occluded blind areas, it takes the spatiotemporal joint features of the surrounding non-occluded areas as input, performs convolution operations on the effective non-occluded pixels through partial convolution, and propagates the effective information to the center of the blind area layer by layer through a mask update mechanism. The output is the infilled precipitation intensity field and raindrop spectrum distribution DSD parameter field. The occlusion infilling subnetwork adopts a U-Net structure decoder, which receives the feature maps after occlusion perception gating in each layer of the encoder through skip connections. The raindrop spectral distribution DSD parameters output by the second output head are described using a normalized Gamma distribution, which is N(D) = N_w·f(μ)·(D / D_m)^μ·exp[-(4+μ)·D / D_m], where N_w is the normalized intercept parameter, D is the raindrop diameter, D_m is the mass-weighted average diameter, μ is the shape parameter, and f(μ) is the normalization function; the second output head outputs respectively The spatial distribution of the three parameters, D_m and μ.

6. The method for quantitative precipitation estimation using mountain topographic radar according to claim 1, characterized in that, The terrain occlusion consistency constraint L_terrain mentioned in step S5 includes the following three sub-items: Sub-item 1 is the occlusion gradient consistency constraint: calculate the cosine similarity between the spatial gradient magnitude of precipitation intensity R and the normalized terrain elevation gradient magnitude, and apply a penalty term L_terrain_grad to areas where the two are inconsistent. Where w(BBR) is a weighting function that increases with the beam obstruction rate, R is the estimated precipitation intensity output by the network, and H is the terrain elevation; Sub-item 2 is the joint consistency constraint for occlusion attenuation: In partially occluded areas, the polarization parameter after occlusion correction and precipitation attenuation correction by the digital elevation model (DEM) is consistent with the equivalent polarization parameter calculated from the precipitation intensity R output by the physical constraint deep learning network using the ZR relationship. Sub-item 3 is the spatial continuity constraint for the shading blind zone: It constrains the minimization of the weighted spatial second-order difference between the precipitation intensity of the shading-filled area and the adjacent unshading area. The minimized value of the weighted spatial second-order difference is L_terrain_continuity=Σw_i·|R_i-2·R_j+R_k|², where i is a grid point inside the blind zone, j and k are two adjacent grid points along the terrain gradient direction, w_i is the weight that decreases as the distance between grid point i and the edge of the blind zone increases, R_i is the estimated precipitation intensity of grid point i inside the blind zone, R_j is the estimated precipitation intensity of adjacent grid point j, and R_k is the estimated precipitation intensity of adjacent grid point k. The weight w_i in sub-item 3 adopts a distance decay function, which is at least one of the following: exponential decay exp(-d_i / σ), Gaussian decay, or a decay function based on terrain connectivity, where d_i is the distance of grid point i from the nearest unshading grid point at the edge of the blind zone, and σ is the decay scale. The calculation method for the vertical profile consistency constraint term L_VPR is as follows: For the same horizontal position, forward inference is performed using low-elevation and high-elevation data obtained from the adaptive elevation angle selection results. The difference in precipitation intensity between the two inferences, L_VPR, is calculated as a constraint: L_VPR=(1 / N_pair)·Σ(R_low-R_high)²; where R_low is the precipitation intensity inferred using low-elevation data, R_high is the precipitation intensity inferred using high-elevation data, and N_pair is the effective logarithm of elevation angles. In unobstructed areas where the beam obstruction rate of both elevation angle data is less than or equal to 0.1, the differential L_VPR is applied. The physical consistency constraint L_DSD is calculated as follows: Based on the raindrop spectral distribution DSD parameters output by the network, including the normalized intercept parameter N_w, the mass-weighted average diameter D_m, and the shape parameter μ, the equivalent polarization parameter prediction values ​​Z_H^pred, Z_DR^pred, and K_DP^pred are calculated forward using the T-matrix scattering model. These prediction values ​​are then compared with the radar-measured polarization parameters Z_H^obs, Z_DR^obs, and K_DP^obs. The differences are used as constraints. Constraints are applied to grid points in unobstructed regions where the beam obstruction rate is less than or equal to 0.1 to avoid polarization parameter distortion caused by obstruction affecting the reliability of the constraints. The weight coefficients λ1 to λ6 are selected using an adaptive scheduling strategy: in the early stage of training, λ5 and λ6 are selected with smaller values, ranging from 0.01 to 0.05, while λ1 to λ4 are selected with medium values, ranging from 0.1 to 0.5, so that data-driven terms and basic physical constraints dominate the fast convergence; in the middle and later stages of training, λ5 and λ6 are gradually increased to 0.5-1.5, so that terrain-related physical constraints play a dominant role.

7. The method for quantitative precipitation estimation using mountain topographic radar according to claim 1, characterized in that, The target mountainous terrain region has a dual structure of front mountain and back mountain. With the regional fault zone as the boundary, the front mountain region is a hilly, low mountain, and medium mountain landform with an altitude below 2,500 meters, while the back mountain region is a medium mountain and high mountain landform with an altitude above 3,500 meters. The ridge of the front mountain partially blocks the radar deployed on the edge of the basin, while the back mountain region forms a severe block or a complete blind spot. Adaptive elevation angle selection and terrain occlusion consistency constraints employ differentiated threshold parameters and constraint weights for partially occluded areas in the front mountain and severely occluded areas in the back mountain. When S-band or C-band operational weather radars are deployed simultaneously in the study area, the polarization parameters of the S-band or C-band radars are added as an additional input channel to the multi-channel input data block; a cross-band feature alignment module is added to the network to learn the feature correspondence between S-band, C-band and X-band data at corresponding spatial locations through a spatial attention mechanism; a cross-band consistency constraint term L_crossband is added to the total loss function to constrain the consistency of precipitation intensity retrieved from different bands at the same horizontal location.

8. A mountain terrain radar quantitative precipitation estimation system, which employs any one of the mountain terrain radar quantitative precipitation estimation methods described in claims 1 to 7, characterized in that, The mountain terrain radar quantitative precipitation estimation system includes the following modules: The data acquisition and preprocessing module is used to acquire X-band dual-polarization radar base data, ground rain gauge data, and digital elevation model (DEM), extract four polarization parameters: horizontal reflectivity factor Z_H, differential reflectivity factor Z_DR, differential propagation phase shift rate K_DP, and correlation coefficient CC, and perform non-meteorological echo identification and removal, differential propagation phase shift Φ_DP defolding and quality control, and system deviation correction. The terrain obstruction pre-analysis and adaptive elevation angle selection module is used for the digital elevation model (DEM). It uses the ray tracing method to calculate the beam obstruction rate (BBR) of the radar at each azimuth-radial range-elevation layer, and performs adaptive elevation angle selection based on the beam obstruction rate to construct hybrid elevation polarization parameter data. The multi-channel input data construction module is used to construct standardized multi-channel input data blocks from the polarization parameter channel, beam obstruction rate channel, beam height above ground channel, and terrain elevation gradient channel, and stack them in the time dimension to form a four-dimensional spatiotemporal input sample. The Physically Constrained Deep Learning Mountain QPE Network Module includes a terrain occlusion awareness coding submodule, a spatial feature extraction submodule, a temporal context aggregation submodule, a dual-head output submodule, an occlusion filling subnetwork, and a physical consistency verification submodule. An occlusion perception gating mechanism is introduced into the spatial feature extraction submodule to adaptively adjust the feature contribution weight of regions with different occlusion levels based on the beam occlusion rate. The dual-head output submodule includes a parallel precipitation intensity output head and a DSD parameter output head; The loss function calculation module is used to construct the total loss function, which includes the data-driven loss term L_data, the physical consistency constraint term L_DSD, the attenuation equation constraint term L_attenuation, the nonnegation constraint term L_nonneg, the spatial smoothness constraint term L_smooth, the terrain occlusion consistency constraint term L_terrain, and the vertical profile consistency constraint term L_VPR. Among them, L_terrain includes the occlusion gradient consistency constraint sub-term, the occlusion attenuation joint consistency constraint sub-term, and the occlusion blind zone spatial continuity constraint sub-term. The model training module is used to optimize network parameters based on training data with the total loss function as the objective. The QPE product generation module is used to receive radar data in real time and output quantitative precipitation estimation grid products. The quantitative precipitation estimation grid products include occlusion probability products, which indicate the beam occlusion rate (BBR) value of each grid point, and occlusion filling confidence products, which indicate the filling confidence score of each grid point for forecasters to refer to and judge.

9. A storage medium for a mountain terrain radar quantitative precipitation estimation system, used to store the calculation program of the mountain terrain radar quantitative precipitation estimation method according to any one of claims 1 to 7, characterized in that, The storage medium of the mountain terrain radar quantitative precipitation estimation system is installed on a computer, and the calculation program of the mountain terrain radar quantitative precipitation estimation method is run through the computer.

10. A mountain terrain radar quantitative precipitation estimation device, used to install the mountain terrain radar quantitative precipitation estimation system of claim 8, characterized in that, The aforementioned mountain terrain radar quantitative precipitation estimation device includes a memory and a processor. The memory stores the calculation program for the mountain terrain radar quantitative precipitation estimation method, and the processor executes the calculation program for the mountain terrain radar quantitative precipitation estimation method for field testing.

Images