Rapid simulation method and system for micro-scale meteorological elements in complex terrain area

Abstract

The invention relates to a quick simulation method and system for micro-scale meteorological elements in a complex terrain area. The method comprises the following steps: calculating underlying surface information of a to-be-simulated area; based on the underlying surface information, utilizing a terrain parameter weighted regression model to generate an air temperature, air pressure and relative humidity grid point meteorological element background field; generating a wind direction and wind speed grid point background field by using inverse distance weight interpolation; based on the grid point meteorological element background field, performing preliminary error correction by using STMAS multi-grid variation to obtain a grid point meteorological element initial value; generating a three-dimensional meteorological initial value according to the grid point meteorological element initial value; inputting the three-dimensional wind field initial value into the wind field rapid diagnosis model to obtain a corrected three-dimensional fine grid point wind field; and on the basis of the corrected three-dimensional fine grid point wind field and other three-dimensional meteorological element initial values, simulating to obtain fine grid point air temperature and humidity through an integral micro-scale mode. The problems that a current meteorological element rapid simulation method is not suitable for a complex terrain region and the simulation error is large are solved.

Description

Technical Field

[0001] The present application relates to the field of meteorological monitoring technology, and in particular to a method and system for rapid simulation of micro-scale meteorological elements in complex terrain areas. Background Art

[0002] Compared with plain areas, complex terrain areas (such as mountains and hills) are characterized by large terrain undulations, strong localization and suddenness of disasters and secondary disasters, which puts higher demands on the refinement and accuracy of local meteorological services and disaster prevention and mitigation work.

[0003] Meteorological stations in complex terrain areas are unevenly distributed and local climate characteristics are difficult to distinguish. To improve the quality of meteorological services and disaster prevention and mitigation, the only way is to convert station data into fine grid data based on mathematical models. There are many methods for making grid meteorological elements based on station meteorological observation data, among which spatial interpolation of station meteorological observations, traditional numerical simulation and fast simulation are three commonly used methods, and all three methods have their own advantages and disadvantages.

[0004] Spatial interpolation has the advantage of high computational efficiency, and some models that take terrain factors into account are also well applicable in complex terrain areas. However, there is still great uncertainty in the application of interpolation algorithms in meteorological services and disaster prevention and mitigation work in complex terrain areas: for meteorological elements such as wind fields and humidity that have strong spatial discontinuities and are highly dependent on atmospheric physical dynamic processes, it is difficult to obtain ideal results, whether it is traditional interpolation algorithms or interpolation algorithms that take terrain into account; at the same time, affected by the characteristics of the algorithm, it is also difficult to reflect the fine land use types and their changing characteristics through interpolation, which further limits the application scenarios of interpolation algorithms.

[0005] Traditional numerical simulation is based on mature atmospheric dynamics, thermodynamic equations and physical parameterization schemes, and can be used with data assimilation and other means to obtain grid meteorological elements with high temporal and spatial resolution and low errors. However, the technical threshold of traditional numerical simulation is high and the calculation time is long. In actual disaster prevention and mitigation or service scenarios, business departments are often more concerned about how to produce refined data products as quickly and conveniently as possible while ensuring a certain accuracy rate, so the application of traditional numerical models is restricted. At the same time, the calculation stability of the model will face huge challenges when conducting high temporal and spatial resolution simulations in complex terrain areas. In addition, the systematic errors of the model will also bring a lot of uncertainty-these factors limit the application of traditional numerical simulations 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: 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. Use inverse distance weighted interpolation to generate the wind direction and wind speed grid background field; Based on the background field of grid meteorological elements, the STMAS multi-grid variation constrained by terrain parameters is used to make preliminary error corrections; the grid wind direction and wind speed background field 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.

[0009] By adopting the above technical solutions, a set of rapid simulation technologies for micro-scale meteorological elements that are more suitable for complex terrain areas was designed. It is used to quickly generate fine grid data of common meteorological elements such as temperature, humidity and wind field in complex terrain areas, solving the problem that the current rapid simulation methods of meteorological elements are not applicable to complex terrain areas and have large simulation errors.

[0010] Optionally, the terrain parameter weighted regression model includes establishing a univariate 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 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 substituting it into the regression equation to calculate the final interpolation result.

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

[0012] Optionally, 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.

[0013] By adopting the above technical solution and using different total weight coefficients for different meteorological elements, the interpolation accuracy of the weighted regression model of terrain parameters can be improved and the simulation effect can be improved.

[0014] Optionally, the generating of the background field of the temperature, air pressure and relative humidity grid meteorological elements by using the weighted regression model of terrain parameters based on the underlying surface information 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.

[0015] Optionally, the method of generating a background field of temperature, air pressure and relative humidity grid meteorological elements based on underlying surface information by using a terrain parameter weighted regression model further 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.

[0016] Optionally, the generating of the wind direction and wind speed grid point background field by using inverse distance weighted interpolation includes: 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 (m / s), 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) swdir Msin = sin(tan -1 (svv,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).

[0017] Optionally, the STMAS multi-grid variation includes using 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, which is usually between 0 and 1. T represents the transposition. J k is the name of the cost function of the kth layer.

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

[0019] 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 formula3d ; 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.

[0020] In the second aspect, the present application provides a system for rapid simulation of micro-scale meteorological elements in complex terrain areas, including: 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, 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 by 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; and the initial value of the grid meteorological element is obtained; 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 corrected three-dimensional fine grid wind field and other three-dimensional meteorological element initial values ​​through the integrated microscale model.

[0021] In summary, this application at least includes the following beneficial technical effects: 1. It solves the problem 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

[0022] Figure 1A schematic diagram of a process for rapid simulation of micro-scale meteorological elements in a complex terrain area in an embodiment of the present application; Figure 2 A flowchart of a method for rapid simulation of micro-scale meteorological elements in a complex terrain area in an embodiment of the present application; Figure 3 This is a schematic diagram of the multi-year average 10m height wind field of the Tainong Garden simulated in the embodiments of the present application; Figure 4 This is a schematic diagram of the average temperature at 2m above the Tainong Garden over many years simulated in the embodiments of this application; Figure 5 This is a schematic diagram of the structure of a system for rapid simulation of micro-scale meteorological elements in complex terrain areas in an embodiment of the present application. DETAILED DESCRIPTION

[0023] In order to make the purpose, technical solution and advantages of the present application more clearly understood, the present application is further described in detail below in conjunction with 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.

[0024] 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 be used as limitations to the present application. As used in the specification and the appended claims of the present application, the singular expressions "one", "a kind of", "said", "above", "the" and "this" are intended to also include plural expressions, unless there is a clear contrary indication in its context. 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 listed items. The term "exemplary" means "used as an example, embodiment or illustrative", and any embodiment described as "exemplary" here does not need to be interpreted as being superior to or better than other embodiments. The terms "first" and "second" are only used for descriptive purposes, and cannot be understood as implying or implying relative importance or implicitly indicating the number of indicated technical features. Thus, the features defined as "first" and "second" can expressly or implicitly include one or more of these features, and in the description of the embodiments of the present application, unless otherwise specified, the meaning of "multiple" is two or more.

[0025] refer to Figure 1 , a method for rapid simulation of micro-scale meteorological elements in complex terrain areas, comprising the following steps: S101: Calculate the underlying surface information of the area to be simulated; S102: Based on the underlying surface information, a terrain parameter weighted regression model is used to generate a background field of meteorological elements of air temperature, air pressure, and relative humidity grid points; S103: Generate a wind direction and wind speed grid background field using inverse distance weighted interpolation; S104: Based on the grid meteorological element background field, preliminary error correction is performed using STMAS multi-grid variation constrained by terrain parameters; S105: Correct the grid wind direction and speed background field using conventional STMAS; S106: generating initial values ​​of meteorological elements of three-dimensional temperature, three-dimensional relative humidity, three-dimensional air pressure and three-dimensional wind field for the grid meteorological element initial values ​​obtained after correction; S107: inputting the initial value of the three-dimensional wind field into the wind field rapid diagnosis model to obtain a corrected three-dimensional fine grid wind field; S108: Based on the corrected three-dimensional fine grid wind field and the three-dimensional temperature, three-dimensional relative humidity, and three-dimensional pressure meteorological initial values, the fine grid temperature and humidity are simulated through the integrated micro-scale model.

[0026] The embodiment of the present application designs a method for rapid simulation of micro-scale meteorological elements that is more suitable for application in complex terrain areas. The 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 rapid simulation of common meteorological elements such as micro-scale temperature, humidity, and wind field in complex terrain areas. It is not only a perfection of rapid simulation methods and methods for developing 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.

[0027] See also Figure 2 The following is a detailed description of the method for rapid simulation of micro-scale meteorological elements in complex terrain areas provided by this embodiment: (1) Collect and organize multi-source data such as site 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.

[0028] Among them, the station meteorological data includes 2-meter temperature, 2-meter relative humidity, ground pressure, and 10-meter wind direction and speed; geographic information data includes digital terrain elevation (hereinafter referred to as DEM) and land use type data. Land use type refers to the classification of land into different categories based on the natural and economic characteristics of the land for reasonable use and management. Common land use types include cultivated land, gardens, forests, grasslands, waters, urban and rural settlements, industrial and mining land, and transportation land.

[0029] Slope aspect refers to the inclination direction of a point on the surface of the earth. 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 north direction is 0°, and increases clockwise. The maximum value is 359°59'59', and the unit is degree. The calculation formula is as follows (1), 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.

[0030] The effective terrain index is a standardized index based on DEM calculation to determine the complexity of a region's terrain. The effective terrain index value ranges from 0 to 1. The larger the value, the more complex the terrain at the grid point. The calculation process of the effective terrain is as follows: ① Count the minimum elevation value in the 15km radius circle around each grid point in the DEM to obtain the minimum elevation data, and calculate the average value of the 15km radius circle for each grid point on the minimum elevation data to achieve data smoothing; ② Subtract the smoothed minimum elevation data from the original DEM to obtain the initial effective terrain height; ③ For each grid point of the initial effective terrain height, the grid point average value is calculated according to a circular domain with a radius of 6 km and smoothed to obtain the final smoothed effective terrain height.

[0031] ④ Calculate the index I according to formula (2) 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.

[0032] ⑤ Calculate the index I according to formula (3) 3a . where h a It is the average effective terrain height within 8km around the grid point, calculated with the inverse distance as the weight. h2 and h3 are the same as ④.

[0033] ⑥Finally, calculate the effective terrain index I according to formula (4): 3d .

[0034] I 3d =max(I 3c ,I 3a ) (4) Since inversion layers (i.e., the temperature increases with increasing altitude) often exist in complex terrain areas, it is necessary to consider the differences in temperature variability at different altitudes when subsequently generating the temperature background field, so the vertical stratification data must be calculated in advance.

[0035] Methods for calculating vertical stratification include: ① Count the minimum altitude within a horizontal radius of 10km at a certain grid point; ② 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; ③ Repeat ① and ② at each grid point to obtain vertically stratified grid data.

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

[0037] The terrain index is calculated as follows: ① Count the minimum terrain altitude within a horizontal radius of 10km at a certain grid point; ② Perform Gaussian smoothing on the minimum terrain elevation grid point with a radius of 10 km; ③ Subtract the minimum terrain altitude from the original terrain altitude of the grid point to obtain the terrain index data.

[0038] (2) Based on the various underlying surface information obtained from the above calculations, the terrain parameter weighted regression model is used to generate the grid meteorological element background field of temperature, air pressure, and relative humidity; the inverse distance weighted interpolation is used to generate the grid background field of wind direction and wind speed.

[0039] The terrain parameter weighted regression model used in the present invention is improved based on the PRISM (Parameter-elevationRegressions 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 altitude, and the meteorological elements are approximately linearly related under the combined influence of various geographical factors such as altitude. The PRISM model uses terrain altitude as the independent variable and meteorological element values ​​as the dependent variable to establish a univariate linear regression equation; calculates the comprehensive weight coefficient of each sample based on the difference in underlying surface information between the sample and the grid point; uses the weighted least squares method to solve the regression equation coefficient and brings it into the regression equation to calculate the final interpolation result.

[0040] The PRISM regression equation can be specifically formulated as follows: Y=β1X+β0 (5) 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; The regression coefficients β1 and β0 can be calculated by the following equations: In the formula, I 3d is the effective terrain index, x i Represents the altitude of each sample location, y i Represents the meteorological element value of each sample, wi Represents the comprehensive weight coefficient of each sample, represents the weighted average of sample altitudes, Represents the weighted average of sample meteorological element values.

[0041] Weighted average of sample altitude It can be calculated by the following formula: Weighted average of sample meteorological elements It can be calculated by the following formula: Comprehensive weight coefficient w i Sample-based distance weight W d , height weight W z , aspect weight W f , vertical layer weight W l , effective terrain weight W e , terrain index weight W t Specifically, the comprehensive weight coefficient w i The calculation formula is: In the formula, F d and F z They are distance weights W d and height weight W z The respective proportions, F d With F z The sum is equal to 1. Among them, F d With F z The value of can be flexibly determined based on experience according to the type of meteorological elements and local terrain characteristics; the commonly used setting is F d =0.7, F z =0.3.

[0042] The distance weight W of the sample d It can be calculated by the following formula: 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.

[0043] When the distance d from the sample point to the interpolation grid point is less than or equal to the influence radius r m , then the similarity 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 larger than the influence radius r mAs the distance increases, the similarity of meteorological elements between two points will gradually decrease. This form can better avoid the "bull's eye" phenomenon in the interpolation results (referring to the circle phenomenon with the interpolation point as the center formed by some large or small data during the interpolation process).

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

[0045] The height weight W of the sample z Calculated by the following formula: Where Δz is the absolute value of the height difference between the sample point and the interpolation grid point; Δz m and Δz x They represent the minimum height difference and the maximum height difference between each sample point and the interpolation grid point respectively; b is the gain coefficient, which is usually taken as 1, that is, the linear effect of height change on feature similarity is considered; I 3d is the effective terrain index, that is, as the terrain of the interpolation target grid point gradually becomes flat, the height weights of each sample will tend to be consistent. The calculation of the effective terrain index is as shown in (4). The similarity of meteorological elements between two points is also closely related to the height difference. When the height difference is the minimum, the similarity of meteorological elements between the two points is the highest; when the height difference exceeds a certain level, the correlation gradually decreases or even becomes completely irrelevant.

[0046] The aspect weight W of the sample f Calculated by the following formula: In the formula, Δf represents the slope difference, c is the gain coefficient, which is usually 1; I 3d is the effective terrain index, that is, as the terrain of the interpolation target grid point becomes gradually flat, the slope weights of each sample will tend to be consistent. The calculation of the effective terrain index is as shown in (4). The meteorological elements of two spatial points with similar slopes should be more similar.

[0047] Considering that inversion layers are easily formed in complex terrain areas, different vertical temperature gradients may exist at different heights above the ground. Therefore, when interpolating a specific interpolation grid point, samples in the same vertical layer as the grid point should be given a higher weight.

[0048] The vertical stratification weight W of the sample l Calculated by the following formula: Where Δl is the absolute value of the stratified difference between the sample point and the interpolation grid point, c is the gain coefficient; Δz m and Δz x They represent the minimum and maximum height differences between each sample point and the interpolation grid point, respectively.

[0049] 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 valleys are more susceptible to cold air, and samples in similar areas should be given higher weights. The terrain index weight of the sample W t Calculated by the following formula: Where Δt represents the absolute value of the difference between the terrain index of 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 Take 0.1km, Δt x Take 0.5km; z represents the gain index, usually 2.

[0050] Taking into account the representativeness of samples in complex terrain areas and flat terrain areas, the effective terrain weight W is used for each grid point. e Weight or reduce the weight of different types of samples. The effective terrain weight formula is as follows: Among them, I 3ds and I 3dg 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 shown in the previous text.

[0051] In order to better consider the urban climate effect in the background field, the present invention further adds the urban land weight W in the comprehensive weight coefficient. u , at this time, the comprehensive weight coefficient is calculated using the following formula: Urban land weight W u , calculated by the following formula: 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 weight reduction coefficient set based on 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, the proportion of urban meteorological observation stations in the sample is increased. By introducing the urban land weight W u, and controlling the sample ratio of urban sites, the impact of urban land use on meteorological elements can be better considered in interpolation.

[0052] In addition, since relative humidity and ground pressure are not linearly related to altitude, directly using a linear regression model is not applicable. Therefore, the terrain weighted regression model used in the present invention performs linear processing on the samples.

[0053] For relative humidity, the interpolation process is as follows: ① Use the station temperature observation to convert relative humidity into dew point temperature. The conversion formula is shown in equations (21)-(23). In equation (21), a′, b′, and c′ are constants, with values ​​of 611.21, 17.502, and 240.97, respectively; T a represents temperature (℃), RH represents relative humidity (%), e s represents the saturated water vapor pressure, e represents the water vapor pressure, T d Represents the dew point temperature (℃).

[0054] e=RH / 100·es (22) ② Taking the station dew point temperature as the sample, the terrain weighted regression method is used to interpolate the grid dew point temperature; taking the station air temperature as the sample, the terrain weighted regression method is used to interpolate the grid air temperature.

[0055] ③ Using the grid dew point temperature data and grid air temperature data, refer to equations (21)-(23) to calculate the grid relative humidity.

[0056] For surface pressure, the interpolation process is as follows: ① Count the maximum and minimum ground pressure in the sample, and standardize it using formula (24) to control the pressure value in the range of 0 to 1. norm Represents the standardized air pressure sample, y min Represents the minimum sample pressure, y max Represents the maximum sample pressure.

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

[0058] ③Based on the logarithmic samples, the terrain weighted regression method is used to interpolate the grid data.

[0059] ④ According to equations (24) and (25), the grid point data obtained in ③ are converted back to the surface pressure through back calculation to complete the surface pressure interpolation.

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

[0061] Table 1 Weight coefficients used when interpolating background fields of different element grids Since the relationship between wind direction and wind speed and altitude is relatively complex, it is not possible to use terrain parameter weighted regression interpolation to obtain the grid background field. The embodiment of the present application uses the inverse distance weighted interpolation method. The inverse distance weighted interpolation formula is shown in (26), 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; p represents the gain coefficient, which is usually 2.

[0062] Before interpolating the site wind direction to the grid point, equations (27)-(30) must be used to convert the site wind direction into standardized east-west and north-south components. Where swspd represents the site wind speed (m / s), swdir represents the site wind direction (radians), su represents the site wind speed east-west component (m / s), and sv represents the site wind speed north-south component (m / s). swdir Msin ,swdir Mcos represent the site-normalized east-west and north-south components, respectively.

[0063] su=-swspd·sin(swdir) (27) sv=-swspd·cos(swdir) (28) swdir Msin = sin(tan -1 (sv,su)) (29) swdir Mcos = cos(tan -1 (sv,su)) (30) Calculate swdir Msin ,swdir Mcos After that, change swspd and swdir Msin ,swdir Mcos Substitute x into equation (26) and interpolate to obtain the grid-standardized east-west and north-south components. Then, the grid wind direction is obtained by inverse calculation of equations (31)-(34), 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, gridded east-west wind speed, and gridded north-south wind speed, respectively.

[0064] wdir M =tan -1 (wdir Msin ,wdir Mcos ) (31) u=wspd·sin(wdir M ) (32) v = wspd.cos(wdir M ) (33) wdir=tan -1 (-v,-u) (34) Through the above algorithm steps, the present invention can interpolate the gridded temperature, air pressure, relative humidity, wind direction and wind speed background field of the complex terrain area based on the site meteorological data and geographic information data.

[0065] (3) Based on the background field of grid meteorological elements, the STMAS multi-grid variational method constrained by terrain parameters is used to perform preliminary error correction to obtain the initial values ​​of grid meteorological elements.

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

[0067] 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 to be set, which is determined according to the actual situation of the grid dimension. O represents the empirical ratio of the observation error covariance to the background field error covariance, which is usually between 0 and 1. T represents the transpose. J k is the name of the cost function of the kth layer.

[0068] When solving the STMAS equation, the grid background field that needs to be corrected is first interpolated to a 3×3 dimensional grid. Then, the variables in equation (35) are calculated using equation (36). The quasi-Newton method is used to solve equation (37) to obtain x 1 Then, the background field is interpolated to a 6×6 dimensional grid, and the variables in equation (35) are calculated using equation (36), and equation (37) is also solved to obtain x 2 In this way, the grid dimension of the background field interpolation is increased step by step by a factor of 2, and the corresponding level background field increment x is obtained by solving k , until the dimension is the same as the original dimension of the background field. Finally, the corrected background field is calculated using equation (38).

[0069] It should be understood that due to the observation operator h used in STMAS operation k Bilinear interpolation is used, so it is easy to cause local "bull's eye" linearity in complex terrain areas. k It is replaced by the terrain parameter weighted regression model mentioned above, thus realizing the STMAS multi-grid variational method constrained by terrain parameters.

[0070] In this way, the grid temperature, relative humidity, and air pressure background fields obtained by interpolating the terrain parameter weighted regression model are corrected using the terrain parameter weight constraints STMAS and the corresponding station observation data; the grid wind direction and wind speed background fields are corrected using conventional STMAS; and the initial values ​​of the ground grid meteorological elements required subsequently are obtained.

[0071] (4) Using land use types and empirical parameterization schemes, generate initial meteorological values ​​of three-dimensional temperature, three-dimensional humidity, three-dimensional air pressure, and three-dimensional wind field.

[0072] First, according to actual needs, set the vertical grid and calculate the altitude of each grid layer. The vertical grid used in the present invention adopts the sigma terrain following coordinate system, and the height of each grid point of the three-dimensional grid is calculated as follows: 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 can be set according to actual needs (for example, the lowest layer can be set to 5, the altitude of the highest layer can be set to hmax, and the values ​​of the intermediate layers can 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 represent the grid point indexes in the east-west, north-south, and vertical directions of the three-dimensional grid, respectively.

[0073] For the temperature, the vertical temperature lapse rate (i.e., β1) obtained by weighted regression of the initial grid point temperature with the terrain parameters is first interpolated to the bottom layer of the grid, and then the initial three-dimensional temperature of each grid point is obtained by equation (40): t i,j,k =gamma i，j，k ·(zp i,j,k -zp i，j，k-1 )+t i,j,k-1 (40) Among them, gamma represents the vertical decrease rate of the current grid temperature, and t represents the temperature.

[0074] The vertical temperature lapse rate gamma is calculated as follows: Where h is the height set according to experience, which can usually be set to 500 meters; gamma0 is the artificially set empirical temperature lapse rate, which is usually -0.006. Formula (41) actually means that when the height above the ground is less than h, the temperature lapse rate of each layer is linearly transformed to gamma0 with the height; when the height above the ground is greater than h, the temperature lapse rate remains unchanged at gamma0.

[0075] For relative humidity, first refer to equations (21)-(23) to convert it into dew point temperature, and then generate the three-dimensional grid dew point temperature in a similar way to the air temperature, except that gamma0 is usually set to -0.003. After obtaining the three-dimensional dew point temperature, inversely calculate equations (21)-(23) to convert the three-dimensional dew point temperature into three-dimensional relative humidity.

[0076] For air pressure, the initial three-dimensional air pressure is generated according to equations (42)-(44):

[0077] Among them, the values ​​of a1, a2, a3 and a4 are -3.9082017e-2, -1.1526465e-3, 3.2891937e-5, and -2.0494958e-7 respectively; the values ​​of b1, b2, b3 and b4 are -4.9244637e-3, -1.2984142e-6, -1.5701595e-6, and 1.5535974e-8 respectively. i,j,k , B i,j,k is an intermediate variable.

[0078] After obtaining the initial three-dimensional air pressure, replace zp in (42)-(44) with the terrain altitude zg to calculate the estimated ground pressure Formula (45) uses terrain weighted regression interpolation to calculate the surface air pressure Initial 3D pressure Make corrections.

[0079] In the formula, k represents the vertical grid index. For the horizontal wind speed, it is calculated according to formula (46) according to the type of underlying surface.

[0080] Among them, wspd i,j,k Represents the final generated three-dimensional wind speed initial value; represents the initial value of the two-dimensional wind speed near the ground; h0 represents the height of the initial wind speed from the ground (usually 10 meters); r0 represents the rough length of the ground type; h mix represents the distance from the top of the mixed layer to the ground, which can be empirically set to 800 meters; is the altitude of the sounding data (the altitude of the sounding data at 700hPa can be used); is the wind speed at the height of the mixing layer; represents the wind speed at the sounding observation height (the 700hPa sounding observation wind speed can be used); the parameters of linear regression fitting. Formula (46) represents that the wind speed changes at a logarithmic rate within the mixed layer, and the wind speed changes linearly to the wind speed at the sounding observation height above the mixed layer; when the height exceeds the sounding height, the wind speed remains unchanged.

[0081] The specific land use type classification and corresponding rough length are explained later.

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

[0083] Among them, wdir i,j,k Represents the final generated initial value of the three-dimensional wind direction; Represents the initial value of the near-surface two-dimensional wind speed; represents the wind direction at the sounding height, and the other variables are the same as (46). Formula (47) means that the wind direction remains unchanged within the mixed layer, and the wind direction gradually changes to the sounding height wind direction as the mixed layer goes up; when the height exceeds the sounding height, the wind direction remains unchanged.

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

[0085] After generating the initial three-dimensional temperature, three-dimensional air pressure, three-dimensional relative humidity, horizontal wind speed, horizontal wind direction and vertical wind speed initial values ​​according to the above steps, a parameterization scheme is used to superimpose the correction coefficient in the initial wind field, taking into account the influence of the complex terrain slope flow (i.e., the thermal circulation of the valley wind). Since the parameterization scheme is relatively complex and is not the focus of the present invention, it will not be repeated here. For example, reference can be made to "Scire JS, Robe FRFine-scale application of the CALMET meteorological model to a complex terrain site. 1997. In'Air&WasteManagement Association's 90th Annual Meeting&Exhibition'. pp. 1-16. (Toronto, Ontario, Canada)".

[0086] (5) The initial value of the three-dimensional wind field is input into the wind field rapid diagnosis model, and the corrected three-dimensional fine grid wind field is obtained after adjustment based on the mass conservation law.

[0087] The calculation process of the wind field rapid diagnosis model used in the present invention is as follows: ① Use equations (48)-(50) to superimpose the terrain following coordinate system correction coefficient on the initial three-dimensional wind field. are the components of the initial three-dimensional wind field in the east-west, north-south and vertical directions; is the wind field component after coordinate system transformation.

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

[0089] In the formula, x, y, and z represent the east-west, north-south, and vertical grids respectively; and are the horizontal and vertical scaling factors, respectively.

[0090] It should be noted that Controls the strength of the model's adjustment of the vertical and horizontal wind components during diagnosis: When it approaches infinity, the model only adjusts the vertical wind component; when When it approaches 0, only the horizontal wind component is adjusted.

[0091] To adaptively calculate A stability parameterization scheme is used here. The equations of the parameterization scheme are as follows (52) and (53): Among them, pt i,j,k represents the three-dimensional potential temperature, and its calculation formula is shown in the following formula (57); wspd i,j,k represents the three-dimensional wind speed, and max represents the maximum value of the two, that is, the wind speed during calculation is not less than 0.2m / s; represents the absolute difference between the grid altitude and the highest terrain altitude within a 10km radius; g represents the acceleration of gravity, which is 9.8; F i,j,k is the Strouhal number. Finally, we calculate The maximum value does not exceed 30, and the minimum value does not fall below 0.033. Automatically calculated After that, just give Just assign a value, usually a simple value of 1.

[0092] ③ Solve to get λ i,j,k Then, we substitute into equations (54)-(56) to calculate the wind component corrected by the terrain following coordinate system and satisfying the mass conservation constraint: Finally, Substituting the left side into equations (48)-(50) and inversely calculating the corresponding wind component on the right side, we can obtain the corrected three-dimensional fine grid wind field.

[0093] (6) Based on the three-dimensional fine grid wind field, the initial values ​​of potential temperature, specific humidity and air pressure are input, and by integrating a microscale model, the three-dimensional air temperature and humidity that takes into account the influence of fine land use are simulated.

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

[0095] The initial value of relative humidity is calculated using the initial values ​​of air temperature, air pressure and relative humidity, as shown in formulas (58)-(59). i,j,k is the initial value of three-dimensional relative humidity; esat i,j,k is the saturated water vapor pressure; q i,j,k For specific humidity.

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

[0097] The land surface process parameterization model used in the present invention considers a thin layer of soil on the ground and uses the forced recovery method to calculate the surface temperature. The energy balance equation of the model is: Among them, T g is the surface temperature (the initial field generation scheme is described later), t is the integration step of the land surface process model, R n is the net radiation, H e and H s are ground latent heat flux and ground sensible heat flux, respectively, H m is the heat flux 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 heat capacity, K e is the heat transfer coefficient. K e and C s The size of is determined by the type of surface vegetation cover (see Table 2).

[0098] Table 2 Land use types and corresponding thermodynamic parameters used in the present invention In formula (60), the net radiation R n It can be calculated as follows (61): R n =R s +I ↑ +I ↓ (61) Where R s is the net shortwave radiation flux, calculated as follows: R s =S0(1-A)D se p m sin(s h ) (62) D se =1.000110+0.034221*cos(d0)+0.001280*sin(d0)+0.000719*cos(2*d0)+ 0.000077*sin(2*d0)(64) d0=(noday-1)*2*3.14159 / 365 (65) Among them, S0 is the solar constant (valued at 1370.0); A is the surface shortwave albedo (see Table 2); D se is the distance factor between the sun and the earth, p is the atmospheric transparency coefficient, m is the atmospheric optical quality, d0 represents the arc coordinate of the earth relative to the sun calculated based on the current date, noday represents the total number of days in the current date in the whole year, int represents integer, and glat represents latitude. h and θ are the solar altitude angle and the zenith angle respectively. The specific calculation method can adopt the existing method, for example, refer to "Xu Yumao, Liu Hongnian, Xu Guiyu. 2000. Introduction to Atmospheric Science. Nanjing: Nanjing University Press.".

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

[0100] In formula (61), I ↑ is the upward longwave radiation flux from the ground, I ↓ is the downward longwave radiation flux of the atmosphere, calculated as follows (67)-(68): Where, T a is the surface air temperature, T g is the surface temperature; ε a , ε g are the infrared emissivity of the atmosphere (valued at 0.725) and the surface (see Table 2). σ is the Boltzmann constant, valued at 5.67e －8 .

[0101] Ground sensible heat flux H s According to formula (69), H s =ρ a C p C H |V a |(T g -T a ) (69) Among them, ρ a is the surface air density (calculated using ground pressure, relative humidity and temperature), Cp is the constant pressure specific heat of air (valued at 1004.67), C H is the thermal drag coefficient (can be 2.5e-3), V a is the surface wind speed.

[0102] Ground latent heat flux H e Calculate according to formula (70): H e =ηρ a LC H |V a |(q sg -q a ) (70) Where η represents the surface saturation (see Table 2), L is the heat of vaporization of liquid water (valued at 2.5e6), and q sg Based on the surface temperature (T g ), the surface air pressure (PSFC0) is calculated as the ground saturated specific humidity, as shown in formula (71). Where a″, b″ are state coefficients. When T g When the temperature is less than 0 degrees Celsius (i.e. less than 273.15 degrees Kelvin), the two values ​​are 21.87 and 7.66 respectively, and at other times they are 17.269 and 35.86. PSFC0 is obtained by terrain weighted regression interpolation and STMAS correction.

[0103] q a is based on the specific humidity of the air at the lowest level of the model grid, i.e., q when k = 1 i,j,k ;q i,j,k The calculation of is as shown in equations (58)-(59).

[0104] Heat flux H flowing into the deep soil layer m Calculated by formula (72); where T m is the deep soil temperature, which remains unchanged during the integration process. The initial field generation scheme is described later.

[0105] H m =K e C g (T g -T m ) (72) In order to better reflect the anthropogenic heat factor of urban land, the present invention adds anthropogenic heat flux to the energy balance equation of the model. That is, equation (60) is changed to equation (73): Among them, H a represents the anthropogenic heat flux.

[0106] H a According to formula (74), we can get: H a =Anth*(1-0.6cos(3.1415*(hour-5.5) / 12) (74) Among them, Anth represents the empirical coefficient of anthropogenic heat flux for different land uses (see Table 2), and hour represents the current simulation time, which describes the daily variation characteristics of anthropogenic heat flux.

[0107] Table 2 shows the land use classification and corresponding parameters used in the present invention. The specific parameters can be adjusted according to local conditions in actual application. In actual simulation, the surface temperature, sensible heat flux and latent heat flux of each land use type are calculated for each grid point, 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.

[0108] After the surface temperature, surface sensible heat and latent heat fluxes are calculated through the land surface model, the air temperature and humidity advection simulation is carried out through the integrated temperature and humidity advection equation, and the impact of land surface processes is transmitted to the atmosphere in the form of flux.

[0109] The temperature and humidity advection equations used in the present invention are as follows (75) and (76): In formulas (75)-(76), pt i,j,k represents potential temperature, calculated based on air temperature and air pressure, as shown in equation (57); q i,j,k represents specific humidity, which is calculated based on air temperature, relative humidity and air pressure, as shown in equations (58) and (59); u i,j,k 、v i,j,k Represents the north-south and east-west components of horizontal wind speed; ω i,j,k represents the vertical wind speed component in the terrain following coordinate system and is calculated according to formula (77): w i,j,k represents the vertical wind speed component; represents the terrain following coordinate system correction coefficient, calculated according to equations (78)-(80): In formulas (77)-(80), sigma k 、hmax、zg i,j They represent the model sigma layer value, the model highest layer altitude and the terrain altitude respectively.

[0110] In formulas (75)-(76), represents the horizontal diffusion coefficient and is calculated according to equations (81) and (82): In formulas (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, it is impossible to calculate the high-altitude turbulent energy and dissipation rate. Therefore, the following equations (83)-(84) are used to calculate and For model layers above 200 m above the ground, the vertical diffusion coefficients are set to 0; for vertical layers below 200 m above the ground, they are calculated according to equations (83)-(84). The superscript k represents the number of vertical layers; and They are the ground sensible heat flux and latent heat flux respectively. In this way, the land surface model simulation results can be transferred to the temperature and humidity advection process in the form of lower boundary flux.

[0111] In addition, when integrating the temperature and humidity advection equation, the mode space difference format uses the central difference, the time integration uses the Euler format, and the boundary conditions except the lower boundary use open boundary conditions, that is, the variable gradient at the boundary is 0. The integration step size and grid resolution are determined according to the actual situation.

[0112] Since time integration and spatial differentiation are inevitably accompanied by nonlinear error growth, the embodiment of the present application uses a data assimilation method to constrain the model error. The present invention uses Nudging as an assimilation algorithm.

[0113] Nudging adds approximation terms to the model integral equation so that each step of the integral can approach the observed field, and the longer the integration time, the better the assimilation effect. Nudging adds relaxation forcing terms to the model integral equation so that the error correction acts on each integration step. Its assimilation process will not cause serious discontinuity problems; since the model integral is used for dynamic constraints, the physical balance between the model variables will not be destroyed. Taking the temperature advection equation as an example, the temperature advection equation with the addition of Nudging relaxation forcing terms is shown in equation (85): Where, τ is the Nudging assimilation time window; G i,j,k is the observation increment. G of each grid point of the model is calculated as shown in formula (86), where the subscript i represents the observation data number, N obs represents the number of observations at the site; represents the site potential temperature observation; h represents the observation operator for interpolating the grid data to the site (linear interpolation is used in the present invention); pt represents the model grid potential temperature.

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

[0115] In formula (86), W t,i is the time weight coefficient, calculated according to formula (88). Where t represents the current mode integration time; t o Represents the observation time of the observation data; τ is the Nudging assimilation time window, in seconds.

[0116] The above are all the technical points of the microscale model. In actual simulation, follow the steps below: ① Based on the fine grid wind field and initial values ​​of three-dimensional air temperature, relative humidity, air pressure, and ground air pressure, model variables such as potential temperature and specific humidity are calculated. Among them, the initial value of surface temperature is replaced by the 2-meter air temperature interpolated by PRISM; soil temperature is calculated based on existing research. For example, see Liang LL, DARiveros-Iregui., REEmanuel., et al. A simple framework to estimate distributed soil temperature from discrete airtemperature measurements in data scarce regions. 2013. J. Geophys. Res. Atmos., 119, 407–417.

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

[0118] ③Simultaneously integrate the land surface model and temperature and humidity advection equation, and assimilate the meteorological station observation data.

[0119] ④ Based on the difference in potential temperature and specific humidity obtained from two consecutive integrations, determine whether a steady state has been reached as the criterion for ending the integration; or use the integration time determined by experience as the criterion for ending.

[0120] ⑤ Convert the simulated potential temperature and specific humidity into air temperature and relative humidity, and output the final three-dimensional air temperature, three-dimensional relative humidity, and three-dimensional wind field simulation results. In addition, the simulation results of ground pressure, surface temperature and other factors can also be output.

[0121] It should be noted that in order to achieve better simulation results, it is recommended to use 08:00 local time as the starting simulation time when performing hourly continuous simulation: at this time, the values ​​of near-surface air temperature, surface temperature and soil temperature are similar, and the model initialization effect is better. In other simulation times after the first time, the initial values ​​of surface temperature and soil temperature can use the simulation results of the previous time.

[0122] The present invention can quickly simulate micro-scale temperature, humidity and wind field elements with a resolution of less than 100m in complex terrain areas. Taking the simulation case of Chongqing Tainongyuan area as an example, the simulation grid spatial resolution is 30m, the horizontal dimension is 267×333, the vertical dimension is 21, and a total of 15,000 times are simulated. The case parameter setting, on a desktop computer with Windows operating system, uses 20 threads for calculation, and the total time consumption does not exceed 2 days. Figure 3-4 The multi-year average near-surface wind field and temperature field of Tainongyuan obtained by simulating the model constructed based on the present invention are given. It can be seen that the model can not only realize the rapid production of micro-scale meteorological element grid data, but also better reflect the impact of different land use types and complex terrain on the elements.

[0123] This embodiment also provides a system for rapid simulation of micro-scale meteorological elements in complex terrain areas.

[0124] refer to Figure 5 , a rapid simulation system for micro-scale meteorological elements in complex terrain areas, including: An underlying surface information calculation module 51 is used to calculate the underlying surface information of the area to be simulated; The complex terrain two-dimensional grid point background field generation module 52 is used to generate the temperature, pressure, relative humidity grid point meteorological element background field based on the underlying surface information by using the terrain parameter weighted regression model; and to generate the wind direction and wind speed grid point background field by using the inverse distance weighted interpolation; The error correction module 53 is used to perform preliminary error correction based on the grid meteorological element background field by using the STMAS multi-grid variation constrained by terrain parameters; to correct the grid wind direction and wind speed background field by using conventional STMAS; and to obtain the initial value of the grid meteorological element; The three-dimensional meteorological element initial field generation module 54 is used to generate three-dimensional meteorological initial values ​​according to the grid meteorological element initial values ​​by using land use types; 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; A wind field diagnosis module 55, used to input the initial value of the three-dimensional wind field into the wind field rapid diagnosis model to obtain a corrected three-dimensional fine grid wind field; The micro-scale model calculation module 56 simulates the fine grid air temperature and humidity based on the corrected three-dimensional fine grid wind field and other three-dimensional meteorological element initial values ​​through the integrated micro-scale model.

[0125] The various variations and specific examples of the methods provided in the above embodiments are also applicable to the system of this embodiment. Through the above detailed description of the method, those skilled in the art can clearly understand the implementation method of the system in this embodiment. For the sake of brevity, it will not be described in detail here.

[0126] In order to better execute the program of the above method, an embodiment of the present application also provides a computer device, which includes a processor and a memory.

[0127] Computer devices may be implemented in various forms, including computers, servers, and other devices.

[0128] The memory may be used to store instructions, programs, codes, code sets or instruction sets. The memory may include a program storage area and a data storage area, wherein the program storage area may store instructions for implementing an operating system, instructions for at least one function, and instructions for implementing the method provided in the above embodiment, etc.; the data storage area may store data involved in the method provided in the above embodiment, etc.

[0129] The processor may include one or more processing cores. The processor calls the data stored in the memory by running or executing the instructions, programs, code sets or instruction sets stored in the memory, performs various functions of the present application and processes data. The processor may 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 is understandable that for different devices, the electronic device for realizing the above-mentioned processor function can also be other, and the embodiments of the present application are not specifically limited.

[0130] The present application embodiment provides a computer-readable storage medium, for example, including: a USB flash drive, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk, and other media that can store program codes. The computer-readable storage medium stores a computer program that can be loaded by a processor and execute the method of the above embodiment.

[0131] An embodiment of the present application also provides a computer program product, which includes a computer program tangibly contained in a readable medium, wherein the computer program includes program code for executing any method in the embodiments of the present application, and the computer program can be downloaded and installed via a network, and / or installed from a removable medium (such as a disk, an optical disk, a magneto-optical disk, a semiconductor memory, etc.).

[0132] As described above, the above embodiments are only used to introduce the technical solution of the present application in detail, but the description of the above embodiments is only used to help understand the method and core idea of ​​the present application and should not be understood as limiting the present application. Any changes or substitutions that can be easily thought of by technicians in this technical field within the technical scope disclosed in the present application should be included in 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, the STMAS multi-grid variation constrained by terrain parameters is used to make preliminary error corrections; the grid wind direction and wind speed background field 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 site 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 transposition. 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