A small-scale meteorological numerical simulation method based on large eddy simulation

Abstract

The application belongs to the field of weather forecast, and discloses a small-scale weather numerical simulation method based on large eddy simulation, which comprises the following steps: obtaining driving field data of a target area; performing data preprocessing to obtain initial weather field data which can be interpolated to a three-dimensional grid system with a preset resolution; performing vortex separation on the initial weather field data based on a spatial filtering function to obtain large-scale vortices and small-scale vortices, wherein a filtering scale threshold is determined by the grid resolution of the three-dimensional grid system; performing explicit numerical solution on the large-scale vortices by using a filtered Navier-Stokes equation, and performing closed solution on the small-scale vortices by introducing a sub-grid stress term, and obtaining three-dimensional weather element prediction data of the target area in a space below a preset height by simulation. The method can provide multi-element and high-resolution three-dimensional weather prediction products.

Description

Technical Field

[0001] This application belongs to the field of meteorological forecasting, and specifically relates to a small-scale meteorological numerical simulation method based on large eddy simulation. Background Technology

[0002] With the rapid development of the low-altitude economy, unmanned aerial vehicles (UAVs) are increasingly widely used in logistics, agriculture, surveying and mapping, and emergency rescue. Meteorological conditions are a crucial factor affecting the safety of UAVs flying at low altitudes, an important basis for the development and utilization of suitable airspace and temporal resources, and a key factor to consider in low-altitude flight route planning and real-time support. Within the UAV flight path area, the interaction between the complex urban underlying surface and the local micro-meteorological environment poses a significant challenge to real-time, high-resolution meteorological forecasting. Currently, forecasting and early warning technologies for low-altitude meteorological elements are still immature both domestically and internationally, especially the refined meteorological services specifically for UAV flights, which are still in the exploratory stage.

[0003] Currently used numerical weather prediction models are typically built upon mesoscale atmospheric models, with grid resolutions usually greater than 4 km. All turbulent motions are considered subgrid processes and can only be approximated using Reynolds-averaged parameterization methods. To improve the precision of numerical weather prediction and services, increasing the model grid resolution as much as possible, within the limits of computational capabilities, is an important development direction. However, since most mesoscale numerical models' subgrid schemes are based on coarse-resolution grids, their theoretical assumptions and framework designs are not applicable to finer grid scales. Therefore, simply increasing the model resolution does not necessarily improve the simulation accuracy of mesoscale models; it may even lead to the anomalous result of increased resolution but decreased accuracy.

[0004] Existing mesoscale numerical simulation techniques have the following drawbacks: First, the grid resolution cannot resolve small-scale core meteorological processes. Traditional mesoscale models typically have a grid resolution greater than 4 km, making it difficult to capture small-scale turbulent vortex structures below 100 meters that are crucial for low-altitude flight safety.

[0005] Second, to compensate for insufficient resolution, mesoscale numerical simulation systems rely on numerous parameterization schemes to characterize unanalyzed small-scale processes, such as turbulent mixing, cloud microphysics, and surface-atmosphere energy exchange. These parameterization schemes are mostly derived from fitting observational data under ideal conditions, exhibiting strong empirical limitations and failing to adapt to the variability of small-scale meteorological processes under complex underlying surface conditions.

[0006] Third, the errors in the parameterization scheme have a cumulative effect, which will continue to amplify as the simulation time increases and the spatial range expands, further reducing the reliability of small-scale meteorological simulations.

[0007] Therefore, how to provide a method that can achieve small-scale, high-resolution meteorological numerical simulation and accurately capture the spatial heterogeneity of small-scale meteorological fields under complex underlying surface conditions is a technical problem that urgently needs to be solved by those skilled in the art. Summary of the Invention

[0008] To address the aforementioned issues, this application provides a small-scale meteorological numerical simulation method based on large eddy simulation, which can provide multi-element, high-resolution three-dimensional meteorological forecast products.

[0009] A small-scale meteorological numerical simulation method based on large eddy simulation includes: Acquire driving field data for the target area, wherein the driving field data is at least one of global weather forecast data and the previous simulation forecast data; The driving field data is preprocessed to obtain initial meteorological field data for a three-dimensional grid system that can be interpolated to a preset resolution; Vortex separation is performed on the initial meteorological field data based on a spatial filtering function. Vortexes with a scale larger than the filtering scale threshold are classified as large-scale vortices that can be resolved by the grid, while vortices with a scale smaller than or equal to the filtering scale threshold are classified as small-scale vortices that cannot be resolved by the grid. The filtering scale threshold is determined by the grid resolution of the three-dimensional grid system. The large-scale vortex is solved explicitly using the filtered Navier-Stokes equations, and the small-scale vortex is solved by introducing a subgrid stress term to simulate and obtain three-dimensional meteorological element forecast data for the target area in the airspace below the preset height.

[0010] Compared with the prior art, this application has the following advantages: (1) The large eddy simulation (LES) method is used to replace the parameterization scheme of the traditional mesoscale simulation. The large-scale eddies and small-scale eddies in the three-dimensional atmosphere are separated by the spatial filtering function. The large-scale eddies are solved explicitly by the filtered Navier-Stokes equations, and the small-scale eddies are solved by the subgrid stress terms. The grid spatial resolution is improved to the level of several meters to tens of meters, realizing the high-resolution numerical weather prediction function for the airspace below 1000 meters altitude, and significantly improving the precision of the weather forecast. (2) The large eddy simulation method is used to directly analyze large-scale turbulent motion, which reduces the intervention of empirical parameterization schemes and avoids the cumulative effect of parameterization scheme errors being amplified as the simulation time and spatial range increase in traditional mesoscale simulation, thus improving the reliability and stability of small-scale meteorological simulation. (3) To meet the meteorological forecasting needs of low-altitude airspace below 1000 meters, high-resolution meteorological element forecast products such as wind field, temperature field, humidity field, and air pressure field are provided. These products can be widely used in low-altitude economic fields such as UAV logistics, agriculture, surveying and mapping, and emergency rescue. They provide scientific basis for the development and utilization of suitable airspace and time-domain resources, low-altitude flight route planning, and real-time support, and have significant social and economic benefits.

[0011] Other features and advantages of this application will be set forth in the description which follows, and will be apparent in part from the description, or may be learned by practicing the application. The objectives and other advantages of this application may be realized and obtained by means of the structures pointed out in the description, claims and drawings. Attached Figure Description

[0012] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are some embodiments of this application. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0013] Figure 1 A flowchart of a small-scale meteorological numerical simulation method based on large eddy simulation is shown. Figure 2 The flowchart of the driving field data preprocessing is shown; Figure 3 A schematic diagram of a multi-layered nested mesh is shown; Figure 4 A daily variation graph showing the average deviation between the 2-meter temperature forecast and the actual observation at the station; Figure 5 A graph showing the daily variation of the average deviation between the 2-meter relative humidity forecast and the actual observation at the station; Figure 6 A daily variation diagram of the average deviation between the predicted U-component of 10m wind speed and actual station observations; Figure 7 The graph shows the daily variation of the average deviation between the predicted V component of 10m wind speed and the actual observation at the station. Detailed Implementation

[0014] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0015] like Figure 1 As shown in the embodiments of this application, a small-scale meteorological numerical simulation method based on large eddy simulation is provided, including: S1. Obtain driving field data for the target area, wherein the driving field data is at least one of global weather forecast data and the previous simulation forecast data; S2. Preprocess the driving field data to obtain initial meteorological field data of a three-dimensional grid system that can be interpolated to a preset resolution; S3. Based on the spatial filtering function, vortex separation is performed on the initial meteorological field data. Vortices with a scale greater than the filtering scale threshold are regarded as large-scale vortices that can be resolved by the grid, and vortices with a scale less than or equal to the filtering scale threshold are regarded as small-scale vortices that cannot be resolved by the grid. The filtering scale threshold is determined by the grid resolution of the three-dimensional grid system. S4. The large-scale vortex is solved explicitly using the filtered Navier-Stokes equations, and the small-scale vortex is solved by introducing a subgrid stress term to simulate and obtain three-dimensional meteorological element forecast data of the target area in the airspace below the preset height.

[0016] The launch of a numerical weather prediction model requires driving field data, namely the model's initial field and boundary conditions. An automatic collection and management system for driving field data monitors the real-time update status of the driving data, collects and stores the necessary driving data, and provides data support for the launch of the numerical weather prediction model.

[0017] In this embodiment, the driving field data mainly comes from large-scale global model forecast fields and the model's previous forecast field, specifically including: (1) Global Weather Forecasting System Data The forecast field is collected from the Global Forecast System (GFS), for example, with a horizontal resolution of 0.5°×0.5° (the basic horizontal resolution between grid points is about 28 kilometers), a temporal resolution of 3 hours, an update frequency of four times a day, and a forecast lead time of 16 days (384 hours). Obtain the GFS dataset, which contains atmospheric and land-soil variables including, but not limited to, atmospheric variables (temperature, wind field, precipitation, atmospheric ozone concentration, etc.) and land-soil variables (soil moisture, soil temperature, etc.), with the data type GRIB2. The GFS model is a coupled model consisting of four separate models (atmospheric model, ocean model, land / soil model, and sea ice model), which work together to provide accurate weather condition forecasts.

[0018] (2) The forecast field at one time in the model When no analysis field is available at the current moment, the forecast field from the previous timeframe of the model can be used as the initial field for the current forecast. This initial field includes the following meteorological elements: model layer temperature, wind field, water vapor mixing ratio, cloud-water mixing ratio, rainwater mixing ratio, soil temperature, and soil moisture, with data type NetCDF. The horizontal resolution of the forecast field used for the hot start is consistent with the resolution set in the current model, ensuring seamless data transfer between different forecast batches.

[0019] Furthermore, Shell scripts are used to automatically download and categorize the data from the driving field, ensuring the timeliness and reliability of data acquisition. The specific steps are as follows: 1) Set system environment variables, including: driver farm storage directory (used to specify the local storage path of driver farm data), download point link (FTP address of data sources such as GFS), download software directory (installation path of download tools such as wget and curl), and common command path (path of Linux system commands). 2) Read the current system time and split the time variable, breaking down "year, month, day, hour, minute" into independent year, month, day, and hour / minute time variables; The split time variables are used to construct the URL path of the data file, generate the directory structure of local storage, and format the naming of the downloaded files.

[0020] 3) Create a drive field storage subdirectory named after the time variable. The directory structure is usually organized in the hierarchy of "year / month / day" or "year-month-day". Before creation, the script checks if the directory already exists. If the directory already exists, it further checks if the corresponding driver data has been downloaded and decides whether it needs to be downloaded again based on the check results (e.g., if the data file is corrupted or missing). If the directory does not exist, it automatically creates a new directory.

[0021] 4) Convert the split time variables into the standard naming format for GFS data; 5) Construct a complete download URL according to the standard naming format, start the download software (such as wget) to download the driver data to the corresponding storage subdirectory, and output the download information to the download log file.

[0022] The log content includes the download start time, the name and URL of the downloaded file, the download completion status, and error messages. The logging function facilitates subsequent review of download records and troubleshooting of download failures and other anomalies.

[0023] In this embodiment, the driving field data collection follows a short-term forecast initiation method. Specifically, when the system is first started (cold start mode) or needs re-initialization, the GFS global forecast field is used as the initial field. Cold starts are typically performed every 6 or 12 hours to obtain an initial state independent of previous forecasts. In contrast, in the continuous rolling forecast mode (hot start mode), the forecast field from the previous time period is used as the initial value for the current forecast. Hot starts can be performed hourly, maintaining forecast continuity and computational efficiency. The combined use of these two initiation methods ensures both the accuracy of the initial forecast field (by periodically introducing global observation assimilation data through cold starts) and the real-time nature of high-frequency rolling forecasts (by updating hourly through hot starts).

[0024] Data preprocessing is a crucial preliminary step in numerical weather prediction systems. Its core task is to process raw driving field data and static geographic data into initial field files recognizable by numerical models. This application's embodiment employs the WRF Preprocessing System (WPS) to implement data preprocessing functions. WPS consists of three core programs: geogrid is used to determine the scope and projection method of the simulated area, and to interpolate static geographic data to grid points; ungrib is used to decompress and reanalyze three-dimensional meteorological elements in data. Metgrid is used to horizontally interpolate decompressed meteorological elements onto static data processed by geogrid.

[0025] After the above preprocessing steps, the met_em file is finally obtained, which can be directly used for WRF numerical mode input.

[0026] In this embodiment of the application, the driving field data is preprocessed to obtain initial meteorological field data of a three-dimensional grid system that can be interpolated to a preset resolution, such as... Figure 2 As shown, the preprocessing process includes the following steps: S21. Determine the range of the target simulation area, the map projection type, and the grid parameters of the three-dimensional grid system. The grid parameters include the number of grid points, the grid resolution, the number of grid nesting layers, and the nesting ratio between adjacent grid levels. S22. Interpolate static geographic data to each layer of the three-dimensional grid system. The static geographic data includes terrain height, land use type, soil type, and vegetation cover. The aforementioned static geographic data can reflect the real geographic and surface characteristics of the underlying surface of the target simulation area. Interpolating it to grid points can enable the simulation grid to have real terrain and surface attributes, thereby more accurately characterizing the dynamic and thermal influence of the underlying surface on the near-surface meteorological field and improving the consistency between the initial field and the actual environment.

[0027] S23. Decode the driving field data, extract meteorological elements, including temperature, wind field, humidity, air pressure, soil temperature and soil moisture, and interpolate the meteorological elements horizontally and vertically onto the static geographic data of the corresponding grid points to obtain the initial meteorological field data.

[0028] Horizontal interpolation is used to map coarse-resolution driving field data to the horizontal spatial location of a high-resolution 3D grid, while vertical interpolation is used to assign meteorological elements to grid nodes at different altitudes. Through bidirectional horizontal and vertical interpolation, the initial meteorological field is made complete, continuous, and physically consistent on the spatial grid, providing reliable initial conditions for subsequent eddy separation and numerical solutions based on large eddy simulation.

[0029] The data preprocessing process is described in detail below: I. Configuring a 3D Mesh System The three-dimensional mesh system is a multi-layer nested mesh structure, which includes a root domain mesh and at least one level of sub-domain mesh arranged sequentially from the outside to the inside. Each level of mesh works together. The resolution of the root domain mesh is at the kilometer level, and the resolution of the sub-domain mesh is progressively refined to the meter level or the hundred-meter level. The nesting ratio of adjacent levels of mesh is an odd number, and the starting coordinates of the sub-domain mesh and the starting coordinates of the root domain mesh satisfy an integer multiple matching relationship.

[0030] Specifically, the root domain grid, serving as the outermost grid, has a resolution set to the kilometer level, covering the entire target simulation area to capture large-scale meteorological field evolution characteristics. Based on the root domain grid, one or more sub-domain grids are set according to forecast requirements. The grid resolution of the sub-domain grids is progressively finer with each nested level, eventually reaching the hundred-meter or even meter level, thus achieving fine analysis of small-scale meteorological processes in key areas. To ensure the stability of numerical calculations and the accuracy of interpolation, the following constraints must be met between adjacent grid levels: (1) The nesting ratio is odd. In numerical models, the variable values ​​of subdomain grids need to be interpolated through the root domain grid. Therefore, the nesting ratio of adjacent grid levels (i.e., the ratio of the root domain grid resolution to the subdomain grid resolution) must be set to an odd number to ensure that the center point of the subdomain grid coincides with a point of the root domain grid, thereby avoiding phase shift during the interpolation process and ensuring the consistency of the calculation. Commonly used nesting ratios are 1:3 and 1:5.

[0031] like Figure 3 As shown, the root domain relationship `parent_id` is used to specify the parent domain (supervisor domain) of each nested region. For example, `parent_id=1, 1, 2, 3` respectively represent: Region 1 (d01) is the root region, and its parent region is itself; Region 2 (d02) is nested within Region 1 (d01), and its parent region is Region 1; Region 3 (d03) is nested within Region 2 (d02), with Region 2 being the parent region; Region 4 (d04) is nested within Region 3 (d03), with Region 3 being the parent region.

[0032] The nesting ratio `parent_grid_ratio` determines the mesh spacing refinement ratio of the subgrid grid relative to its root grid. For example, `parent_grid_ratio=1, 3, 3, 3` represent: The ratio of the root domain (d01) is 1 (itself); The mesh refinement ratio of region 2 (d02) relative to the root region (d01) is 3; The mesh refinement ratio of region 3 (d03) relative to region 2 (d02) is 3; The mesh refinement ratio of region 4 (d04) to region 3 (d03) is 3.

[0033] If the root domain grid size dx = 3000 meters, then the resolution of d02 is 1000 meters (3000 ÷ 3), the resolution of d03 is approximately 333 meters (1000 ÷ 3, rounded down), and the resolution of d04 is approximately 111 meters (333 ÷ 3), thus achieving progressive encryption from kilometer level to hundred-meter level.

[0034] When using real data for simulation, the nesting ratio `parent_grid_ratio` must be set to an odd number to ensure that the center point of the subgrid coincides with a point of the root grid, thereby avoiding phase shift during interpolation and ensuring consistency of calculations.

[0035] (2) Matching relationship of integer multiples of starting coordinates Assuming (i_parent_start, j_parent_start) are the starting coordinates of the lower left corner of the subdomain grid within its root domain grid, when i_parent_start = 1, 65, 65, 65 and j_parent_start = 1, 42, 42, 42, respectively, they represent: The starting position of d01 is (1, 1); The starting position of d02 in d01 is (65, 42); The starting position of d03 in d02 is (65, 42); The starting position of d04 in d03 is (65, 42).

[0036] The parameters e_we and e_sn specify the number of grid points in each region in the east-west and north-south directions, respectively. The starting grid point value in the north-south direction (s_sn) and the east-west direction (s_we) must be set to 1, and the ending grid point value (e_sn and e_we) determines the nesting size.

[0037] The parameters e_we and e_sn must satisfy the following constraints: (e_we-1) / (parent_grid_ratio) and (e_sn-1) / (parent_grid_ratio) are both integers. The significance of the above constraints is that the upper right boundary of the subdomain grid can coincide with a grid point in the parent domain grid, ensuring complete spatial matching of nested grids and avoiding grid misalignment or boundary discontinuity.

[0038] By rationally dividing the hierarchical resolution, standardizing the nesting ratio and coordinate matching relationship to configure a multi-layer nested grid structure, the three-dimensional grid system constructed in this application can realize multi-scale meteorological simulation from kilometer level to meter level or hundred-meter level while ensuring computational efficiency. The outer root domain grid covers the entire target simulation area and captures the evolution characteristics of large-scale meteorological fields. The inner subdomain grid is gradually densified to perform fine analysis of small-scale meteorological processes in key areas. The coordinated cooperation of each level of grid provides accurate spatial support for high-resolution meteorological forecasting in low-altitude airspace (below 1000 meters).

[0039] II. Determining the Map Projection Method Since numerical simulations require data conversion between the three-dimensional Earth surface and a two-dimensional planar grid, a suitable map projection method must be selected to map the spatial information of the Earth's surface onto the planar grid of the simulation area. The essence of map projection is to establish a mathematical correspondence between points on the Earth's surface and points on the planar grid. In step S21, the map projection type is selected from one of the following: Lambert projection, polar projection, Mercator projection, or lat-lon projection.

[0040] The four projection methods described above produce maps with different shapes and are suitable for different geographical latitudes and simulated areas. Their specific characteristics are as follows: Lambert projection: a conformal conical projection that exhibits less distortion in mid-latitude regions and is suitable for simulating large areas of mid-latitudes; Polar projection: a type of conformal azimuth projection that exhibits less distortion in high-latitude regions and is suitable for simulating polar or high-latitude areas; Mercator projection: a conformal cylindrical projection that exhibits minimal deformation near the equator and is suitable for simulating low-latitude regions; Latitude and longitude projection (lat-lon): directly uses latitude and longitude as grid coordinates, suitable for simulations on a global scale.

[0041] The location and extent of the total simulation region on the planar grid are determined by the following parameters: reference longitude ref_lon, reference latitude ref_lat, number of grid points in the east-west direction e_we, number of grid points in the north-south direction e_sn, grid size dx (meters) in the east-west direction of the root region, and grid size dy (meters) in the north-south direction of the root region. Among them, the reference longitude and latitude (ref_lat and ref_lon) determine the center position or reference point of the simulation region, and the number of grid points and grid size (e_we, e_sn, dx, dy) determine the actual spatial coverage of the simulation region.

[0042] For example, when using the Lambert projection, the projection parameters that need to be configured include the reference latitude ref_lat, the reference longitude ref_lon, the first true latitude truelat1, the second true latitude truelat2, and the standard longitude stand_lon.

[0043] The ref_lat and ref_lon parameters specify the location of a reference point on the Earth's surface for the simulation region. This reference point is usually located near the center of the simulation region and is used to determine the reference position for the projection. The numerical model uses this reference point as a reference to convert the latitude and longitude coordinates of other points on the Earth's surface into coordinates of planar grid points.

[0044] `truelat1` and `truelat2` are two standard parallels in the Lambert projection. Along these two parallels, the projected map remains undistorted, meaning distances on the map maintain a precise proportional relationship to actual distances on the Earth's surface. Distortion is minimal between the two standard parallels, gradually increasing outwards. When `truelat1` equals `truelat2`, it's a standard Lambert projection with only one standard parallel. When `truelat1` ≠ `truelat2`, it's a double-standard-parallel Lambert projection, where two standard parallels better reduce projection distortion in mid-latitude zonal regions. For example, `truelat1` = 30 and `truelat2` = 60 cover most of the mid-latitude region, effectively controlling projection distortion within the simulated area.

[0045] `stand_lon` specifies the meridian that, after the Lambert projection, is parallel (i.e., perpendicular) to the y-axis of the planar grid coordinate system. This meridian maintains its north-south orientation after projection and does not undergo rotation.

[0046] stand_lon is used to determine the orientation of the grid after projection, so that the grid coordinate system of the simulated area maintains a reasonable correspondence with the actual geographic orientation. The stand_lon is configured to be 116.688, which is consistent with the reference longitude ref_lon, to ensure that the north direction of the grid after projection is basically consistent with the geographic north direction.

[0047] For other projection types, polar projection only requires defining truelat1 and stand_lon, lat-lon projection requires defining pole_lat, pole_lon, and stand_lon, and mercator projection requires defining truelat1. This parameter determines the deformation control near the equator.

[0048] By selecting the map projection type and configuring the projection parameters appropriately, the optimal projection method is chosen based on the geographical location and shape characteristics of the simulated target area, minimizing the impact of projection distortion on simulation accuracy. Simultaneously, the reasonable projection parameter settings ensure consistency between the spatial extent and grid orientation of the simulated area and the actual geographical location, providing an accurate spatial reference for subsequent meteorological element interpolation and numerical simulation.

[0049] III. Configuring Dynamic Simulation Time The start and end times of the numerical simulation are set to determine the forecast time window. For forecast systems that start frequently every hour, the start and end times of the simulation need to be dynamically adjusted continuously. In this embodiment, the start and end times of the simulation are set through the following parameters in the namelist.wps file: `start_date` is the simulated start time, in the format "YYYY-MM-DD_HH:00:00" end_date, the simulated end time, formatted the same as above. interval_seconds, the time interval between meteorological data files It should be noted that the three-dimensional mesh system constructed in this application is a multi-layer nested mesh structure. start_date and end_date need to specify the start and end times for each nested region, that is, each region corresponds to a time value.

[0050] This application also employs hourly rolling forecasts to meet the real-time requirements of low-altitude airspace for weather forecasts. For forecast systems that start frequently every hour, the start and end times of the simulation need to be constantly adjusted; this can be achieved using shell scripts to dynamically update the time parameters.

[0051] The specific mechanism is as follows: each time the forecast is started, the Shell script reads the current system time; uses the current time as the new start time (start_date); calculates the end time (end_date) based on the forecast duration (e.g., the next 24 hours); automatically updates the time parameters in the namelist.wps file; and calls the preprocessor to perform data processing.

[0052] IV. Processing Static Geographic Data The geogrid program is responsible for interpolating static geographic data onto the grid points of each layer of a 3D grid system. The static geographic data includes, but is not limited to: Terrain data, such as digital elevation models (DEMs); Land use data, such as land use types. Soil data, such as soil type and soil texture; Vegetation data, such as vegetation cover (e.g., Normalized Difference Vegetation Index NDVI) and Leaf Area Index (LAI). Other surface parameters, such as surface roughness, albedo, and surface emissivity.

[0053] The main process of geogrid performing geographic data interpolation is as follows: 1. Read the simulation region parameters (including projection method, number of grid points, grid size, etc.) defined in the namelist.wps file; 2. Calculate the geographical location (latitude and longitude) of each grid point based on the regional parameters; 3. Read the static geographic data corresponding to each grid point from the geographic data source file; 4. Use appropriate interpolation methods (such as bilinear interpolation, nearest neighbor interpolation, etc.) to map geographic data onto grid points; 5. Output a file containing static geographic data for use by metgrid.

[0054] V. Extracting meteorological elements The driving field data used when starting a numerical weather prediction model typically employs a special storage format (such as GRIB) and contains various meteorological variables, making it unsuitable for directly driving the model. Therefore, the ungrib program is used to decode the raw driving field data and extract the meteorological elements required by the model, converting it into an intermediate format recognizable by the model. The driving field data decoding process mainly includes: selecting a suitable Vtable encoding table, running the ungrib program to decode the GRIB file, and outputting the intermediate format meteorological data file.

[0055] Because different data sources (such as NCEP's GFS data and ECMWF data) use their own different encoding tables for GRIB format data, a suitable Vtable file needs to be selected before driving field extraction. The Vtable file is an encoding lookup table used to guide the decoding program in identifying and extracting specific meteorological elements from the GRIB file.

[0056] Here is an example of a Vtable file:

[0057] The meanings of each field are as follows:

[0058] For data released by ECMWF, upper-air and surface feature fields are usually stored separately. Therefore, during decoding and extraction, the upper-air and surface features need to be processed separately. That is, the ungrib program is executed once for each of the upper-air and surface feature field files, and the data in the two files are decoded and then merged to form a complete meteorological feature field.

[0059] After the driving field data is decoded, the extracted meteorological elements are written into an intermediate format data file. In the namelist.wps file, the out_format parameter is used to specify the format of the output intermediate data, which is set to WPS. The prefix parameter is used to describe the path and prefix of the intermediate file. For example, when prefix is ​​set to FILE, the created intermediate file will be named in the format "FILE:YYYY-MM-DD_HH", indicating the valid time of the data, such as FILE:2025-05-03_00, FILE:2025-05-03_03, etc.

[0060] The aforementioned driving field decoding and meteorological element extraction methods offer the following advantages: They are compatible with multiple data sources; the Vtable mechanism provides flexible data mapping capabilities, supporting access to various data sources such as GFS and ECMWF, demonstrating strong versatility and scalability; they uniformly convert raw driving field data from different sources and in different formats into an intermediate format recognizable by numerical models, ensuring data input compatibility and consistency; through precise mapping using the Vtable lookup table, they accurately extract specific meteorological elements required by the model from the original GRIB file containing hundreds or thousands of variables, avoiding data redundancy; and they support the complete extraction of vertically stratified data such as multi-layered soil temperature and soil moisture, providing comprehensive data support for subsequent vertical interpolation.

[0061] VI. Interpolated meteorological elements (1) Horizontal interpolation The metgrid program horizontally interpolates the decompressed meteorological elements onto the geogrid-processed static data. The key control file for horizontal interpolation is METGRID.TBL, which specifies the stratification intervals and corresponding interpolation methods for each meteorological element.

[0062] The purpose of the METGRID.TBL file is to define the vertical stratification intervals for each meteorological element, specify the interpolation method (such as linear interpolation, logarithmic interpolation, etc.) for each stratification interval, and ensure that meteorological elements from different sources can be correctly interpolated onto the model grid.

[0063] Taking soil temperature interpolation configuration as an example, the METGRID.TBL file is as follows: name=ST z_dim_name=num_st_layers derived=yes IF fill_lev=10 : ST000010(200100) fill_lev=40: ST010040(200100) fill_lev=100: ST040100(200100) fill_lev=200 : ST100200(200100) In the above configuration, soil temperature (ST) is vertically stratified into 4 layers (0-10cm, 10-40cm, 40-100cm, 100-200cm), with each layer provided by the corresponding source data variable.

[0064] (2) Vertical interpolation Vertical interpolation is performed after horizontal interpolation. Its control parameters are set in the domains and dynamics name list records of the namelist.input file. The main parameters include: Model top pressure, which is used to define the upper boundary pressure value of the model; Vertical layer number (e_vert), which defines the number of layers in the pattern in the vertical direction; Vertical coordinate type, which is used to select the vertical coordinate system.

[0065] This application employs a hybrid vertical coordinate (HVC) system for vertical interpolation. The HVC coordinate system is characterized by: using terrain-following coordinates (TF) near the ground to ensure the model's bottom-layer grid accurately reflects actual terrain undulations; using isobaric coordinates in higher regions, meaning the coordinate plane remains horizontal above the user-defined pressure layer and does not follow terrain variations; and a smooth transition between the two coordinate systems in transitional regions. The main function of the hybrid vertical coordinate system is to reduce the artificial influence of terrain on the top of the model and avoid computational instability in steep terrain areas.

[0066] Formula for calculating dry pressure in mixed vertical coordinates: PDRY(i,j,k)=B(k)(PDRY SFC (i,j)–PTOP)+(λ(k)–B(k))(P0–PTOP)+PTOP Where PDRY(i,j,k) is the dry air pressure at grid point (i,j,k), PDRY SFC (i,j) is the surface dry air pressure; B(k) is a 1D weighted array used to control the transition of terrain following coordinates to isobaric coordinates; PTOP stands for mode top pressure. λ(k) represents the original eta coordinates; P0 is the reference pressure (usually 1000 hPa). The function of B(k) is as follows: when B(k) ≡ λ(k), the formula simplifies to pure terrain-following coordinates; when B(k) ≡ 0, the formula simplifies to pure isobaric coordinates; when B(k) takes an intermediate value, it achieves a smooth transition from terrain-following coordinates near the ground to isobaric coordinates at high altitudes.

[0067] The critical value for B(k) transitioning to isobaric layers determines how many eta layers (counting downwards from the top of the model layer) remain isobaric surfaces. The adjustable range of the critical value is 0 to 1. Increasing the critical value results in more eta layers being affected by the HVC, enhancing the "flattening" effect of the coordinate surface. However, if the critical value is too large, in higher terrain areas, the vertical eta layers may be excessively compressed, affecting the uniformity of vertical resolution. In practical applications, the critical value needs to be set appropriately based on the topographic relief characteristics of the simulation area.

[0068] VII. Lateral Boundary Constraints Since regional models only cover a limited geographical area, atmospheric motion will constantly cross the regional boundaries. Therefore, when generating or updating side boundaries for model forecasts, the side boundary conditions should be as close as possible to the actual atmospheric conditions, and should also be consistent in forecasts and simulations.

[0069] In this embodiment of the application, a lateral boundary constraint is provided for the three-dimensional grid system based on global weather forecast data, and a lateral boundary file containing the current value of meteorological elements and the trend term to the next time moment is generated, wherein the trend term is the rate of change of meteorological elements between two adjacent boundary time points.

[0070] Each lateral boundary field is defined along the north, south, east, and west sides of a rectangular grid, and includes boundary values ​​and trend terms for the following elements: east-west wind speed components. North-South Wind Speed ​​Components Vertical wind speed component w Temperature Water vapor mixing ratio Disturbance potential Mass of disturbed dry air column These variable prediction models constrain the lateral boundaries.

[0071] Each variable in the lateral boundary file has a valid value at the initial time of the lateral boundary, as well as a trend term for the next boundary time period. The time dimension of the lateral boundary file is set as follows: assuming large-scale data is available every 3 hours and the forecast lead time is 12 hours, the lateral boundary file will contain data at the following time points: the initial values ​​of each variable at hour 0, and the trend terms of each variable at hours 3, 6, 9, and 12. Since the trend term for the next time period does not need to be calculated for the last time point, the time period (number of time points) of the boundary file is one less than the time period after interpolation of the large-scale data.

[0072] For example, for the east-west wind speed component The first time period in the horizontal boundary file contains the 0 hour point. Data U 0hThe calculation formula is as follows:

[0073] U 0h This represents the relevant quantities used to constrain the horizontal wind speed U-prediction results of the regional forecast model in the lateral boundary conditions at the initial time (0-hour point). Representative x Mass of the dry air column after directional averaging u Represents horizontal wind speed u Quantity, Representative x Map projection factor after directional averaging.

[0074] Trend value is defined as:

[0075] This trend value represents the value that allows the grid points to be advanced from their initial values ​​to a 3-hour simulation time.

[0076] The coupling relationship is: horizontal momentum field ( u , v )and Coupled with inverse map factors; other three-dimensional fields ( , and Only with coupling.

[0077] In this embodiment, the boundary values ​​and tendencies of vertical velocity and non-vapor moisture species are included in the outer lateral boundary file, but serve as placeholders for nested boundary data within the fine mesh. The width of the lateral boundary regions along the four sides is a user-selectable parameter used to control the thickness of the boundary transition regions.

[0078] In the Large Eddy Simulation (LES) method, the basic principle of spatial filtering is that atmospheric turbulent motion contains vortex structures of different scales. Large-scale vortices contain most of the turbulent kinetic energy, while small-scale vortices mainly dissipate it. Using a spatial filtering function, large-scale and small-scale vortices in the three-dimensional atmosphere can be separated. Large-scale vortices are solved explicitly using the turbulence equations, while small-scale vortices are solved by adding stress terms to the turbulence equations for closed-loop solutions.

[0079] In this embodiment, a spatial filtering function is used to perform spatial filtering on the instantaneous physical quantities of the initial meteorological field. Vortices with a scale larger than the filtering scale threshold are regarded as large-scale vortices that can be resolved by the grid, while vortices with a scale smaller than or equal to the filtering scale threshold are regarded as small-scale vortices that cannot be resolved by the grid.

[0080] Specifically, the instantaneous physical quantities of the initial meteorological field are obtained by using a spatial filtering function. After spatial filtering, the large-scale physical quantity components that can be resolved by the mesh are obtained. for:

[0081] Among them, the instantaneous physical quantities of the initial meteorological field It can be any one of wind speed, temperature, humidity, or air pressure; D Indicates a region of atmospheric flow; Represents the filtered grid space coordinates. These are the original spatial coordinates of the actual flow field; This represents the filtering kernel function.

[0082] Filter kernel function The form and scale parameters determine the cutoff scale of the filtering operation, i.e., the filtering scale threshold, which is the dividing line between large-scale vortices and small-scale vortices: vortices with a scale larger than the filtering scale threshold are retained in the filtered field. In the filter, large-scale vortices are analyzed by the mesh; vortices with a scale smaller than or equal to the filtering scale threshold are filtered out, and their effects are parameterized through subgrid stress terms. Commonly used filtering kernel functions include: box filtering, Gaussian filtering, and spectral truncation filtering. The appropriate type of filtering kernel function should be selected based on the actual simulation requirements and computational resources.

[0083] In this embodiment, the filtering scale threshold is adaptively set according to the mesh resolution of the 3D mesh system. Specifically, the filtering scale threshold is positively correlated with the mesh resolution; the higher the mesh resolution (the smaller the mesh spacing), the smaller the filtering scale threshold, enabling the resolution of smaller-scale vortices; conversely, the lower the mesh resolution (the larger the mesh spacing), the larger the filtering scale threshold, resulting in the filtering out of more relatively small-scale vortices, which are then processed by the sub-mesh model. The filtering scale is automatically adjusted according to different simulation requirements (such as nested meshes of different resolutions) to ensure the applicability and computational efficiency of the large eddy simulation method at different resolutions.

[0084] In this embodiment of the application, the filtered Navier-Stokes equations used for explicit numerical solutions of the large-scale vortices and closed-loop solutions of the small-scale vortices are as follows: Variables marked with a superscript "~" are the solvable scale components obtained through the filtering function, where, , They represent along , The velocity components of the solvable scale in the direction, subscript iWhen j=1,2,3, they represent the flow direction, spanwise direction, and vertical direction, respectively.

[0085] This represents the pressure at a solvable scale.

[0086] This indicates air density.

[0087] This represents the subgrid stress term, used to characterize the momentum dissipation effect of small-scale vortices on large-scale flow fields.

[0088] To close the above equations, the subgrid stress terms need to be... Perform parameterized characterization. Use... Smagorinsky The model is parameterized. Smagorinsky The mathematical expression of the model is as follows:

[0089] in, For the Kronecker function, For subgrid stress isotropic components, is the subgrid eddy viscosity coefficient.

[0090] The filtered velocity-strain rate tensor This describes the deformation rate of fluid micro-elements.

[0091] The basic assumption of the Smagorinsky model is that the deviatoric stress component of the subgrid stress ( ) With the solvable scale velocity strain rate tensor A linear relationship exists, with a proportionality constant of 1. .

[0092] In large eddy simulation, subgrid stress terms The role of subgrid stress is twofold: small-scale vortices extract energy from large-scale flows and dissipate this energy as heat (in some cases, subgrid stress can also transfer energy from small-scale to large-scale, i.e., energy backflow), effectively suppressing numerical oscillations and improving the stability of numerical simulations. Secondly, through the Smagorinsky model, the subgrid stress term... It is represented as a function of the solvable scale strain rate tensor, thus making the filtered Navier-Stokes equations a closed system of equations that can be solved numerically.

[0093] Based on this, the embodiments of this application use the filtered Navier-Stokes equations for explicit solution of large-scale vortices, which fully preserves the dynamic characteristics of the main energy carriers (large-scale vortices) in atmospheric turbulence, and achieves accurate characterization of large-scale vortices in atmospheric turbulence; for small-scale vortices, a simple Smagorinsky eddy viscosity model is used for parameterization, which simplifies the modeling of small-scale vortices. The Smagorinsky model is based on the turbulence dissipation mechanism, has a clear physical picture, low parameter sensitivity, and improves the reliability of the simulation.

[0094] In this embodiment of the application, three-dimensional meteorological element forecast data of the target area in the airspace below a preset altitude is obtained by large eddy simulation. Preferably, the airspace below the preset altitude is the low-altitude airspace below 1000 meters, and its three-dimensional meteorological element forecast data includes wind field, temperature field, humidity field and pressure field data.

[0095] To verify the effectiveness of the technical solution in this application, the forecast results of the large eddy simulation were compared and verified with the actual observations at the site. The verification results are presented below from four aspects: 2-meter temperature, 2-meter relative humidity, and 10-meter wind speed (U component and V component).

[0096] Figure 4 This is a graph showing the daily variation of the average deviation between the 2-meter temperature forecast and the actual station observations. The horizontal axis represents time (UTC, from 0:00 to 24:00), and the vertical axis represents the average temperature deviation (°C). Figure 4 As can be seen, during the period from 0:00 to 11:00, the predicted temperature is relatively close to the observed temperature, with an average temperature deviation within the range of 0-1.31℃, indicating that this application has good forecast accuracy in the early stages of atmospheric boundary layer development. During the period from 12:00 to 18:00, the average temperature error gradually increases, reaching a maximum of -3.83℃ (the predicted value is too low). This period corresponds to the afternoon to evening, which is the period when atmospheric boundary layer turbulence is most intense. During the period from 19:00 to 24:00, the average temperature deviation decreases to some extent, indicating that as solar radiation weakens, the model's ability to simulate temperature changes recovers to some extent.

[0097] Figure 5 This is a graph showing the daily variation of the average deviation between the 2-meter relative humidity (RH) forecast and the actual observation at the station. The horizontal axis represents time (UTC, from 0:00 to 24:00), and the vertical axis represents the average deviation of relative humidity (unit: %). Figure 5 As can be seen, the relative humidity forecast results generally show a negative deviation trend, that is, the forecast humidity is lower than the observed humidity, and the average deviation of relative humidity fluctuates in the range of -10.59% to 4.75%.

[0098] Figure 6This is a daily variation graph showing the average deviation between the predicted U-component wind speed at 10m and the actual observations at the stations. The horizontal axis represents time (UTC, from 0:00 to 24:00), and the vertical axis represents the average deviation of the U-component wind speed (unit: m / s). Figure 6 It can be seen that within the forecast period of 0:00-07:00, the forecast effect of near-surface U-wind (east-west wind) is better, with an average error within 0.37 m / s, indicating that the simulation of horizontal wind field in this application has high accuracy from night to early morning. During the period of 8:00-12:00, the average error gradually increases, reaching a maximum of -1.52 m / s (the forecast U-wind is weak). This period corresponds to the increased atmospheric instability after sunrise, the intensified turbulent activity, and the increased complexity of wind field changes.

[0099] Figure 7 This is a daily variation graph showing the average deviation between the predicted V component of 10m wind speed and the actual observations at the stations. The vertical axis represents the average deviation of the V component (unit: m / s). Figure 7 It can be seen that the average forecast error of the 10-meter wind speed V component (north-south wind) is in the range of -0.37 to 1.86 m / s. Compared with the U component, the error fluctuation range of the V component is slightly larger.

[0100] The above verification results show that the 2-meter temperature forecast can simulate the diurnal variation trend well, with an average error within 1.31℃ during the night to early morning period, meeting the meteorological support requirements for low-altitude flights. The average error of the 10-meter wind speed U-component during the night to early morning period is within 0.37 m / s, demonstrating high forecast accuracy. Combined with spatial distribution characteristic analysis, this application can clearly simulate the spatial distribution characteristics of regional meteorological elements, such as the temperature and wind speed differences between the south and north, as well as the meteorological gradient caused by the north-south topography or underlying surface. This application can also simulate the temporal variation trend of meteorological elements, such as the temperature rise caused by increased solar radiation after sunrise and the near-surface wind speed increase caused by increased atmospheric instability. The synergistic analysis of the temperature and wind fields in this application can provide multi-element, high-resolution three-dimensional meteorological forecast products, providing comprehensive data support for meteorological support for low-altitude flights.

[0101] It should be noted that, for the sake of simplicity, the foregoing method embodiments are all described as a series of actions. However, those skilled in the art should understand that the present invention is not limited to the described order of actions, because according to the present invention, some steps can be performed in other orders or simultaneously. Furthermore, those skilled in the art should also understand that the embodiments described in the specification are preferred embodiments, and the actions and modules involved are not necessarily essential to the present invention.

[0102] In the above embodiments, the descriptions of each embodiment have their own emphasis. Parts not described in detail in a particular embodiment can be found in the relevant descriptions of other embodiments. 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 foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A small-scale meteorological numerical simulation method based on large eddy simulation, characterized in that, The method includes: Acquire driving field data for the target area, wherein the driving field data is at least one of global weather forecast data and the previous simulation forecast data; The driving field data is preprocessed to obtain initial meteorological field data for a three-dimensional grid system that can be interpolated to a preset resolution; Vortex separation is performed on the initial meteorological field data based on a spatial filtering function. Vortexes with a scale larger than the filtering scale threshold are classified as large-scale vortices that can be resolved by the grid, while vortices with a scale smaller than or equal to the filtering scale threshold are classified as small-scale vortices that cannot be resolved by the grid. The filtering scale threshold is determined by the grid resolution of the three-dimensional grid system. The large-scale vortex is solved explicitly using the filtered Navier-Stokes equations, and the small-scale vortex is solved by introducing a subgrid stress term to simulate and obtain three-dimensional meteorological element forecast data for the target area in the airspace below the preset height.

2. The method according to claim 1, characterized in that, Preprocessing the driving field data includes: Determine the extent of the target simulation area, the map projection type, and the grid parameters of the three-dimensional grid system, including the number of grid points, grid resolution, number of grid nesting layers, and nesting ratio; The static geographic data is interpolated to the grid points of each layer of the three-dimensional grid system. The static geographic data includes terrain height, land use type, soil type and vegetation cover. The driving field data is decoded, and meteorological elements are extracted, including temperature, wind field, humidity, air pressure, soil temperature and soil moisture. The meteorological elements are then horizontally and vertically interpolated onto the static geographic data of the corresponding grid points to obtain the initial meteorological field data.

3. The method according to claim 2, characterized in that, The three-dimensional mesh system is a multi-layer nested mesh structure, which includes a root domain mesh and at least one level subdomain mesh arranged sequentially from the outside to the inside. The root domain grid has a resolution of kilometers, the sub-domain grids have resolutions progressively reduced to meters or hundreds of meters, the nesting ratio of adjacent grid levels is odd, and the starting coordinates of the sub-domain grids and the starting coordinates of the root domain grids satisfy an integer multiple matching relationship.

4. The method according to claim 2, characterized in that, The map projection type is selected from one of the following: Lambert projection, Polar projection, Mercator projection, and Latt-Lon projection.

5. The method according to claim 2, characterized in that, When performing horizontal interpolation on the meteorological elements, a stratified interval and corresponding interpolation method are specified for each meteorological element. When performing vertical interpolation on the meteorological elements, a hybrid vertical coordinate system is used, and the hybrid vertical coordinate system adopts terrain-following coordinates near the ground and isobaric coordinates above the preset pressure layer.

6. The method according to claim 1, characterized in that, After preprocessing the driving field data, the process further includes: Based on global weather forecast data, the 3D grid system is provided with lateral boundary constraints, generating a lateral boundary file containing the current values ​​of meteorological elements and the trend terms to the next time step. The trend term is the rate of change of meteorological elements between two adjacent boundary time points; Each lateral boundary field is defined along the north, south, east, and west sides of a rectangular grid, and includes boundary values ​​and trend terms for the following elements: east-west wind speed components. North-South Wind Speed ​​Components Vertical wind speed component w Temperature Water vapor mixing ratio Disturbance potential Mass of disturbed dry air column .

7. The method according to claim 1, characterized in that, Using spatial filtering functions to analyze the instantaneous physical quantities of the initial meteorological field After spatial filtering, the large-scale physical quantity components that can be resolved by the mesh are obtained. for: Among them, the instantaneous physical quantities of the initial meteorological field It can be any one of wind speed, temperature, humidity, or air pressure; D Indicates a region of atmospheric flow; Represents the filtered grid space coordinates. These are the original spatial coordinates of the actual flow field; The filter kernel function determines the filter scale threshold, which is adaptively set according to the grid resolution of the three-dimensional mesh system, so that large-scale vortices larger than the filter scale threshold can be directly solved.

8. The method according to claim 1 or 7, characterized in that, The filtered Navier-Stokes equations used for explicit numerical solutions of the large-scale vortices and closed-loop solutions of the small-scale vortices are as follows: Variables marked with a superscript "~" are the solvable scale components obtained through the filtering function, where, , They represent along , The velocity components of the solvable scale in the direction, subscript i When j=1,2,3, they represent the flow direction, spanwise direction, and vertical direction, respectively. Represents the pressure at a solvable scale; Indicates air density; This represents the subgrid stress term, used to characterize the momentum dissipation effect of small-scale vortices on large-scale flow fields.

9. The method according to claim 8, characterized in that, The subgrid stress term use Smagorinsky Model parameterization: in, For the Kronecker function, For subgrid stress isotropic components, The subgrid eddy viscosity coefficient; The filtered velocity-strain rate tensor .

10. The method according to claim 1, characterized in that, The airspace below the preset altitude is the low-altitude airspace below 1000 meters, and the three-dimensional meteorological element forecast data includes wind field, temperature field, humidity field and pressure field data.

Images