Quantitative model-based radar combined reflectivity correction method

Abstract

The invention discloses a radar combined reflectivity correction method based on a quantitative model, and the method comprises the steps: collecting the observation data of a preset radar network, building a grid search box with an S-band radar as a benchmark, calculating the deviation between a networking reflectivity mean value and a reference mean value according to the observation data, and carrying out the correction when the difference value exceeds a deviation threshold value, performing error decomposition on the corrected deviation to obtain a system deviation and a random error; for a specific resolution reflectivity field, matching to a radar grid by adopting bilinear interpolation and a space-time sliding window, synchronously extracting a terrain height and a beam blocking rate, and constructing a three-dimensional feature matrix; performing two-stage attenuation correction on the system deviation and the random error, introducing a terrain correction factor to construct a time-varying coefficient equation, and constructing a space-time adaptive correction quantitative model according to the time-varying coefficient equation; and extracting features to generate a spatial weight map to obtain an enhanced feature map, and outputting a correction result based on the space-time adaptive correction quantization model of the enhanced feature map.

Description

Technical Field

[0001] This invention relates to the field of meteorology, and more particularly to a method for correcting radar composite reflectivity based on a quantization model. Background Technology

[0002] In the field of meteorological observation and forecasting, radar combined reflectivity is a core data point reflecting the intensity, distribution, and evolution of precipitation. Its accuracy directly affects the timeliness and accuracy of early warnings for severe weather events such as short-term severe convection and rainstorms. With the promotion of radar network observation technology, the application of multi-source radar data fusion has become a trend, but many technical bottlenecks still exist in actual observations.

[0003] Radar observations are susceptible to system biases and propagation attenuation, such as radome attenuation, atmospheric attenuation, and beam obstruction caused by terrain, resulting in discrepancies between raw reflectivity data and actual precipitation. Performance differences between different radar models (such as S-band and other band radars) and spatiotemporal inconsistencies in observation further exacerbate the heterogeneity of network data.

[0004] Existing correction methods mostly focus on correcting single error sources, lacking adaptability to spatiotemporal dynamic changes and failing to take into account multi-dimensional influencing factors such as terrain and time lag. Traditional interpolation and correction models are insufficient in characterizing the topographic features of complex underlying surfaces and lack specificity in error correction for strong and weak echo areas, resulting in poor spatial continuity and temporal consistency of the correction results.

[0005] Meanwhile, the rapid evolution of weather systems requires correction models to have real-time response capabilities. However, existing methods are insufficient in terms of the comprehensiveness of feature extraction and model convergence efficiency, making it difficult to meet the demand for high-resolution, high-precision reflectance data in refined forecasting. Therefore, there is an urgent need to construct a quantitative correction method that integrates multi-source features and adapts to spatiotemporal changes to improve the quality of radar composite reflectance data and provide reliable support for meteorological early warning and disaster prevention and mitigation. Summary of the Invention

[0006] The purpose of this invention is to provide a method for correcting radar composite reflectivity based on a quantization model.

[0007] To achieve the above objectives, the present invention is implemented according to the following technical solution: This invention includes the following steps: Collect observation data from a pre-defined radar network and preprocess the observation data; the observation data includes model data, radar data, and auxiliary data. A regional detection algorithm is used to eliminate the effects of system bias and attenuation. A grid search box is established based on the S-band radar. The deviation between the mean network reflectivity and the reference mean is calculated based on the observation data. When the difference exceeds the deviation threshold, correction is performed. The deviation that has been corrected is decomposed to obtain the system bias and random error. For a specific resolution reflectivity field, bilinear interpolation and spatiotemporal sliding window matching are used to match the radar grid, simultaneously extracting terrain height and beam blocking rate, and constructing a three-dimensional feature matrix based on mode reflectivity, terrain, and time lag. A two-stage attenuation correction method is used to correct the system deviation and the random error by employing a network attenuation algorithm and gradient descent method. A terrain correction factor is introduced to construct a time-varying coefficient equation, and a spatiotemporal adaptive correction quantization model is constructed based on the time-varying coefficient equation. Multi-scale attention gating is introduced, and features are extracted through multi-scale dilated convolution to generate spatial weight maps to obtain enhanced feature maps. The data to be corrected is input into a spatiotemporal adaptive correction quantization model based on the enhanced feature maps, and the correction results are output.

[0008] In this embodiment, the method for calculating the deviation between the mean reflectance of the network and the reference mean includes: Spatial scanning of the networked radar reflectivity data is performed using a grid sliding window, and the deviation between the observed mean within the window and the S-band radar reference mean is calculated: ; in The mean reflectivity of the network is given. This is the S-band reference mean. This is the original observed reflectance value. This represents the deviation between the observed mean and the S-band radar reference mean.

[0009] Furthermore, the method of matching radar grids using bilinear interpolation and spatiotemporal sliding windows includes: The latitude and longitude coordinates of the model data are uniformly converted into the Albers equal-area projection used by the radar data. The projection distortion error is eliminated by the seven-parameter Bursa model. A 20-minute assimilation time window is set, and time interpolation is performed on the model data to generate high-frequency data sequences with equal time intervals. For any target point in the radar grid Calculate the row and column indices of the target point in the pattern grid. Determine the coordinates of the four surrounding pattern grid points, with the top left corner as the coordinate point. Top right corner Bottom left corner bottom right corner ; Bilinear interpolation is performed based on the weighting coefficients in the horizontal and vertical directions, expressed as follows: ; ; ; in These are the weighting coefficients in the horizontal direction. The weighting coefficient is in the vertical direction. The x-coordinate of the target point The ordinate of the target point. The x-coordinate of the top-left grid point. The x-coordinate of the top right grid point is... The ordinate of the top-left grid point is... The ordinate of the bottom left grid point is... The reflectance value of the pattern grid point in the i-th row and j-th column is... For preliminary interpolation results, The reflectance value of the pattern grid point in the (i+1)th row and jth column; For radar grid points within a 100km radius of the model grid boundary, range-weighted averaging is used instead of pure bilinear interpolation, with weighting coefficients... ,in This represents the distance from the target point to the pattern grid point. When the terrain height of the target point exceeds the height corresponding to the lowest layer of the pattern, vertical interpolation correction is automatically enabled, with the correction coefficient... ; A window size is defined, a spatial weight matrix is ​​generated based on the Gaussian kernel function, and a time weight matrix is ​​generated based on the exponential decay function. A spatiotemporal sliding window is constructed based on the window size, spatial weight matrix, and time weight matrix. A three-dimensional spatiotemporal sliding window is adopted, with a spatial dimension grid and three time dimensions, forming a cubic data structure. The Gaussian kernel function is... The exponential decay function is ; The mean, standard deviation, and gradient magnitude of the model reflectivity within the window are calculated as matching evaluation indicators. The similarity between the model forecast sequence and the radar observation sequence is calculated using the dynamic time warping algorithm, and the forecast time corresponding to the minimum distance is selected as the benchmark. The spatial weights, temporal weights, and gradient magnitudes are weighted and fused together, expressed as follows: ; in The reflectance value of the mode within the window at the k-th spatial grid and the t-th time interval is... The total number of spatial grids, For the number of times in time, To match the fused reflectivity values ​​to radar grid points, The spatial weight of the k-th spatial grid. For spatial grid indexing, For time-time index, Let t be the time weight for the t-th time period.

[0010] Furthermore, the method for constructing the three-dimensional feature matrix based on pattern reflectivity, terrain, and time lag includes: Acquire pattern reflection field data, terrain data, and beam blocking rate data, and perform standardization processing to extract terrain features, including terrain factors and radar beam influence factors. Terrain factors include regional terrain type and relative height. Regional terrain type is classified based on slope threshold and encoded as a one-hot vector. Relative height is calculated by determining the standard deviation of terrain height within a neighborhood window. Radar beam influence factors include effective beam height and blocking correction coefficient. The expression is: ; in The elevation of the radar station. The height of the beam center. The effective beam height; ,in This is the blocking correction factor. To determine the beam blocking rate, a Gaussian distribution mode is used to simulate beam energy attenuation. The reflectivity values ​​of the mode at three time lags are extracted from the radar grid points. The reflectivity change rate is calculated to capture the evolution trend of echo intensity and obtain time-varying characteristics. A sinusoidal function is introduced to encode the time periodicity to obtain the diurnal periodicity variation factor. Time feature alignment is performed based on the time-varying characteristics and the diurnal periodicity variation factor. The expression is: ; ; in The current mode reflectivity value. The reflectance values ​​for the first 6 minutes of the model. The reflectance values ​​for the first 12 minutes of the model. The rate of change of reflectance is 1. The rate of change of reflectance is 2; ,in This is a daily cycle variation factor; The three-dimensional feature matrix includes a time dimension, a spatial dimension, and a feature dimension. The spatial dimension is the latitude and longitude range of the radar grid. The time dimension is the 0-12 minute lead time window of the model forecast corresponding to three time lags. The feature dimension contains 12 feature variables for each spatiotemporal region. The matrix is ​​concatenated in the order of space, time, and features and output as a three-dimensional feature matrix. / / The 12 feature variables are respectively... , , , , Slope threshold, terrain height, terrain height standard deviation, beam blocking rate, blocking correction factor, diurnal variation factor, and terrain type coding; Furthermore, the method for two-stage attenuation correction of the systematic bias and the random error includes: After marking clutter areas based on the skewness-kurtosis joint identification method, effective echoes from the upper elevation angle are used to fill the clutter area and the scanning time of radars within the network is synchronously calibrated; among them, meteorological echo skewness > 0 and kurtosis > 3, and ground clutter skewness < 0 and kurtosis < 2. Using a stable S-band radar within the network as a reference, the mean reflectivity deviation between the reference radar and other radars is calculated. During clear, precipitation-free periods, the reflectivity difference between the reference radar and the target radar is statistically analyzed. The radome attenuation coefficient is fitted, and a first-order correction is applied to the target radar reflectivity. ,in The attenuation coefficient of the radome is... The target radar reflectivity after first-order correction. The original reflectivity is given; the expression for the radome attenuation coefficient is: ; in , Atmospheric attenuation coefficient, To detect distance, This represents the average reflectance deviation. Beam blocking rate is calculated based on terrain data, and a lookup table based on terrain height and blocking rate is constructed. Attenuation is calculated segment by segment. When the beam blocking rate is less than or equal to 0.2: When the beam blocking ratio is greater than 0.2 and less than or equal to 0.5: When the beam blocking ratio is greater than 0.5: ;in This is the attenuation amount. Beam blocking ratio; For regions with abrupt changes in blocking rate, gridded Gaussian smoothing is applied. Based on the first-order corrected target radar reflectivity and attenuation, a first-order corrected reflectivity field is output, expressed as:

[0011] in This is the first-order corrected reflectivity field; Based on the first-order corrected reflectivity field, the objective function for random error optimization is defined as follows:

[0012] in Let the first-order corrected reflectivity field be the k-th spatial grid point. Let k be the reference radar observation reflectivity of the k-th spatial grid point. This is the proportionality coefficient. This is the offset. For parameters to be optimized, For L2 regularization parameters, Let the objective function be the parameter to be optimized. This represents the total number of spatial grids. The gradient is calculated using mini-batch gradient descent, and the coefficients to be optimized are dynamically updated using a cosine annealing strategy. The process is iterated until the change in the objective function is less than the error threshold for five consecutive iterations or the maximum number of iterations reaches 500. Then the iteration stops and the optimized parameters are output. Second-order correction is performed based on the optimized parameters, expressed as follows:

[0013] in The optimized scaling factor, This is the optimized offset. This is the second-order corrected reflectivity field; For the strong echo region, nonlinear fine-tuning is performed, and an intensity correction factor is introduced, expressed as follows:

[0014] in For the second-order corrected reflectivity field in the strong echo region, This is the intensity correction factor.

[0015] Furthermore, the method for constructing time-varying coefficient equations by introducing terrain correction factors includes:

[0016]

[0017] in For terrain correction factors, For terrain height, For final reflectivity, For the second-order corrected reflectivity field, For beam blocking variables, For the residual term, Let be the model prediction weight coefficient for the t-th time period. The terrain factor weight coefficient is the value at time t. The blocking factor weight coefficient is the value at time t.

[0018] Furthermore, the method for constructing a spatiotemporal adaptive correction quantization model based on the time-varying coefficient equation includes: Using the final reflectance in the time-varying coefficient equation as the objective function of the spatiotemporal adaptive correction quantization model, the terrain height is compressed to the [0,1] interval using the hyperbolic tangent function. Nonlinear weights are set according to the beam-blocking variable: when the beam-blocking variable is less than 0.2, the nonlinear weight is set to 0.1; when the beam-blocking variable is greater than or equal to 0.2 and less than 0.5, the nonlinear weight is set to 0.4; when the beam-blocking variable is less than or equal to 0.5, the nonlinear weight is set to 0.8. The expression for compressed terrain height is: ; in This is the compressed terrain height, which is... ; Introducing the zero-order optimization concept of the Zebra Optimization algorithm, the gradient is estimated through two forward propagations. The first forward propagation calculates the loss based on the current coefficients. The second forward propagation adds random perturbations to the model forecast weight coefficients, and the calculation... ; We construct a weather type and coefficient matrix mapping library to store the optimal coefficient combinations for typical weather events such as rainstorms, typhoons, and cold waves. We identify the current weather type through KL divergence and call on historical knowledge to accelerate convergence.

[0019] Furthermore, methods for obtaining enhanced feature maps by extracting features through multi-scale dilated convolution to generate spatial weight maps include: Multi-scale dilated convolution comprises an input layer, a feature extraction module, a weight generation module, and a feature enhancement module. The feature extraction module consists of three parallel dilated convolution branches: branch 1 extracts detail feature maps, branch 2 extracts mesoscale feature maps, and branch 3 extracts macroscale feature maps. The input layer receives the spatiotemporally matched pattern reflectivity field and normalizes its size. The feature extraction module consists of three parallel dilated convolution branches with different dilation rates (d=1,3,5), corresponding to receptive field sizes of 3×3, 7×7, and 11×11. The weight generation module calculates the feature importance of each branch using a channel attention mechanism, generating a 512×512×1 spatial weight map. The feature enhancement module performs element-wise multiplication of the weight map with the original feature map to highlight the feature responses of strongly weighted regions. A fusion coefficient is introduced, and the weights of each branch are dynamically adjusted through backpropagation. The feature maps of the three branches are then concatenated along the channel dimension using the branch weights to reduce the dimensionality and obtain the fused feature map. Global average pooling is applied to the fused feature map to generate channel descriptors. Channel weights are learned through a two-layer fully connected network. The channel weights are multiplied by the fused feature map, summed along the channel dimension, and then normalized using the S-function to generate a spatial weight map, expressed as: ; ; ; in For the spatial weight map of row a and column u, For spatial row index, For spatial column indexes, Let c be the weight of the channel. For channel index, For the c-th channel descriptor, The weights are for the first fully connected layer. These are the weights for the second fully connected layer. It is the ReLU activation function. It is the Sigmoid activation function. For the space is high, For the space is wide; The spatial weight map is multiplied pixel-by-pixel with the original reflectivity feature map to generate the enhanced feature map, expressed as: ; in This is the original reflectivity feature map. To enhance the feature map.

[0020] The beneficial effects of this invention are: This invention is a correction method for radar combined reflectivity based on a quantization model. Compared with the prior art, this invention has the following technical advantages: This invention eliminates the effects of system bias, attenuation, and terrain occlusion simultaneously through preprocessing, bilinear interpolation and spatiotemporal sliding window matching, construction of a three-dimensional feature matrix, two-stage attenuation correction, model building, and acquisition of enhanced feature maps. Compared to single error correction, this is more comprehensive and improves data accuracy. Combining time-varying coefficient equations and multi-scale attention, it dynamically adapts to weather evolution and terrain differences, outperforming traditional static models. Two-stage correction and historical knowledge accelerate convergence, and nonlinear fine-tuning in strong echo areas enhances targeting while ensuring real-time performance. Adaptable to different radar models and weather types, it can directly serve meteorological early warning, making it more practical. Attached Figure Description

[0021] Figure 1 This is a flowchart illustrating the steps of the radar combined reflectivity correction method based on a quantization model according to the present invention. Detailed Implementation

[0022] The present invention will be further described below through specific embodiments. The illustrative embodiments and descriptions herein are used to explain the present invention, but are not intended to limit the present invention.

[0023] The radar combined reflectivity correction method based on a quantization model of the present invention includes the following steps: like Figure 1 As shown, this embodiment includes the following steps: Collect observation data from a pre-defined radar network and preprocess the observation data; the observation data includes model data, radar data, and auxiliary data. In the actual evaluation, the radar network in a certain eastern region was used as the application scenario. The network includes 8 S-band radars and 12 C-band radars, covering a range of 115°-125° east longitude and 28°-38° north latitude, with a horizontal resolution of 1km×1km and a time resolution of 6 minutes. Data types acquired: Model data consists of 0-12 hour forecast data output from WRF mode (spatial resolution 3km), radar data consists of the original reflectivity of the networked radar (ZDR, KDP and other polarization parameters are acquired synchronously), and auxiliary data includes 30-second resolution terrain data and ERA5 atmospheric profile data. Preprocessing: Quality control is performed on the radar data, removing noise data with reflectivity < -10dBZ; after Albers equal-area projection transformation, the model data is processed using a seven-parameter Bursa model (translation parameters ΔX=123.4m, ΔY=456.7m, ΔZ=78.9m, rotation parameters...). 0.0001 rad 0.0002 rad To eliminate projection errors, a 20-minute assimilation time window was set, and linear interpolation was used to generate high-frequency sequences with 6-minute intervals. A regional detection algorithm is used to eliminate the effects of system bias and attenuation. A grid search box is established based on the S-band radar. The deviation between the mean network reflectivity and the reference mean is calculated based on the observation data. When the difference exceeds the deviation threshold, correction is performed. The deviation that has been corrected is decomposed to obtain the system bias and random error. For a specific resolution reflectivity field, bilinear interpolation and spatiotemporal sliding window matching are used to match the radar grid, simultaneously extracting terrain height and beam blocking rate, and constructing a three-dimensional feature matrix based on mode reflectivity, terrain, and time lag. A two-stage attenuation correction method is used to correct the system deviation and the random error by employing a network attenuation algorithm and gradient descent method. A terrain correction factor is introduced to construct a time-varying coefficient equation, and a spatiotemporal adaptive correction quantization model is constructed based on the time-varying coefficient equation. Multi-scale attention gating is introduced, and features are extracted through multi-scale dilated convolution to generate spatial weight maps to obtain enhanced feature maps. The data to be corrected is input into a spatiotemporal adaptive correction quantization model based on the enhanced feature maps, and the correction results are output.

[0024] In this embodiment, the method for calculating the deviation between the mean reflectance of the network and the reference mean includes: Spatial scanning of the networked radar reflectivity data is performed using a grid sliding window, and the deviation between the observed mean within the window and the S-band radar reference mean is calculated: ; in The mean reflectivity of the network is given. This is the S-band reference mean. This is the original observed reflectance value. This represents the deviation between the observed mean and the S-band radar reference mean. In the actual assessment, the grid sliding window was 5×5 and the horizontal resolution was 1km×1km; Grid search box settings: A 5×5 grid sliding window is used to perform spatial scanning of the network data. The Nanjing radar with stable performance is selected as the reference benchmark for the S-band radar. Deviation calculation and correction triggering: via formula Calculate the deviation within the window, setting the deviation threshold to 3 dBZ. Initiate the correction process immediately; Error decomposition: The least squares method is used to decompose the deviation into systematic deviation (fixed deviation caused by radome attenuation and beam blocking) and random error (instantaneous deviation caused by atmospheric turbulence). After decomposition, the systematic deviation accounts for about 65%-75% and the random error accounts for about 25%-35%.

[0025] In this embodiment, the method of matching radar grids using bilinear interpolation and spatiotemporal sliding windows includes: The latitude and longitude coordinates of the model data are uniformly converted into the Albers equal-area projection used by the radar data. The projection distortion error is eliminated by the seven-parameter Bursa model. A 20-minute assimilation time window is set, and time interpolation is performed on the model data to generate high-frequency data sequences with equal time intervals. For any target point in the radar grid Calculate the row and column indices of the target point in the pattern grid. Determine the coordinates of the four surrounding pattern grid points, with the top left corner as the coordinate point. Top right corner Bottom left corner bottom right corner ; Bilinear interpolation is performed based on the weighting coefficients in the horizontal and vertical directions, expressed as follows: ; ; ; in These are the weighting coefficients in the horizontal direction. The weighting coefficient is in the vertical direction. The x-coordinate of the target point The ordinate of the target point. The x-coordinate of the top-left grid point. The x-coordinate of the top right grid point is... The ordinate of the top-left grid point is... The ordinate of the bottom left grid point is... The reflectance value of the pattern grid point in the i-th row and j-th column is... For preliminary interpolation results, The reflectance value of the pattern grid point in the (i+1)th row and jth column; For radar grid points within a 100km radius of the model grid boundary, range-weighted averaging is used instead of pure bilinear interpolation, with weighting coefficients... ,in This represents the distance from the target point to the pattern grid point. When the terrain height of the target point exceeds the height corresponding to the lowest layer of the pattern, vertical interpolation correction is automatically enabled, with the correction coefficient... ; A window size is defined, a spatial weight matrix is ​​generated based on the Gaussian kernel function, and a time weight matrix is ​​generated based on the exponential decay function. A spatiotemporal sliding window is constructed based on the window size, spatial weight matrix, and time weight matrix. A three-dimensional spatiotemporal sliding window is adopted, with a spatial dimension grid and three time dimensions, forming a cubic data structure. The Gaussian kernel function is... The exponential decay function is ; The mean, standard deviation, and gradient magnitude of the model reflectivity within the window are calculated as matching evaluation indicators. The similarity between the model forecast sequence and the radar observation sequence is calculated using the dynamic time warping algorithm, and the forecast time corresponding to the minimum distance is selected as the benchmark. The spatial weights, temporal weights, and gradient magnitudes are weighted and fused together, expressed as follows: ; in The reflectance value of the mode within the window at the k-th spatial grid and the t-th time interval is... The total number of spatial grids, For the number of times in time, To match the fused reflectivity values ​​to radar grid points, The spatial weight of the k-th spatial grid. For spatial grid indexing, For time-time index, The time weight for the t-th time period; In the actual evaluation, the time interval was 6 minutes, with three time dimensions: the current time, the previous 6 minutes, and the next 6 minutes. The cube data structure was 15×15×3, and the spatial dimension of the grid was 5×5. Bilinear interpolation implementation: For the radar grid target point (x,y), after determining the four grid points around the model grid, the horizontal weighting coefficient and the vertical weighting coefficient are calculated according to the formula to complete the preliminary interpolation; for the 236 radar grid points within 100km of the model grid boundary, range-weighted average is used to replace pure bilinear interpolation. Spatiotemporal sliding window construction: The spatial dimension adopts a 5×5 grid, and the time dimension includes three time periods: the current time, the previous 6 minutes, and the next 6 minutes, forming a 15×15×3 cubic data structure; the spatial weight matrix is ​​generated by the Gaussian kernel function (σ=2.0), and the time weight matrix is ​​generated by the exponential decay function (τ=30 minutes).

[0026] In this embodiment, the method for constructing the three-dimensional feature matrix based on mode reflectivity, terrain, and time lag includes: Acquire pattern reflection field data, terrain data, and beam blocking rate data, and perform standardization processing to extract terrain features, including terrain factors and radar beam influence factors. Terrain factors include regional terrain type and relative height. Regional terrain type is classified based on slope threshold and encoded as a one-hot vector. Relative height is calculated by determining the standard deviation of terrain height within a neighborhood window. Radar beam influence factors include effective beam height and blocking correction coefficient. The expression is: ; in The elevation of the radar station. The height of the beam center. The effective beam height; ,in This is the blocking correction factor. To determine the beam blocking rate, a Gaussian distribution mode is used to simulate beam energy attenuation. The reflectivity values ​​of the mode at three time lags are extracted from the radar grid points. The reflectivity change rate is calculated to capture the evolution trend of echo intensity and obtain time-varying characteristics. A sinusoidal function is introduced to encode the time periodicity to obtain the diurnal periodicity variation factor. Time feature alignment is performed based on the time-varying characteristics and the diurnal periodicity variation factor. The expression is: ; ; in The current mode reflectivity value. The reflectance values ​​for the first 6 minutes of the model. The reflectance values ​​for the first 12 minutes of the model. The rate of change of reflectance is 1. The rate of change of reflectance is 2; ,in This is a daily cycle variation factor; The three-dimensional feature matrix includes a time dimension, a spatial dimension, and a feature dimension. The spatial dimension is the latitude and longitude range of the radar grid. The time dimension is the 0-12 minute lead time window of the model forecast corresponding to three time lags. The feature dimension contains 12 feature variables for each spatiotemporal region. The matrix is ​​concatenated in the order of space, time, and features and output as a three-dimensional feature matrix. / / The 12 feature variables are respectively... , , , , Slope threshold, terrain height, terrain height standard deviation, beam blocking rate, blocking correction factor, diurnal variation factor, and terrain type coding; In the actual assessment, the three-dimensional feature matrix was generated by extracting 12 feature variables. Terrain type was encoded into a 3D one-hot vector based on slope thresholds (<5° for plains, 5°-15° for hills, >15° for mountains). Temporal features were... The daily cycle changes are encoded, and a three-dimensional feature matrix with dimensions of 800×1000×3×12 is finally generated (800×1000 is the latitude and longitude range of the radar grid).

[0027] In this embodiment, the method for performing two-stage attenuation correction on the systematic bias and the random error includes: After marking clutter areas based on the skewness-kurtosis joint identification method, effective echoes from the upper elevation angle are used to fill the clutter area and the scanning time of radars within the network is synchronously calibrated; among them, meteorological echo skewness > 0 and kurtosis > 3, and ground clutter skewness < 0 and kurtosis < 2. Using a stable S-band radar within the network as a reference, the mean reflectivity deviation between the reference radar and other radars is calculated. During clear, precipitation-free periods, the reflectivity difference between the reference radar and the target radar is statistically analyzed. The radome attenuation coefficient is fitted, and a first-order correction is applied to the target radar reflectivity. ,in The attenuation coefficient of the radome is... The target radar reflectivity after first-order correction. The original reflectivity is given; the expression for the radome attenuation coefficient is: ; in , Atmospheric attenuation coefficient, To detect distance, This represents the average reflectance deviation. Beam blocking rate is calculated based on terrain data, and a lookup table based on terrain height and blocking rate is constructed. Attenuation is calculated segment by segment. When the beam blocking rate is less than or equal to 0.2: When the beam blocking ratio is greater than 0.2 and less than or equal to 0.5: When the beam blocking ratio is greater than 0.5: ;in This is the attenuation amount. Beam blocking ratio; For regions with abrupt changes in blocking rate, gridded Gaussian smoothing is applied. Based on the first-order corrected target radar reflectivity and attenuation, a first-order corrected reflectivity field is output, expressed as:

[0028] in This is the first-order corrected reflectivity field; Based on the first-order corrected reflectivity field, the objective function for random error optimization is defined as follows:

[0029] in Let the first-order corrected reflectivity field be the k-th spatial grid point. Let k be the reference radar observation reflectivity of the k-th spatial grid point. This is the proportionality coefficient. This is the offset. For parameters to be optimized, For L2 regularization parameters, Let the objective function be the parameter to be optimized. This represents the total number of spatial grids. The gradient is calculated using mini-batch gradient descent, and the coefficients to be optimized are dynamically updated using a cosine annealing strategy. The process is iterated until the change in the objective function is less than the error threshold for five consecutive iterations or the maximum number of iterations reaches 500. Then the iteration stops and the optimized parameters are output. Second-order correction is performed based on the optimized parameters, expressed as follows:

[0030] in The optimized scaling factor, This is the optimized offset. This is the second-order corrected reflectivity field; For the strong echo region, nonlinear fine-tuning is performed, and an intensity correction factor is introduced, expressed as follows:

[0031] in For the second-order corrected reflectivity field in the strong echo region, This is the intensity correction factor; In practical assessments, clear skies with no precipitation are defined as reflectivity < 5 dBZ; beam blocking ratio is calculated using 30-second resolution terrain data; the beam blocking ratio ranges from 0 to 1, with 1 indicating complete blocking; the grid in the gridded Gaussian smoothing is 5×5, σ = 1.5; the error threshold is 10. -4 , First-order correction: Clutter Removal and Scan Calibration: Ground clutter areas (skewness < 0, kurtosis < 2) are marked based on the skewness-kurtosis joint identification method, and effective echoes at a 3° upper elevation angle are used to fill the clutter area; the scanning time of the networked radar is synchronized, and the maximum time difference is controlled within 30 seconds; Antenna attenuation correction: During clear, precipitation-free periods (reflectivity < 5 dBZ), the reflectivity difference between the reference radar and the target radar is statistically analyzed, and the radome attenuation coefficient is fitted to obtain the coefficient. The average value for C-band radar is 1.2 dB, and the average value for S-band radar is 0.8 dB. The coefficient is then calculated using the formula... Complete the correction; Beam blocking attenuation correction: The beam blocking ratio (BBR) is calculated based on terrain data, and the attenuation is calculated piecewise using a lookup table. For radar grid points in mountainous areas with a BBR > 0.5, the attenuation can reach a maximum of 8.5 dB. A 5×5 grid Gaussian smoothing (σ = 1.5) is used to process regions with abrupt changes in the blocking ratio, and a first-order corrected reflectivity field is output. Second-order correction: Objective function optimization: The L2 regularization parameter is set to 0.001, and the gradient is calculated using mini-batch gradient descent (batch size = 64). The parameters are dynamically updated using a cosine annealing strategy. After 420 iterations, the change in the objective function is less than 10. -4 Stop iteration; Fine-tuning in strong echo regions: For strong echo regions with reflectivity > 45 dBZ, an intensity correction factor is introduced. The final output is a second-order corrected reflectivity field, with a correction amplitude of 1.5-3.0 dB in the strong echo region.

[0032] In this embodiment, the method for constructing time-varying coefficient equations by introducing terrain correction factors includes:

[0033]

[0034] in For terrain correction factors, For terrain height, For final reflectivity, For the second-order corrected reflectivity field, For beam blocking variables, For the residual term, Let be the model prediction weight coefficient for the t-th time period. The terrain factor weight coefficient is the value at time t. The blocking factor weight coefficient is the value of the time-th time step. In practical assessments, the time-varying coefficient equation is constructed using the terrain correction factor. Where H is the terrain height (unit: m); nonlinear weights are set according to the beam blocking variable Q, with weights of 0.1 when Q < 0.2, 0.4 when 0.2 ≤ Q < 0.5, and 0.8 when Q ≥ 0.5.

[0035] In this embodiment, the method for constructing a spatiotemporal adaptive correction quantization model based on the time-varying coefficient equation includes: Using the final reflectance in the time-varying coefficient equation as the objective function of the spatiotemporal adaptive correction quantization model, the terrain height is compressed to the [0,1] interval using the hyperbolic tangent function. Nonlinear weights are set according to the beam-blocking variable: when the beam-blocking variable is less than 0.2, the nonlinear weight is set to 0.1; when the beam-blocking variable is greater than or equal to 0.2 and less than 0.5, the nonlinear weight is set to 0.4; when the beam-blocking variable is less than or equal to 0.5, the nonlinear weight is set to 0.8. The expression for compressed terrain height is: ; in This is the compressed terrain height, which is... ; Introducing the zero-order optimization concept of the Zebra Optimization algorithm, the gradient is estimated through two forward propagations. The first forward propagation calculates the loss based on the current coefficients. The second forward propagation adds random perturbations to the model forecast weight coefficients, and the calculation... ; Construct a weather type and coefficient matrix mapping library to store the optimal coefficient combinations for typical weather events such as rainstorms, typhoons, and cold waves. Identify the current weather type through KL divergence and call on historical knowledge to accelerate convergence. In the actual evaluation, the zero-order optimization idea of ​​the Zebra optimization algorithm was introduced, with a random perturbation amplitude of 0.001. A weather type coefficient mapping library was constructed to store the optimal coefficient combinations for rainstorms (Z>50dBZ), typhoons (wind speed>12), and cold waves (24-hour temperature drop>10℃). The current weather type was identified by KL divergence, and the convergence speed was improved by 40% after calling historical coefficients.

[0036] In this embodiment, the method for obtaining an enhanced feature map by extracting features through multi-scale dilated convolution to generate a spatial weight map includes: Multi-scale dilated convolution comprises an input layer, a feature extraction module, a weight generation module, and a feature enhancement module. The feature extraction module consists of three parallel dilated convolution branches: branch 1 extracts detail feature maps, branch 2 extracts mesoscale feature maps, and branch 3 extracts macroscale feature maps. The input layer receives the spatiotemporally matched pattern reflectivity field and normalizes its size. The feature extraction module consists of three parallel dilated convolution branches with different dilation rates (d=1,3,5), corresponding to receptive field sizes of 3×3, 7×7, and 11×11. The weight generation module calculates the feature importance of each branch using a channel attention mechanism, generating a 512×512×1 spatial weight map. The feature enhancement module performs element-wise multiplication of the weight map with the original feature map to highlight the feature responses of strongly weighted regions. A fusion coefficient is introduced, and the weights of each branch are dynamically adjusted through backpropagation. The feature maps of the three branches are then concatenated along the channel dimension using the branch weights to reduce the dimensionality and obtain the fused feature map. Global average pooling is applied to the fused feature map to generate channel descriptors. Channel weights are learned through a two-layer fully connected network. The channel weights are multiplied by the fused feature map, summed along the channel dimension, and then normalized using the S-function to generate a spatial weight map, expressed as: ; ; ; in For the spatial weight map of row a and column u, For spatial row index, For spatial column indexes, Let c be the weight of the channel. For channel index, For the c-th channel descriptor, The weights are for the first fully connected layer. These are the weights for the second fully connected layer. It is the ReLU activation function. It is the Sigmoid activation function. For the space is high, For the space is wide; The spatial weight map is multiplied pixel-by-pixel with the original reflectivity feature map to generate the enhanced feature map, expressed as: ; in This is the original reflectivity feature map. To enhance the feature map; In the actual evaluation, the feature map was 64×512×512×3. The stitched feature map was then reduced to 64 channels using a 1×1 convolution. The multi-scale dilated convolution branches have dilation rates of d=1, 3, and 5, corresponding to receptive fields of 3×3, 7×7, and 11×11, respectively. The fused feature map is reduced to 64 channels by 1×1 convolution, and channel weights are learned through a two-layer fully connected network (128 hidden nodes) to generate a 512×512×1 spatial weight map. Finally, the weights are calculated according to the formula... An enhanced feature map is obtained.

[0037] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method for the correction of radar composite reflectivity based on a quantitative model, characterized in that, The method comprises the following steps: Collecting observation data of a preset radar network, and preprocessing the observation data; the observation data comprises mode data, radar data and auxiliary data; A regional detection algorithm is used to eliminate system bias and attenuation effects, a grid search frame is established based on an S-band radar, a bias between network reflectivity mean and reference mean is calculated according to the observation data, when the difference value exceeds a bias threshold, a correction is performed, error decomposition is performed on the bias to obtain system bias and random error; For a specific resolution reflectivity field, a bilinear interpolation and a time-space sliding window are matched to a radar grid, terrain height and beam blocking rate are synchronously extracted, and a three-dimensional feature matrix based on mode reflectivity, terrain and time lag is constructed; A network attenuation algorithm and a gradient descent method are used to perform two-stage attenuation correction on the system bias and the random error, a terrain correction factor is introduced to construct a time-varying coefficient equation, and a space-time adaptive correction quantization model is constructed according to the time-varying coefficient equation; A multi-scale attention gate is introduced, a spatial weight map is generated by multi-scale hollow convolution to obtain an enhanced feature map, and the to-be-corrected data is input into the space-time adaptive correction quantization model based on the enhanced feature map, and a correction result is output.

2. The method of claim 1, wherein the radar composite reflectivity is based on a quantitative model. The method for calculating the bias between the network reflectivity mean and the reference mean comprises: The network radar reflectivity data is scanned in space by a grid sliding window, and the bias value between the observation mean in the window and the reference mean of the S-band radar is calculated: ； wherein is the mean of the networked reflectivity, is the mean of the S-band reference, is the raw observed reflectivity value, is the bias value of the observed mean from the S-band radar reference mean.

3. The method of claim 1, wherein the radar composite reflectivity is based on a quantitative model. The method for matching the bilinear interpolation and the time-space sliding window to the radar grid comprises: The latitude and longitude coordinates of the mode data are uniformly converted into Albers equal-area projection adopted by the radar data, the projection distortion error is eliminated through a seven-parameter Bursa model, a 20-minute assimilation time window is set, time interpolation is performed on the mode data, and a high-frequency data sequence with equal time intervals is generated; For any target point in the radar grid , calculate the row and column index of the target point in the pattern grid , determine the coordinates of the surrounding 4 pattern grid points, the 4 pattern grid points are the upper left corner , the upper right corner , the lower left corner , the lower right corner ; Bilinear interpolation is performed according to the weight coefficients in the horizontal direction and the weight coefficients in the vertical direction, and the expression is: ； ； ； wherein is a weight coefficient in the horizontal direction, is a weight coefficient in the vertical direction, is the horizontal coordinate of the target point, is the vertical coordinate of the target point, is the horizontal coordinate of the upper-left corner mode grid point, is the horizontal coordinate of the upper-right corner mode grid point, is the vertical coordinate of the upper-left corner mode grid point, is the vertical coordinate of the lower-left corner mode grid point, is the reflectivity value of the mode grid point in the i-th row and j-th column, is the preliminary interpolation result, is the reflectivity value of the mode grid point in the i+1-th row and j-th column; For radar grid points within the range of the pattern grid boundary, distance weighted average is used instead of pure bilinear interpolation, and the weight coefficient wherein is the distance from the target point to the pattern grid point; when the terrain height of the target point exceeds the corresponding height of the lowest layer of the pattern, vertical interpolation result correction is automatically enabled, and the correction coefficient ; The window size is set, a spatial weight matrix is generated according to a Gaussian kernel function, a time weight matrix is generated according to an exponential decay function, and a time-space sliding window is constructed according to the window size, the spatial weight matrix and the time weight matrix; The mean, standard deviation and gradient amplitude of the mode reflectivity in the window are calculated as matching degree evaluation indexes, the similarity between the mode prediction sequence and the radar observation sequence is calculated through a dynamic time warping algorithm, and the prediction time corresponding to the minimum distance is selected as the reference; The spatial weight, the time weight and the gradient amplitude are weighted and fused, and the expression is: ； wherein is the windowed mode reflectivity value for the kth spatial grid, the tth time instance, is the total number of spatial grids, is the number of time instances, is the fused reflectivity value matched to the radar grid point, is the spatial weight for the kth spatial grid, is the spatial grid index, is the time instance index, is the temporal weight for the tth time instance.

4. The method of claim 1, wherein the radar composite reflectivity is based on a quantitative model. The method for constructing the three-dimensional feature matrix based on the mode reflectivity, the terrain and the time lag comprises: Mode reflectivity field data, terrain data and beam blocking rate data are obtained and standardized, terrain features are extracted, and the terrain features include terrain factors and radar beam influence factors; the terrain factors include regional terrain types and relative heights; the regional terrain types are divided based on a slope threshold and are coded into a one-hot vector; the relative heights calculate the standard deviation of the terrain height in the neighborhood window; the radar beam influence factors include beam effective height and blocking correction coefficient; The radar grid point is extracted with three time lag mode reflectivity values, the reflectivity change rate is calculated, the echo intensity evolution trend is captured, the time-varying characteristics are obtained, the daily cycle change factor is obtained by introducing a sine function to code the time periodicity, and the time characteristics are aligned according to the time-varying characteristics and the daily cycle change factor; The three-dimensional feature matrix includes a time dimension, a space dimension and a feature dimension, the space dimension is the latitude and longitude range of the radar grid, the time dimension is a 0 to 12 minute time window corresponding to the mode prediction of the three time lags, and the feature dimension is 12 feature variables contained in each space-time, and the matrix is spliced in the order of space, time and feature and output as a three-dimensional feature matrix.

5. The method of claim 1, wherein the radar composite reflectivity is based on a quantitative model. The method for two-stage attenuation correction of the system bias and the random error comprises: After the clutter area is marked based on the joint skewness-kurtosis identification method, the effective echo of the upper elevation is filled, and the scanning time of the radars in the network is calibrated synchronously; The S-band radar with stable performance in networking is taken as a reference, the reflectivity deviation mean of the reference and other radars is calculated, the reflectivity difference between the reference radar and the target radar is counted in the clear sky and non-rainfall period, the radome attenuation coefficient is fitted, and the first-order correction is made on the reflectivity of the target radar: wherein is the radome attenuation coefficient, is the first-order corrected target radar reflectivity, is the original reflectivity; wherein , is the atmospheric attenuation coefficient, is the detection distance, is the reflectivity deviation mean; According to the terrain data, the beam block rate is calculated, a lookup table based on the terrain height and the block rate is constructed, and the attenuation amount is calculated in sections; when the beam block rate is less than or equal to 0.2: ; when the beam block rate is greater than 0.2 and less than or equal to 0.5: ; when the beam block rate is greater than 0.5: ; wherein is the attenuation amount, is the beam block rate; The first-order corrected target radar reflectivity and attenuation are used to output the first-order corrected reflectivity field in the region with sudden change of the blockage rate; The target function for random error optimization is defined according to the first-order corrected reflectivity field, and the expression is: ； wherein is the first order corrected reflectivity field at the kth spatial grid point, is the reference radar observed reflectivity at the kth spatial grid point, is a proportionality coefficient, is a bias, is a parameter to be optimized, is an L2 regularization parameter, is an objective function of the parameter to be optimized, is the total number of spatial grids; The gradient is calculated by using the small batch gradient descent, the to-be-optimized coefficients are dynamically updated by using the cosine annealing strategy, and the iteration is continuously performed until the variation of the target function is less than the error threshold for five consecutive iterations or the maximum number of iterations reaches 500, and then the iteration is stopped and the optimized parameters are output; The second-order correction is performed based on the optimized parameters, and the expression is: ； wherein is the optimized scale factor, is the optimized offset, is the second-order corrected reflectivity field, is the first-order corrected reflectivity field; The nonlinear fine tuning is performed on the strong echo area, and the intensity correction factor is introduced, and the expression is: ； wherein is the reflectivity field for the second order correction of the strong echo region, is the intensity correction factor.

6. The method of claim 1, wherein the radar composite reflectivity is based on a quantitative model. The method for constructing the time-varying coefficient equation by introducing the terrain correction factor comprises: ； ； wherein is a terrain correction factor, is a terrain height, is a final reflectivity, is a second order corrected reflectivity field, is a beam blockage variable, is a residual term, is a model prediction weight coefficient for the tth time instance, is a terrain factor weight coefficient for the tth time instance, is a blockage factor weight coefficient for the tth time instance.

7. The method of claim 1, wherein the radar composite reflectivity is based on a quantitative model. The method for constructing the spatio-temporal adaptive correction quantization model according to the time-varying coefficient equation comprises: The terminal corrected reflectivity in the time-varying coefficient equation is taken as the target function of the spatio-temporal adaptive correction quantization model, the terrain height is compressed to the [0, 1] interval through the hyperbolic tangent function, and the nonlinear weight is set according to the size of the beam blockage variable: when the beam blockage variable is less than 0.2, the nonlinear weight is set to 0.1; when the beam blockage variable is greater than or equal to 0.2 and less than 0.5, the nonlinear weight is set to 0.4; and when the beam blockage variable is less than or equal to 0.5, the nonlinear weight is set to 0.8; The zero-order optimization idea of the zebra optimization algorithm is introduced, and the gradient is estimated by twice forward propagation. The first forward propagation calculates the loss according to the current coefficient , and the second forward propagation adds random disturbance to the mode prediction weight coefficient and calculates ; The weather type and coefficient matrix mapping library is constructed, the optimal coefficient combination of typical weather such as rainstorm, typhoon and cold wave is stored, the current weather type is identified by KL divergence, and historical knowledge is called to accelerate convergence.

8. The method of claim 1, wherein the radar composite reflectivity is based on a quantitative model. The method for obtaining an enhanced feature map by extracting features through multi-scale hollow convolution comprises: The multi-scale hollow convolution comprises an input layer, a feature extraction module, a weight generation module and a feature enhancement module, the feature extraction module is composed of three parallel hollow convolution branches, branch 1 extracts a detailed feature map, branch 2 extracts a medium-scale feature map, and branch 3 extracts a macro feature map; The fusion coefficient is introduced, the weights of the branches are dynamically adjusted through back propagation, the feature maps of the three branches are spliced along the channel dimension through the branch weights, and the fusion feature map is obtained by dimension reduction. The global average pooling is performed on the fused feature map to generate a channel descriptor, a two-layer fully connected network is used to learn a channel weight, the channel weight is multiplied by the fused feature map, summation is performed along the channel dimension, and then a spatial weight map is normalized by using an S function, and the expression is as follows: ； ； ； wherein is a spatial weight map for the a-th row and u-th column, is a spatial row index, is a spatial column index, is a c-th channel weight, is a channel index, is a c-th channel descriptor, is a first layer fully connected layer weight, is a second layer fully connected layer weight, is a ReLU activation function, is a Sigmoid activation function, is a spatial height, is a spatial width; The spatial weight map is pixel-by-pixel multiplied by the original reflectivity feature map to generate an enhanced feature map, and the expression is as follows: The spatial weight map is pixel-by-pixel multiplied by the original reflectivity feature map to generate an enhanced feature map, and the expression is as follows: ； wherein is the original reflectance feature map, is the enhanced feature map.

Images