A method and system for rapid simulation of microscale meteorological elements in complex terrain regions

Abstract

The present application relates to a method and system for rapidly simulating microscale meteorological elements in complex terrain areas. The method includes calculating the underlying surface information of the area to be simulated; based on the underlying surface information, generating a background field of grid meteorological elements such as air temperature, air pressure, and relative humidity by using a terrain parameter weighted regression model; generating a background field of grid wind direction and wind speed by using inverse distance weighted interpolation; based on the background field of grid meteorological elements, performing preliminary error correction by using STMAS multi-grid variational method to obtain the initial values of grid meteorological elements; generating three-dimensional meteorological initial values for the initial values of grid meteorological elements; inputting the three-dimensional wind field initial values into a wind field rapid diagnosis model to obtain a corrected three-dimensional fine grid wind field; based on the corrected three-dimensional fine grid wind field and other three-dimensional meteorological element initial values, simulating and obtaining fine grid air temperature and humidity through an integrated microscale model. The problems that the current rapid simulation method of meteorological elements is not applicable to complex terrain areas and has large simulation errors are solved.

Description

Technical Field

[0001] The present application relates to the technical field of meteorological monitoring, and in particular, to a method and system for rapidly simulating microscale meteorological elements in complex terrain areas. Background Art

[0002] Compared with plain areas, complex terrain areas (such as mountains, hills, etc.) have characteristics such as large terrain undulations, strong local nature of disasters and secondary disasters, and strong suddenness, which pose higher requirements for the refinement and accuracy of local meteorological services and disaster prevention and mitigation work.

[0003] In complex terrain areas, the spatial distribution of meteorological stations is uneven and it is difficult to distinguish local climate characteristics. To improve the quality of meteorological services and disaster prevention and mitigation work, it is only possible to convert site data into fine grid data based on a mathematical model. There are many methods for producing grid meteorological elements based on site meteorological observation data. Among them, spatial interpolation of site meteorological observations, traditional numerical simulation, and rapid simulation are three commonly used methods, and each method has its own advantages and disadvantages.

[0004] Spatial interpolation has the advantage of high calculation efficiency, and some models considering terrain factors also have good applicability in complex terrain areas. However, there are still great uncertainties in the application of interpolation algorithms in meteorological services and disaster prevention and mitigation work in complex terrain areas: for meteorological elements with strong spatial discontinuity such as wind fields and humidity, and strong dependence on atmospheric physical dynamic processes, neither traditional interpolation algorithms nor terrain-considering interpolation algorithms can obtain ideal results; at the same time, affected by the characteristics of the algorithm, it is also very difficult to reflect fine land use types and their change characteristics through interpolation, which further limits the application scenarios of interpolation algorithms.

[0005] Traditional numerical simulation is based on mature atmospheric dynamic, thermodynamic equations and physical parameterization schemes. With the help of data assimilation and other means, grid meteorological elements with high spatio-temporal resolution and low error can be obtained. However, the technical threshold of traditional numerical simulation is high and the operation time is long. In actual disaster prevention and mitigation or service scenarios, the business department often pays more attention to how to produce refined data products as quickly and conveniently as possible while ensuring a certain accuracy rate. Therefore, the application of traditional numerical models is restricted. At the same time, when carrying out high spatio-temporal resolution simulation in complex terrain areas, the computational stability of the model will face huge challenges. In addition, the systematic error of the model will also bring a lot of uncertainties - these factors all limit the application of traditional numerical simulation in meteorological services in complex terrain areas.

[0006] Due to the above problems of spatial interpolation and traditional numerical simulation methods, some scholars at home and abroad have begun to develop some lightweight meteorological element simulation methods. There is no unified name for this type of method, but its common feature is the comprehensive use of spatial interpolation, parameterization schemes and empirical or semi-empirical meteorological element diagnostic relationships to simulate grid meteorological data that take into account the characteristics of the fine underlying surface. This type of method reduces the reliance on the integration of atmospheric dynamic equations and has a faster calculation speed, so it is called a "fast simulation model" here. Guo Xiaoran et al. published "Establishment and Application of Rapid Simulation Methods for Urban Microscale Meteorological Elements" in Chinese Science. For wind, temperature, humidity and radiation elements, they developed a set of urban microscale meteorological element rapid simulation technology by comprehensively applying wind field diagnostic models, land surface parameterization models and temperature and humidity advection equations, providing new ideas for the rapid response of urban fine meteorological services and disaster prevention and mitigation work. However, this technical solution cannot be directly applied to complex terrain areas. The main reasons are: ① Due to the lack of fluid dynamic equation integration, the rapid diagnostic simulation method needs to superimpose the influence of some underlying surface factors on the initial grid meteorological elements in advance. The main influencing factors of micro-scale meteorological elements in plain areas come from land use types. Therefore, this known technology first constructs uniformly distributed wind fields, temperature and humidity fields, and then superimposes the influence of buildings, integral land surface parameterization models and temperature and humidity advection equations in the wind field to better simulate the distribution of urban micro-scale meteorological elements. However, for complex terrain areas, terrain factors cannot be ignored, and how to further introduce terrain effects into the initial wind field, temperature and humidity remains to be studied. ② Although the influence of land use types can be added to temperature and humidity through integral land surface parameterization models and temperature and humidity advection equations, time integration is inevitably accompanied by an increase in nonlinear errors, resulting in the final simulated temperature and humidity deviating from the actual site meteorological observations. How to reduce this error and obtain grid temperature and humidity fields with higher accuracy remains to be further studied. Summary of the invention

[0007] In order to solve the problem that the current rapid simulation method of meteorological elements is not applicable to complex terrain areas and the simulation error is large, the present application provides a rapid simulation method and system for micro-scale meteorological elements in complex terrain areas.

[0008] In the first aspect, the present application provides a method for rapid simulation of micro-scale meteorological elements in complex terrain areas:

[0009] Calculate the underlying surface information of the area to be simulated;

[0010] Based on the underlying surface information, the terrain parameter weighted regression model is used to generate the background field of temperature, pressure and relative humidity grid meteorological elements.

[0011] Generate wind direction and wind speed grid background field using inverse distance weighted interpolation;

[0012] Based on the gridded meteorological element background field, preliminary error correction is carried out using STMAS multigrid variational constrained by terrain parameters; the background field of gridded wind direction and wind speed is corrected using conventional STMAS; the initial values of gridded meteorological elements are obtained;

[0013] For the initial values of gridded meteorological elements, three-dimensional meteorological initial values are generated; the three-dimensional meteorological initial values include the initial value of the three-dimensional wind field and the initial values of other three-dimensional meteorological elements; the initial values of other three-dimensional meteorological elements include the initial value of three-dimensional temperature, the initial value of three-dimensional relative humidity, and the initial value of three-dimensional air pressure;

[0014] The initial value of the three-dimensional wind field is input into the wind field rapid diagnosis model to obtain the corrected three-dimensional fine gridded wind field;

[0015] Based on the corrected three-dimensional fine gridded wind field and the initial values of other three-dimensional meteorological elements, the fine gridded temperature and humidity are simulated through integrating the microscale model.

[0016] By adopting the above technical solutions, a set of rapid simulation technology for microscale meteorological elements more suitable for complex terrain areas is designed, which is used to quickly generate fine gridded data of common meteorological elements such as temperature, humidity, and wind field in complex terrain areas, and solves the problems that the current rapid simulation methods of meteorological elements are not applicable to complex terrain areas and have large simulation errors.

[0017] Optionally, the terrain parameter weighted regression model includes establishing a unary linear regression equation with terrain altitude as the independent variable and meteorological element value as the dependent variable; calculating the comprehensive weight coefficient of each sample according to the difference in underlying surface information between the sample and the grid point; using the weighted least squares method to solve the regression equation coefficients and substituting them into the regression equation to calculate the final interpolation result.

[0018] By adopting the above technical solutions, compared with the spatial interpolation of station data, the terrain parameter weighted regression model interpolation can better consider the influence of fine underlying surface characteristics, and can obtain more realistic wind field and temperature and humidity advection movement results; compared with the traditional numerical simulation method, the calculation amount is smaller, and the three-dimensional data of microscale temperature, humidity, and wind field grid points can be obtained quickly; compared with the conventional rapid simulation model, it is more suitable for complex terrain area conditions and has lower simulation errors.

[0019] Optionally, the method further includes:

[0020] When the meteorological element is temperature, the comprehensive weight coefficient is calculated using the distance weight W d of the sample, the height weight W z of the sample, the slope aspect weight W f of the sample, the vertical stratification weight W l of the sample, the terrain index weight W t of the sample, and the urban land weight W u of the sample;

[0021] When the meteorological element is relative humidity, the comprehensive weight coefficient is calculated using the distance weight W based on samples d , altitude weight W z , slope aspect weight W f , effective terrain weight W e and urban land use weight W u ;

[0022] When the meteorological element is air pressure, the comprehensive weight coefficient is calculated using the distance weight W of the samples d , altitude weight W z and effective terrain weight W e .

[0023] By adopting the above technical solution, different total weight coefficients are used for different meteorological elements, which can improve the interpolation accuracy of the terrain parameter weighted regression model and enhance the simulation effect.

[0024] Optionally, the generation of the grid meteorological element background fields of air temperature, air pressure, and relative humidity using the terrain parameter weighted regression model based on the underlying surface information includes:

[0025] For air temperature, using the station air temperature as a sample, the grid air temperature is interpolated using the terrain parameter weighted regression model;

[0026] For relative humidity, the station relative humidity is converted into the station dew point temperature using the following conversion formula with the station air temperature observation:

[0027]

[0028] where a, b, and c are constants, T a represents the station air temperature (°C), RH represents the station relative humidity (%), e s represents the saturated water vapor pressure, e represents the water vapor pressure, and T d represents the station dew point temperature (°C);

[0029] Using the station dew point temperature as a sample, the grid dew point temperature is interpolated using the terrain parameter weighted regression model;

[0030] Using the grid dew point temperature and grid air temperature data, the grid relative humidity is calculated by inverse calculation using this conversion formula.

[0031] Optionally, the generation of the grid meteorological element background fields of air temperature, air pressure, and relative humidity using the terrain parameter weighted regression model based on the underlying surface information further includes:

[0032] For air pressure, based on the maximum and minimum ground air pressures in the samples, the station air pressure data is standardized to control the air pressure value within the range of 0 to 1;

[0033] Convert the standardized barometric data into logarithmic form;

[0034] Use the barometric data in logarithmic form as samples, and interpolate the grid barometric pressure by using the terrain parameter weighted regression model.

[0035] Optionally, the generation of the grid background field of wind direction and wind speed by using inverse distance weighting interpolation includes:

[0036] Adopt the inverse distance weighting interpolation formula, and the grid wind speed can be directly interpolated based on the station wind speed; where x represents the sample value, the subscript i represents the sample serial number; y represents the interpolation result; d represents the distance from the sample to the prediction grid point; p represents the gain coefficient;

[0037]

[0038] When interpolating the wind direction, convert the station wind direction into standardized east-west component and north-south component according to the following formula; where swspd represents the station wind speed (m / s), swdir represents the station wind direction, su represents the east-west component of the station wind speed, sv represents the north-south component of the station wind speed, swdir Msin 、swdir Mcos represent the standardized east-west component of the station wind direction and the standardized north-south component of the station wind direction respectively;

[0039] su = -swspd·sin(swdir)

[0040] sv = -swspd·cos(swdir)

[0041] swdir Msin =sin(tan -1 (sv,su))

[0042] swdir Mcos =cos(tan -1 (sv,su))

[0043] After calculating swdir Msin 、swdir Mcos , use the inverse distance weighting interpolation formula to interpolate the grid standardized east-west and north-south components, and then calculate the grid wind direction through the following formula, where wdir M is an intermediate variable, wspd, wdir Msin 、wdir Mcos 、u、v represent the grid wind speed, the grid standardized east-west component, the grid standardized north-south component, the grid standardized east-west wind speed, and the grid standardized north-south wind speed respectively;

[0044] wdir M =tan -1(wdir Msin , wdir Mcos )

[0045] u = wspd·sin(wdir M )

[0046] v = wspd·cos(wdir M )

[0047] wdir = tan -1 (-v, -u).

[0048] Optionally, the STMAS multigrid variation includes using the following formula:

[0049]

[0050]

[0051] where x b represents the grid point background field; y k represents the site observation increment used in multigrid solution, y o represents the site observation value, x k represents the grid point background field correction increment to be calculated for each layer in multigrid solution, x represents the finally obtained corrected background field, h k represents the observation operator that interpolates the grid point background field to the site coordinates, the superscript k represents different layers of the multigrid; kmax represents the set number of multigrid levels; O represents the ratio of the empirical observation error covariance to the background field error covariance, usually taking values between 0 and 1; T represents taking the transpose; J k is the name of the cost function for the k-th layer.

[0052] Optionally, the underlying surface information includes at least one of slope aspect, effective terrain index, vertical stratification, and terrain index.

[0053] Calculating the underlying surface information of the area to be simulated includes:

[0054] For the slope aspect, it is calculated using the following formula:

[0055]

[0056] where asp represents the slope aspect, dz represents the elevation difference between two grid points, dy represents the vertical distance between two grid points, and dx represents the horizontal distance between two grid points;

[0057] For the effective terrain index, the minimum elevation value within a circular area with a radius of 15 km around each grid point is statistically obtained to get the minimum elevation data, and the average value is calculated for each grid point within the circular area with a radius of 15 km on the minimum elevation data to achieve data smoothing;

[0058] Subtract the smoothed minimum elevation data from the original elevation data to obtain the initial effective terrain height;

[0059] For each grid point of the initial effective terrain height, calculate the grid point average value within a circular area with a radius of 6 km for smoothing to obtain the finally smoothed effective terrain height;

[0060] Calculate the index I according to the following formula 3c ;

[0061]

[0062] where h c is the finally smoothed effective terrain height; h 2 and h 3 are the 2D terrain and 3D terrain thresholds respectively, which are set to 100 m and 600 m respectively;

[0063] Calculate the index I according to the following formula 3a ;

[0064]

[0065] where h a is the average value of the effective terrain height calculated with the distance inverse as the weight within a range of 8 km around the grid point;

[0066] Calculate the effective terrain index I according to the following formula 3d ;

[0067] I 3d = max(I 3c , I 3a )

[0068] For vertical stratification, count the minimum altitude within a horizontal radius of 10 km of the grid point;

[0069] Subtract the minimum altitude from the altitude of the current grid point. If the difference exceeds 150 m, the layer where the grid point is located is 2; otherwise, the layer where the grid point is located is 1;

[0070] For the terrain index, count the minimum terrain altitude within a horizontal radius of 10 km of the grid point;

[0071] Perform Gaussian smoothing on the grid point with the minimum terrain altitude with a radius of 10 km;

[0072] Subtract the minimum terrain altitude from the original terrain altitude of the grid point to obtain the terrain index data.

[0073] In a second aspect, the present application provides a rapid simulation system for microscale meteorological elements in complex terrain areas, including:

[0074] The underlying surface information calculation module is used to calculate the underlying surface information of the area to be simulated;

[0075] The complex terrain two-dimensional grid background field generation module is used to generate the grid meteorological element background fields of air temperature, air pressure, and relative humidity based on the underlying surface information by using the terrain parameter weighted regression model; generate the grid background field of wind direction and wind speed by using inverse distance weight interpolation; the error correction module is used to perform preliminary error correction based on the grid meteorological element background fields by using the STMAS multi-grid variational method constrained by terrain parameters; correct the grid background field of wind direction and wind speed by using the conventional STMAS; obtain the initial values of grid meteorological elements;

[0076] The three-dimensional meteorological element initial field generation module is used to generate three-dimensional meteorological initial values for the initial values of grid meteorological elements; the three-dimensional meteorological initial values include the initial value of the three-dimensional wind field and the initial values of other three-dimensional meteorological elements; the initial values of other three-dimensional meteorological elements include the initial value of three-dimensional air temperature, the initial value of three-dimensional relative humidity, and the initial value of three-dimensional air pressure;

[0077] The wind field diagnosis module is used to input the initial value of the three-dimensional wind field into the wind field rapid diagnosis model to obtain the corrected three-dimensional fine grid wind field;

[0078] The microscale model calculation module is used to simulate and obtain the fine grid air temperature and humidity by integrating the microscale model based on the corrected three-dimensional fine grid wind field and the initial values of other three-dimensional meteorological elements.

[0079] In summary, the present application at least includes the following beneficial technical effects:

[0080] 1. Solve the problems that the current rapid simulation method of meteorological elements is not applicable to complex terrain areas and has large simulation errors. Brief Description of the Drawings

[0081] Figure 1 It is a schematic flow chart of a rapid simulation method for microscale meteorological elements in a complex terrain area in an embodiment of the present application;

[0082] Figure 2 It is a flow block diagram of a rapid simulation method for microscale meteorological elements in a complex terrain area in an embodiment of the present application;

[0083] Figure 3 It is a schematic diagram of the average annual wind field at a height of 10 m in the Tainong Orchard simulated in an embodiment of the present application;

[0084] Figure 4 It is a schematic diagram of the average annual air temperature at a height of 2 m in the Tainong Orchard simulated in an embodiment of the present application;

[0085] Figure 5 It is a schematic structural diagram of a rapid simulation system for microscale meteorological elements in a complex terrain area in an embodiment of the present application. Detailed implementation manners

[0086] In order to make the objectives, technical solutions and advantages of the present application more clear and understandable, the present application will be further described in detail below with reference to the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain the present application and are not used to limit the present application.

[0087] The terms used in the following embodiments of the present application are only for the purpose of describing specific embodiments and are not intended to limit the present application. As used in the specification and appended claims of the present application, the singular forms "a", "an", "the", "above-mentioned", "said", and "this" are also intended to include the plural forms unless the context clearly indicates otherwise. It should also be understood that the term "and / or" used in the present application refers to and includes any or all possible combinations of one or more of the listed items. The term "exemplary" means "serving as an example, embodiment, or illustration", and any embodiment described as "exemplary" here does not have to be construed as superior or better than other embodiments. The terms "first" and "second" are only used for descriptive purposes and cannot be understood as implying or suggesting relative importance or implicitly indicating the quantity of the indicated technical features. Thus, features defined with "first" and "second" may explicitly or implicitly include one or more of such features. In the description of the embodiments of the present application, unless otherwise specified, the meaning of "a plurality" is two or more.

[0088] Reference Figure 1 , a rapid simulation method for microscale meteorological elements in complex terrain areas, comprising the following steps:

[0089] S101: Calculate the underlying surface information of the area to be simulated;

[0090] S102: Based on the underlying surface information, use a terrain parameter weighted regression model to generate grid meteorological element background fields of air temperature, air pressure, and relative humidity;

[0091] S103: Use inverse distance weighted interpolation to generate a grid background field of wind direction and wind speed;

[0092] S104: Based on the grid meteorological element background fields, use terrain parameter-constrained STMAS multi-grid variational to perform preliminary error correction;

[0093] S105: Use conventional STMAS to correct the grid background field of wind direction and wind speed;

[0094] S106: For the initial values of the grid meteorological elements obtained after correction, generate initial values of meteorological elements of three-dimensional air temperature, three-dimensional relative humidity, three-dimensional air pressure, and three-dimensional wind field;

[0095] S107: Input the initial three-dimensional wind field into the fast wind field diagnostic model to obtain the corrected three-dimensional fine grid wind field;

[0096] S108: Based on the corrected three-dimensional fine grid wind field and the initial meteorological values of three-dimensional air temperature, three-dimensional relative humidity, and three-dimensional air pressure, simulate the fine grid air temperature and humidity through integrating the microscale model.

[0097] The embodiment of this application designs a fast simulation method for microscale meteorological elements more suitable for application in complex terrain areas. This method comprehensively applies terrain parameter weighted regression, wind field diagnostic model, land surface process model, temperature and humidity advection equation, and data assimilation method to achieve the fast simulation of common meteorological elements such as microscale air temperature, humidity, and wind field in complex terrain areas. It is not only a perfection of the fast simulation method and the development method of fine meteorological data products, but also can provide technical support for fine meteorological services and meteorological disaster prevention and mitigation work in complex terrain areas.

[0098] See Figure 2 , and the following provides a detailed description of the fast simulation method for microscale meteorological elements in complex terrain areas provided in this embodiment:

[0099] (1) Collect and organize multi-source data such as station meteorological observations and geographic information data, carry out data quality control and standardization processing, and calculate underlying surface information such as slope aspect, effective terrain index, vertical stratification, and terrain index required for subsequent operations.

[0100] Among them, station meteorological data includes 2-meter air temperature, 2-meter relative humidity, surface air pressure, 10-meter wind direction and wind speed; geographic information data includes digital terrain elevation (hereinafter referred to as DEM) and land use type data. Land use type refers to dividing land into different categories according to the natural and economic characteristics of the land for reasonable utilization and management. Common land use types include cultivated land, garden land, forest land, grassland, water area, urban and rural residential and industrial and mining land, transportation land, etc.

[0101] Slope aspect refers to the inclination direction of a certain point on the ground surface; it can be determined by measuring the slope and direction of the terrain. The measurement method of slope aspect is based on the north, that is, the due north direction is 0°, increasing clockwise, with the maximum value of 359°59'59'', and the unit is degree. The calculation formula is as shown in formula (1), where asp represents slope aspect, dz represents the elevation difference between two grid points, dy represents the vertical distance between two grid points, and dx represents the horizontal distance between two grid points.

[0102]

[0103] The effective terrain index is a standardized index based on DEM calculation to determine the terrain complexity of a region. The value of the effective terrain index ranges from 0 to 1. The larger the value, the more complex the terrain at the grid point location. The calculation process of the effective terrain is as follows:

[0104] ① Statistically calculate the minimum elevation value within a circular area with a radius of 15 km around each grid point in the DEM to obtain the minimum elevation data, and calculate the average value for each grid point according to the circular area with a radius of 15 km on the minimum elevation data to achieve data smoothing;

[0105] ② Subtract the smoothed minimum elevation data from the original DEM to obtain the initial effective terrain height;

[0106] ③ For each grid point of the initial effective terrain height, calculate the grid point average value according to a circular area with a radius of 6 km for smoothing to obtain the finally smoothed effective terrain height.

[0107] ④ Calculate the index I according to equation (2) 3c . Where h c is the finally smoothed effective terrain height; h 2 and h 3 are the 2D terrain and 3D terrain thresholds respectively, which are set to 100 m and 600 m respectively.

[0108]

[0109] ⑤ Calculate the index I according to equation (3) 3a . Where h a is the average value of the effective terrain height calculated with the distance inverse as the weight within 8 km around the grid point, and h 2 and h 3 are the same as ④.

[0110]

[0111] ⑥ Finally, calculate the effective terrain index I according to equation (4) 3d .

[0112] I 3d = max(I 3c , I 3a ) (4)

[0113] In complex terrain areas, due to the frequent existence of inversion layers (i.e., the temperature increases with height), when generating the temperature background field subsequently, it is necessary to consider the variability differences of temperature in different height levels, so it is necessary to calculate the vertical stratification data in advance.

[0114] The methods for calculating vertical stratification include:

[0115] ① Statistically calculate the minimum altitude within a horizontal radius of 10 km of a certain grid point;

[0116] ② Subtract the minimum elevation of the current grid point from the elevation of the grid point. If the difference exceeds 150 m, the level where the grid point is located is 2; otherwise, the level where the grid point is located is 1.

[0117] ③ Repeat ① and ② for each grid point to obtain the vertical stratification grid point data.

[0118] Considering the influence of cold air deposition in the low-altitude area in the complex terrain area, it is necessary to calculate the terrain index to highlight the characteristics of local depressions and valley bottoms.

[0119] The calculation method of the terrain index is as follows:

[0120] ① Statistically analyze the minimum terrain elevation within a horizontal radius of 10 km of a certain grid point;

[0121] ② Perform Gaussian smoothing on the grid point with the minimum terrain elevation with a radius of 10 km;

[0122] ③ Subtract the minimum terrain elevation from the original terrain elevation of the grid point to obtain the terrain index data.

[0123] (2) Based on the various underlying surface information calculated above, use the terrain parameter weighted regression model to generate the grid point meteorological element background fields of temperature, pressure, and relative humidity; use the inverse distance weight interpolation to generate the grid point background fields of wind direction and wind speed.

[0124] The terrain parameter weighted regression model used in the present invention is improved based on the PRISM (Parameter-elevation Regressions on Independent Slopes Model) model. The PRISM model believes that the most important factor affecting the spatial distribution of mountain meteorological elements is the terrain elevation, and the meteorological elements show an approximately linear relationship under the combined influence of elevation and other geographical factors. The PRISM model takes the terrain elevation as the independent variable and the meteorological element value as the dependent variable to establish a unary linear regression equation; according to the difference in the underlying surface information between the sample and the grid point, calculate the comprehensive weight coefficient of each sample; use the weighted least squares method to solve the regression equation coefficients and substitute them into the regression equation to calculate the final interpolation result.

[0125] The specific PRISM regression equation can adopt the following formula:

[0126] Y = β 1 X + β 0 (5)

[0127] In the formula, X represents the grid point elevation, Y represents the interpolation result, and β 1 and β 0 are the regression coefficients of the equation;

[0128] The regression coefficient β of the equation1 and β 0 can be calculated through the following system of equations:

[0129]

[0130] In the formula, I 3d is the effective terrain index, x i represents the altitude of the location where each sample is located, y i represents the meteorological element value of each sample, w i represents the comprehensive weight coefficient of each sample, represents the weighted average of sample altitudes, represents the weighted average of sample meteorological element values.

[0131] Weighted average of sample altitudes can be calculated through the following formula:

[0132]

[0133] Weighted average of sample meteorological elements can be calculated through the following formula:

[0134]

[0135] Comprehensive weight coefficient w i is calculated based on the distance weight W d of the sample, height weight W z slope aspect weight W f vertical stratification weight W l effective terrain weight W e terrain index weight W t Specifically, the comprehensive weight coefficient w i The calculation formula is:

[0136]

[0137] In the formula, F d and F z are the respective proportions of the distance weight W d and the height weight W z The sum of F d and F z is equal to 1. Among them, F d and F z can be flexibly determined according to the type of meteorological elements and local terrain characteristics by experience; the commonly used setting is F d = 0.7, F z = 0.3.

[0138] Distance weight W of the sample dIt can be calculated by the following formula:

[0139]

[0140] In the formula, d represents the distance from the sample point to the interpolation grid point, r m is the influence radius, and a is the gain coefficient.

[0141] When the distance d from the sample point to the interpolation grid point is less than or equal to the influence radius r m , the similarity degree of meteorological elements between the two points is the highest; when the distance d from the sample point to the interpolation grid point gradually becomes greater than the influence radius r m , as the distance increases, the similarity degree of meteorological elements between the two points will gradually decrease with the increase of the distance. Using such a form can better avoid the "bull's-eye" phenomenon in the interpolation result (referring to the circular phenomenon centered on the interpolation point formed by some too large or too small data during the interpolation process).

[0142] The specific numerical range of the influence radius r m is determined according to the actual situation (such as the number of samples, spatial density, etc.) and combined with manual experience. In the optional embodiment of the present application, the influence radius r m is, for example, set to 1 km. The gain coefficient a is determined according to the actual situation and combined with manual experience, and usually takes a value of 2.

[0143] The height weight W of the sample z is calculated by the following formula:

[0144]

[0145] In the formula, Δz is the absolute value of the height difference between the sample point and the interpolation grid point; Δz m and Δz x respectively represent the minimum height difference and the maximum height difference between each sample point and the interpolation grid point; b is the gain coefficient, usually taken as 1, that is, considering the linear influence of height change on the element similarity; I 3d is the effective terrain index, that is, as the terrain of the interpolation target grid point gradually flattens, the height weights of each sample will tend to be consistent, and the calculation of the effective terrain index is as in (4). The similarity degree of meteorological elements between two points also has a great relationship with the height difference. When it is the minimum height difference, the similarity degree of their meteorological elements is the highest; when it exceeds a certain height difference, its correlation gradually decreases or even becomes completely irrelevant.

[0146] The slope aspect weight W of the sample f is calculated by the following formula:

[0147]

[0148] In the formula, Δf represents the aspect difference, c is the gain coefficient, usually taking the value of 1; I 3d is the effective terrain index. That is, as the terrain of the interpolation target grid point gradually flattens, the aspect weights of each sample will tend to be consistent. The calculation of the effective terrain index is as shown in (4). For two spatial points with similar aspects, the similarity of meteorological elements should be higher.

[0149] Considering that an inversion layer is likely to form in complex terrain areas and there may be different vertical temperature lapse rates at different heights from the ground, when interpolating for a specific interpolation grid point, samples within the same vertical layer as this grid point should be given higher weights.

[0150] The vertical layer weight W of the sample l is calculated through the following formula:

[0151]

[0152] In the formula, Δl is the absolute value of the layer difference between the sample point and the interpolation grid point, c is the gain coefficient; Δz m and Δz x respectively represent the minimum height difference and the maximum height difference between each sample point and the interpolation grid point.

[0153] W t represents the terrain index weight of the sample. The purpose of setting this weight is to capture the influence of cold air. Grid points located in depressions and valley bottoms are more vulnerable to the influence of cold air, and higher weights should be given to samples in similar areas. The terrain index weight W of the sample t is calculated through the following formula:

[0154]

[0155] In the formula, Δt represents the absolute value of the difference in terrain index between the sample point and the interpolation grid point, Δt n and Δt x respectively represent the minimum set threshold and the maximum set threshold set according to experience; optionally, Δt n takes 0.1 km, and Δt x takes 0.5 km; z represents the gain index, usually taking 2.

[0156] Considering the representativeness of samples in complex terrain areas and flat terrain areas, for each grid point, the effective terrain weight W e is used to weight or down-weight different types of samples. The formula for the effective terrain weight is as follows:

[0157]

[0158] Among them, I 3ds and I 3dgThey are the effective terrain indices of the meteorological observation point and the current interpolation grid point respectively. The calculation of the effective terrain index is described in the previous text.

[0159] In order to better consider the urban climate effect in the background field, the present invention further adds the urban land use weight W to the comprehensive weight coefficient. u At this time, the comprehensive weight coefficient is calculated by the following formula:

[0160]

[0161] The urban land use weight W u is calculated by the following formula:

[0162]

[0163] In the formula, Δurb represents whether the sample point and the interpolation grid point are of the same land use type; if so, Δurb = 0, if not, Δurb ≠ 0. d is a downweighting coefficient set according to experience and can be set to 2. At the same time, when performing weighted least squares fitting, when the interpolation grid point is urban land use, the proportion of the number of urban meteorological observation stations in the sample is increased. By introducing the urban land use weight W u and controlling the proportion of urban site samples, the influence of urban land use on meteorological elements can be better considered in the interpolation.

[0164] In addition, since the relationship between relative humidity and surface air pressure and altitude is not linear, directly using a linear regression model is not applicable. Therefore, the terrain weighted regression model used in the present invention linearizes the samples.

[0165] For relative humidity, the interpolation process is as follows:

[0166] ① Use the site air temperature observation to convert the relative humidity into the dew point temperature. The conversion formula is shown in formulas (21)-(23). In formula (21), a', b', and c' are constants, taking values of 611.21, 17.502, and 240.97 respectively; T a represents the air temperature (°C), RH represents the relative humidity (%), e s represents the saturated water vapor pressure, e represents the water vapor pressure, and T d represents the dew point temperature (°C).

[0167]

[0168] e = RH / 100·es (22)

[0169]

[0170] ②Using the station dew point temperature as a sample, the grid dew point temperature is interpolated by the terrain weighted regression method; using the station air temperature as a sample, the grid air temperature is interpolated by the terrain weighted regression method.

[0171] ③Using the grid dew point temperature data and the grid air temperature data, referring to equations (21)-(23), the grid relative humidity is calculated by back-calculation.

[0172] For the surface air pressure, the interpolation process is as follows:

[0173] ①Statistically analyze the maximum and minimum values of the surface air pressure in the sample, and standardize them using formula (24) to control the air pressure value within the range of 0-1. Among them, y norm represents the standardized air pressure sample, y min represents the minimum value of the sample air pressure, and y max represents the maximum value of the sample air pressure.

[0174]

[0175] ②Convert the standardized sample data into logarithmic form according to formula (25), where y log represents the sample in logarithmic form.

[0176]

[0177] ③Based on the sample in logarithmic form, the grid data is interpolated by the terrain weighted regression method.

[0178] ④According to equations (24) and (25), the grid data obtained in ③ is converted back to the surface air pressure by back-calculation to complete the interpolation of the surface air pressure.

[0179] Due to the different variable attributes of air temperature, relative humidity, and air pressure, not all parameters in the comprehensive weight coefficient are applicable. Therefore, different combinations of weight coefficients are used when interpolating different meteorological elements. The specific combinations are shown in Table 1.

[0180] Table 1 Combinations of weight coefficients used when interpolating grid background fields of different elements

[0181]

[0182]

[0183] Since the relationship between wind direction and wind speed and altitude is relatively complex, the terrain parameter weighted regression interpolation cannot be used to obtain the grid background field. In the embodiment of the present application, the inverse distance weight interpolation method is used. The inverse distance weight interpolation formula is as shown in (26), where x represents the sample value, the subscript i represents the sample serial number; y represents the interpolation result; d represents the distance from the sample to the prediction grid point; p represents the gain coefficient, usually taking 2.

[0184]

[0185] Before interpolating the station wind direction to a grid, equations (27)-(30) need to be used to convert the station wind direction into standardized east-west and north-south components. Here, swspd represents the station wind speed (m / s), swdir represents the station wind direction (radians), su represents the east-west component of the station wind speed (m / s), and sv represents the north-south component of the station wind speed (m / s). swdir Msin and swdir Mcos represent the standardized east-west and north-south components of the station respectively.

[0186] su = -swspd·sin(swdir) (27)

[0187] sv = -swspd·cos(swdir) (28)

[0188] swdir Msin = sin(tan -1 (sv,su)) (29)

[0189] swdir Mcos = cos(tan -1 (sv,su)) (30)

[0190] After calculating swdir Msin and swdir Mcos , substitute swspd, swdir Msin and swdir Mcos into x in equation (26), and interpolate to obtain the grid standardized east-west and north-south components. Then, through the inverse equations (31)-(34), obtain the grid wind direction, where wdir M is an intermediate variable, and wspd, wdir Msin , wdir Mcos , u, and v represent the grid wind speed, grid standardized east-west component, grid standardized north-south component, grid east-west wind speed, and grid north-south wind speed respectively.

[0191] wdir M = tan -1 (wdir Msin ,wdir Mcos ) (31)

[0192] u = wspd·sin(wdir M ) (32)

[0193] v = wspd·cos(wdir M ) (33)

[0194] wdir = tan-1 (-v, -u) (34)

[0195] Through the above algorithm steps, the present invention can interpolate the grid-based air temperature, air pressure, relative humidity, wind direction, and wind speed background fields in complex terrain areas based on site meteorological data and geographic information data.

[0196] (3) Based on the grid-based meteorological element background field, use the STMAS multi-grid variational constrained by terrain parameters for preliminary error correction to obtain the initial values of grid-based meteorological elements.

[0197] Considering that there are errors in the terrain parameter weighted regression interpolation model, in order to obtain a more accurate initial field, it is necessary to further correct its errors. STMAS is a variational data assimilation algorithm based on the multi-grid idea, and its equations are shown in Equations (35)-(38).

[0198]

[0199] Among them, x b represents the grid-based background field; y k represents the site observation increment used in multi-grid solution, y o represents the site observation value, x k represents the grid-based background field correction increment that needs to be calculated for each layer in multi-grid solution, x represents the finally obtained corrected background field, h k represents the observation operator that interpolates the grid-based background field to the site coordinates. The superscript k represents different layers of the multi-grid, kmax represents the set number of multi-grid levels, which is determined according to the actual situation of the grid dimension; O represents the ratio of the empirical observation error covariance to the background field error covariance, usually taking values between 0 and 1; T represents taking the transpose; J k is the name of the cost function for the kth layer.

[0200] When solving the STMAS equation, first interpolate the grid-based background field whose error needs to be corrected to a 3×3 dimension grid, then calculate the variables in Equation (35) using Equation (36), and use the quasi-Newton method to solve Equation (37) to obtain x 1 . Then interpolate the background field to a 6×6 dimension grid, calculate the variables in Equation (35) using Equation (36), and also solve Equation (37) to obtain x 2 . In this way, the grid dimension of the interpolated background field is increased step by step in a 2-fold ratio, and the corresponding level background field increment x k is solved until the dimension is the same as the original dimension of the background field. Finally, the corrected background field is calculated using Equation (38).

[0201] It should be understood that since the observation operator h used in STMAS operationk Bilinear interpolation is adopted, so it is easy to cause local "bull's-eye" linearity in complex terrain areas. In the present invention, the observation operator h k is replaced by the terrain parameter weighted regression model described above, thereby realizing STMAS multigrid variational with terrain parameter constraints.

[0202] In this way, the grid point air temperature, relative humidity, and barometric background fields interpolated from the terrain parameter weighted regression model are corrected using the terrain parameter weight to constrain STMAS and the corresponding station observation data; the grid point wind direction and wind speed background fields are corrected using conventional STMAS; and the initial values of the required ground grid meteorological elements are obtained.

[0203] (4) Using the land use type and empirical parameterization scheme, generate the initial meteorological values of three-dimensional air temperature, three-dimensional humidity, three-dimensional barometric pressure, and three-dimensional wind field.

[0204] First, according to actual needs, set the vertical grid and calculate the altitude of each layer of the grid. The vertical grid used in the present invention adopts the sigma terrain-following coordinate system, and the height of each grid point in the three-dimensional grid is calculated as follows:

[0205]

[0206] Among them, zg represents the terrain altitude, hmax represents the altitude of the highest layer of the three-dimensional grid; sigma represents the coordinate value of each layer of the sigma terrain-following coordinate system, which is set according to actual needs (for example, the lowest layer can be set to 5, the highest layer altitude can be set to hmax, and the values of the intermediate layers decrease monotonically); zp is the actual altitude of each grid point in the three-dimensional space of the sigma coordinate system. The subscripts i, j, and k respectively represent the grid point indices in the east-west, north-south, and vertical directions of the three-dimensional grid.

[0207] For air temperature, first interpolate the initial grid point air temperature to the bottom layer of the grid according to the vertical lapse rate of air temperature obtained by terrain parameter weighted regression (i.e., β 1 ), and then obtain the initial three-dimensional air temperature of each grid point in the grid through equation (40):

[0208] t i,j,k = gamma i,j,k ·(zp i,j,k - zp i,j,k-1 ) + t i,j,k-1 (40)

[0209] Among them, gamma represents the vertical lapse rate of air temperature at the current grid point, and t represents air temperature.

[0210] The vertical lapse rate of air temperature gamma is calculated as follows:

[0211]

[0212] where h is the height set according to experience, usually set to 500 meters; gamma 0 is the empirically set lapse rate of air temperature, usually -0.006. Equation (41) actually means that when the height above the ground is less than h, the lapse rate of air temperature in each layer linearly changes with height to gamma 0 ; when the height above the ground is greater than h, the lapse rate of air temperature remains gamma 0 unchanged.

[0213] For relative humidity, first convert it to dew point temperature with reference to equations (21)-(23), and then generate the three-dimensional grid dew point temperature in a similar way to air temperature, except that gamma 0 is usually set to -0.003. After obtaining the three-dimensional dew point temperature, then invert equations (21)-(23) to convert the three-dimensional dew point temperature to three-dimensional relative humidity.

[0214] For air pressure, generate the initial three-dimensional air pressure according to equations (42)-(44)

[0215]

[0216] where a 1 、a 2 、a 3 and a 4 take values of -3.9082017e-2, -1.1526465e-3, 3.2891937e-5, -2.0494958e-7 respectively; b 1 、b 2 、b 3 and b 4 take values of -4.9244637e-3, -1.2984142e-6, -1.5701595e-6, 1.5535974e-8 respectively. A i,j,k 、B i,j,k are intermediate variables.

[0217] After obtaining the initial three-dimensional air pressure, replace zp in (42)-(44) with the terrain elevation zg, and calculate the estimated ground air pressure Equation (45) uses terrain-weighted regression interpolation for the ground air pressure to correct the initial three-dimensional air pressure

[0218]

[0219] ​In the formula, k represents the grid point index in the vertical direction. For the horizontal wind speed, it is calculated according to Equation (46) based on different underlying surface types.

[0220]

[0221] Among them, wspd i,j,k represents the initial value of the three-dimensional wind speed finally generated; represents the initial value of the two-dimensional wind speed near the ground; h 0 represents the height of the wind speed initial value from the ground (usually 10 meters); r 0 represents the roughness length of the ground land use type; h mix represents the distance from the ground to the top of the mixed layer, which can be empirically set to 800 meters; is the altitude of the sounding data (the height of the 700 hPa sounding data can be used); is the wind speed at the mixed layer height; represents the wind speed at the sounding observation height (the 700 hPa sounding observation wind speed can be used); the parameters obtained by linear regression fitting. Equation (46) represents that the wind speed changes logarithmically within the mixed layer, and linearly above the mixed layer to the wind speed at the sounding observation height; when the height exceeds the sounding height, the wind speed remains unchanged.

[0222] The specific classification of land use types and the corresponding roughness lengths are described later.

[0223] For the horizontal wind direction, the initial value of the three-dimensional grid wind field is calculated according to Equation (47).

[0224]

[0225] Among them, wdir i,j,k represents the initial value of the three-dimensional wind direction finally generated; represents the initial value of the two-dimensional wind speed near the ground; represents the wind direction at the sounding height, and the other variables are the same as (46). Equation (47) represents that the wind direction remains unchanged within the mixed layer, and gradually changes to the wind direction at the sounding height above the mixed layer; when the height exceeds the sounding height, the wind direction remains unchanged.

[0226] For the vertical wind speed, the initial value is set to 0.

[0227] After generating the initial values of three-dimensional air temperature, three-dimensional air pressure, three-dimensional relative humidity, horizontal wind speed, horizontal wind direction, and vertical wind speed according to the above steps, considering the influence of the complex terrain slope flow effect (i.e., the valley wind thermal circulation effect), a parameterization scheme is used to superimpose a correction coefficient on the initial wind field. Since this parameterization scheme is relatively complex and not the key content of the present invention, it will not be elaborated here. For example, reference can be made to "Scire J.S., Robe F.R. Fine-scale application of the CALMET meteorological model to a complex terrain site. 1997. In 'Air & Waste Management Association's 90th Annual Meeting & Exhibition'. pp. 1-16. (Toronto, Ontario, Canada)".

[0228] (5) Input the initial value of the three-dimensional wind field into the fast diagnostic model of the wind field, and after being adjusted by the law of mass conservation, the revised three-dimensional fine grid wind field is obtained.

[0229] The calculation process of the fast diagnostic model of the wind field used in the present invention is as follows:

[0230] ① Superimpose the terrain-following coordinate system correction coefficient on the initial three-dimensional wind field using equations (48)-(50). Where are the components of the initial three-dimensional wind field in the east-west, north-south, and vertical directions; are the components of the wind field after coordinate transformation.

[0231]

[0232]

[0233] ② Substitute the wind field components calculated in ① into the Poisson equation, equation (51), and use over-relaxation iteration to solve for the Lagrangian operator λ i,j,k .

[0234]

[0235] In the formula, x, y, and z represent the grid in the east-west, north-south, and vertical directions respectively; and are the horizontal and vertical scale coefficients respectively.

[0236] It should be noted that controls the adjustment intensity of the model for the vertical and horizontal wind components during diagnosis: when approaches infinity, the model only adjusts the vertical wind component; when When approaching 0, only the horizontal wind component is adjusted.

[0237] To calculate adaptively A stability parameterization scheme is used here. The equations of the parameterization scheme are as shown in Eqs. (52) and (53):

[0238]

[0239]

[0240] where \(p_t\) i,j,k represents the three-dimensional potential temperature, and its calculation formula is shown in Eq. (57) later; \(wspd\) i,j,k represents the three-dimensional wind speed, and \(max\) represents taking the maximum of the two, that is, the wind speed is not less than 0.2 m / s during calculation; represents the absolute difference between the grid elevation and the highest terrain elevation within a 10 km radius; \(g\) represents the acceleration due to gravity, with a value of 9.8; \(F\) i,j,k is the Strouhal number. Finally, the calculated The maximum value does not exceed 30, and the minimum value is not less than 0.033. After automatically calculating only need to assign a value to according to experience, and usually it can be simply assigned a value of 1.

[0241] ③ After solving to obtain \(\lambda\) i,j,k substitute it into Eqs. (54)-(56) to calculate the wind components corrected by the terrain-following coordinate system and satisfying the mass conservation constraint Finally, substitute into the left side of Eqs. (48)-(50) and back-calculate the corresponding wind components on the right side, which is the corrected three-dimensional fine grid wind field.

[0242]

[0243]

[0244] (6) Based on the three-dimensional fine grid wind field, input the initial values of potential temperature, specific humidity, and air pressure, and through integrating a microscale model, simulate the three-dimensional air temperature and humidity considering the influence of fine land use.

[0245] The initial value of potential temperature is calculated using the initial values of air temperature and air pressure, and the formula is as shown in (57). Where \(t\) i,j,k is the air temperature, and \(P\) i,j,k is the air pressure.

[0246]

[0247] The initial value of specific humidity is calculated using the initial values of air temperature, air pressure, and relative humidity, and the formulas are as shown in (58)-(59). Where \(rh\)i,j,k is the initial value of three-dimensional relative humidity; esat i,j,k is the saturation water vapor pressure; q i,j,k is the specific humidity.

[0248]

[0249] The microscale model used in the present invention includes three main technical points: the land surface process parameterization model, the temperature and humidity advection equation, and the Nudging data assimilation.

[0250] The land surface process parameterization model used in the present invention considers a thin layer of soil on the ground and calculates the surface temperature using the forced recovery method. The energy balance equation of the model is:

[0251]

[0252] where, T g is the surface temperature (the generation scheme of its initial field is shown later), t is the integration time step of the land surface process model, R n is the net radiation, H e and H s are the ground latent heat flux and the ground sensible heat flux respectively, H m is the heat flux flowing into the deep soil layer. C g is the heat capacity per unit volume, C g = C s *(86400K e ) 0.5 , C s represents the heat capacity, K e is the heat transfer coefficient. K e and C s are determined by the surface vegetation cover type (see Table 2).

[0253] Table 2 Land use types and corresponding thermodynamic parameters used in the present invention

[0254]

[0255] In formula (60), the net radiation R n can be calculated according to the following formula (61):

[0256] R n = R s + I ↑ + I ↓ (61)

[0257] In the formula, R s is the net shortwave radiation flux, and the calculation method is formula (62)-(63):

[0258] Rs = S 0 (1 - A)D se p m sin(s h ) (62)

[0259]

[0260] D se = 1.000110 + 0.034221 * cos(d0) + 0.001280 * sin(d0) + 0.000719 * cos(2 * d0) + 0.000077 * sin(2 * d0) (64)

[0261] d0 = (noday - 1) * 2 * 3.14159 / 365 (65)

[0262]

[0263] Among them, S 0 is the solar constant (with a value of 1370.0); A is the surface short - wave albedo (see Table 2); D se is the sun - earth distance factor, p is the atmospheric transparency coefficient, m is the atmospheric optical mass, d0 represents the angular coordinate of the Earth relative to the Sun calculated according to the current date, noday represents the total number of days of the current date in the whole year, int represents taking the integer, and glat represents the latitude. s h and s θ are the solar altitude angle and the zenith angle respectively. The specific calculation method can adopt the existing method. For example, it can refer to "Xu Yumao, Liu Hongnian, Xu Guiyu. 2000. Introduction to Atmospheric Science. Nanjing: Nanjing University Press.".

[0264] In addition, according to the solar altitude angle and the terrain gradient, the sun - shaded area can be calculated. The specific calculation method can adopt the existing method. For example, it can refer to "Xu Yumao, Liu Hongnian, Xu Guiyu. 2000. Introduction to Atmospheric Science. Nanjing: Nanjing University Press.". The solar radiation amount in these areas is taken as 0.6 times the solar short - wave radiation reaching the ground in the unshaded area, that is, it is considered that the intensity of the scattered radiation and the reflected radiation received is equivalent to 0.6 times the solar short - wave radiation.

[0265] In Equation (61), I ↑ is the upward long - wave radiation flux from the ground, I ↓ is the downward long - wave radiation flux from the atmosphere, and the calculation is as follows in Equations (67) - (68):

[0266]

[0267] In the formula, T a is the near - surface air temperature, T gis the surface temperature; ε a , ε g are the infrared emissivities of the atmosphere (taking a value of 0.725) and the surface (see Table 2), respectively. σ is the Boltzmann constant, taking a value of 5.67e －8 .

[0268] The surface sensible heat flux H s is calculated according to Equation (69):

[0269] H s = ρ a C p C H |V a |(T g - T a )(69)

[0270] where ρ a is the surface air density (calculated using surface pressure, relative humidity, and air temperature), C p is the specific heat at constant pressure of air (taking a value of 1004.67), C H is the drag coefficient (which can take a value of 2.5e-3), and V a is the near-surface wind speed.

[0271] The surface latent heat flux H e is calculated according to Equation (70):

[0272] H e = ηρ a LC H |V a |(q sg - q a )(70)

[0273] where η represents the surface saturation (see Table 2), L is the latent heat of vaporization of liquid water (taking a value of 2.5e6), q sg is the surface saturated specific humidity calculated based on the surface temperature (T g ), surface pressure (PSFC 0 ), as shown in Equation (71). Here, a” and b” are state coefficients. When T g is less than 0 degrees Celsius (i.e., less than 273.15 Kelvin), the two take values of 21.87 and 7.66 respectively, and at other times take 17.269 and 35.86. PSFC 0 is obtained by terrain-weighted regression interpolation and STMAS correction from the previous text.

[0274]

[0275] q a is the specific humidity of the air at the lowest layer of the model grid, that is, q i,j,k when k = 1.; q i,j,k is calculated as shown in Equations (58)-(59).

[0276] The heat flux H flowing into the deep soil m is calculated through Equation (72); where T m is the deep soil temperature, which remains constant during the integration process. The initial field generation scheme is described later.

[0277] H m = K e C g (T g - T m ) (72)

[0278] To better reflect the anthropogenic heat factor of urban land use, the present invention adds an anthropogenic heat flux to the energy balance equation of the model. That is, Equation (60) is changed to the form of Equation (73):

[0279]

[0280] where H a represents the anthropogenic heat flux.

[0281] H a is calculated according to Equation (74):

[0282] H a = Anth * (1 - 0.6 cos(3.1415 * (hour - 5.5) / 12) (74)

[0283] where Anth represents the empirical coefficient of anthropogenic heat flux for different land uses (see Table 2), and hour represents the current simulation time, that is, to describe the daily variation characteristics of the anthropogenic heat flux.

[0284] Table 2 gives the land use classification and corresponding parameters used in the present invention. The specific parameters can be adjusted according to local conditions in actual applications. During actual simulations, for each grid point, the surface temperature, sensible heat flux, and latent heat flux of each land use type are calculated, and then the weighted average is calculated according to the proportion of different land use types in the grid point to obtain the final simulated surface temperature and other elements.

[0285] After calculating the surface temperature, surface sensible heat, and latent heat flux through the land surface model, the air temperature and humidity advection are simulated by integrating the temperature and humidity advection equations, and the land surface process effects are transferred to the atmosphere in the form of fluxes.

[0286] The temperature and humidity advection equation sets used in the present invention are as shown in Equations (75) and (76):

[0287]

[0288] In Equations (75)-(76), pti,j,k represents the potential temperature, calculated based on air temperature and atmospheric pressure, as shown in Equation (57); q i,j,k represents the specific humidity, calculated based on air temperature, relative humidity, and atmospheric pressure, as shown in Equations (58) and (59); u i,j,k , v i,j,k represent the north-south and east-west components of the horizontal wind speed; ω i,j,k represents the vertical wind speed component in the terrain-following coordinate system, calculated according to Equation (77):

[0289]

[0290] w i,j,k represents the vertical wind speed component; represents the correction coefficient in the terrain-following coordinate system, calculated according to Equations (78)-(80):

[0291]

[0292] In Equations (77)-(80), sigma k , hmax, zg i,j represent the model sigma layer value, the altitude of the highest layer of the model, and the terrain altitude, respectively.

[0293] In Equations (75)-(76), represents the horizontal diffusion coefficient, calculated according to Equations (81) and (82):

[0294]

[0295] In Equations (75)-(76), and represent the vertical diffusion coefficient of heat and the vertical diffusion coefficient of humidity, respectively. Since the model used in the present invention does not include a turbulence parameterization scheme and cannot calculate the turbulent energy and dissipation rate at high altitudes, the following Equations (83)-(84) are used for calculation and

[0296]

[0297] For model levels more than 200 meters above the ground, the vertical diffusion coefficient is set to 0; for vertical levels less than 200 meters above the ground, it is calculated according to Equations (83)-(84). The superscript k represents the vertical layer number; and are the surface sensible heat flux and latent heat flux, respectively. In this way, the simulation results of the land surface model can be input into the temperature and humidity advection process in the form of lower boundary fluxes.

[0298] In addition, when integrating the integrated temperature and humidity advection equation, the central difference is adopted for the spatial difference format of the model, the Euler format is used for time integration, and open boundary conditions are adopted for the boundary conditions except the lower boundary, that is, the variable gradient is 0 at the boundary. The integration step size and grid resolution are determined according to the actual situation.

[0299] Since time integration and spatial difference are inevitably accompanied by the growth of non-linear errors, the data assimilation method is used in the embodiments of the present application to constrain the model errors. The present invention uses the relaxation approximation (Nudging) as the assimilation algorithm.

[0300] The relaxation approximation (Nudging) makes each integration step approximate to the observation field by adding an approximation term to the model integration equation, and the longer the integration time, the better the assimilation effect. Nudging adds a relaxation forcing term to the model integration equation, so that the error correction acts on each integration step, and its assimilation process will not cause serious discontinuity problems; since the dynamic constraint is carried out by using the model integration, the physical balance between the variables of the model will not be destroyed. Taking the temperature advection equation as an example, the temperature advection equation with the Nudging relaxation forcing term added is shown in Equation (85):

[0301]

[0302] In the formula, τ is the Nudging assimilation time window; G i,j,k is the observation increment. G at each grid point of the model is calculated as shown in Equation (86), where the subscript i represents the observation data serial number, and N obs represents the number of site observations; represents the site potential temperature observation; h represents the observation operator for interpolating the grid point data to the site (linear interpolation is used in the present invention); pt represents the model grid point potential temperature.

[0303]

[0304] In Equation (86), W xyz,i is the spatial weight coefficient, which is calculated using Equation (87), where R represents the influence radius, which is set according to the type of observation data and time situation; d represents the distance from the grid point to the observation point.

[0305]

[0306] In Equation (86), W t,i is the time weight coefficient, which is calculated according to Equation (88). Where t represents the current model integration time; t o represents the observation time of the observation data; τ is the Nudging assimilation time window, and the unit is seconds.

[0307]

[0308] The above are all the technical points of the microscale model. During actual simulation, proceed as follows:

[0309] ① Based on the initial values of the fine-grid wind field and three-dimensional air temperature, relative humidity, air pressure, surface air pressure, etc., calculate model variables such as potential temperature and specific humidity. Among them, the initial value of the surface temperature is replaced by the 2-meter air temperature interpolated by PRISM; the soil temperature is calculated according to existing research. For example, see Liang L.L., D.A. Riveros-Iregui., R.E. Emanuel., et al. A simple framework to estimate distributed soil temperature from discrete air temperature measurements in data-scarce regions. 2013. J. Geophys. Res. Atmos., 119, 407–417.

[0310] ② Integrate the integrated land surface model for 5 minutes to initialize the surface temperature.

[0311] ③ Integrate the land surface model and the temperature and humidity advection equations simultaneously, and assimilate the meteorological station observation data.

[0312] ④ Based on the differences in potential temperature and specific humidity obtained from two adjacent integrations, determine whether the steady state is reached as the criterion for ending the integration; or use the integration duration determined empirically as the ending criterion.

[0313] ⑤ Convert the potential temperature and specific humidity obtained from the simulation into air temperature and relative humidity, and output the simulation results of the final three-dimensional air temperature, three-dimensional relative humidity, and three-dimensional wind field. Additionally, the simulation results of elements such as surface air pressure and surface temperature can also be output.

[0314] It should be noted that to achieve better simulation results, it is recommended to use 08:00 in the local time as the starting simulation time when conducting continuous hourly simulations: at this time, the near-surface air temperature, surface temperature, and soil temperature are close in value, and the model initialization effect is better. For other simulation time steps after the first time step, the initial values of the surface temperature and soil temperature can use the simulation results of the previous time step.

[0315] The present invention can quickly simulate microscale air temperature, humidity, and wind field elements with a resolution below 100 m in complex terrain areas. Taking the simulation case of the Chongqing Taicongyuan area as an example, the simulation grid spatial resolution is 30 m, the horizontal dimension is 267×333, the vertical dimension is 21, and a total of 15,000 time steps are simulated. For the parameter settings of this case, on a desktop computer with a windows operating system, using 20 threads for calculation, the total time consumption does not exceed 2 days. Figure 3-4The figure shows the multi-year average near-surface wind field and air temperature field simulated by building a model based on the present invention. It can be seen that the model can not only quickly generate grid data of microscale meteorological elements, but also better reflect the influence of different land use types and complex terrain on the elements.

[0316] This embodiment also provides a rapid simulation system for microscale meteorological elements in complex terrain areas.

[0317] Reference Figure 5 , a rapid simulation system for microscale meteorological elements in complex terrain areas, comprising:

[0318] A underlying surface information calculation module 51, configured to calculate the underlying surface information of the area to be simulated;

[0319] A complex terrain two-dimensional grid background field generation module 52, configured to generate grid meteorological element background fields of air temperature, air pressure, and relative humidity by using a terrain parameter weighted regression model based on the underlying surface information; generate a grid background field of wind direction and wind speed by using inverse distance weighted interpolation;

[0320] An error correction module 53, configured to perform preliminary error correction by using STMAS multi-grid variational constrained by terrain parameters based on the grid meteorological element background fields; correct the grid wind direction and wind speed background fields by using conventional STMAS; obtain the initial values of grid meteorological elements;

[0321] A three-dimensional meteorological element initial field generation module 54, configured to generate three-dimensional meteorological initial values by using the land use type for the initial values of grid meteorological elements; the three-dimensional meteorological initial values include an initial value of a three-dimensional wind field and initial values of other three-dimensional meteorological elements; the initial values of other three-dimensional meteorological elements include an initial value of three-dimensional air temperature, an initial value of three-dimensional relative humidity, and an initial value of three-dimensional air pressure;

[0322] A wind field diagnosis module 55, configured to input the initial value of the three-dimensional wind field into a wind field rapid diagnosis model to obtain a revised three-dimensional fine grid wind field;

[0323] A microscale model calculation module 56, based on the revised three-dimensional fine grid wind field and the initial values of other three-dimensional meteorological elements, simulates and obtains fine grid air temperature and humidity by integrating the microscale model.

[0324] The various change modes and specific examples in the method provided in the above embodiment are equally applicable to the system in this embodiment. Through the foregoing detailed description of the method, those skilled in the art can clearly know the implementation method of the system in this embodiment. For the sake of brevity of the specification, it will not be elaborated herein.

[0325] To better execute the program of the above method, an embodiment of the present application further provides a computer device, which includes a processor and a memory.

[0326] A computer device can be implemented in various forms, including devices such as computers and servers.

[0327] Among them, the memory can be used to store instructions, programs, codes, code sets or instruction sets. The memory can include a program storage area and a data storage area. Among them, the program storage area can store instructions for implementing the operating system, instructions for at least one function, and instructions for implementing the methods provided in the above embodiments, etc.; the data storage area can store data involved in the methods provided in the above embodiments, etc.

[0328] The processor can include one or more processing cores. By running or executing instructions, programs, code sets or instruction sets stored in the memory, the processor calls the data stored in the memory and executes various functions of this application and processes the data. The processor can be at least one of an Application Specific Integrated Circuit (ASIC), a Digital Signal Processor (DSP), a Digital Signal Processing Device (DSPD), a Programmable Logic Device (PLD), a Field Programmable Gate Array (FPGA), a Central Processing Unit (CPU), a controller, a microcontroller, and a microprocessor. It can be understood that for different devices, the electronic devices for implementing the above processor functions can also be others, and the embodiments of this application do not make specific limitations.

[0329] The embodiments of this application provide a computer-readable storage medium, for example, including various media that can store program codes such as USB flash drives, mobile hard disks, Read Only Memory (ROM), Random Access Memory (RAM), magnetic disks, or optical discs. The computer-readable storage medium stores a computer program that can be loaded and executed by the processor to execute the methods in the above embodiments.

[0330] The embodiments of this application also provide a computer program product, which includes a computer program tangibly contained on a computer-readable medium. The computer program includes program codes for executing any method in the embodiments of this application. The computer program can be downloaded and installed over a network and / or installed from a removable medium (such as a disk, an optical disc, a magneto-optical disc, a semiconductor memory, etc.).

[0331] As described above, the above embodiments are only used to introduce the technical solutions of the present application in detail. However, the description of the above embodiments is only used to help understand the method and its core idea of the present application, and should not be construed as a limitation of the present application. Any changes or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present application should be covered within the protection scope of the present application.

Claims

1. A method for rapid simulation of micro-scale meteorological elements in complex terrain areas, characterized in that: The method comprises: Calculate the underlying surface information of the area to be simulated; Based on the underlying surface information, the terrain parameter weighted regression model is used to generate the background field of temperature, pressure and relative humidity grid meteorological elements. Generate wind direction and wind speed grid background field using inverse distance weighted interpolation; Based on the background field of grid meteorological elements, STMAS multi-grid variation constrained by terrain parameters is used to make preliminary error correction; the background field of grid wind direction and speed is corrected using conventional STMAS; the initial values ​​of grid meteorological elements are obtained; Generate three-dimensional meteorological initial values ​​for the initial values ​​of the grid meteorological elements; the three-dimensional meteorological initial values ​​include three-dimensional wind field initial values ​​and other three-dimensional meteorological element initial values; the other three-dimensional meteorological element initial values ​​include three-dimensional temperature initial values, three-dimensional relative humidity initial values ​​and three-dimensional air pressure initial values; Input the initial value of the three-dimensional wind field into the wind field rapid diagnosis model to obtain the corrected three-dimensional fine grid wind field; Based on the corrected three-dimensional fine grid wind field and other three-dimensional meteorological element initial values, the fine grid air temperature and humidity are simulated through the integrated microscale model.

2. The method according to claim 1, characterized in that The terrain parameter weighted regression model includes establishing a univariate linear regression equation with terrain altitude as the independent variable and meteorological element values ​​as the dependent variable; calculating the comprehensive weight coefficient of each sample based on the difference in underlying surface information between the sample and the grid point; and solving the regression equation coefficient using the weighted least squares method and bringing it into the regression equation to calculate the final interpolation result.

3. The method according to claim 2, characterized in that The method further comprises: When the meteorological element is temperature, the comprehensive weight coefficient adopts the distance weight W based on the sample. d , height weight W z , aspect weight W f , vertical layer weight W l , terrain index weight W t and urban land weight W u Calculated; When the meteorological element is relative humidity, the comprehensive weight coefficient adopts the distance weight W based on the sample. d , height weight W z , aspect weight W f , effective terrain weight W e and urban land weight W u Calculated; When the meteorological element is air pressure, the comprehensive weight coefficient adopts the distance weight W of the sample. d , height weight W z and the effective terrain weight W e Calculated.

4. The method according to claim 3, characterized in that The method of generating the background field of the temperature, air pressure and relative humidity grid meteorological elements based on the underlying surface information and using the terrain parameter weighted regression model includes: For the temperature, the station temperature is used as a sample, and the grid temperature is interpolated using the terrain parameter weighted regression model; For relative humidity, the site temperature observation is used to convert the site relative humidity into the site dew point temperature using the following conversion formula: e=RH / 100·es Among them, a, b, c are constants, T a represents the station temperature (℃), RH represents the station relative humidity (%), e s represents the saturated water vapor pressure, e represents the water vapor pressure, T d Representative site dew point temperature (℃); Taking the station dew point temperature as a sample, the grid dew point temperature is interpolated using the terrain parameter weighted regression model; The grid dew point temperature and grid air temperature data are used to back-calculate the grid relative humidity using the conversion formula.

5. The method according to claim 4, characterized in that The method of generating the background field of the temperature, air pressure and relative humidity grid meteorological elements based on the underlying surface information and using the terrain parameter weighted regression model also includes: For air pressure, the site air pressure data is standardized based on the maximum and minimum ground air pressure in the sample to control the air pressure value in the range of 0 to 1; Convert the normalized air pressure data into logarithmic form; The logarithmic air pressure data are used as samples, and the grid pressure is interpolated using the terrain parameter weighted regression model.

6. The method according to claim 1, characterized in that The method of generating a wind direction and wind speed grid point background field by using inverse distance weighted interpolation comprises: The inverse distance weighted interpolation formula can be used to directly interpolate the grid wind speed based on the site wind speed; where x represents the sample value, the subscript i represents the sample number; y represents the interpolation result; d represents the distance from the sample to the predicted grid point; and p represents the gain coefficient. When interpolating wind direction, convert the site wind direction into standardized east-west and north-south components according to the following formula; where swspd represents the site wind speed, swdir represents the site wind direction, su represents the east-west component of the site wind speed, sv represents the north-south component of the site wind speed, and swdir Msin ,swdir Mcos They represent the standardized east-west component of the site wind direction and the standardized north-south component of the site wind direction, respectively; su=-swspd·sin(swdir) sv=-swspd·cos(swdir) <h2 style=";text-align:left;direction:ltr">swdir<h2 style=";text-align:left;direction:ltr"> Msin <h2 style=";text-align:left;direction:ltr"> (sin(tan)<h2 style=";text-align:left;direction:ltr"> -1 <h2 style=";text-align:left;direction:ltr"> (sv,su)) swdir Mcos =cos(tan -1 (sv,su)) Calculate swdir Msin ,swdir Mcos Then, the inverse distance weighted interpolation formula is used to interpolate the standardized east-west and north-south components of the grid points, and then the grid wind direction is calculated by the following formula, where wdir M For intermediate variables, wspd, wdir Msin 、wdir Mcos , u, and v represent the grid wind speed, grid standardized east-west component, grid standardized north-south component, grid standardized east-west wind speed, and grid standardized north-south wind speed, respectively; wdir M =tan -1 (wdir Msin ,wdir Mcos ) u=wspd·sin(wdir M ) v=wspd·cos(wdir M ) wdir=tan -1 (-v,-u)。 7. The method according to claim 1, characterized in that The STMAS multigrid variational method includes the following formula: Among them, x b represents the grid background field; y k Represents the site observation increment used in multi-grid solution, y o represents the site observation value, x k represents the grid point background field correction increment that needs to be calculated at each layer during multi-grid solution, x represents the final corrected background field, and h k represents the observation operator that interpolates the grid background field to the station coordinates. The superscript k represents the different levels of the multigrid. kmax represents the number of multigrid levels. o represents the empirical ratio of the observation error covariance to the background field error covariance. T represents the transpose. J k is the name of the cost function of the kth layer.

8. The method according to any one of claims 1 to 7, characterized in that: The underlying surface information includes at least one of slope aspect, effective terrain index, vertical stratification, and terrain index.

9. The method according to claim 8, characterized in that The calculation of the underlying surface information of the area to be simulated includes: For the slope aspect, the following formula is used to calculate: Where asp represents the slope aspect, dz represents the elevation difference between two grid points, dy represents the vertical distance between two grid points, and dx represents the horizontal distance between two grid points. For the effective terrain index, the minimum elevation value in the circular domain with a radius of 15km around each grid point is counted to obtain the minimum elevation data, and the average value of the circular domain with a radius of 15km is calculated for each grid point on the minimum elevation data to achieve data smoothing; Subtract the smoothed minimum elevation data from the original elevation data to obtain the initial effective terrain height; For each grid point of the initial effective terrain height, the grid point average value is calculated and smoothed according to a circular domain with a radius of 6 km to obtain the final smoothed effective terrain height; The index I was calculated according to the following formula: 3c ; where h c is the effective terrain height after final smoothing; h2 and h3 are the 2D terrain and 3D terrain thresholds, which are set to 100m and 600m respectively; The index I was calculated according to the following formula: 3a ; where h a It is the average effective terrain height within 8 km around the grid point, calculated with the inverse distance as the weight; The effective terrain index I is calculated according to the following formula 3d ; I 3d =max(I 3c ,I 3a ) For vertical stratification, the minimum altitude within a horizontal radius of 10 km of the grid points is counted; Subtract the minimum altitude from the current grid point altitude. If the difference is more than 150m, the grid point is at level 2; otherwise, the grid point is at level 1. For the terrain index, the minimum terrain altitude within a horizontal radius of 10 km of the grid point is counted; Gaussian smoothing is performed on the minimum terrain elevation grid point with a radius of 10 km; The terrain index data is obtained by subtracting the minimum terrain altitude from the original terrain altitude of the grid point.

10. A rapid simulation system for micro-scale meteorological elements in complex terrain areas, characterized by: include: An underlying surface information calculation module is used to calculate the underlying surface information of the area to be simulated; The complex terrain two-dimensional grid background field generation module is used to generate the grid meteorological element background field of temperature, air pressure and relative humidity based on the underlying surface information and the terrain parameter weighted regression model; the wind direction and wind speed grid background field is generated using the inverse distance weighted interpolation; the error correction module is used to perform preliminary error correction based on the grid meteorological element background field using the STMAS multi-grid variation constrained by terrain parameters; the grid wind direction and wind speed background field is corrected using conventional STMAS; Get the initial values ​​of meteorological elements at two-dimensional grid points; A three-dimensional meteorological element initial field generation module is used to generate three-dimensional meteorological initial values ​​for grid meteorological element initial values; the three-dimensional meteorological initial values ​​include three-dimensional wind field initial values ​​and other three-dimensional meteorological element initial values; the other three-dimensional meteorological element initial values ​​include three-dimensional temperature initial values, three-dimensional relative humidity initial values ​​and three-dimensional air pressure initial values; The wind field diagnosis module is used to input the initial value of the three-dimensional wind field into the wind field rapid diagnosis model to obtain the corrected three-dimensional fine grid wind field; The microscale model calculation module is used to simulate the fine grid air temperature and humidity based on the revised initial values ​​of the three-dimensional fine grid wind field and other three-dimensional meteorological elements through the integrated microscale model.

Images