A dynamic surface data updating method based on numerical weather forecast model

Abstract

The present invention discloses a method for updating dynamic surface data based on a numerical weather forecast model. The present invention collects the static geographic data and the corresponding atmospheric initial field and boundary field data required for the study area and performs spatial interpolation and preprocessing on the static geographic data; decodes and processes the meteorological data based on the geographical environment characteristics and climate background conditions of the study area to generate boundary conditions that meet the simulation requirements; determines the appropriate physical process parameterization scheme and numerical calculation method based on the hydrological and meteorological characteristics of the target area and the research purpose, and then couples the dynamic surface update module with the WRF model. During the operation of the WRF model, the dynamic surface update module is called to continuously update the changes in the water surface range and water level caused by reservoir flooding during the simulation period. Based on the parallel simulation principle of the WRF model, the present invention improves the simulation accuracy and enhances the numerical stability by considering the impact of the dynamic surface on the atmospheric process.

Description

Technical Field

[0001] The present invention relates to the field of land-atmosphere coupling process and the technical field of numerical weather forecasting and climate simulation, and in particular to a dynamic surface data updating method based on a numerical weather forecast model. Background Art

[0002] Land surface processes and their interactions with climate have long been an important area of ​​scientific research. With the continuous development and improvement of observational and satellite remote sensing data, a growing body of research has demonstrated the importance of land surface processes in predicting climate change and in determining adaptation measures. However, current mainstream numerical weather prediction models, such as the WRF model, typically load static surface data during initialization, primarily consisting of surface elevation and land cover information, and maintain this static data throughout the simulation.

[0003] However, in practical applications, scenarios involving rapid, dynamic surface evolution often arise. These include, but are not limited to, changes in local topography and water cover caused by large-scale water conservancy projects, such as the Three Gorges Reservoir's periodic storage and release of water, which can lead to significant fluctuations in regional water levels and shorelines. Land use shifts caused by rapid urbanization, including the rapid conversion of agricultural land to urban construction, can significantly alter surface energy and water flux characteristics. Extreme natural disasters (such as landslides, debris flows, and floods) can cause short-term, drastic changes in topography, roughness, and surface physical properties. In these scenarios, traditional static surface processing methods are unable to promptly reflect the feedback effects of surface changes on the atmospheric system, resulting in significant deviations from actual simulation results. This is particularly evident in key areas such as local precipitation distribution, boundary layer thermal structure, and surface energy and water exchange. These deviations are particularly pronounced in specific application areas such as coupled hydrological simulation, disaster warning, and urban climate analysis, significantly limiting the effectiveness of model applications.

[0004] Currently, numerical weather models, such as WRF, lack an effective mechanism for dynamically adjusting surface properties in real time during runtime. In particular, in a distributed parallel computing environment, there are still significant technical difficulties and challenges in efficiently and securely implementing local refreshes of dynamic surface properties, parallel data synchronization, and consistent updates of physically derived variables. Therefore, there is an urgent need for a method that can dynamically update underlying surface characteristics on demand during the numerical model integration process while ensuring the consistency and numerical stability of the physical processes. This is intended to improve the accuracy of numerical weather simulations under dynamically changing surface scenarios, thereby broadening the application potential of numerical models in multiple fields such as environmental change, disaster risk assessment, and urban sustainable development. Summary of the Invention

[0005] The purpose of the present invention is to provide a dynamic surface data updating method based on a numerical weather forecast model. Through progressive updating, adaptive mixing factors, distributed data exchange and a complete error handling mechanism, the method can achieve efficient and stable dynamic updating of terrain and land use data, improve the accuracy of meteorological numerical simulation, and provide more accurate surface boundary conditions for weather forecast and climate simulation.

[0006] In order to solve the above technical problems, the present invention provides the following technical solutions:

[0007] A method for updating dynamic surface data based on a numerical weather forecast model, comprising the following steps:

[0008] During the data preparation phase, the static geographic data required for the study area and the corresponding atmospheric initial field and boundary field data were collected; the static geographic data included elevation data, land cover data, soil data, vegetation leaf area index, albedo, wind speed, geopotential height, temperature, pressure, sea level pressure, relative humidity, 2m air temperature, surface air temperature, four-layer soil temperature, four-layer soil moisture, etc.

[0009] In the static geographic data preprocessing stage, the spatial interpolation and preprocessing of static geographic data are performed, and a geographic background field containing information such as terrain height and land use type is generated by setting the spatial resolution, determining the location of the nested area, and the land use classification method;

[0010] In the meteorological data decoding stage, the meteorological data is decoded based on the generated geographic background field;

[0011] In the boundary generation phase, global-scale meteorological data are interpolated into regional simulation grids, and multi-time data are interpolated one by one to ensure smooth boundary transition and avoid computational instability. The geographic background field and decoded meteorological field data are integrated, and the meteorological data are converted from their original coordinate system to the map projection coordinate system used by the WRF model simulation to ensure spatial resolution matching. The vertical hierarchy is adjusted to the sigma coordinate system required by the WRF simulation, thereby generating the horizontal grid dataset required for WRF model operation initialization. The horizontal grid dataset provides a complete three-dimensional atmospheric state at the start of the simulation for each nested simulation area, and contains the boundary condition time series of each nested simulation area throughout the simulation period to support the subsequent time integration operation of the WRF model, ensuring the continuity of boundary conditions and the consistency of physical state during the simulation process.

[0012] In the physical process and numerical simulation setting stage, according to the research objectives and application scenario analysis, the physical parameterization process scheme is optimized, and the numerical calculation method is selected for model configuration;

[0013] In the coupled dynamic surface update phase, the dynamic surface update module is coupled with the WRF model. The dynamic surface update module is called during the operation of the WRF model to continuously update the changes in water surface range and water level caused by reservoir inundation during the simulation period.

[0014] According to the above technical solution, the steps executed by the dynamic surface update module include:

[0015] Parse scheduling information and construct dynamic elevation data;

[0016] Construct dynamic land cover maps based on dynamic elevation data.

[0017] According to the above technical solution, the step of constructing dynamic elevation data includes:

[0018] Obtaining actual operational data of the reservoir over many years and using the Lagrange interpolation method to supplement missing data; the actual operational data includes: water level-capacity relationship curve data of the reservoir, daily water level process line of the reservoir, normal water level of the reservoir, dead water level, and flood control limit water level;

[0019] Extract DEM grid flooded areas based on actual dispatching operation data;

[0020] For each grid point in the flooded area, calculate its water depth value and update the DEM elevation data;

[0021] A reservoir underlying surface topography dataset at a daily scale is integrated, wherein the reservoir underlying surface topography dataset includes a basic topography elevation layer, an updated topography elevation layer, a water depth distribution layer, and a quality control mark layer.

[0022] According to the above technical solution, the DEM grid flooded area extraction step includes:

[0023] Perform pothole filling and filtering on the DEM elevation data to remove noise and outliers, and demarcate the reservoir dam site in the DEM elevation data;

[0024] Based on the water level-storage capacity relationship curve of the reservoir, a corresponding relationship table between water level and flooded area is established to determine the water catchment range and boundary threshold of the reservoir;

[0025] Based on the reservoir's daily water level process data, GIS spatial analysis capabilities were used to extract all pixels below the set water level in the DEM elevation data as inundated areas. Connectivity analysis was then applied to ensure that the extracted areas were connected to the reservoir theme. Boundary and range constraints were imposed based on the reservoir water level-capacity-area relationship curve, and a daily-scale reservoir boundary polygon dataset was established.

[0026] Gaussian smoothing is performed on the daily flooded area boundaries to ensure that there are no jagged edges in the flooded area boundaries;

[0027] Based on the daily-scale reservoir boundary polygon dataset, the water depth value is calculated for each grid point within the reservoir boundary and a daily-scale water depth dataset is generated.

[0028] According to the above technical solution, the water depth calculation formula is:

[0029] ;

[0030] Where, represents the water depth at the (i, j)th grid point that is submerged by water after the reservoir is built and put into operation. It indicates the water level at the current time after the reservoir is constructed and put into operation. Indicates the original terrain elevation value at the (i, j)th grid point when it is not submerged by water.

[0031] According to the above technical solution, the DEM elevation data updating step includes:

[0032] when When , it means that the grid point is not actually submerged, and it is set to the minimum water depth, while the original DEM elevation data is still maintained;

[0033] when When , it means that the grid point is flooded. For the flooded area, the original DEM elevation data is replaced with the water level value at the current moment to achieve a flat water surface effect;

[0034] In the transition area at the edge of the reservoir, the distance-weighted average method is used to achieve a smooth transition of the DEM elevation data:

[0035] ;

[0036] in, Represents the surface elevation data of the (i, j)th grid point after update, represents the original terrain elevation value at the (i, j)th grid point when it is not submerged by water; w is a distance-based weight coefficient (value range: 0≤w≤1), indicating the degree to which the current grid point is affected by water.

[0037] ;

[0038] Where d is the distance from the grid point to the water body boundary, reflecting the spatial relationship between the grid point and the water body boundary; σ is a smoothing parameter, which is determined according to the actual application scenario and is usually a few to dozens of grid units.

[0039] According to the above technical solution, the step of constructing the dynamic land cover map includes:

[0040] Based on MODIS land cover data, the land use type data of the study area were clipped according to the simulation area set in the WRF model;

[0041] Water mask data is constructed based on the daily-scale reservoir boundary polygon dataset created in the dynamic elevation data construction step. Specifically, grid points within the reservoir boundary polygon are considered water grid points and their value is set to 2. The remaining grid points outside the reservoir boundary polygon but within the study area are considered land grid points and their value is set to 1, thus obtaining daily-scale water mask data. The water mask data refers to the extent of the reservoir's submerged waters, as extracted previously.

[0042] The original land use data were then recoded, and the water mask data with relevant spatial range and resolution were fused with the land use data. Specifically, the spatial location of the water body grid point with a median value of 2 in the water mask data was determined, and the corresponding grid point with the same spatial location in the land use data was recoded to 17 (the grid point with a land use type of 17 is considered a water body in the WRF model). The other grid points (i.e., the grid points with the same location as the grid point with a median value of 1 in the water mask data) retained the original land use type code.

[0043] By fusing the daily-scale water mask dataset and the land use data of the study area, a dynamic land cover dataset containing the daily-scale flooded water area information of the reservoir is generated.

[0044] According to the above technical solution, the steps of executing the coupled dynamic surface update phase include:

[0045] Data preparation and preprocessing: producing elevation and land cover data corresponding to the time point and coverage range according to the simulation area, time period and update frequency set in the WRF model;

[0046] Module initialization: Initialize in WRF mode and build update indexes, read / write channels and data structures, and read initial terrain and land cover status;

[0047] Dynamic update trigger judgment: in each major time step of the WRF mode, it is determined whether the current time triggers the surface update task, and whether the current area has been updated;

[0048] File detection and reading: automatically read external elevation and land cover change files at the trigger time point. The reading range supports global or local changes, and detects, fills in and corrects the data validity.

[0049] Local grid update and distributed data exchange: perform local surface variable updates for the changed area and synchronize the updated local surface variables to other processing nodes;

[0050] The surface physical quantity refresh mechanism uses adaptive mixing factor technology to dynamically adjust the mixing coefficient according to the magnitude of terrain changes, and update the surface data in a progressive manner;

[0051] Physical consistency adjustment and verification: the adjusted physical variables are checked hierarchically based on the principles of thermodynamic and static equilibrium. If an anomaly is found, a rollback operation will be triggered to restore the physical variables to a safe state. Specifically, the physical consistency adjustment is achieved through a phased and conservative physical field update strategy: based on the physical field calculation rules in the WRF mode, the physical field under the new underlying surface conditions is updated through the following stages: basic pressure field adjustment: the surface pressure is calculated according to the new terrain to ensure that the vertical distribution of the pressure is monotonically decreasing; temperature and inverse density field adjustment: using segmented defined temperature profiles to consider the influence of water vapor; potential height field adjustment: based on the static equilibrium equation, to ensure that the potential height is monotonically increasing; pressure perturbation and total pressure calculation: the gas state equation is applied to calculate the new pressure perturbation; final consistency check: verify the rationality and monotonicity of each physical field.

[0052] Error handling and rollback mechanism. The program has a complete built-in error monitoring, counting and rollback mechanism. Once the physical field exceeds the limit or boundary overflow occurs, the system will call the backup data for recovery to ensure that the simulation is not interrupted. The specific steps are to establish a structured error classification system, including physical errors (situations where the physical field violates the laws of physics, such as negative air pressure and negative temperature), boundary errors (array index out of range) and monotonicity errors (vertical distribution is not monotonic, such as air pressure not decreasing with height); set multi-level error thresholds to track physical errors, boundary errors and monotonicity errors; before updating the physical field, back up the key physical field in advance, and when the error exceeds the set threshold during the update process, roll back to the initial state; use segmented recovery to avoid large-scale data operations, and synchronize the rollback status of all processes through the file system or MPI to ensure global consistency.

[0053] Module resource cleanup: After the simulation is completed, all dynamic memory space requested by the update module is released, temporary variables and communication cache are cleared, and update log information is recorded.

[0054] According to the above technical solution, the mixing factor is determined by the terrain change amplitude and the vertical layer position. , to achieve a smooth transition of the physical field, the mixing factor is calculated by the following formula:

[0055] ;

[0056] Where i,j are horizontal grid positions, k is the vertical layer position, is the change in terrain height, is a change in physical quantity, represents the basic mixing factor function, represents the vertical layer correction function, Represents the correction function for physical quantity changes. The terrain height change is calculated using the following formula:

[0057] ;

[0058] The basic mixing factor function implements a three-level dynamic mixing strategy based on the magnitude of terrain change, dynamically adjusting the basic mixing factor coefficient based on the maximum terrain change according to the following rules:

[0059] ;

[0060] The vertical layer correction function adjusts the blending factor according to the vertical height to achieve differentiated processing in the vertical direction, which is determined by the following formula:

[0061] ;

[0062] Where, is the attenuation coefficient, which has different values ​​according to different physical quantities (0.3-0.4), is the bottom index in the vertical direction, is the top index in the vertical direction;

[0063] Changes in physical quantities Calculated by the following formula:

[0064] ;

[0065] Where, represents the physical quantity calculated based on the updated underlying surface data, Indicates the original physical quantity when the underlying surface data is not updated. It can be physical quantities such as temperature (T), air pressure (P), potential height (Φ), etc.

[0066] It is a physical quantity change correction function, which is used to further adjust the mixing factor according to the actual change amplitude of the physical quantity. Its value follows the following rules:

[0067] ;

[0068] Where, is the threshold value of physical quantity change;

[0069] Final value of the physical quantity after updating the underlying surface data Calculated by the following formula:

[0070] .

[0071] Compared with the prior art, the present invention has the following beneficial effects:

[0072] 1. The full-process dynamic surface update capability enables, for the first time, dynamic support for adjusting underlying surface information at any time, in any region, and at any number of times during the WRF time step, making the simulation more closely aligned with real-world geographic changes. This significantly improves the forecast accuracy of surface-sensitive weather systems (such as local circulation and precipitation).

[0073] 2. An adaptive hybrid interpolation strategy controls the slope of surface variable changes through multiple factor weights, achieving progressive updates and adaptive hybrid strategies. This effectively avoids numerical fluctuations caused by sudden surface changes and maintains the computational stability of the model.

[0074] 3. A complete parallel communication guarantee mechanism, compatible with WRF's parallel operation architecture, ensures data boundary consistency and communication efficiency. The improved distributed data exchange mechanism significantly reduces communication overhead, reduces deadlock risks, and improves the efficiency of large-scale parallel computing;

[0075] 4. Automatic physical variable refresh link: automatically adjusts related physical variables when surface information changes, maintains energy conservation, momentum continuity, and water balance, ensures overall simulation physical consistency, and avoids human modification and omissions;

[0076] 5. Robust rollback fault tolerance mechanism, which can roll back and restart in the event of data anomalies or calculation failures, avoiding simulation termination or crash;

[0077] 6. System Applicability: This method is applicable to WRF simulations of different resolutions and in different regions. It has broad application prospects and can provide important technical support for fields such as weather forecasting, climate simulation, and environmental assessment. BRIEF DESCRIPTION OF THE DRAWINGS

[0078] The accompanying drawings are used to provide a further understanding of the present invention and constitute a part of the specification. Together with the embodiments of the present invention, they are used to explain the present invention and do not constitute a limitation of the present invention. In the accompanying drawings:

[0079] Figure 1 It is a flowchart of the steps of coupling operation of the dynamic surface update module in the WRF mode in the present invention. DETAILED DESCRIPTION

[0080] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0081] The study area for this proposal was the Three Gorges Reservoir, using ERA-5 reanalysis data provided by the European Centre for Medium-Range Weather Forecasts (ECMWF). ERA5 reanalysis data have a spatial resolution of 0.25° × 0.25° and a temporal resolution of 1 hour. The data spans from 1979 to the present, including surface boundary field data and 37 vertical layers of standard isobaric surface altitude data from 1000 hPa to 1 hPa. The static data uses MODIS land cover data, which includes 20 land use types. Elevation data with a spatial resolution of 30 meters were acquired from the Shuttle Radar and Topography Mission (SRTM30).

[0082] The surface data is updated by a dynamic surface data updating method based on a numerical weather forecast model. The specific steps include:

[0083] S1. Convert high-resolution terrain / land cover data into static field files consistent with the model's horizontal grid, projection, and nested structure. First, determine the geographic boundaries of the simulation area (the latitude and longitude range, encompassing the entire Three Gorges Reservoir area and the surrounding mountainous area), set the simulation grid resolution and number of grid points, and determine the map projection (e.g., Lambert or Mercator). Then, interpolate the multi-source static geographic data onto the simulation grid points to ensure smooth terrain data and avoid numerical instability.

[0084] S2. For the generated geographic background field, the meteorological data from the ERA5 surface and upper-air fields are converted into an intermediate format that can be directly read by the model. The meteorological elements and layers required for the simulation are extracted, the temporal continuity and spatial integrity of the data are checked, missing data and outliers are processed, and finally the meteorological data are interpolated onto the regional simulation grid to generate the horizontal grid dataset required for WRF model initialization. The horizontal grid dataset provides the complete three-dimensional atmospheric state at the simulation start time for each nested simulation region and contains the time series of boundary conditions for each nested simulation region throughout the simulation period. This supports the subsequent time integration of the WRF model and ensures the continuity of boundary conditions and the consistency of physical states during the simulation.

[0085] S3, the physical process and numerical simulation setup phase, specifically sets the simulation start and end times and total simulation duration, determines the model integration time step (usually based on the grid resolution), sets the model output frequency and saved variables, and also sets the number of vertical layers and the inter-layer distance distribution. Based on this, physical parameterization schemes are determined, including microphysical process schemes, radiative transfer schemes, planetary boundary layer schemes, land surface process schemes, and cumulus convection parameterization schemes. Appropriate numerical calculation methods are then selected for model configuration, including the selection of a discrete format for the advection term, configuration of the time integration scheme, treatment of lateral boundary conditions, setting of the top absorbing layer of the computational domain, and parallel computation partitioning strategies.

[0086] S4. Couple the dynamic surface update module with the WRF model and start the dynamic surface update module synchronously when the WRF model is initialized. During the operation of the WRF model, the dynamic surface update module is called to continuously update the changes in the water surface range and water level caused by the reservoir flooding during the simulation period. The specific steps include:

[0087] S401: Data preparation and preprocessing: According to the simulation area, time period and update frequency set in the WRF model, the elevation and land cover data corresponding to the time point and coverage range are produced.

[0088] S402: Module initialization: Initialize in WRF mode and build update indexes, read / write channels, and data structures to read initial terrain and land cover status.

[0089] S403, dynamic update trigger judgment, in each main time step of the WRF mode, determine whether the current time triggers the surface update task, and check whether the current area has been updated; specifically, in the main time loop, each time step enters the dynamic update process to determine whether the specified update time has been reached, allowing a time error of ±1 minute to ensure correct update even if the time step is irregular. Considering seasonal factors, the inspection frequency can be increased during specific periods.

[0090] S404, file detection and reading. When an update is required, the data reading process is executed, the timeout period is set to 60 seconds, the number of retries is set to 5, and the external elevation and land cover change files are automatically read at the trigger time point. The reading range supports global or local changes, and the data validity is detected, filled and corrected.

[0091] S405. Execute local surface variable updates for the changed area and synchronize them to other processing nodes through subroutines such as secure parallel block broadcast. Specifically, a batch broadcast strategy is used in a parallel environment: each batch processes no more than 1,000 rows of data, uses non-blocking communication to reduce waiting time, and sets communication timeout and retry mechanisms to ensure parallel consistency.

[0092] MPI communication is optimized through the following hierarchical batching strategy:

[0093] Define global constants and limit the maximum amount of data for a single MPI communication to prevent communication congestion. Create a data broadcast function that receives the data array, number of elements, root process ID, communicator, and success flag as parameters. When the amount of data to be transmitted exceeds one-quarter of the global constant, start the block broadcast strategy: divide the data into multiple blocks of one-quarter the global constant in size, perform a separate broadcast operation for each data block, add a short delay to the inter-block broadcast, reduce network pressure, and improve the success rate of large-scale communication. Continuously monitor the error status during the broadcast process. When a batch fails to be broadcast, record the error but continue to process the remaining batches, and fill the failed area with reasonable default values ​​to ensure that the program will not be interrupted by communication errors.

[0094] S406, the surface physical quantity refresh mechanism uses adaptive mixing factor technology to dynamically adjust the mixing coefficient according to the magnitude of terrain changes, and updates the surface data in a progressive manner to avoid numerical oscillations.

[0095] Among them, the mixing factor is determined by the terrain change amplitude and the vertical layer position. , to achieve a smooth transition of the physical field, the mixing factor is calculated by the following formula:

[0096] ;

[0097] Where i,j are horizontal grid positions, k is the vertical layer position, is the change in terrain height, is a change in physical quantity, represents the basic mixing factor function, represents the vertical layer correction function, Represents the correction function for physical quantity changes. The terrain height change is calculated using the following formula:

[0098] ;

[0099] The basic mixing factor function implements a three-level dynamic mixing strategy based on the magnitude of terrain change, dynamically adjusting the basic mixing factor coefficient based on the maximum terrain change according to the following rules:

[0100] ;

[0101] The vertical layer correction function adjusts the blending factor according to the vertical height to achieve differentiated processing in the vertical direction, which is determined by the following formula:

[0102] ;

[0103] Where, is the attenuation coefficient, which has different values ​​according to different physical quantities (0.3-0.4), is the bottom index in the vertical direction, is the top index in the vertical direction;

[0104] Changes in physical quantities Calculated by the following formula:

[0105] ;

[0106] Where, represents the physical quantity calculated based on the updated underlying surface data, Indicates the original physical quantity when the underlying surface data is not updated. It can be physical quantities such as temperature (T), air pressure (P), potential height (Φ), etc.

[0107] It is a physical quantity change correction function, which is used to further adjust the mixing factor according to the actual change amplitude of the physical quantity. Its value follows the following rules:

[0108] ;

[0109] Where, is the threshold value of physical quantity change;

[0110] The final value of the physical quantity after the underlying surface data is updated ( ) is calculated using the following formula:

[0111] .

[0112] Physical consistency adjustment and verification: The adjusted physical variables are checked at different levels based on the thermodynamic and static equilibrium principles. If an anomaly is found, a rollback operation will be triggered to restore the physical variables to a safe state.

[0113] At the same time, WRF automatically triggers the update of surface-related physical variables and calls the physical field update subroutine. Physical variables include roughness, height, temperature, albedo, surface flux, etc. to ensure thermodynamic-dynamic consistency. Specifically, the physical consistency adjustment is achieved through a phased and conservative physical field update strategy:

[0114] Based on the physical field calculation rules in the WRF model, the physical fields under the underlying surface conditions are updated through the following stages: basic pressure field adjustment: surface pressure is calculated according to the new terrain to ensure that the vertical distribution of air pressure is monotonically decreasing; temperature and inverse density field adjustment: using segmented temperature profiles to consider the influence of water vapor; potential height field adjustment: based on the static equilibrium equation, ensuring that the potential height is monotonically increasing; pressure perturbation and total pressure calculation: applying the gas state equation to calculate the new pressure perturbation; final consistency check: verifying the rationality and monotonicity of each physical field.

[0115] S408, error handling and rollback mechanism. The program has a complete built-in error monitoring, counting and rollback mechanism. Once the physical limit is exceeded or the boundary overflow occurs, the system will call the backup data for recovery to ensure that the simulation is not interrupted.

[0116] The error detection and rollback mechanism is implemented through the following steps:

[0117] Establish a structured error classification system, including physical errors (situations where the physical field violates physical laws, such as negative air pressure and negative temperature), boundary errors (array index out of range), and monotonicity errors (vertical distribution is not monotonic, such as air pressure not decreasing with altitude); set multi-level error thresholds to track physical errors, boundary errors, and monotonicity errors. For example, perform boundary and physical rationality checks on updated variables (such as altitude monotonicity and air pressure non-negativity). Set multi-level error thresholds:

[0118] Abnormal temperature: T<-50°C or T>60°C

[0119] Abnormal air pressure: P<500hPa or P>1100hPa

[0120] Monotonicity check: air pressure should decrease with altitude, and geopotential height should increase with altitude

[0121] Physical quantity rationality: humidity range [0,1], wind speed less than 100m / s,

[0122] If the conditions are not met, the rollback function is automatically executed to restore all physical fields to the state before the update, record error information, and generate a diagnostic report for analysis to a safe state; segmented recovery is used to avoid large amounts of data operations, and the rollback status of all processes is synchronized through the file system or MPI to ensure global consistency.

[0123] S409, module resource cleanup: after the simulation is complete, all dynamic memory space requested by the update module is released, temporary variables and communication cache are cleared, and update log information is recorded. The specific steps include:

[0124] (1) Release dynamically allocated memory space, including temporary arrays, cache space, and communication buffers

[0125] (2) Close all open data files to ensure data integrity

[0126] (3) Generate update statistics report, including: update times and time distribution, success rate and failure reasons, computing efficiency statistics, resource usage

[0127] (4) Clean up temporary files and synchronization directories, free up disk space, release dynamically allocated memory, and record update logs to ensure system stability.

[0128] After the model is completed, the simulation results are tested using observation data (meteorological stations, radar, and satellite products), and error indicators such as RMSE, NSE, and Bias are calculated to verify the improvement effect of the dynamic update scheme on simulation accuracy and spatial response.

[0129] The implementation of the dynamic surface update module is mainly achieved through the following steps:

[0130] Step 1: Scheduling information analysis and dynamic elevation data (DEM) construction

[0131] 1.1 Preparation of water level data: Obtain the actual operation data of the reservoir over many years, including daily (or hourly) water level process lines. During the data preparation process, Lagrange interpolation method is used to supplement the missing data as the input basis for dynamic changes;

[0132] ;

[0133] Where L(x) is the water level at the interpolation point, x represents the time point to be interpolated (the time corresponding to the missing data), and n represents the number of known data points minus 1 (if there are n+1 known points, then n is the highest order term). is the i-th known time point (i=0,1,2,...,n), is the water level value corresponding to the i-th known time point, is the jth known time point, which is used to construct the Lagrangian basis function.

[0134] 1.2 DEM grid inundation extraction: Based on SRTM 30-meter resolution DEM elevation data and reservoir capacity curves, the inundation area is extracted according to the following steps:

[0135] (1) Perform pothole filling processing on DEM elevation data and use the Wang-Liu algorithm to eliminate false potholes in DEM elevation data;

[0136] (2) Determine the possible inundation area based on the historical flood range and establish the maximum range polygon of the reservoir boundary;

[0137] (3) Based on the given water level, the seed filling algorithm is used to extract the flooded area: starting from the center of the reservoir and expanding to the surrounding areas, all grids with ground elevations lower than the given water level are marked as flooded areas;

[0138] (4) Gaussian smoothing is performed on the boundary of the submerged area, and the smoothing radius is set to 3 pixels to ensure that the boundary of the submerged area has no jagged edges;

[0139] (5) Based on the reservoir water level-capacity-area relationship curve and the maximum range polygon of the reservoir boundary, boundary and range constraints are imposed, and a daily-scale reservoir boundary polygon dataset is established based on the extracted flooded area range.

[0140] 1.3 Water depth calculation and DEM elevation data update: For each grid point in the flooded area, calculate its water depth value:

[0141] ;

[0142] Where, represents the water depth at the (i, j)th grid point that is submerged by water after the reservoir is built and put into operation. It indicates the water level at the current time after the reservoir is constructed and put into operation. Indicates the original terrain elevation value at the (i, j)th grid point when it is not submerged by water.

[0143] The steps to update DEM elevation data include:

[0144] when When , it means that the grid point is not actually submerged, and it is set to the minimum water depth, while the original DEM elevation value is still maintained;

[0145] when When , it means that the grid point is flooded. For the flooded area, the original DEM elevation data is replaced with the water level value at the current moment to achieve a flat water surface effect;

[0146] In the transition area at the edge of the reservoir, distance-weighted averaging is used to achieve a smooth transition:

[0147] ;

[0148] in, Represents the surface elevation data of the (i, j)th grid point after update, represents the original terrain elevation value at the (i, j)th grid point when it is not submerged by water; w is a distance-based weight coefficient (value range: 0≤w≤1), indicating the degree to which the current grid point is affected by water.

[0149] ;

[0150] Where d is the distance from the grid point to the water body boundary, reflecting the spatial relationship between the grid point and the water body boundary; σ is a smoothing parameter, which is determined according to the actual application scenario and is usually a few to dozens of grid units.

[0151] 1.4 Generate daily-scale DEM elevation dataset: Assemble the reservoir underlying surface topography dataset at daily scale, covering the entire simulation process, including: basic terrain elevation layer (original_elevation), water depth distribution layer (water_depth), updated terrain elevation layer (updated_elevation), quality control flag layer (quality_flag), and mark the data source and processing method.

[0152] Step 2: Dynamic land cover map construction

[0153] 2.1 Original land cover processing: Use MODIS MCD12Q1 product (500 m resolution) or high-resolution satellite data (such as GF-1) to crop the land types in the study area.

[0154] 2.2 Land type recoding: Based on the daily-scale reservoir boundary polygon dataset extracted in step 1, water body mask data are constructed. The grid points within the reservoir boundary polygon are regarded as water body grid points, and their values ​​are set to 2. The remaining grid points outside the reservoir boundary polygon and within the study area are regarded as land grid points, and their values ​​are set to 1, thus obtaining the daily-scale water body mask data.

[0155] The original land use data were then recoded, and the water mask data with relevant spatial range and resolution were fused with the land use data. Specifically, the spatial location of the water body grid point with a median value of 2 in the water mask data was determined, and the corresponding grid point with the same spatial location in the land use data was recoded to 17 (the grid point with a land use type of 17 is considered a water body in the WRF model). The other grid points (i.e., the grid points with the same location as the grid point with a median value of 1 in the water mask data) retained the original land use type code.

[0156] By fusing the daily-scale water mask dataset and the land use data of the study area, a dynamic land cover dataset containing the daily-scale flooded water area information of the reservoir is generated.

[0157] 2.3 Generate daily-scale land cover dataset: Generate an independent land use data file for each simulation date to cover the entire simulation process.

[0158] It should be noted that, in this document, relational terms such as first and second, etc., are used only to distinguish one entity or operation from another entity or operation, and do not necessarily require or imply any actual relationship or order between these entities or operations. Moreover, the terms "comprises," "comprising," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that includes a list of elements includes not only those elements but also other elements not explicitly listed, or elements inherent to such process, method, article, or apparatus.

[0159] Finally, it should be noted that the above descriptions are merely preferred embodiments of the present invention and are not intended to limit the present invention. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art will be able to modify the technical solutions described in the aforementioned embodiments or substitute equivalents for some of the technical features. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention shall be included within the scope of protection of the present invention.

Claims

1. A dynamic surface data updating method based on a numerical weather forecast model, characterized in that: The steps include: In the data preparation stage, the static geographic data and corresponding atmospheric initial field and boundary field data required for the study area are collected; In the static geographic data preprocessing stage, spatial interpolation and preprocessing of static geographic data are performed based on the underlying surface and climate background of the study area to generate a geographic background field; In the meteorological data decoding stage, the meteorological data is decoded based on the generated geographic background field; In the boundary generation phase, global-scale meteorological data are interpolated onto the regional simulation grid. Multi-time data are interpolated one by one, and the geographic background field and decoded meteorological field data are integrated and interpolated onto the WRF model operation grid. This generates a horizontal grid dataset that integrates the initial conditions and boundary conditions required for WRF model operation, realizing the real data initialization required for WRF model simulation. In the physical process and numerical simulation setting stage, according to the research objectives and application scenario analysis, the physical parameterization process scheme is optimized, and the numerical calculation method is selected for model configuration; In the coupled dynamic surface update phase, the dynamic surface update module is coupled with the WRF model. The dynamic surface update module is called during the WRF model operation to continuously update the changes in the water surface range and water level caused by reservoir inundation during the simulation period. The steps of executing the coupled dynamic surface update phase include: Data preparation and preprocessing: producing elevation and land cover data corresponding to the time point and coverage range according to the simulation area, time period and update frequency set in the WRF model; Module initialization: Initialize in WRF mode and build update indexes, read / write channels and data structures, and read initial terrain and land cover status; Dynamic update trigger judgment: in each major time step of the WRF mode, it is determined whether the current time triggers the surface update task, and whether the current area has been updated; File detection and reading: automatically read external elevation and land cover change files at the trigger time point. The reading range supports global or local changes, and detects, fills in and corrects the data validity. Local grid update and distributed data exchange: perform local surface variable updates for the changed area and synchronize the updated local surface variables to other processing nodes; The surface physical quantity refresh mechanism determines the mixing factor based on the magnitude of terrain change and the vertical layer position. Using adaptive mixing factor technology, WRF dynamically adjusts the mixing coefficient according to the magnitude of terrain change and updates the surface data in a progressive manner. Physical consistency adjustment and verification: The adjusted physical variables are checked at different levels based on the thermodynamic and static equilibrium principles. If an anomaly is found, a rollback operation will be triggered to restore the physical variables to a safe state. Error handling and rollback mechanism: The program has a complete built-in error monitoring, counting and rollback mechanism. Once the physical limit is exceeded or the boundary overflow occurs, the system will call the backup data for recovery to ensure that the simulation is not interrupted. Module resource cleanup: After the simulation is completed, all dynamic memory space requested by the update module is released, temporary variables and communication cache are cleared, and update log information is recorded; The mixing factor Calculation formula: ； Where i and j are horizontal grid positions, k is the vertical layer position, is the change in terrain height, is a change in physical quantity, represents the basic mixing factor function, represents the vertical layer correction function, Represents the correction function for changes in physical quantities.

2. A method for updating dynamic surface data based on a numerical weather forecast model according to claim 1, characterized in that: The steps of executing the dynamic surface update module include: Parse scheduling information and construct dynamic elevation data; Construct dynamic land cover maps based on dynamic elevation data.

3. A method for updating dynamic surface data based on a numerical weather forecast model according to claim 2, characterized in that: The steps of constructing dynamic elevation data include: Obtaining multiple years of actual reservoir operation data and using Lagrange interpolation to complete missing data; the actual operation data includes: reservoir water level-storage capacity relationship curve data, reservoir daily water level process line, reservoir normal water level, dead water level, and flood control limit water level; Extract DEM grid flooded areas based on actual dispatching operation data; For each grid point in the flooded area, calculate its water depth value and update the DEM elevation data; The reservoir underlying surface topography data set at a daily scale is integrated, and the reservoir underlying surface topography data includes a basic topography elevation layer, an updated topography elevation layer, a water depth distribution layer, and a quality control mark layer.

4. A method for updating dynamic surface data based on a numerical weather forecast model according to claim 3, characterized in that: The DEM grid flooded area extraction step includes: Perform pothole filling and filtering on the DEM elevation data to remove noise and outliers, and demarcate the reservoir dam site in the DEM elevation data; Based on the water level-storage capacity relationship curve of the reservoir, a corresponding relationship table between water level and flooded area is established to determine the water catchment range and boundary threshold of the reservoir; Based on the reservoir's daily water level process data, GIS spatial analysis capabilities were used to extract all pixels below the set water level in the DEM elevation data as inundated areas. Connectivity analysis was then applied to ensure that the extracted areas were connected to the reservoir theme. Boundary and range constraints were imposed based on the reservoir water level-capacity-area relationship curve, and a daily-scale reservoir boundary polygon dataset was established. Gaussian smoothing is performed on the daily flooded area boundaries to ensure that there are no jagged edges in the flooded area boundaries; Based on the daily-scale reservoir boundary polygon dataset, the water depth value is calculated for each grid point within the reservoir boundary and a daily-scale water depth dataset is generated.

5. The method for updating dynamic surface data based on a numerical weather forecast model according to claim 3, characterized in that: The calculation formula of the water depth value corresponding to the grid point is: ； Where, represents the water depth at the (i, j)th grid point that is submerged by water after the reservoir is built and put into operation. It indicates the water level at the current time after the reservoir is constructed and put into operation. Indicates the original terrain elevation value at the (i, j)th grid point when it is not submerged by water.

6. A method for updating dynamic surface data based on a numerical weather forecast model according to claim 3, characterized in that: The DEM elevation data updating step includes: when When , it means that the grid point is not actually submerged, and it is set to the minimum water depth, while the original DEM elevation data is still maintained; when When , it means that the grid point is flooded. For the flooded area, the original DEM elevation data is replaced with the water level value at the current moment to achieve a flat water surface effect; In the transition area at the edge of the reservoir, the distance-weighted average method is used to achieve a smooth transition of the DEM elevation data: ； in, Represents the elevation data of the (i, j)th grid point after update, represents the original elevation value at the (i, j)th grid point when it is not submerged by water; w is the distance-based weight coefficient, which indicates the degree to which the current grid point is affected by the water body; ； Where d is the distance from the grid point to the water body boundary, reflecting the spatial relationship between the grid point and the water body boundary; σ is the smoothing parameter.

7. The method for updating dynamic surface data based on a numerical weather forecast model according to claim 2, characterized in that: The dynamic land cover map construction step includes: Based on MODIS land cover data, the land use type data of the study area were clipped according to the simulation area set in the WRF model; Construct water body mask data based on daily-scale reservoir boundary polygon dataset; Recode the original land use data and fuse the water mask data and land use data with relevant spatial range and resolution; By fusing the daily-scale water mask dataset and the land use data of the study area, a dynamic land cover dataset containing the daily-scale flooded water area information of the reservoir is generated.

Images