A method for interpolating meteorological elements in three-dimensional space

Abstract

The application discloses a meteorological element three-dimensional space interpolation method based on an exponential base function, obtains an initial field and a terrain field of a target station or a grid point and carries out data unit conversion; a plurality of auxiliary stations or grid points closest to the target station or grid point in the three-dimensional space are selected, distance relations between each point in the auxiliary stations or grid points and the target station or grid point are obtained as first weight coefficients; the first weight coefficients are adjusted according to the latitude and longitude direction relation to obtain second weight coefficients; the first weight coefficients are adjusted according to the height direction relation to obtain third weight coefficients; a weighted summation equation is established based on the first, second and third weight coefficients to obtain a three-dimensional space interpolation result of the meteorological element of the target station or grid point. The application solves the shortcoming that the dependent variable amplitude decreases with the independent variable high value area when the cosine base function is used to calculate the direction weight, avoids the limitation of the regression equation, has a wide application range, and reduces the complexity of operation and running time consumption.

Description

Technical Field

[0001] This invention belongs to the field of meteorological engineering technology, specifically relating to a three-dimensional spatial interpolation method for meteorological elements. Background Technology

[0002] Interpolation of meteorological elements is widely used in the field of meteorological engineering technology. For example, in order to achieve point-to-point weather forecasting, the gridded forecast results of numerical forecasting models are interpolated to a specific meteorological observation station, or the data of meteorological observation stations are further interpolated to certain specific geographic coordinates to obtain more refined observation results that are easier to process later. All of these require the use of spatial interpolation algorithms.

[0003] Spatial interpolation has a long history of development and application in geographic information processing. Traditional spatial interpolation methods mainly include Thiessen polygon method, Kriging interpolation method, inverse distance weighting method (IDW), distance direction weighting method (ADW), etc. However, after introducing these methods into the meteorological field, an inevitable problem arises: the use of four-dimensional data. This is because the interpolation target of geographic information is elevation, and the data used includes longitude, latitude, and elevation in three dimensions. However, the interpolation target in meteorological engineering is a certain meteorological element, and the data used includes longitude, latitude, elevation, and meteorological element in four dimensions. Therefore, the existing processing method is to consider longitude, latitude, and elevation separately. Among them, longitude and latitude are processed by distance direction weighting method to obtain a weight, and then elevation is fitted with a regression equation to obtain another weight. Finally, the weights are combined to obtain the final interpolation result (Comparative documents: [1] Spatial interpolation method and system for precipitation, electronic equipment and readable storage medium, patent application publication number CN 113553696 A; [2] Spatial interpolation method for precipitation data in complex terrain area and computer storage medium, patent application publication number CN 113779505 A). However, the problem with this method is that: (1) the cosine basis function for calculating the directional weights fails to efficiently reflect the principle that the fewer auxiliary stations or grid points in a certain direction around the target station or grid point, the greater the weight of that station; (2) the use of regression equations to establish the elevation gradient relationship of meteorological elements is limited by the linear nature of the equations, and can only fit a few meteorological elements such as precipitation, and also requires the subjective division of regional grids. Summary of the Invention

[0004] To address the technical problems mentioned in the background section, this invention proposes a three-dimensional spatial interpolation method for meteorological elements.

[0005] To achieve the above-mentioned technical objectives, the technical solution of the present invention is as follows:

[0006] A three-dimensional spatial interpolation method for meteorological elements includes the following steps:

[0007] Acquire the initial field and terrain field of the target site or grid point and perform data unit conversion;

[0008] Select multiple auxiliary stations or grid points in three-dimensional space that are closest to the interpolation target station or grid point, and obtain the distance relationship between each point in each auxiliary station or grid point and the target station or grid point as the first weight coefficient;

[0009] Based on the latitude and longitude relationship between each auxiliary station or grid point and the target station or grid point, the first weighting coefficient is adjusted to obtain the second weighting coefficient.

[0010] Based on the elevation and directional relationship between each auxiliary station or grid point and the target station or grid point, the first weighting coefficient is adjusted to obtain the third weighting coefficient.

[0011] Based on the first, second, and third weighting coefficients mentioned above, a weighted summation equation is established to obtain the three-dimensional spatial interpolation results of meteorological elements at the target station or grid point.

[0012] Preferably, the data unit conversion specifically refers to converting latitude and longitude coordinates into rectangular coordinates, that is, converting geographic coordinates WGS84 to the WebMercator projection plane, as expressed by the following formula:

[0013] X = lon × equ / 180 (1)

[0014] Y=log(tan((90+lat)×π / 360)) / (π / 180)×equ / 180 (2)

[0015] Z=Ge / 9.8×9.8 / g+2≈Ge / 9.8+2 (3)

[0016] In the formula, X and Y are the planar coordinates after WebMercator projection, in meters; lon and lat are the latitude and longitude coordinates, in degrees; equ is half the circumference of the Earth's equator; Z is the sea level height, i.e., elevation, in meters; Ge is the potential; and g is the gravitational acceleration.

[0017] Preferably, the formula for the first weighting coefficient w1i is as follows:

[0018]

[0019] In the formula, d is the Euclidean distance in three-dimensional space between the target station or grid point and its nearest auxiliary station or grid point i, i∈{1,2...N}, and N is the number of auxiliary stations or grid points. c The critical distance is the average distance between all auxiliary stations or grid points in three-dimensional space; CF is the correction coefficient.

[0020] Preferably, the correction coefficient CF is set to 4.

[0021] Preferably, the number of auxiliary stations or grid points selected in three-dimensional space that are closest to the interpolation target station or grid point does not exceed 9.

[0022] Preferably, the formula for the second weighting coefficient w2i is as follows:

[0023]

[0024] In the formula, w1i is the first weight coefficient; N is the number of auxiliary stations or grid points; j is the complement of auxiliary station or grid point i relative to the universal set {1,2...N}, with a quantity of N-1; θ is the angle in radians between the lines connecting auxiliary station or grid point i and j to the target station or grid point in the two-dimensional latitude and longitude space and the target station or grid point. The calculation formula is expressed as follows:

[0025]

[0026] In the formula, doi is the distance between the target station or grid point and the auxiliary station or grid point i, doj is the distance between the target station or grid point and the auxiliary station or grid point j, and dij is the distance between the auxiliary station or grid point i and j.

[0027] Preferably, the third weighting coefficient w3i is expressed by the following formula:

[0028] w3 i =w1 i (2-2 2φ / π (7)

[0029] In the formula, w1i is the first weighting coefficient, and φ is the vertical elevation or depression angle in radians of the auxiliary station or grid point i relative to the target station or grid point. The calculation formula is expressed as follows:

[0030] φ=arctan(|Z o -Z i | / d oi (8)

[0031] In the formula, Zo is the elevation of the target station or grid point, Zi is the elevation of the auxiliary station or grid point i, and doi is the distance between the target station or grid point and the auxiliary station or grid point i in the two-dimensional space of latitude and longitude.

[0032] Preferably, a weighted summation equation is established to obtain the three-dimensional spatial interpolation result of a certain meteorological element at the target station or grid point. The MEo formula is expressed as follows:

[0033]

[0034] In the formula, MEi is the meteorological element value of the auxiliary station or grid point i, w1i is the first weight coefficient; N is the number of auxiliary stations or grid points, w2i is the second weight coefficient; and w3i is the third weight coefficient.

[0035] The beneficial effects of adopting the above technical solution are as follows:

[0036] The significance of the first, second, and third weighting coefficients:

[0037] The first weight coefficient describes the reverse distance relationship between the target station or grid point and the auxiliary station or grid point. That is, the closer the distance, the higher the weight coefficient, and the greater the influence of the auxiliary station or grid point on the target station or grid point. According to the comparison documents [1][2], the Gaussian kernel function is used in the form of equation (4). Since the elevation is integrated into the direction weight calculation, the distance calculation also takes into account the elevation factor. It is based on the Euclidean distance in the three-dimensional space of longitude, latitude and elevation.

[0038] The second and third weighting coefficients describe the latitude-longitude direction relationship and the elevation direction relationship, where the exponential basis function 2 for the latitude-longitude direction relationship... θ / π -1 As shown in equation (5), it is a function of the angle θ between the auxiliary station or grid point and the remaining N-1 auxiliary stations or grid points at the target station or grid point. If there are many auxiliary stations or grid points clustered together in a certain direction centered on the target station or grid point, the angle between these auxiliary stations or grid points will be relatively small after being mapped by the exponential basis function. Therefore, the weight coefficients of their latitude and longitude directions will also be relatively small. This is in line with the principle that the influence of each point in the clustered auxiliary stations or grid points on the target station or grid point is relatively small. The principle that the elevation direction relationship needs to satisfy is that the smaller the vertical elevation or depression angle of the auxiliary station or grid point relative to the target station or grid point, the greater the influence of the auxiliary station or grid point on the target station or grid point. Therefore, the exponential basis function 2-2 is established. 2φ / π As shown in equation (7), the smaller the pitch angle φ, the larger the function value.

[0039] The reason why the second and third weighting coefficients were not directly merged into the directional relationship in three-dimensional space (that is, the formula (5) is extended to three-dimensional space) is that when the elevation differences between auxiliary stations or grid points are significant, the angle between the auxiliary stations or grid points and the target station or grid point in three-dimensional space cannot correctly reflect the above two principles: the more auxiliary stations or grid points that are clustered in the horizontal direction, the smaller their weight; the smaller the elevation angle of the auxiliary stations or grid points in the vertical direction, the greater their weight. In this case, such as Figure 1 As shown;

[0040] Figure 1In a two-dimensional latitude and longitude space, auxiliary stations or grid points 1, 2, and 3 are evenly spaced and equidistant from the target station or grid point. Therefore, the second weight coefficients of auxiliary stations or grid points 1 and 3 are equal. However, the elevation angle of auxiliary station or grid point 1 relative to the target station or grid point is greater than that of auxiliary station or grid point 3, so the third weight coefficient of auxiliary station or grid point 1 is smaller than that of auxiliary station or grid point 3. This makes the influence of auxiliary station or grid point 1 on the target station or grid point less than that of auxiliary station or grid point 3. However, if the second and third weight coefficients are directly combined into a directional relationship in three-dimensional space, the weight of auxiliary station or grid point 1 will be the largest. This is because, with significant elevation differences, the angle between auxiliary station or grid point 1 and auxiliary stations or grid points 2 and 3 at the target station or grid point will be greatly increased, thus raising its weight to exceed that of auxiliary stations or grid points 2 and 3.

[0041] 2. Exponential basis functions

[0042] The graphs of the exponential basis functions used for the second and third weighting coefficients are shown below. Figure 2 (b) and Figure 2 (c), as can be seen in the figure, the 1-cosθ basis functions used in the comparison documents [1][2] are as follows: Figure 2 (a) One problem is that the increase in the dependent variable y decreases in the high-value region of the independent variable θ, which contradicts the principle that the larger the θ, the greater the increase in its weight. Therefore, we change to using 2. θ / π This problem can be solved by using a -1 basis function. Therefore, it is evident that using such a method... Figure 2 (b) and Figure 2 (c) exponential basis function ratio Figure 2 (a) The cosine basis function is more suitable for calculating the direction weight coefficient.

[0043] 3. The interpolation method is independent of the characteristics of meteorological elements.

[0044] Because existing spatial interpolation methods for meteorological elements use regression equations to fit the relationship between meteorological elements and elevation gradients to obtain weights when dealing with elevation problems, the establishment of regression equations should conform to the relationship between meteorological elements and elevation (such as linear relationship). This results in a significant limitation on the meteorological elements that the interpolation method can be applied to (such as precipitation elements that can establish linear regression relationships in the comparative document [1]). However, this invention does not use regression equations and only considers the distance and direction relationship between auxiliary stations or grid points and target stations or grid points in three-dimensional space, thus making it more universally applicable. At the same time, because it reduces a large number of algorithm steps for fitting the relationship between meteorological elements and elevation gradients, the computational complexity and running time of this invention are reduced. Attached Figure Description

[0045] Figure 1 This is a schematic diagram of auxiliary sites or grid points in the three-dimensional space of this patent;

[0046] Figure 2 This is a graph of the basis functions of this patent, in which... Figure 2 (a) shows the graph of the 1-cosθ basis functions; Figure 2 (b) is 2 θ / π -1 basis function graph; Figure 2 (c) is 2-2 2φ / π Graphs of basis functions;

[0047] Figure 3 This is a flowchart of the three-dimensional spatial interpolation method for meteorological elements in this patent;

[0048] Figure 4 This is the interpolation result of the 2-meter dew point temperature at the national meteorological observation station at 00:00 (UTC) on November 5, 2021. Detailed Implementation

[0049] The technical solution of the present invention will be described in detail below with reference to the accompanying drawings.

[0050] This invention patent can be applied to any kind of meteorological element observed or reanalyzed on the ground, and there are no restrictions on whether the interpolation target station or grid point and the auxiliary station or grid point are station data or grid point data. Here, for the purpose of illustrating a specific implementation, we take the 2-meter dew point temperature on the ground as an example to perform three-dimensional spatial interpolation from grid point to station.

[0051] The flowchart of this invention patent is as follows Figure 3 As shown, the data source is the initial grid field and topographic field of the ECMWF fine-grid numerical weather prediction model (0.125° x 0.125° grid spacing). The meteorological element selected is the surface 2-meter dew point temperature at 00:00 (UTC) daily in China (3°51'-53°33' N, 73°33'-135°05' E) in 2021. The interpolation target is the surface 2-meter dew point temperature at 00:00 (UTC) daily from 2165 national meteorological observation stations across the country.

[0052] First, the EC model initial field and national station data are converted to rectangular coordinates, that is, the geographic coordinates WGS84 are converted to the WebMercator projection plane, as shown in formulas (10) and (11):

[0053] X=lon×20037508.342789 / 180 (10)

[0054] Y=log(tan((90+lat)×π / 360)) / (π / 180)×20037508.342789 / 180 (11)

[0055] In formulas (10) and (11), X and Y are the planar coordinates after WebMercator projection, in meters (m); lon and lat are the latitude and longitude coordinates, in degrees (°).

[0056] Secondly, the unit of EC topographic field data is potential (m). 2 / s 2 Therefore, it also needs to be converted to meters (m), as shown in formula (12):

[0057] Z=Ge / 9.8×9.8 / g+2≈Ge / 9.8+2 (12)

[0058] In formula (12), Z is the sea level height, i.e., elevation, in meters (m); Ge is the potential; g is the gravitational acceleration. Considering that g does not change much with latitude at near-ground height, g≈9.8.

[0059] Finally, after obtaining gridded dew point temperature data with unified units of X (longitude), Y (latitude), and Z (elevation) and interpolated target station location data (station elevation and topographic field elevation should both be +2m), the first, second, and third weighting coefficients are calculated using the method proposed in this patent. The first weighting coefficient is w1. i =(e -d / d c ) 4 Where i∈{1,2...4}, d and d c The meaning is as shown in equation (4); the second weighting coefficient Where j has the meaning as shown in equation (5), and θ has the meaning as shown in equation (6); the third weighting coefficient w3 i =w1 i (2-2 2φ / π ), where φ has the meaning as shown in equation (8), and finally the dew point temperature interpolation result of the target station is obtained. ME i The meaning is as shown in equation (9), and the example effect is as follows. Figure 4 As shown.

[0060] This invention patent's embodiments compare the errors of interpolation results considering elevation and not considering elevation. Specifically, the ground 2-meter dew point temperature observations at 00:00 (UTC) daily in 2021 from 2165 national meteorological observation stations across China were used as the verification standard. According to the above-described implementation method, weighted summation equations were established using the first and second weighting coefficients (Scheme 1) and the first, second, and third weighting coefficients (Scheme 2) to obtain the interpolation results. Then, the monthly mean absolute error (MMAE) of the two schemes and the observations was compared, and the ratio of MMAE2 (Scheme 2) to MMAE1 (Scheme 1) under different elevation intervals was statistically analyzed. The results are shown in Table 1.

[0061] Table 1. Improvement of error after considering elevation element interpolation in different elevation intervals.

[0062]

[0063]

[0064] As shown in Table 1, the MMAE2 / MMAE1 ratio does not exceed 100% in each elevation interval. Therefore, the error of Scheme 2 is smaller than that of Scheme 1. The ratio decreases overall as the elevation increases, indicating that the error improvement effect of Scheme 2 is better at higher elevations. This is because the terrain around stations with high elevations is often more complex (such as mountainous areas), resulting in significant differences in elevation between nearby points. At this time, the advantages of the third weighting coefficient and exponential basis function of Scheme 2 are better reflected, achieving better results than Scheme 1.

[0065] Those skilled in the art will understand that embodiments of this application can be provided as methods, systems, or computer program products. Therefore, this application can take the form of a completely hardware embodiment, a completely software embodiment, or an embodiment combining software and hardware aspects. Furthermore, this application can take the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) containing computer-usable program code. The solutions in the embodiments of this application can be implemented in various computer languages, such as the object-oriented programming language Java and the interpreted scripting language JavaScript.

[0066] This application is described with reference to flowchart illustrations and / or block diagrams of methods, apparatus (systems), and computer program products according to embodiments of this application. It will be understood that each block of the flowchart illustrations and / or block diagrams, and combinations of blocks in the flowchart illustrations and / or block diagrams, can be implemented by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, special-purpose computer, embedded processor, or other programmable data processing apparatus to produce a machine, such that the instructions, which execute via the processor of the computer or other programmable data processing apparatus, generate instructions for implementing the flowchart... Figure 1 One or more processes and / or boxes Figure 1 A device that provides the functions specified in one or more boxes.

[0067] These computer program instructions may also be stored in a computer-readable storage medium that can direct a computer or other programmable data processing device to function in a particular manner, such that the instructions stored in the computer-readable storage medium produce an article of manufacture including instruction means, which are implemented in a process Figure 1 One or more processes and / or boxes Figure 1 The function specified in one or more boxes.

[0068] These computer program instructions may also be loaded onto a computer or other programmable data processing equipment to cause a series of operational steps to be performed on the computer or other programmable equipment to produce a computer-implemented process, thereby providing instructions that execute on the computer or other programmable equipment for implementing the process. Figure 1 One or more processes and / or boxes Figure 1 The steps of the function specified in one or more boxes.

[0069] Although preferred embodiments of this application have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including the preferred embodiments as well as all changes and modifications falling within the scope of this application.

[0070] Obviously, those skilled in the art can make various modifications and variations to this application without departing from the spirit and scope of this application. Therefore, if such modifications and variations fall within the scope of the claims of this application and their equivalents, this application also intends to include such modifications and variations.

Claims

1. A three-dimensional spatial interpolation method for meteorological elements, characterized in that, Includes the following steps: Obtain the initial field and topographic field of the target site and perform data unit conversion; Select multiple auxiliary stations in three-dimensional space that are closest to the interpolation target station, and obtain the distance relationship between each point in each auxiliary station and the target station as the first weight coefficient; Based on the latitude and longitude relationship between each auxiliary station and the target station, the first weighting coefficient is adjusted to obtain the second weighting coefficient; Based on the elevation and directional relationship between each auxiliary station and the target station, the first weighting coefficient is adjusted to obtain the third weighting coefficient; A weighted summation equation is established based on the first weighting coefficient, the second weighting coefficient, and the third weighting coefficient to obtain the three-dimensional spatial interpolation results of meteorological elements at the target station. Among them, the first weighting coefficient w1 i The formula is expressed as follows: ； （1） In the formula, d is the Euclidean distance between the target station and its nearest auxiliary station i in three-dimensional space. N is the number of auxiliary stations, d c The critical distance is the average distance between all auxiliary stations in three-dimensional space; CF is the correction coefficient. Second weighting coefficient w2 i The formula is expressed as follows: ； （2） In the formula, The formula for calculating the angle in radians between the lines connecting auxiliary stations i and j to the target station in two-dimensional latitude and longitude space and the target station is as follows: ； （3） In the formula, d oi Let d be the distance between the target station and the auxiliary station i in a two-dimensional space of latitude and longitude. oj Let d be the distance between the target site and the auxiliary site j. ij The distance between auxiliary stations i and j; Third weighting coefficient w3 i The formula is expressed as follows: ； （4） In the formula, The vertical elevation or depression angle (in radians) of auxiliary station i relative to the target station is calculated using the following formula: ； （5） In the formula, Z o Z is the elevation of the target station. i The elevation of auxiliary station i; A weighted summation equation is established to obtain the three-dimensional spatial interpolation result (ME) of a certain meteorological element at the target station. o The formula is expressed as follows: ； （6） In the formula, ME i For auxiliary station i, the meteorological element values ​​are provided. Data unit conversion specifically refers to converting latitude and longitude coordinates into rectangular coordinates, that is, converting geographic coordinates WGS84 to the WebMercator projection plane. The formula is expressed as follows: ； （7） ； （8） ； （9） In the formula, X and Y are the planar coordinates after WebMercator projection, in meters; lon, 1 represents latitude and longitude coordinates in degrees; equ is half the circumference of the Earth's equator; Z is the sea level (elevation) in meters; Ge is the potential; g is the acceleration due to gravity.

2. The three-dimensional spatial interpolation method for meteorological elements according to claim 1, characterized in that, The correction factor CF is set to 4.

3. The three-dimensional spatial interpolation method for meteorological elements according to claim 1, characterized in that, The number of auxiliary stations selected in three-dimensional space that are closest to the interpolation target station shall not exceed 9.

Images