A high-resolution weather downscaling method fusing dynamic modes, terrain physical correction and station residual learning

Abstract

The application belongs to the technical field of high-resolution meteorological downscaling, and particularly relates to a high-resolution meteorological downscaling method fusing a dynamic model, terrain physical correction and station residual learning, a data preprocessing module, a high-resolution target grid construction module, a multi-layer terrain background and multi-scale terrain derived factor construction module, a dynamic background interpolation and water vapor budget decomposition module, a terrain-driven temperature pressure humidity wind field physical correction module, a first multivariate physical consistency closure module, a station residual model training module, a residual inference correction module, a second physical consistency closure module, a meteorological grid product output module; each module is executed in sequence along a directed acyclic call chain.

Description

Technical Field

[0001] This invention belongs to the field of high-resolution meteorological downscaling technology, and in particular relates to a high-resolution meteorological downscaling method that integrates dynamic models, topographic physical correction and station residual learning. Background Technology

[0002] Currently, high-resolution meteorological downscaling mainly falls into two categories: dynamic downscaling and statistical downscaling. Traditional dynamic downscaling relies on numerical models to directly refine the grid, which is limited by computational resources and makes it difficult to achieve fine output at the hundred-meter level. Furthermore, the model's built-in terrain data is only the average elevation of the coarse grid, which cannot characterize micro-topographic undulations such as canyons, ridges, and slopes. Near-surface temperature, pressure, humidity, and wind fields are difficult to match the real local topographic thermal and dynamic characteristics. Conventional statistical downscaling often uses linear regression and ordinary machine learning to fit the relationship between meteorological observations and the large-scale background. It simply introduces a single topographic factor of elevation and lacks physical priors such as multi-scale topographic morphology and dynamic forcing of water vapor budget. The model is prone to learning spurious correlations without physical meaning, and the extrapolation error is significant in areas with drastic topographic undulations and sparse stations. At the same time, both methods lack physical consistency verification of multi-variable linkages, and the output results often show non-physical numerical anomalies such as oversaturation, static imbalance, and wind speed exceeding limits.

[0003] Most existing fusion-based downscaling schemes simply overlay dynamic simulation results with statistical correction residuals, failing to construct a complete hierarchical physical correction link. They lack refined physical diagnostic modules such as segmented iterative correction of elevation differences, multi-scale topographic derived factors, and six-component decomposition of water vapor budget. The feature system in the residual modeling stage is singular, failing to distinguish between multiple feature clusters such as topography, underlying surface, dynamics, water vapor, time, and reliability. Furthermore, they do not set correction weight attenuation mechanisms for topographic extrapolation areas and sparse station areas, and only perform simple range pruning once before and after statistical correction. They lack two hierarchical physical closure constraints to coordinate the balance of multiple variables such as thermodynamics, dynamics, and water vapor. The final product has insufficient spatiotemporal continuity and topographic adaptability, making it difficult to meet the high-precision gridded meteorological data requirements of complex terrain scenarios such as mountain hydrology, wind energy assessment, and refined meteorological services. Summary of the Invention

[0004] In view of the aforementioned problems, and in conjunction with the first aspect of the present invention, embodiments of the present invention provide a high-resolution meteorological downscaling method that integrates dynamic models, topographic physical correction, and station residual learning. The method includes: a data preprocessing module, a high-resolution target grid construction module, a multi-layer topographic background and multi-scale topographic derived factor construction module, a dynamic background interpolation and water vapor budget decomposition module, a topographic-driven physical correction module for temperature, pressure, humidity, and wind fields, a first multivariate physical consistency closure module, a station residual model training module, a residual inference correction module, a second physical consistency closure module, and a meteorological grid product output module. Each module is executed sequentially along a directed acyclic call chain, and the implementation flow is as follows: Step 1: Multi-source meteorological and geographic data reading, coordinate unification, spatiotemporal matching and basic quality control. Input data includes regional dynamic model output, high-resolution DEM, underlying static auxiliary data, and optional ground station observation data. Step 2: Generate a standard high-resolution target mesh based on the dynamic mode spatial range and the preset target resolution; Step 3: Based on the DEM, construct a high-resolution real topographic field and a dynamic model equivalent coarse topographic field, calculate the elevation difference between the two, and generate single-scale and multi-scale topographic morphological features and regional underlying surface auxiliary features. Step 4: Interpolate the near-surface meteorological variables of the dynamic model to the target grid to obtain the dynamic background field; use the three-dimensional temperature, pressure, wind, and water vapor fields to complete the six-component physical decomposition of the entire layer's water vapor budget and generate dynamic meteorological feature clusters; Step 5: Based on the terrain elevation difference, perform stepwise thermodynamic and physical correction of near-surface air temperature, surface air pressure, and specific humidity. Combine slope, aspect, curvature, and terrain blocking effect to complete the 10m wind field diagnosis and correction, and obtain a purely physical downscaling meteorological field. Step 6: Perform the first global multivariate physical consistency closure constraint on the purely physical downscaling field to eliminate non-physical states such as numerical anomalies, violations of hydrostatic equilibrium, water vapor supersaturation, and wind speed exceeding limits, and output a physically consistent meteorological state field that can be used for machine learning training; Step 7: If valid station observation data exists, match the station observations with the corresponding grid point physical downscaling results to construct meteorological residual labels, and train the variate residual regression model by integrating multiple feature clusters such as topography, region, dynamics, thermodynamics, water vapor, time, and confidence; if there is no station data, skip this step and directly output the pure physical downscaling product. Step 8: In the inference stage, the trained residual model is called to predict the grid residuals of each meteorological variable, limit the amplitude, and then superimpose them onto the physical downscaling field to obtain the hybrid-statistical correction meteorological field. Step 9: Apply a second multivariate physical consistency closure constraint to the corrected meteorological field to slightly correct the physical imbalance problem introduced by the statistical residuals; Step 10: Output the high-resolution temperature, pressure, humidity, and style point meteorological fields, as well as topography and quality control markers, which have undergone two physical closure constraints, as standardized meteorological files.

[0005] Preferably, the multi-source data preprocessing in step 1 specifically includes: reading the output files of WRF regional dynamic models, reanalysis, or climate models, and extracting grid latitude and longitude, three-dimensional geopotential height, disturbed temperature and pressure, water vapor, near-surface 2m temperature and humidity, 10m wind, surface air pressure, precipitation, and surface flux variables; reading the high-resolution DEM digital elevation model; reading static auxiliary data such as land use, vegetation cover, snow cover, coastline, and administrative boundaries; optionally reading station table data containing time, latitude and longitude, elevation, temperature, pressure, humidity, and wind observations; unifying the coordinate reference system of all data, using conservative averaging / bilinear resampling for continuous geographic variables, and using mode resampling for categorical underlying surface variables; completing the unit conversion, time alignment, and missing value identification of meteorological variables, and simultaneously generating data credibility marker clusters; when performing water vapor budget decomposition, verifying the completeness of three-dimensional wind, vertical velocity, water vapor mixing ratio, and pressure layer thickness, calculating only solvable decomposition items for missing variables and updating credibility markers, and prohibiting zero-value filling of missing physical quantities.

[0006] Preferably, the high-resolution target mesh construction method in step 2 is as follows: the spatial boundary is determined by the extreme values ​​of latitude and longitude of the original mesh in the dynamic mode, and the latitude cosine correction of the regional center is introduced in combination with the geometric relationship of the earth. The latitudinal and longitudinal mesh intervals are calculated respectively. Based on the mesh interval and latitude and longitude boundaries, the total number of north-south and east-west grid points of the target mesh is calculated to generate a one-dimensional latitude and longitude array and a two-dimensional target latitude and longitude mesh. A block-based lazy loading strategy is adopted to divide the target mesh into 256×256 or 512×512 blocks, and an overlap buffer is set between the blocks. After the calculation is completed, the overlapping area is cut off and the entire field is stitched together to reduce memory usage.

[0007] Preferably, the multi-layer terrain background construction process in step 3 includes: cropping the DEM, transforming its coordinates, and interpolating it to the target mesh to obtain a high-resolution terrain field. Using a coarse grid of dynamic model as the unit, the average elevation of all high-resolution grid points within the unit is aggregated to obtain the equivalent coarse terrain. Then, the coarse terrain is interpolated to the target grid; the elevation difference is calculated point by point. ; Set a fixed elevation integration step size for... Temperature, virtual temperature, and air pressure are updated iteratively in segments for grid points exceeding the threshold. When the absolute value of the elevation difference exceeds the 99th percentile of the training samples or 1500m, the terrain extrapolation area is marked, and the correction weight of the stations in that area is reduced in the subsequent residual learning stage.

[0008] Preferably, the terrain and regional auxiliary features generated synchronously in step 3 include: single-scale basic terrain factors: Slope, aspect, and topographic curvature are calculated using the central difference / Sobel operator; multi-scale topographic morphology factors: multiple physical distance neighborhoods are set up to calculate topographic relief, topographic location index (TPI), topographic roughness (TRI), sky visibility factor (SVF), and windward shading index along the incoming flow direction at each scale; spatial distance factors: Calculate the shortest planar distance to the coastline, water body, ridge, valley, and study area boundary point by grid point; Underlying surface and station representativeness factors: land use unique thermal coding, vegetation / impermeable surface / snow cover ratio, and roughness length; construct station spatial representativeness weights to characterize the horizontal and vertical representativeness of observation stations around grid points.

[0009] Preferably, the six-element decomposition of water vapor budget in step 4 is implemented by performing area low-pass aggregation on the coarse grid field of the dynamic mode to obtain a large-scale background field, and the difference between the original field and the background field is defined as the local disturbance. The entire layer water vapor flux convergence (MFC) is solved by integrating three-dimensional specific humidity, horizontal wind field, and pressure layer thickness, and decomposed into dynamic advection term, dynamic convergence term, thermal advection term, thermal convergence term, nonlinear advection term, and nonlinear convergence term. The sum of the six decomposition terms equals the water vapor convergence increment of the high-resolution field relative to the coarse background field; Simultaneously calculate the total precipitable water and water vapor collection pillar residuals, and use all decomposed quantities as dynamic input features for residual learning.

[0010] Preferably, the physical correction process for temperature, air pressure, and humidity in step 5 is as follows: The height, actual air temperature, and virtual temperature of each layer are obtained from the multi-layer profile of the dynamic model. The vertical lapse rate Γ of the virtual temperature in the near-surface layer is extracted and a reasonable range is defined. 2m air temperature is corrected based on piecewise iterative correction of elevation difference; The average virtual temperature is calculated by taking the average temperature before and after the correction, and the surface pressure topography is corrected by combining the hydrostatic index equation. The original relative humidity is calculated from the interpolated background temperature, pressure and specific humidity. Under the corrected temperature and pressure conditions, the saturated vapor pressure and actual vapor pressure are recalculated using the Bolton saturated vapor pressure formula. The specific humidity after topographic correction is solved in reverse and a saturation constraint is applied to limit the specific humidity to not exceed the saturated specific humidity and eliminate the supersaturated state.

[0011] Preferably, step 5, the 10m wind field correction step, includes: Calculate wind speed and meteorological wind direction from interpolated wind components; calculate slope acceleration / deceleration factor based on the angle between wind direction and slope aspect. Valley passage areas are identified by slope and curvature, and valley wind deflection is calculated. A simplified Froude number is introduced to characterize the topographic airflow blocking effect, and the topographic blocking factor is solved. After the wind speed is corrected by coupling the slope factor and the blocking factor, the wind direction deflection is superimposed to re-decompose the east-west and north-south wind components, and a global upper limit constraint on the wind speed is set. Optionally, a weakly weighted ground-to-wind-trend hybrid correction can be introduced to modify the large-scale wind field background.

[0012] Preferably, the first physical consistency closure in step 6 includes three layers of constraints: Hard constraints: Eliminate infinite and missing values, limit the global physical range of temperature, pressure, humidity, and wind speed, and forcibly constrain relative humidity to 0-100% to eliminate oversaturation; Differentiate between formulas for saturated vapor pressure on ice / water surfaces, and empirical formulas for switching ice surfaces at sub-zero temperatures; Soft constraints: Construct cost functions for static equilibrium residuals and column water vapor balance residuals, and introduce multivariate weights to balance various physical constraints; Iterative projection optimization: The temperature, pressure, humidity and wind state vectors are updated iteratively using the projection gradient method until the change in variables is less than the convergence threshold or the maximum number of iterations is reached; Unconverged grid points retain their valid physical states and are marked as closure anomalies.

[0013] Preferably, the residual model training method in step 7 includes: The site observations are matched to the corresponding target grid points, and the observations minus the physical downscaling results of the same grid point are used as the residual labels. Seven feature clusters are constructed: topographic morphology cluster, regional underlying surface cluster, dynamic background cluster, thermal and water vapor cluster, water vapor budget decomposition cluster, time period cluster, and data credibility cluster, which are then concatenated to form a complete input feature vector. Continuous features are processed using robust median-interquartile range standardization. Training, validation, and test sets are divided by site / time block, and random splitting of consecutive time blocks at the same site is prohibited. Regression models were trained independently for 2m air temperature, ground air pressure, specific humidity, and 10m wind component. Random forest, gradient boosting tree, or neural network were selected for the models. After training, the feature order, standardized parameters, and residual amplitude threshold were saved.

[0014] Based on the above, we first rely on multi-layer topographic fields and multi-scale topographic derived factors to complete the stepwise thermodynamic and kinetic physical correction of temperature, pressure, humidity and wind. Combined with the six-component decomposition of the whole layer water vapor budget to supplement the dynamic meteorological forcing characteristics, we restore the modulation effect of micro-topography on local meteorological elements from the physical mechanism level, which greatly alleviates the distortion problem of local meteorological fields caused by traditional dynamic downscaling of coarse topography. The first physical consistency closure eliminates non-physical states such as oversaturation, static imbalance, and wind speed exceeding limits, providing a physically reasonable reference field for subsequent residual learning, avoiding the machine learning to fit false residual laws. When there are no station observations, reliable pure physical downscaling products can be directly output, which is suitable for remote mountainous and plateau areas without observation data.

[0015] A discrete-variable residual regression model is introduced, and seven complete feature term clusters are constructed to carry out station residual correction. At the same time, the correction weights are automatically attenuated for areas with excessive topographic elevation differences and insufficient station representativeness to reduce model extrapolation bias. After superimposing residual predictions, a second physical closure constraint is performed to slightly correct the multivariate imbalance problem introduced by statistical correction. It takes into account the advantages of observation data correction and meteorological physical conservation laws. The output 100-meter grid meteorological field has significantly improved the consistency between temperature, air pressure, humidity, 10m wind field and station observations in complex landforms such as mountains, canyons and basins, and has stronger spatiotemporal continuity. It can directly support a variety of business applications such as hydrological driving, wind energy resource assessment, refined weather forecasting, and disaster risk simulation. Attached Figure Description

[0016] Figure 1 This is a diagram illustrating the overall technical architecture of a high-resolution meteorological downscaling method that integrates dynamical models, topographic and regional auxiliary information, station residual learning, and physical closure constraints, as provided in this embodiment of the invention. Figure 2 This is a schematic diagram of a high-resolution meteorological downscaling method that integrates dynamic models, topographic physical correction, and station residual learning, provided by an embodiment of the present invention, which includes topographic and regional auxiliary information preprocessing and the composition of seven types of item cluster features; Figure 3 This is a schematic diagram of a high-resolution meteorological downscaling method that integrates dynamic models, topographic physical correction, and station residual learning, which includes six-term decomposition of water vapor budget, hierarchical integration, and column water vapor diagnosis. Figure 4 This invention provides a high-resolution meteorological downscaling method that integrates dynamic models, topographic physical correction, and station residual learning. The method involves the synergistic relationship between the first physical closure, station residual learning, statistical correction, and the second physical closure. Figure 5 This is a specific application embodiment of the high-resolution meteorological downscaling method that integrates dynamic models, topographic physical correction, and station residual learning provided by the present invention. Figure 6 This is a comparison of the results before and after the patented process in a specific application example of a high-resolution meteorological downscaling method that integrates dynamic models, terrain physical correction, and station residual learning provided by an embodiment of the present invention. Figure 7 This invention provides a specific application example of a high-resolution meteorological downscaling method that integrates dynamic models, topographic physical correction, and station residual learning, showing the spatial distribution of topographic forcing and temperature and pressure responses. Figure 8This invention provides a specific application example of a high-resolution meteorological downscaling method that integrates dynamic models, terrain physical correction, and station residual learning, including physical closure and quality control diagnosis. Figure 9 This is a specific application example of the 24-hour daily variation of physical quantities in a high-resolution meteorological downscaling method that integrates dynamic models, topographic physical correction, and station residual learning, provided by an embodiment of the present invention. Figure 10 This is a specific application example of the high-resolution meteorological downscaling method that integrates dynamic models, topographic physical correction, and station residual learning provided in this invention, which describes the elevation zone physical response. Detailed Implementation

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

[0018] like Figures 1 to 10 As shown, this invention provides a high-resolution meteorological downscaling method that integrates dynamic models, topographic physical correction, and station residual learning. This method can be executed by an electronic computing device, which can be a server, workstation, high-performance computing node, or business platform with data processing capabilities. The device includes at least a processor, a memory, and program instructions stored in the memory. These program instructions execute the method of this invention when run on the processor.

[0019] In system implementation, each module executes according to a directed acyclic call chain. The data acquisition and configuration module first completes file indexing, unit identification, time alignment, and quality control; the target grid and terrain background module generates a unified target grid, high-resolution terrain, and equivalent coarse terrain; the terrain and regional auxiliary information module generates static term families; the dynamic background interpolation and water vapor budget decomposition module generates hourly dynamic term families; the temperature, pressure, and humidity correction and wind field diagnosis module forms a purely physical downscaling field; the first physical closure constraint module outputs a physically consistent state that can be used for training; the site residual learning module trains the model by variables and saves the feature order, standardized parameters, and applicable range; during inference, the statistical residual correction module reads the same feature pattern for prediction, the second physical closure constraint module reprojects the final state, and finally, the output module writes out the grid product and quality control markers. If any optional module is missing, the system maintains the purely physical result without interrupting the main process according to preset rollback rules.

[0020] The system can be deployed on a single server, workstation, high-performance computing node, or distributed computing platform. The CPU is responsible for reading and writing NetCDF, GeoTIFF, and site tables, resampling, and physical diagnostics. The GPU is only used as an optional acceleration device when a neural network residual model is selected. To control memory usage, 3D variables are read lazily by time interval, by vertical layer, and by spatial block. The target mesh is processed as 256×256 or 512×512 pixel blocks with overlapping boundaries. The overlap width between blocks is not less than the maximum neighborhood radius or the difference template radius. After closure, the overlapping area is trimmed and stitched together, thus avoiding loading all spatiotemporal data at once.

[0021] Data acquisition and configuration module; Obtain output files from dynamical models such as WRF, high-resolution DEM files, and optional site observation files. The dynamical model output files preferably include latitude and longitude grids XLAT and XLONG, near-surface variables U10, V10, T2, Q2, PSFC, and three-dimensional geopotential height, perturbation geopotential height, perturbation temperature, air pressure, and water vapor variables. Site observation files may include fields such as time, lat, lon, elevation, T2, Q2, PSFC, U10, V10, and RH2.

[0022] In one implementation, the system determines the input / output paths, target date, target hours, target resolution, original resolution of the dynamic mode, statistical residual model path, output NetCDF path, and physical constraint thresholds through a configuration file or default configuration. The target resolution can be set to 300m, the original resolution of the dynamic mode can be set to approximately 3000m, and the processing period can be set to 24 hours of the target date.

[0023] When performing the six-item decomposition of the water vapor budget, the dynamic model file should preferably also include three-dimensional horizontal wind U and V, vertical velocity or vertical pressure velocity, three-dimensional water vapor mixing ratio or specific humidity, pressure thickness or recoverable pressure P and PB of each layer, geopotential height PH and PHB, as well as cumulative precipitation and surface latent heat flux or evaporation. When some variables are missing, only the decomposition items with available data can be calculated, and the missing indicators should be added to the data confidence item cluster. It is prohibited to replace missing physical quantities with zero values.

[0024] Surface and regional ancillary data may consist of DEM, land use or land cover, vegetation cover, surface roughness, albedo, snow cover, water bodies, coastlines, watershed boundaries, administrative boundaries, or climate zone data. All static data are converted to a uniform coordinate reference system before being entered into the model. Continuous variables are resampled using conservative averaging, bilinear, or spline methods, while categorical variables are resampled using mode or nearest neighbor methods. Source resolution, interpolation method, and effective pixel ratio are retained as quality control metadata.

[0025] High-resolution target mesh building module; The system reads the latitude and longitude range from the dynamic model reference file, converts the metric distance into latitude and longitude intervals based on the target resolution, and generates a one-dimensional array of target latitude, a one-dimensional array of target longitude, and a two-dimensional target latitude and longitude grid. The longitude intervals can be cosine corrected based on the center latitude of the study area.

[0026] (1); In equation (1), Indicates the target horizontal resolution. This represents the approximate distance corresponding to each degree in the latitudinal direction. This represents the approximate distance per degree of longitude near the equator. Indicates the latitude of the center of the study area. and These are the latitudinal and longitudinal intervals of the target grid, respectively.

[0027] (2); In equation (2), and These represent the number of north-south grid points and the number of east-west grid points in the target's high-resolution mesh, respectively. These represent the latitude and longitude boundaries of the area covered by the dynamic mode, respectively.

[0028] Terrain background construction module; The DEM data was cropped based on the coverage of the dynamic model, and the DEM coordinate system was converted to the geographic coordinate system. Then, linear interpolation was used to resample the DEM elevation to the target high-resolution grid; when linear interpolation produced null values, nearest neighbor interpolation was used to fill in the missing values.

[0029] (3); In equation (3), Represents the target high-resolution mesh. DEM elevation at each grid point This represents the spatial interpolation operator from DEM source data to the target grid.

[0030] After obtaining the high-resolution topographic field, the system spatially aggregates the high-resolution DEM according to the dynamic model grid scale to obtain a coarse topographic background corresponding to the dynamic model grid, and then interpolates the coarse topographic background to the target high-resolution grid.

[0031] (4); In equation (4), The first coarse grid representing the dynamic mode The equivalent coarse terrain elevation corresponding to each grid point This represents the set of high-resolution grid points of the target area covered by the coarse grid. This indicates the number of valid grid points in the set.

[0032] (5); In equation (5), This represents the elevation difference between the high-resolution real terrain and the equivalent coarse terrain in the dynamic model. This represents the result after interpolating the coarse terrain background to the target high-resolution mesh. This terrain interpolation is the core terrain forcing factor for subsequent temperature, air pressure, humidity, and wind field corrections.

[0033] The elevation difference uses the notation convention of "real terrain minus equivalent coarse terrain". When ΔH is greater than 0, the target grid point is located on an unresolved ridge, plateau, or high part of a slope. Temperature correction generally decreases along the local lapse rate, while surface pressure decreases according to the hydrostatic relationship. When ΔH is less than 0, the target grid point is located on an unresolved valley, basin, or lowland. Temperature and pressure corrections are in opposite directions. Humidity does not increase or decrease linearly with the elevation difference. Instead, it maintains the physical meaning of background relative humidity or dew point difference under new temperature and pressure conditions and applies saturation constraints again.

[0034] (5-1); In equation (5-1), NH represents the number of vertical segments, and ΔHstep represents the step size for a single height integration. For grid points with large absolute elevation differences, such as plateaus or deep valleys, the system uses... The temperature, virtual temperature, and air pressure are updated iteratively with a step size, thus preserving the physical effects of all elevation differences and avoiding numerical overflow caused by a single exponential correction. The range can be 200–500 m, with a typical value of 300 m; when When the elevation difference exceeds the 99th percentile of the absolute value of the training sample elevation difference or exceeds 1500m, a terrain extrapolation marker is set, and the station correction weight is reduced during the statistical residual stage.

[0035] Topography-derived factor calculation module; like Figure 2 As shown, slope, aspect, and curvature are calculated based on a high-resolution DEM. Slope is determined by the elevation gradients in both the east-west and north-south directions, aspect is obtained by transforming the gradient direction, and curvature is calculated using the second-order difference of elevation. The topographically derived factors are used for wind field diagnosis and correction, and also serve as important input features for the site residual learning model.

[0036] (6); In equation (6), Indicates the slope angle. and These represent the horizontal elevation gradients in the east-west and north-south directions, respectively. In actual calculations, these gradients can be estimated using the Sobel operator or central difference method.

[0037] (7); In equation (7), Indicates the slope angle. Represents the arctangent function in the four quadrants. This represents the modulo operation, used to limit the aspect angle to a specific value. to Within the range.

[0038] (8); In equation (8), It represents the curvature of the terrain. Curvature can be used to distinguish between concave valleys, convex ridges, and gentle areas, and can help determine the valley passage effect or terrain barrier effect.

[0039] Topography and regional auxiliary information construction module; To enable the model to identify ridges, valleys, open slopes, canyon passages, and coastal transition zones, the system constructs multi-scale terrain morphology quantities for each target grid point i and its neighborhood Ωr(i), in addition to slope, aspect, and curvature. The neighborhood scale r can be set to 1, 3, 9, and 15 target grids according to the target resolution, or to approximately 0.3 km, 1 km, 3 km, and 5 km; maintaining approximately consistent physical distances across different resolutions.

[0040] (8-1); In equation (8-1), Rr is the topographic relief at scale r, which is used to represent the maximum elevation difference in the neighborhood.

[0041] (8-2); In Equation (8-2), TPIr is the terrain location index. A positive value indicates that the grid point is higher than the average terrain of the neighborhood, and a negative value indicates that the grid point is lower than the average terrain of the neighborhood.

[0042] (8-3); In equation (8-3), TRIr is the terrain roughness index, which is used to characterize the root mean square magnitude of the elevation difference between the grid point and the surrounding area.

[0043] (8-4); In Equation (8-4), SVF is the sky visibility factor, βi,n is the horizon elevation angle at the nth azimuth angle, and Nα can be 16, 24 or 36. The smaller the SVF, the stronger the terrain shading.

[0044] (8-5); In equation (8-5), βup represents the direction of the incoming flow. The maximum windward shielding angle, SH is the dimensionless windward shielding index, used to distinguish between windward slope acceleration and leeward side obstruction.

[0045] (8-6); In equation (8-6), A can represent the coastline, major water body, ridge, valley, urban boundary or study area boundary, and dA is the shortest planar distance from the grid point to the corresponding regional feature.

[0046] (8-7); In equation (8-7), xreg,i and yreg,i are the normalized planar coordinates of the region, and Re is the average radius of the Earth. and Let L be the coordinates of the center of the study area, and L0 be the length of the regional feature.

[0047] (8-8); In equation (8-8), Li,k is the unique thermal code of land use category k, and z0,k ​​is the roughness length corresponding to the category. After category resampling, the unique thermal code can be replaced by the sub-pixel area ratio.

[0048] (8-9); In equation (8-9), Csta is the spatial representative weight of the site, dsta is the horizontal distance from the target grid point to the nearest valid site, and Ld and LH are the horizontal and vertical representative scales, respectively.

[0049] Ridge and valley masks can be determined jointly using regional quantile thresholds based on multi-scale TPI, curvature, and slope. For example, grid points with TPI above the 80th percentile and negative curvature are marked as ridge candidates, while grid points with TPI below the 20th percentile and positive curvature are marked as valley candidates. The thresholds are determined using training samples from the study area or global quantiles of the DEM to avoid the fixed thresholds becoming invalid due to changes in curvature units with resolution. Regional auxiliary information may also include vegetation cover, leaf area index, impermeable surface ratio, albedo, snow cover ratio, water body ratio, and climate zone coding, with the source data resolution, effective ratio, and missing data markers synchronously written into the confidence item cluster.

[0050] Dynamic background interpolation module; The near-surface variables U10, V10, T2, Q2, and PSFC of the dynamic model are interpolated to the target high-resolution grid to form a dynamic background interpolation field.

[0051] (9); In Equation (9), Fwrf represents the meteorological variables on the original grid of the dynamic model, Fint represents the dynamic background field after interpolation to the target high-resolution grid, and t represents the time index.

[0052] Six-item decomposition module for water vapor budget; like Figure 3 As shown, the system uses three-dimensional specific humidity q and horizontal wind vector V=(u,v) to calculate the overall water vapor flux convergence. To make this diagnostic applicable to downscaling tasks, the background flux is preferably defined as the area-weighted average within the coarse grid of the dynamic model at the same time, or a low flux field equivalent to the effective resolution of the dynamic model. The local disturbance is defined as the deviation of the target high-resolution field relative to this background field. Positive values ​​uniformly represent water vapor convergence or a positive contribution to local humidification.

[0053] (9-1); In Equation (9-1), Ac represents the coarse grid area aggregation of the dynamic mode or the specified low-pass operator, the overline represents the background quantity, and the apostrophe represents the local disturbance quantity.

[0054] (9-2); In equation (9-2), MFC is the overall water vapor flux convergence, ps and pt are the surface pressure and the integral top pressure, respectively, and g is the gravitational acceleration.

[0055] (9-3); Equation (9-3) is the dynamic advection term, which represents the advection contribution of local wind disturbance to the background water vapor gradient.

[0056] (9-4); Equation (9-4) is the dynamic convergence term, which represents the contribution of local wind field convergence or divergence changes to the background specific humidity conditions.

[0057] (9-5); Equation (9-5) is the thermal advection term, which represents the advection contribution of the background wind to the local specific humidity disturbance gradient.

[0058] (9-6); Equation (9-6) is the thermal convergence term, which represents the contribution of local specific humidity disturbance to the background convergence or divergence field.

[0059] (9-7); Equation (9-7) is a nonlinear advection term, representing the quadratic interaction contribution of wind disturbance and specific humidity disturbance gradient.

[0060] (9-8); Equation (9-8) is a nonlinear convergence term, representing the secondary interaction contribution of the convergence of specific humidity disturbance and wind disturbance.

[0061] (9-9); Equation (9-9) gives the water vapor flux convergence increment of the target high-resolution field relative to the dynamic background field. The sum of the six terms should be consistent with the directly calculated ΔMFC.

[0062] (9-10); Equation (9-10) is the discrete vertical integral form of any decomposition term, Ai,k is the integrand of the term in the k-th layer, Δpi,k is the layer thickness, and Ct is in mm·h-1 when it is 3600 and in mm·d-1 when it is 86400.

[0063] (9-11); In equation (9-11), W is the total precipitable water volume, P is the precipitation, E is the surface evaporation, and Rmb is the column water vapor balance residual.

[0064] During computation, it is preferable to use layer thickness for quality weighting on the model's native terrain-following layer and apply an invalid mask to the subsurface layer; alternatively, interpolation to the common pressure layer can be performed before integration. Horizontal gradients, convergence, and divergence must use the same target grid, the same difference order, and consistent boundary treatment. Lower boundary terms caused by spatial variations in surface pressure, transient eddy terms due to insufficient output time resolution, and numerical difference errors are uniformly included in Rmb and quality control markers, and are not mixed into the aforementioned six physical decomposition terms. When P or E is unavailable, Rmb is not forcibly closed; only the six decomposition terms, W, and computable local variation terms are used as residual learning features.

[0065] Temperature and terrain physical correction module; The near-surface temperature vertical lapse rate is extracted from the three-dimensional profile of the dynamic model. Specifically, the mass layer height can be calculated from the geopotential height, the actual temperature can be recovered from the perturbation temperature and air pressure, and the imaginary temperature can be calculated by combining the water vapor term. Then, the vertical lapse rate is estimated based on the imaginary temperature difference at different heights in the near-surface layer.

[0066] (10); In equation (10), Indicates the first Each quality layer height, and These represent the perturbation potential height and the ground state potential height, respectively. It represents the acceleration due to gravity.

[0067] (11); In equation (11), Indicates the first The actual temperature of the layer, This represents the perturbation potential temperature variable in the dynamic mode. and These represent the disturbed pressure and the ground-state pressure, respectively. Desirable , The gas constant of air. This is the specific heat capacity at constant pressure.

[0068] (12); In equation (12), Indicates a false temperature. Indicates the first The water vapor mixing ratio or specific humidity approximation term. The influence of water vapor on air density and hydrostatic relationship can be incorporated into the calculation through virtual temperature.

[0069] (13); In equation (13), This represents the vertical lapse rate of near-surface temperature. and The lower and upper quality layers near the surface can be sampled separately. This indicates that the rate of decline will be limited to a preset reasonable range. Under typical tropospheric conditions, It can be expressed as a negative value.

[0070] (14); In equation (14), T2,corr represents the 2m air temperature after terrain correction, T2,int represents the interpolated dynamic background air temperature, ΔH represents the terrain elevation difference obtained from equation (5), and Tmin and Tmax are the temperature physical range thresholds. When the actual terrain is higher than the equivalent coarse terrain of the dynamic model, the air temperature decreases according to the vertical lapse rate; otherwise, it increases.

[0071] Temperature correction is preferably applied to numerical weather and regional climate data with near-surface air temperatures ranging from 180 to 330 K. The near-surface lapse rate Γ estimated from the model profile can be limited to -0.012 to 0.003 K·m⁻¹. Under typical neutral or weakly unstable conditions, -0.0065 K·m⁻¹ can be used as a backoff value, while stable inversion conditions allow Γ to be positive. If the original or corrected temperature exceeds 180–330 K, the system does not directly input the outlier into the saturated vapor pressure formula. Instead, it first replaces the outlier with the dynamic background from the same time period, adjacent effective grid points, or the previous effective time period, and sets a temperature anomaly marker; then it proceeds to the phase-sensing saturation calculation.

[0072] For grid points with a large |ΔH|, temperature, virtual temperature, and air pressure are updated piecewise according to equation (5-1). At each step, the temperature is first updated using the average lapse rate of the current layer, then the updated average virtual temperature is integrated with the hydrostatic relationship, and finally, the specific humidity is recalculated while maintaining the background relative humidity or dew point difference. This sequence avoids using a single fixed virtual temperature across excessively large altitude differences in high-altitude or deep valley environments and ensures that the air pressure and water vapor variables are corrected in the same direction with altitude.

[0073] Baroque Physical Correction Module The system corrects surface air pressure based on static relationships. First, it estimates the mean virtual temperature based on the corrected temperature and specific humidity, and then calculates the pressure index correction term based on the elevation difference.

[0074] (15); In Equation (15), the overlined Tv represents the average virtual temperature used for hydrostatic pressure correction, and Q2,int represents the interpolated 2m specific humidity. This equation approximates the average temperature of the altitude range with the average of the background temperature and the terrain-corrected temperature, and takes into account the influence of water vapor on air density with the background specific humidity.

[0075] (16); In equation (16), Pcorr represents the surface air pressure after terrain correction, Pint represents the interpolated dynamic background air pressure, Rd is the air gas constant, and Pmin and Pmax are the air pressure physical range thresholds. This equation reflects the static relationship that air pressure decreases when elevation increases and increases when elevation decreases.

[0076] Humidity thermodynamic consistency correction module; The system calculates the actual water vapor pressure based on the temperature, specific humidity, and air pressure before interpolation, and calculates the relative humidity based on the saturated water vapor pressure formula; then, it recalculates the actual water vapor pressure and specific humidity under the corrected temperature and air pressure conditions, so that the humidity variable remains thermodynamically consistent with the temperature and air pressure.

[0077] (17); In equation (17), Indicates temperature The saturated vapor pressure is expressed in Pa. This formula is an empirical expression for saturated vapor pressure of the Bolton type and can be used to calculate near-surface relative humidity and saturated specific humidity.

[0078] (18); In equation (18), Indicates the specific moisture content and air pressure The calculated actual water vapor pressure.

[0079] (19); In equation (19), RHint represents the relative humidity corresponding to the interpolated background field. By first calculating the relative humidity and then restoring the specific humidity under the corrected temperature and corrected pressure, the physical meaning of the water vapor state can be preserved as much as possible.

[0080] (20); In equation (20), enew represents the actual water vapor pressure recalculated under the corrected temperature conditions, and Q2,corr represents the 2m specific humidity after topographic correction. By imposing an upper limit on enew, the water vapor pressure can be prevented from exceeding the allowable range of atmospheric pressure.

[0081] (twenty one); In equation (21), This indicates the saturated specific humidity. When the corrected specific humidity is greater than the saturated specific humidity, the system will limit the specific humidity to the saturated specific humidity, thereby eliminating the supersaturated state.

[0082] Near-surface wind field topography diagnosis and correction module; The system calculates wind speed and direction from the interpolated U10 and V10, and diagnoses and corrects the near-surface wind field by combining slope, aspect, curvature, and topographic elevation differences. First, it determines the windward or leeward effect of the slope based on the angle between the wind direction and the slope aspect, and applies a topographic acceleration or deceleration factor to the wind speed. Second, it identifies valley topography based on slope and curvature, and slightly deflects the wind direction. Third, it describes the topographic blocking effect based on the Froude number. Finally, it applies an upper limit constraint to the wind speed.

[0083] (twenty two); In equation (22), W represents the interpolated background wind speed, θ represents the airflow direction angle taken from the due north vector clockwise, and Uint and Vint represent the interpolated east-west and north-south 10m wind components, respectively; if the meteorological direction angle is used, then add 180° to the angle and take the modulus of 360°.

[0084] (twenty three); In equation (23), Indicates the slope acceleration or deceleration factor. This is the slope effect coefficient. For slope, Slope direction, This indicates the angle between the wind direction and the slope aspect. Different acceleration or deceleration effects can occur under uphill or downhill conditions.

[0085] (twenty four); In equation (24), Indicates a valley passage mask. Indicates the curvature threshold. This indicates the amount of wind deflection caused by the valley passage. The channel effect coefficient, This is the slope attenuation scale.

[0086] (25); In equation (25), Represents the simplified Froude number. For stability parameters, and To prevent small amounts from being divided by zero, This represents the terrain blocking factor. When... At that time, airflow is more easily blocked by terrain, and wind speed decreases accordingly.

[0087] (26); In equation (26), This indicates the corrected wind speed. and These represent the corrected east-west and north-south wind components, respectively. This indicates the maximum wind speed limit.

[0088] Physical consistency post-processing module; After completing the physical corrections for temperature, air pressure, humidity, and wind field, the system performs the first physical consistency post-processing. This post-processing includes non-finite value replacement, temperature range constraints, air pressure range constraints, specific humidity range constraints, wind speed upper limit constraints, saturated specific humidity constraints, relative humidity recalculation, air density rationality checks, detection of abnormal air pressure related to altitude measurements, and optional weak ground-to-wind constraints.

[0089] (27); Equation (27) provides the objective variable range constraint used in the post-processing of physical consistency. This constraint is used to identify and correct numerical outliers, rather than to artificially adjust local climate characteristics.

[0090] like Figure 4 As shown, the physical closure constraint module is further executed in three levels: hard constraints, soft constraints, and iterative projection. Hard constraints are used to eliminate infinite values, non-numerical values, out-of-bounds states, and oversaturated states; soft constraints are used to satisfy static relationships, column water vapor budget, and consistency with the dynamic background as much as possible; iterative projection is used to coordinate variable updates among multiple constraints.

[0091] (27-1); In equation (27-1), Y is the state vector of a single grid point, Ω is the feasible region of hard constraints, and S is the wind speed. The global range is only used for anomaly identification, and the local upper and lower limits are optimized and then tightened by combining the dynamic background and the training residual quantiles.

[0092] (27-2); In equation (27-2), Tc is the temperature in Celsius, and es is in Pa; the saturated vapor pressure of water surface is used above zero degrees Celsius, and the saturated vapor pressure of ice surface is used below zero degrees Celsius. The two equations can be linearly smoothed within the range of -1 to 1℃ to avoid discontinuities caused by phase switching.

[0093] (27-3); In equation (27-3), ε = Rd / Rv is approximately 0.622. When closing the loop, es and qs are first calculated from T and p, and then q is projected onto [0, (1-εs)qs], where the safety factor εs can be taken as 10⁻⁶ to 10⁻³.

[0094] (27-4); In equation (27-4), rh is the static residual, pbg is the dynamic background air pressure, and the overlined Tv represents the average virtual temperature in the correction height range; rh approaches 0 when the ideal static closure is achieved.

[0095] (27-5); In equation (27-5), rw is the column water vapor budget residual. Non-zero closed weights are assigned only when the time resolution of P, E, W, and MFC meets the calculation conditions.

[0096] (27-6); In equation (27-6), Y0 is the physical downscaling field or statistical residual correction field before entering closure, Wd is the variable retention weight; Dρ is the normalized distance when the density exceeds the reasonable range, Dv is the normalized distance when the wind speed exceeds the limit, and λ is the weight of each soft constraint.

[0097] (27-7); In equation (27-7), For operators projected onto the hard-constrained feasible region, The step size is used. Engineering implementation can employ projection gradient, alternating directions, or alternating projection in the order of temperature—air pressure—humidity—wind field.

[0098] (27-8); Equation (27-8) gives the closure stopping condition. When any variable fails to converge within the maximum number of iterations, the system retains the last state that satisfies the hard constraints and outputs a closure non-convergence flag.

[0099] (27-9); In equation (27-9), rm is the original prediction residual of variable m, am is the residual limit, αi is the confidence scaling factor, and Iqc is the quality control weight from 0 to 1. This equation automatically weakens statistical corrections in cases of sparse stations, terrain extrapolation, or input of missing measurement areas.

[0100] The recommended execution order for closure is as follows: First, unify the units of temperature, air pressure, specific humidity, and relative humidity and replace non-finite values; second, perform wide-range projection according to equation (27-1); then, recalculate saturated vapor pressure, actual vapor pressure, saturated specific humidity, and relative humidity according to equations (27-2) and (27-3); then check the static relationship and air density according to equations (27-4) and (28); after that, check the upper limit of wind speed, wind direction continuity, and optional weak constraints on ground-switching winds; finally, when the data is complete, apply the column water vapor budget soft constraint according to equation (27-5), and iterate according to equation (27-7) until equation (27-8) is satisfied. The first closure uses the purely physical downscaling field as Y0, and the second closure uses the statistical residual correction field as Y0, but the adjustment range of the statistical residual for the second closure must not exceed the variable residual limit.

[0101] (28); In equation (28), This indicates air density. If the air density exceeds a preset reasonable range, the system can calculate the air pressure based on the trimmed air density and virtual temperature to correct obviously non-physical air pressure conditions.

[0102] (29); In equation (29), This represents the reference air pressure calculated based on altimetry. This represents the interpolated background air pressure. This indicates the deviation between the current air pressure and the reference air pressure.

[0103] (30); In equation (30), Indicates a weak adjustment coefficient. This represents the threshold for detecting air pressure anomalies. Only when the air pressure deviation exceeds the threshold will the system make a weak adjustment to the altimeter reference air pressure with a small weight, in order to avoid overcorrecting the already reasonable air pressure structure under complex terrain.

[0104] In an alternative implementation, the system can calculate a weak geostrophic wind trend based on a smoothed pressure field and Coriolis parameters, and mix it into the near-surface wind field with a small weight.

[0105] (31); In equation (31), and These respectively represent the trend of weak ground turning into wind. Represents the Coriolis parameter. Represents the Earth's angular velocity of rotation. Indicates latitude.

[0106] (32); In equation (32), This represents the weak mixing coefficient of geostrophic winds. This coefficient is limited to a small range so that the geostrophic constraint only serves as a correction to large-scale trends, rather than dominating the near-surface wind field in complex terrain.

[0107] The applicable temperature range and anomaly handling rules are as follows: The Bolton surface saturated vapor pressure of Equation (17) has high accuracy in the range of approximately -35 to 35°C; In this invention, the ice surface is expressed using Equation (27-2) at sub-zero temperatures, and the numerical evaluation temperature for saturation calculation is limited to -80 to 50°C. When the temperature to be processed exceeds this range, it is first repaired to a wide physical range of 180 to 330 K according to the dynamic background or neighborhood effective value, and then the saturation calculation is performed. The original out-of-bounds value and the repair method are both written into the quality control mark.

[0108] Site residual learning module; When a CSV file of station observations exists, the system reads the station observation data and standardizes common units, such as converting Celsius temperature to Kelvin (K), hectopascal pressure to Pascal (Pa), and relative humidity in the range of 0 to 1 to a percentage. Then, it matches the station observations with the output timestamps of the dynamic model and matches the station's latitude and longitude to the nearest high-resolution grid point.

[0109] For each site sample, the system extracts features such as site location, site elevation, high-resolution DEM elevation, coarse topographic background, elevation difference, slope, sine of slope aspect, cosine of slope aspect, curvature, hour, day of year, T2, Q2, PSFC, U10, V10, RH2 after physical downscaling, and wind speed after physical downscaling.

[0110] (33); In equation (33), Indicates site At any moment For variables residual labels, Represents the station observation value. This indicates the physical downscaling result at the grid point corresponding to the site.

[0111] (34); In equation (34), This represents the residual model input feature vector representing the target grid points. Indicates hourly characteristics, Indicates the characteristics of accumulated days over a year. Represents a set of physically downscaled meteorological variables. This represents the physical downscaled wind speed.

[0112] Item cluster feature definition and construction module (34-1); Equation (34-1) defines the topographic morphology term family, where β is the slope, α is the aspect, C is the curvature, and Iridge and Ivalley are the markers for ridge and valley, respectively; variables with subscript r are calculated at multiple physical scales.

[0113] (34-2); Equation (34-2) defines the region and the underlying surface cluster, where fveg, fimp, and fsnow are the proportions of vegetation, impermeable surface, and snow cover, respectively.

[0114] (34-3); Equation (34-3) defines the dynamic background term family, where PBLH is the boundary layer height and wlow is either the lower-level vertical velocity or the pressure velocity.

[0115] (34-4a); Equation (34-4a) defines the thermo-vapor term cluster, where Td is the dew point temperature and qs-q is the saturated specific humidity difference; when it is necessary to characterize temperature inversion, stability or boundary layer humidity stratification, the low-level temperature and humidity profile difference can be added.

[0116] (34-4b); Equation (34-4b) gives the imaginary temperature, air density, imaginary potential temperature and the square of the Brunt–Väisälä frequency in the thermo-vapor term.

[0117] (34-5); Equation (34-5) defines a cluster of water vapor budget terms, where the first six terms are the six core decomposition terms of this invention, and the remaining terms are used for closure diagnosis and quality control.

[0118] (34-6); Equation (34-6) defines a time periodic term cluster, where h is the hour and d is the yearly cumulative day; leap years can use 365.2425 or be normalized according to the number of days in the year.

[0119] (34-7); Equation (34-7) defines the multi-scale neighborhood mean, standard deviation and horizontal gradient modulus of continuous features, which are used to express local background, heterogeneity and transition zone.

[0120] (34-8); Equation (34-8) defines robust standardization based on the median and interquartile range of the training samples. The tree model can retain the original dimensions, but it must still preserve the missing test surrogate values ​​and feature order from the training phase.

[0121] (34-8a); Equation (34-8a) defines the data confidence term family, where rvalid is the proportion of valid inputs. Encoding the resampling method: Ioos is the terrain or climate extrapolation marker, Imiss is the key variable missing marker, Iclose is the first closure state, and niter is the number of closure iterations.

[0122] (34-9); Equation (34-9) gives the final term cluster feature vector; each variable residual model can select a subset from each term cluster, but must save and reuse the term cluster order, variable units and missing test handling rules from the training phase.

[0123] To avoid target leakage, site observations are used only for constructing residual labels and are not used as input features for the same time-varying variable. The site information that can be used as input is limited to site location, elevation, quality control level, distance to the target grid point, and historical availability. Training, validation, and testing should preferably be spatially and temporally grouped by site or consecutive time blocks. Randomly splitting adjacent time intervals of the same site into the training and testing sets is prohibited as it could overestimate generalization ability.

[0124] Item clusters are feature sets grouped according to physical origin and spatial scale of input features. The residual model for each variable can use all item clusters, or select subsets based on cross-validation and feature importance; the subsets of item clusters used for different variables and their order must be consistent with those used in the training phase.

[0125] For each variable to be corrected, the system trains a residual model separately. The residual model can be an extreme random tree, random forest, or histogram gradient boosting regression model. When the sample size is insufficient or machine learning libraries are unavailable, the system skips statistical training and maintains a purely physical downscaling mode.

[0126] ; In equation (35), Representing variables The corresponding residual learning model, This represents the residual predicted by the model at the target grid point.

[0127] Statistical residual correction module; In the final inference stage, based on the trained residual model file, the system constructs a feature matrix consistent with that in the training stage for each grid point of the target grid, uses the residual model corresponding to each variable to predict the residual, and superimposes the predicted residual onto the physical downscaling result.

[0128] (36); In equation (36), This represents the prediction residual after amplitude limiting. Representing variables The corresponding maximum allowable residual margin. Residual margin limits can reduce the risk of machine learning models making anomalous corrections in areas with sparse sites or terrain extrapolation.

[0129] (37); In equation (37), This represents the hybrid-statistical downscaling result after statistical residual correction. The result still primarily consists of the physical downscaling field, with only the systematic correction obtained from observational residual learning superimposed.

[0130] (38); In equation (38), This represents a multivariate physical consistency post-processing operator.

[0131] The effect of statistical correction can be evaluated during the model training phase using residual RMSE, MAE, and skill scores relative to the zero residual baseline.

[0132] ; In equation (39), and Representing variables respectively The root mean square error and mean absolute error of the residual prediction. This indicates the number of test samples.

[0133] ; In equation (40), This represents the zero residual baseline error without statistical residual correction. This indicates the degree of improvement of the residual model compared to pure physical downscaling.

[0134] Output module; The system writes the final meteorological field into a NetCDF file, with output variables including XLAT, XLONG, HGT, U10, V10, T2, Q2, PSFC, and RH2. Output dimensions include time, south_north, and west_east.

[0135] ; In equation (41), Indicates the time in the output NetCDF file Target grid points The corresponding dataset. This output can be used for mountain meteorological analysis, hydrological model-driven applications, disaster risk assessment, agricultural meteorological services, wind energy resource analysis, and environmental simulation.

[0136] Key parameter values ​​and adjustment rules; Regarding grid and terrain parameters, the target resolution can be 100–1000 m, with a typical value of 300 m; the terrain neighborhood scale should cover approximately 1–15 target grid points or 0.3–5 km; the Earth's average radius Re can be 6371.0–6378.1 km; the elevation step size ΔHstep can be 200–500 m, with a typical value of 300 m; the horizontal representative scale Ld of the station can be 10–50 km, and the vertical representative scale LH can be 300–1000 m, both of which can be determined by leaving one station for verification. The roughness length z0 should be given by the land use lookup table or observation data, and is usually limited to 10⁻⁴–2 m.

[0137] Regarding temperature, pressure, and humidity parameters, the near-surface lapse rate Desirable When a reliable profile is lacking, it can be reverted to The wide temperature anomaly range can be set from 180 to 330 K, the wide pressure anomaly range can be set from 15 to 110 kPa, and the wide specific humidity range can be set... The relative humidity is set to 0–100%. This range is only for identifying numerical anomalies; for actual closure, the local dynamic range is preferably constructed using the dynamic background plus the 0.5%–99.5% quantile of the training residuals. The saturation safety factor εs can be taken as... Air density can be taken over a wide range .

[0138] Regarding wind field diagnostic parameters, the slope effect coefficient can be taken as 0–0.6, with a typical value of 0.2–0.3; the channel deflection coefficient can be taken as 0–0.5, with a typical value of 0.15; the stability parameter N can be taken as… Typical value A Froude number less than 1 increases terrain blocking, while a number equal to or greater than 1 gradually weakens it; the maximum wind speed Smax can be set to... Typical business value The weak geo-wind mixing coefficient can be set to 0–0.15, with a typical value of 0.05, and should only be activated when the absolute latitude is greater than 5° and the pressure gradient smoothing is effective.

[0139] Regarding statistical residuals and closure parameters, the optimal residual limit *am* is the minimum value among the 99th quantile of the absolute value of the training residuals, three times the interquartile range, and the physical upper limit; T² can be 2–8K, and Q² can be... PSFC is 1000-5000Pa, and the single wind component is... The initial search interval. Closing weights. and You can take 1 to 20. When the water vapor balance data is reliable, the value can be taken as 0.1–10; when P or E is missing, it should be set to 0. The projection step size η can be taken as 0.05–0.5, and the maximum number of iterations can be taken as 3–10, with a typical value of 5. The convergence tolerance can be set as follows: The water vapor balance tolerance is acceptable. Or 5% of the absolute value of the local PE, taking the larger of the two.

[0140] All parameters should be calibrated using an independent validation set for the training area, leave-one-site validation, or weather-specific cross-validation. Parameter adjustments should follow the order of ensuring physical closure first, then minimizing observation errors: when the closure residual increases, prioritize reducing the statistical residual limit or increasing the corresponding physical constraint weights; when the mountain extrapolation error increases, prioritize reducing Hoos or increasing the topographic extrapolation attenuation; when flat areas are overcorrected, reduce the slope, passage, and blocking coefficients.

[0141] Specific application examples and effect comparisons; In a specific application embodiment, such as Figure 5 As shown, the system reads the output of the 3km-level dynamic model, the high-resolution digital elevation model, regional auxiliary data, and optional station observation data to construct a 300m target high-resolution grid, generating real terrain, equivalent coarse terrain, and elevation difference, and constructing slope, aspect, curvature, local undulation, and regional auxiliary term clusters. When three-dimensional wind, water vapor, air pressure, and layer thickness data are available, the system outputs dynamic advection, dynamic convergence, thermal advection, thermal convergence, nonlinear advection, and nonlinear convergence terms according to the aforementioned six-term decomposition method of water vapor budget; when no station observations or residual models are connected, the system outputs high-resolution meteorological grid products according to the pure physical branch.

[0142] Figure 6A comparison of the results from the original procedure before modification and the output of the present invention's process is presented at the same representative time point. The figure shows the spatial distribution of 2m air temperature, 10m wind speed, and surface air pressure before modification, the spatial distribution after modification using the patented process, and the differences between the two. This illustrates the impact of terrain physical correction, complex terrain wind field diagnosis, and physical closure constraints on the final gridded field spatial structure. This comparison is intended to demonstrate the changes the present invention's process brings to the case results and is not intended as an evaluation of observation accuracy or a limitation on the scope of protection.

[0143] like Figure 7 As shown, high-resolution terrain elevation and slope can characterize complex landform structures such as canyons, ridges, slopes, and plateau edges. The 2m air temperature and surface pressure field exhibit significant spatial responses with terrain height. This example illustrates that the difference between the actual terrain and the equivalent coarse terrain, and its derived topographic quantities, can provide local terrain forcing for temperature, pressure, and humidity-linked physical correction.

[0144] like Figure 8 As shown, the system can perform a posteriori physical closure diagnosis on the output field. The diagnostic parameters include the ratio of specific humidity to saturated specific humidity, air density, 10m wind speed, and quality control markers. The quality control markers are used to identify areas that violate wide physical boundaries, saturation constraints, density constraints, or wind speed constraints, thereby supporting the inspection of the effectiveness of the first and second physical closures.

[0145] like Figure 9 As shown, the system can further statistically analyze the daily variations of the regional average 2m air temperature, 10m wind speed, 2m relative humidity, and ground air pressure over a 24-hour period, to check the continuity and physical plausibility of the output results in the time dimension. For example... Figure 10 As shown, the system can also statistically analyze the response relationships of 2m air temperature, ground air pressure, 2m relative humidity and 10m wind speed with elevation changes according to elevation zones, which can be used to verify the overall performance of multi-variable linkage correction of temperature, pressure, humidity and wind in complex terrain areas.

[0146] In the above embodiments, the output grid product includes latitude and longitude, terrain elevation, 10m east-west wind component, 10m north-south wind component, 2m air temperature, 2m specific humidity, surface air pressure, and 2m relative humidity. It can be further extended to store elevation difference, terrain term clusters, water vapor budget decomposition terms, number of closed iterations, and quality control markers. The date, region, spatial resolution, variable set, residual model type, and output format used in the embodiments can all be adjusted according to the input data and application requirements, and do not constitute a limitation on the scope of protection of this invention.

[0147] This invention is not limited to the foregoing embodiments. Without departing from the inventive concept, the following alternative solutions exist: First, the output of the dynamical model is not limited to the WRF model; it can also be used for regionalized output of other regional numerical weather models, reanalysis data, ensemble forecast products, or climate models. Secondly, the target resolution is not limited to 300m, and can be set to 1km, 500m, 100m or other spatial resolutions according to the DEM resolution, application requirements and computing resources; Third, DEM resampling methods are not limited to linear interpolation; nearest neighbor interpolation, cubic interpolation, spline interpolation, kriging interpolation, or multi-scale fusion resampling methods can also be used. Fourth, the temperature lapse rate is not limited to being extracted from WRF profiles; it can also be obtained from regional climate statistics lapse rate, seasonal lapse rate, weather-type lapse rate, or station-derived lapse rate. Fifth, the residual model is not limited to extreme random trees, random forests or histogram gradient boosting models, but can also use XGBoost, LightGBM, CatBoost, neural networks, Gaussian process regression or other machine learning models with physical constraints. Sixth, the post-processing of physical consistency is not limited to the constraints listed in this invention, but may also include radiation balance constraints, boundary layer stability constraints, land surface type constraints, snow cover constraints, or land use modulation constraints. Seventh, the output format is not limited to NetCDF; it can also output GeoTIFF, Zarr, HDF5, GRIB, or tiled grid data for web services.

[0148] It should be noted that, in order to simplify the description of the present invention and thus help to understand one or more embodiments of the invention, multiple features may sometimes be grouped into one embodiment, drawing or description thereof in the foregoing description of the embodiments of the present invention.

[0149] The embodiments of this application have been described above with reference to the accompanying drawings. Unless otherwise specified, the embodiments and features in the embodiments of this application can be combined with each other. This application is not limited to the specific embodiments described above. The specific embodiments described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms under the guidance of this application without departing from the spirit and scope of the claims, and all of these forms are within the protection scope of this application.

Claims

1. A high-resolution meteorological downscaling method integrating dynamic models, topographic physical correction, and station residual learning, comprising: The system comprises the following modules: data preprocessing, high-resolution target mesh construction, multi-layer terrain background and multi-scale terrain derived factor construction, dynamic background interpolation and water vapor budget decomposition, terrain-driven physical correction of temperature, pressure, humidity and wind fields, first multivariate physical consistency closure module, station residual model training module, residual inference correction module, second physical consistency closure module, and meteorological grid product output module. Each module is executed sequentially along a directed acyclic call chain, and the implementation flow is as follows: Step 1: Multi-source meteorological and geographic data reading, coordinate unification, spatiotemporal matching and basic quality control. Input data includes regional dynamic model output, high-resolution DEM, underlying static auxiliary data, and optional ground station observation data. Step 2: Generate a standard high-resolution target mesh based on the dynamic mode spatial range and the preset target resolution; Step 3: Based on the DEM, construct a high-resolution real topographic field and a dynamic model equivalent coarse topographic field, calculate the elevation difference between the two, and generate single-scale and multi-scale topographic morphological features and regional underlying surface auxiliary features. Step 4: Interpolate the near-surface meteorological variables of the dynamic model to the target grid to obtain the dynamic background field; use the three-dimensional temperature, pressure, wind, and water vapor fields to complete the six-component physical decomposition of the entire layer's water vapor budget and generate dynamic meteorological feature clusters; Step 5: Based on the terrain elevation difference, perform stepwise thermodynamic and physical correction of near-surface air temperature, surface air pressure, and specific humidity. Combine slope, aspect, curvature, and terrain blocking effect to complete the 10m wind field diagnosis and correction, and obtain a purely physical downscaling meteorological field. Step 6: Perform the first global multivariate physical consistency closure constraint on the purely physical downscaling field to eliminate non-physical states such as numerical anomalies, violations of hydrostatic equilibrium, water vapor supersaturation, and wind speed exceeding limits, and output a physically consistent meteorological state field that can be used for machine learning training; Step 7: If valid station observation data exists, match the station observations with the corresponding grid point physical downscaling results to construct meteorological residual labels, and train the variate residual regression model by integrating multiple feature clusters such as topography, region, dynamics, thermodynamics, water vapor, time, and confidence; if there is no station data, skip this step and directly output the pure physical downscaling product. Step 8: In the inference stage, the trained residual model is called to predict the grid residuals of each meteorological variable, limit the amplitude, and then superimpose them onto the physical downscaling field to obtain the hybrid-statistical correction meteorological field. Step 9: Apply a second multivariate physical consistency closure constraint to the corrected meteorological field to slightly correct the physical imbalance problem introduced by the statistical residuals; Step 10: Output the high-resolution temperature, pressure, humidity, and style point meteorological fields, as well as topography and quality control markers, which have undergone two physical closure constraints, as standardized meteorological files.

2. The high-resolution meteorological downscaling method according to claim 1, which integrates dynamic models, topographic physical correction, and station residual learning, is characterized in that: Step 1, the multi-source data preprocessing, specifically includes: reading the output files of WRF regional dynamic models, reanalysis, or climate models; extracting grid latitude and longitude, three-dimensional geopotential height, disturbed temperature and pressure, water vapor, near-surface 2m temperature and humidity, 10m wind, surface air pressure, precipitation, and surface flux variables; reading the high-resolution DEM digital elevation model; reading static auxiliary data such as land use, vegetation cover, snow cover, coastline, and administrative boundaries; optionally reading station table data containing time, latitude and longitude, elevation, temperature, pressure, humidity, and wind observations; unifying the coordinate reference system of all data, using conservative averaging / bilinear resampling for continuous geographic variables, and mode resampling for categorical underlying surface variables; completing unit conversion, time alignment, and missing measurement value identification for meteorological variables, and simultaneously generating data credibility marker clusters; when performing water vapor budget decomposition, verifying the completeness of three-dimensional wind, vertical velocity, water vapor mixing ratio, and pressure layer thickness, calculating only solvable decomposition items for missing variables and updating credibility markers, and prohibiting zero-value filling of missing physical quantities.

3. The high-resolution meteorological downscaling method according to claim 2, which integrates dynamic models, topographic physical correction, and station residual learning, is characterized in that: The high-resolution target mesh construction method described in step 2 is as follows: the spatial boundary is determined by the extreme values ​​of latitude and longitude of the original mesh in the dynamic mode, and the latitude cosine correction of the regional center is introduced in combination with the geometric relationship of the earth. The latitudinal and longitudinal mesh intervals are calculated respectively. Based on the mesh interval and latitude and longitude boundary, the total number of north-south and east-west grid points of the target mesh is solved to generate a one-dimensional latitude and longitude array and a two-dimensional target latitude and longitude mesh. A block-based lazy loading strategy is adopted, dividing the target mesh into 256×256 or 512×512 blocks, setting an overlap buffer between blocks, and after the calculation is completed, the overlapping areas are trimmed and the entire field is stitched together, reducing memory usage.

4. The high-resolution meteorological downscaling method according to claim 1, which integrates dynamic models, topographic physical correction, and station residual learning, is characterized in that: Step 3, the multi-layer terrain background construction process, includes: cropping the DEM, transforming coordinates, and interpolating to the target mesh to obtain a high-resolution terrain field. Using a coarse grid of dynamic model as the unit, the average elevation of all high-resolution grid points within the unit is aggregated to obtain the equivalent coarse terrain. Then, the coarse terrain is interpolated to the target grid; the elevation difference is calculated point by point. ; Set a fixed elevation integration step size for... Temperature, virtual temperature, and air pressure are updated iteratively in segments for grid points exceeding the threshold. When the absolute value of the elevation difference exceeds the 99th percentile of the training samples or 1500m, the terrain extrapolation area is marked, and the correction weight of the stations in that area is reduced in the subsequent residual learning stage.

5. The high-resolution meteorological downscaling method according to claim 1, which integrates dynamic models, topographic physical correction, and station residual learning, is characterized in that: Step 3 synchronously generates terrain and regional auxiliary features including: single-scale basic terrain factors: Slope, aspect, and topographic curvature are calculated using the central difference / Sobel operator; multi-scale topographic morphology factors: multiple physical distance neighborhoods are set up to calculate topographic relief, topographic location index (TPI), topographic roughness (TRI), sky visibility factor (SVF), and windward shading index along the incoming flow direction at each scale; spatial distance factors: Calculate the shortest planar distance to the coastline, water body, ridge, valley, and study area boundary point by grid point; Underlying surface and station representativeness factors: land use unique thermal coding, vegetation / impermeable surface / snow cover ratio, and roughness length; construct station spatial representativeness weights to characterize the horizontal and vertical representativeness of observation stations around grid points.

6. The high-resolution meteorological downscaling method according to claim 5, which integrates dynamic models, topographic physical correction, and station residual learning, is characterized in that: The six-element decomposition of water vapor budget described in step 4 is implemented as follows: the coarse grid field of the dynamic model is aggregated by area low-pass to obtain a large-scale background field, and the difference between the original field and the background field is defined as the local disturbance quantity. The entire layer water vapor flux convergence (MFC) is solved by integrating three-dimensional specific humidity, horizontal wind field, and pressure layer thickness, and decomposed into dynamic advection term, dynamic convergence term, thermal advection term, thermal convergence term, nonlinear advection term, and nonlinear convergence term. The sum of the six decomposition terms equals the water vapor convergence increment of the high-resolution field relative to the coarse background field; Simultaneously calculate the total precipitable water and water vapor collection pillar residuals, and use all decomposed quantities as dynamic input features for residual learning.

7. The high-resolution meteorological downscaling method according to claim 5, which integrates dynamic models, topographic physical correction, and station residual learning, is characterized in that: The physical correction process for temperature, air pressure, and humidity in step 5 is as follows: The height, actual air temperature, and virtual temperature of each layer are obtained from the multi-layer profile of the dynamic model. The vertical lapse rate Γ of the virtual temperature in the near-surface layer is extracted and a reasonable range is defined. 2m air temperature is corrected based on piecewise iterative correction of elevation difference; The average virtual temperature is calculated by taking the average temperature before and after the correction, and the surface pressure topography is corrected by combining the hydrostatic index equation. The original relative humidity is calculated from the interpolated background temperature, pressure and specific humidity. Under the corrected temperature and pressure conditions, the saturated vapor pressure and actual vapor pressure are recalculated using the Bolton saturated vapor pressure formula. The specific humidity after topographic correction is solved in reverse and a saturation constraint is applied to limit the specific humidity to not exceed the saturated specific humidity and eliminate the supersaturated state.

8. The high-resolution meteorological downscaling method according to claim 1, which integrates dynamic models, topographic physical correction, and station residual learning, is characterized in that: Step 5, the 10m wind field correction steps, include: Calculate wind speed and meteorological wind direction from interpolated wind components; calculate slope acceleration / deceleration factor based on the angle between wind direction and slope aspect. Valley passage areas are identified by slope and curvature, and valley wind deflection is calculated. A simplified Froude number is introduced to characterize the topographic airflow blocking effect, and the topographic blocking factor is solved. After the wind speed is corrected by coupling the slope factor and the blocking factor, the wind direction deflection is superimposed to re-decompose the east-west and north-south wind components, and a global upper limit constraint on the wind speed is set. Optionally, a weakly weighted ground-to-wind-trend hybrid correction can be introduced to modify the large-scale wind field background.

9. The high-resolution meteorological downscaling method according to claim 1, which integrates dynamic models, topographic physical correction, and station residual learning, is characterized in that: Step 6, the first physical consistency closure, includes three layers of constraints: Hard constraints: Eliminate infinite and missing values, limit the global physical range of temperature, pressure, humidity, and wind speed, and forcibly constrain relative humidity to 0-100% to eliminate oversaturation; Differentiate between formulas for saturated vapor pressure on ice / water surfaces, and empirical formulas for switching ice surfaces at sub-zero temperatures; Soft constraints: Construct cost functions for static equilibrium residuals and column water vapor balance residuals, and introduce multivariate weights to balance various physical constraints; Iterative projection optimization: The temperature, pressure, humidity and wind state vectors are updated iteratively using the projection gradient method until the change in variables is less than the convergence threshold or the maximum number of iterations is reached; Unconverged grid points retain their valid physical states and are marked as closure anomalies.

10. The high-resolution meteorological downscaling method according to claim 1, which integrates dynamic models, topographic physical correction, and station residual learning, is characterized in that: The residual model training method described in step 7 includes: The site observations are matched to the corresponding target grid points, and the observations minus the physical downscaling results of the same grid point are used as the residual labels. Seven feature clusters are constructed: topographic morphology cluster, regional underlying surface cluster, dynamic background cluster, thermal and water vapor cluster, water vapor budget decomposition cluster, time period cluster, and data credibility cluster, which are then concatenated to form a complete input feature vector. Continuous features are processed using robust median-interquartile range standardization. Training, validation, and test sets are divided by site / time block, and random splitting of consecutive time blocks at the same site is prohibited. Regression models were trained independently for 2m air temperature, ground air pressure, specific humidity, and 10m wind component. Random forest, gradient boosting tree, or neural network were selected for the models. After training, the feature order, standardized parameters, and residual amplitude threshold were saved.

Images