A dual-weight spatial interpolation method

Abstract

The present invention discloses a dual-weight spatial interpolation method. The method comprises the following steps: S1, determining the interpolation resolution, spatial interpolation range and interpolation attribute field, and dividing the grid unit; S2, calculating the two-dimensional coordinates of each center point of the grid to be interpolated, and screening the sampling points that meet the interpolation requirements based on the principle of closest distance; S3, calculating the distance weight of the sampling point and the local spatial position interpolation attribute variation weight, using the distance and difference coefficient as dual weights to interpolate the grid points to be interpolated; S4, outputting the obtained interpolation results in the form of a grid. Based on the principle of spatial self-similarity considered in the classical interpolation method, the present invention combines the variation characteristics of the attribute values ​​of the sampling points used for interpolation in space, designs a quantitative expression formula for the variability of the sampling points, and applies it as the second weight in the interpolation calculation, which can effectively reduce the adverse effects of sampling point outliers on the interpolation results.

Description

Technical Field

[0001] The present invention belongs to the technical field of spatial interpolation model construction, and relates to considering distance and local sampling point interpolation attribute variability as weight calculation factors in the spatial interpolation process to improve the accuracy of the interpolation result. Background Art

[0002] Spatial interpolation is a geospatial analysis technique that uses data from a limited number of sampling points to calculate the value of a variable at any location within a given area. Spatial interpolation allows a spatial variable to be extended from a point to a surface, revealing its spatial distribution. Therefore, it has widespread application in fields such as geography, the environment, and ecology. Interpolation accuracy is a crucial factor in determining its application value. Consequently, many scholars have discussed the accuracy characteristics of different spatial interpolation techniques, and are continually proposing and optimizing new ones. In particular, given the cost of data sampling, achieving highly accurate interpolation results based on a limited number of sampling points has become a crucial topic in the field of spatial analysis.

[0003] Spatial interpolation is an important method for constructing DEM, and DEM construction has also become an important way to analyze the accuracy characteristics of spatial interpolation methods. Traditional DEM construction is based on discrete sampling points, and the attribute values ​​of unknown elevation points are estimated through spatial interpolation technology, thereby realizing continuous surface modeling of the terrain. As a key technology for DEM construction, many scholars have not only conducted sufficient and in-depth research on the feasibility and characteristics of various traditional methods (such as polynomials, trend surface fitting, spline functions, Fourier series books) in terrain surface modeling, but also introduced interpolation methods from related disciplines such as Kriging, finite element and other interpolation technologies, and at the same time developed various interpolation methods for terrain data based on the characteristics of the terrain surface. The existing mature interpolation methods mainly include Thiessen polygon method, inverse distance weighted method, gradient distance inverse method, spline function method, trend surface method, area interpolation method, ordinary Kriging method, etc. Most of these methods involve interpolation based on distance as weight, such as the inverse distance weighted method, the inverse gradient distance method, and the ordinary kriging method. Although they make good use of the correlation of the planar distribution of sampling points, they ignore the mutation of interpolation attributes within a small range, resulting in a large impact range of individual sampling interpolation attribute mutation points, such as depressions in natural terrain and buildings in urban terrain, causing distortion of the simulated landforms in small areas.

[0004] In order to solve the above problems, it is necessary to further optimize the interpolation method and design a method to mitigate the abnormal interpolation results caused by the sudden change of attribute values ​​of local sampling points. This is also an urgent problem that needs to be solved in the field of spatial interpolation. Summary of the Invention

[0005] In response to the above problems, the present invention proposes a dual-weight spatial interpolation method, which can reduce the DEM construction error caused by the sudden change of the interpolation attribute values ​​of the sampling points.

[0006] The technical solution adopted in the present invention is as follows:

[0007] A dual-weighted spatial interpolation method includes the following steps:

[0008] S1, determine the interpolation resolution, spatial interpolation range and interpolation attribute fields, and divide the grid units;

[0009] S2, calculate the two-dimensional coordinates of the grid points and filter the neighboring sampling points;

[0010] S3, calculating the distance weight coefficient of the sampling point and the difference weight coefficient of the local spatial position interpolation attribute, and using the distance weight coefficient and the difference weight coefficient as double weights to interpolate the interpolation grid points;

[0011] S4, outputs the obtained interpolated attribute value in raster form.

[0012] Furthermore, the specific process of step S2 is as follows:

[0013] S21, searching for neighboring sampling points by buffering a certain distance with the interpolation grid point as the center, and determining whether the number of neighboring sampling points found is greater than 12. If so, executing step S22; otherwise, doubling the buffer radius to search for neighboring sampling points;

[0014] S22, calculating the distance between the sampling points searched in step S21 and the interpolation grid points, bubble sorting the selected sampling points according to the distance, and selecting the 12 sampling points with the smallest distance as the neighboring sampling points of the interpolation grid.

[0015] Furthermore, the specific process of step S3 is as follows:

[0016] S31, calculate the standard deviation and mean of the interpolation attributes of the grid point's neighboring sampling points;

[0017] S32, compare the interpolation attribute value of the neighboring sampling points with the mean difference of the interpolation attribute, determine the difference weight coefficient, and determine the difference weight coefficient according to the different standard deviation value ranges of the neighboring points:

[0018] Case 1: When Q<1, the formula is as follows:

[0019]

[0020] Case 2: When 1≤Q≤10, the formula is as follows:

[0021]

[0022] Case 3: When Q>10, the formula is as follows:

[0023]

[0024] In the above formula, Q and a represent the standard deviation and mean of the interpolation attribute of the neighboring sampling points of the calculated grid point, respectively; z represents the interpolation attribute value of the sampling point; and y represents the difference weight coefficient to be calculated.

[0025] S33, calculate weights and interpolate: Based on the inverse distance weighted interpolation method, add the difference weight coefficient y when calculating the weight, calculate the weight of each adjacent known elevation point and the interpolation grid point interpolation attribute, and interpolate each grid point.

[0026] When estimating unknown points based on known points, the method of the present invention not only considers the distance between known and unknown points, but also the spatial variation characteristics of each known point, and designs an expression to quantify this variability. This can significantly reduce the adverse impact of outliers in known points on the interpolation results and improve the accuracy of spatial interpolation. The technical features and beneficial effects of the present invention are as follows:

[0027] (1) The dual-weighted spatial interpolation method proposed in this invention adds interpolation attribute variation weights on the basis of distance as weight, better integrates the interpolation attribute value information, makes maximum use of the existing sampling point information, and can improve the interpolation accuracy.

[0028] (2) This method uses the variation of the interpolation attributes of the sampling points in a small range as a weight for interpolation. While reducing the error caused by the sampling points with sudden changes in the interpolation attributes on the local area simulation, it can also reduce the impact of the error points in the sampling points on the interpolation.

[0029] (3) Good scalability. The dual-weight spatial interpolation method proposed in the present invention actually provides a paradigm for constructing new interpolation weights. Existing interpolation methods can adopt the weight design in the method of the present invention to further optimize the interpolation method.

[0030] (4) The dual-weight spatial interpolation method proposed in the present invention is used to construct a digital elevation model as an application example, which shows that the accuracy of the method proposed in the present invention is significantly improved compared with the classical construction method. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 Schematic diagram of the process of the present invention.

[0032] Figure 2 This is a schematic diagram of determining the interpolation range from sampling points according to the method of the present invention.

[0033] Figure 3 Schematic diagram of screening adjacent sampling points by the method of the present invention.

[0034] Figure 4 Comparison of hillshade maps generated by (a) inverse distance weighted interpolation method and (b) DEM constructed by the method of the present invention for region 1 in an embodiment of the present invention.

[0035] Figure 5 Comparison of hillshade maps generated by (a) inverse distance weighted interpolation method and (b) DEM constructed by the method of the present invention for region 2 in an embodiment of the present invention. DETAILED DESCRIPTION

[0036] This embodiment provides a dual-weighted spatial interpolation method, such as Figure 1 As shown, the specific steps of the method are as follows:

[0037] Step 1: Initialization:

[0038] 1.a Determine the spatial interpolation range based on the sampling points and interpolation resolution ( Figure 2 ). According to the interpolation resolution l cell and the minimum x coordinate x of the sampling point min , maximum x-coordinate X max , minimum y coordinate y min and the maximum y-coordinate y max , determine the spatial interpolation range coordinate X min 、X max 、Y min and Y max , the calculation formula is as follows:

[0039] X min =x min -l cell / 2

[0040] X max =x max +l cell / 2

[0041] Y min =y min -l cell / 2

[0042] Y max =y max +l cell / 2

[0043] 1.b Calculate the number of rows and columns of the interpolation area grid based on the range and resolution. Divide the spatial interpolation range into several grid cells and calculate according to the following formula:

[0044] Calculated according to the following formula:

[0045] Number of rows = Int((x max -x min ) / l cell )

[0046] Number of columns = Int((y max -y min ) / l cell )

[0047] where x min 、x max The minimum x coordinate and maximum x coordinate of the interpolation range, y min 、y max is the minimum y coordinate and maximum y coordinate of the interpolation range, l cell is the interpolation resolution, and Int() represents the rounding function.

[0048] 1.c Determine the sampling point coordinate attributes and interpolation attributes, expressed as ternary coordinates (x, y, z).

[0049] Step 2: Calculate the two-dimensional coordinates of the center of each grid point according to the row and column numbers of the grid points. The formula is as follows:

[0050] x=jl cell +X min +0.5l cell

[0051] y=Y max -il cell -0.5l cell

[0052] Where i is the row number, j is the column number, x0 and y0 are the coordinates of the starting point of the interpolation area, and l cell is the resolution.

[0053] Step 3: Filter neighbor sampling points ( Figure 3 ), with the interpolation grid point as the center, buffer a certain distance to search for nearby sampling points to determine whether the searched sampling points are greater than 12. If so, stop searching. Otherwise, double the search radius and continue searching. Calculate the distance d between the sampling point and the grid point. Sort the selected sampling points by distance in a bubble order and select the 12 nearest neighboring sampling points with the smallest distance. The formula is as follows:

[0054]

[0055] where x i 、y i are the x-coordinate and y-coordinate of the nearest sampling point respectively. o 、y o are the x-coordinate and y-coordinate of the interpolation grid point respectively.

[0056] Step 4: Calculate the interpolation attribute standard deviation Q and mean a of the grid point's neighboring sampling points. The formula is as follows:

[0057] a=(z1+z2+…zn ) / n

[0058]

[0059] Where z is the interpolation attribute value of the sampling point, and n is the total number of neighboring sampling points.

[0060] Step 5: Compare the interpolation attribute value of the neighboring sampling points with the interpolation attribute mean a, and determine the difference weight coefficient y for calculating the interpolation attribute. The difference weight coefficient y is determined according to the different value ranges of the standard deviation Q value of the neighboring points. There are three cases:

[0061] Case 1: When Q<1, the formula is as follows:

[0062]

[0063] Case 2: When 1≤Q≤10, the formula is as follows:

[0064]

[0065] Case 3: When Q>10, the formula is as follows:

[0066]

[0067] Step 6: Weight calculation and interpolation. Based on the inverse distance weighted interpolation method, the weight λ of the i-th sampling point is calculated. i Add the difference weight coefficient y i , interpolate each grid point. The calculation formula is as follows:

[0068]

[0069] where y i is the difference weight coefficient, d i is the plane distance between the interpolation grid point and the neighboring sampling point, λ i is the weight coefficient of the i-th sampling point, z i is the elevation coordinate of the neighboring sampling point.

[0070] Step 7: Repeat steps 2 to 6 based on the grid divided in step 1.b, traverse each grid point, and calculate the interpolated attribute value of the grid point.

[0071] Step 8: Output the results. The specific process is as follows:

[0072] 8.a. Output the header information of each file in sequence: including the number of columns, number of rows, vertical coordinate of the southwest corner grid cell, vertical coordinate of the southwest corner grid cell, grid spacing, and invalid data area value.

[0073] 8.b. Output the grid point interpolation attribute values ​​in sequence according to the row and column numbers. When the number of output elements equals the number of columns in the result file, start a new row and restart the count.

[0074] 8.c. Output all vector solution results in TIF format.

[0075] Existing classical interpolation methods only consider distance as a weight during interpolation, ignoring sudden changes in the interpolation attributes of sampling points in local areas, which can cause distortion in the simulated landforms of small areas. To verify the practicality of the method of the present invention, this example selected two areas with different terrains and ranges as study areas. This example compared the interpolation method of the present invention with the inverse distance weighted interpolation method integrated with ArcGIS to verify the accuracy advantage.

[0076] Step 1: Separate the training points and verification points with a separation ratio of 19:1.

[0077] Step 2: interpolation of training points: interpolation of training points is performed using the method of the present invention and the classical interpolation method to obtain a grid DEM.

[0078] Compared with the inverse distance weighted interpolation method, the DEM hill shadow contrast obtained by interpolation (such as Figure 4 、 Figure 5 ), it can be seen that the interpolation method of the present invention has certain advantages in terms of smoothness.

[0079] Step 3: Accuracy test: extract the elevation value of the verification point and take the absolute value of the difference between the grid DEM elevation value obtained by interpolation and the original elevation value. By comparing with the classic inverse distance weighted interpolation method (see Tables 1 and 2), it can be seen that this interpolation method has higher accuracy in terms of interpolation precision.

[0080] Table 1 Accuracy comparison of region 1

[0081] Study Area Inverse distance weighted method Method of the present invention Number of modeling points 959 959 Number of verification points 50 50 Minimum error 0.00162 0 Maximum error 1.01698 1.08585 Error sum 10.804301 10.04698 Mean error 0.216086 0.20094 Error standard deviation 0.250164 0.249854

[0082] Table 2 Accuracy comparison of region 2

[0083] Study Area Inverse distance weighted method Method of the present invention Number of modeling points 1321 1321 Number of verification points 70 70 Minimum error 0.00007 0.00004 Maximum error 2.6619 2.6534 Error sum 25.692591 22.325529 Mean error 0.367037 0.318936 Error standard deviation 0.626674 0.614951

Claims

1. A dual-weighted spatial interpolation method, characterized in that: The steps include: S1, determine the interpolation resolution, spatial interpolation range and interpolation attribute fields, and divide the grid units; S2, calculate the two-dimensional coordinates of the grid points and filter the neighboring sampling points; S3, calculating the distance weight coefficient of the sampling point and the difference weight coefficient of the local spatial position interpolation attribute, and using the distance weight coefficient and the difference weight coefficient as double weights to interpolate the interpolation grid points; S4, outputting the obtained interpolated attribute values ​​in raster form; The specific process of step S2 is as follows: S21, searching for neighboring sampling points by buffering a certain distance with the interpolation grid point as the center, and determining whether the number of neighboring sampling points found is greater than 12. If so, executing step S22; otherwise, doubling the buffer radius to search for neighboring sampling points; S22, calculating the distance between the sampling points searched in step S21 and the interpolation grid points, bubble sorting the selected sampling points according to the distance, and selecting the 12 sampling points with the smallest distance as the neighboring sampling points of the interpolation grid; The specific process of step S3 is as follows: S31, calculate the standard deviation and mean of the interpolation attributes of the grid point's neighboring sampling points; S32, compare the interpolation attribute value of the neighboring sampling points with the mean difference of the interpolation attribute, determine the difference weight coefficient, and determine the difference weight coefficient according to the different standard deviation value ranges of the neighboring points: Case 1: When Q<1, the formula is as follows: Case 2: When 1≤Q≤10, the formula is as follows: Case 3: When Q>10, the formula is as follows: In the above formula, Q and a represent the standard deviation and mean of the interpolation attribute of the neighboring sampling points of the calculated grid point, z represents the interpolation attribute value of the sampling point, and y represents the difference weight coefficient to be calculated; S33, calculate weights and interpolate: Based on the inverse distance weighted interpolation method, add the difference weight coefficient y when calculating the weight, calculate the weight of each adjacent known elevation point and the interpolation grid point interpolation attribute, and interpolate each grid point. The calculation formula is as follows: where y i is the difference weight coefficient, d i is the plane distance between the interpolation grid point and the neighboring sampling point, λ i is the weight coefficient of the i-th sampling point, z i is the elevation coordinate of the neighboring sampling point.

Images