A method for extracting abnormal boundary in interpretation map based on geophysical exploration
By inverting and normalizing geophysical exploration data, the problem of inaccurate identification of geological anomaly boundaries was solved, and the clear display and precise positioning of anomaly boundaries were achieved.
Patent Information
- Authority / Receiving Office
- CN · China
- Patent Type
- Patents(China)
- Current Assignee / Owner
- RES INST OF COAL GEOPHYSICAL EXPLORATION
- Filing Date
- 2022-11-30
- Publication Date
- 2026-06-05
AI Technical Summary
Existing geophysical exploration methods cannot effectively eliminate the covering and volume effects of near-surface low-resistivity anomalies, leading to inaccurate identification of geological anomaly boundaries.
By acquiring resistivity data from the surveyed area and performing inversion, a resistivity spatial grid file is established. The median value and normalized derivative operator of any two adjacent grid nodes are calculated, the resistivity is redefined, and a normalized two-dimensional plot is drawn to display the anomalous boundary.
It improves the ability to identify geological anomalies, clearly showing the boundaries of anomalies and providing a basis for the precise location of underground geological bodies.
Smart Images

Figure CN115797652B_ABST
Abstract
Description
Technical Field
[0001] This invention relates to the field of image data processing technology, and more specifically to a method for extracting anomalous boundaries in interpretation maps based on geophysical exploration. Background Technology
[0002] Geophysical exploration involves detecting and studying the distribution and physical properties of rocks, ores, and other geological formations on Earth, such as density, conductivity, and radioactivity. Geophysical exploration technology primarily utilizes instruments like geophysical probes to monitor physical data within the Earth's geological structure, analyze its changing patterns, and then combine this with modern high-tech electronic information technology, computer technology, communication engineering technology, and relevant knowledge of physics to compare and study the differences between the Earth's geological structure and the monitored physical data. The data obtained through this technology possesses high accuracy, scientific validity, and feasibility. Therefore, it is widely used in various industries, including construction engineering, mineral mining, energy exploration, and geological disaster prevention, playing a crucial role in these sectors and significantly promoting the steady development of various industries in my country, thus accelerating the pace of social development.
[0003] Currently, commonly used geophysical exploration methods include direct current resistivity (DC) methods, transient electromagnetic methods, magnetotellurics (MET) methods, and wide-area electromagnetic methods. After processing the collected data, resistivity one-dimensional curves, two-dimensional cross-sectional maps, planar maps, and three-dimensional maps describing subsurface structures and geological formations can be obtained. These maps represent the direct response of geophysical exploration methods to subsurface geological bodies. However, existing geophysical exploration methods, during use, cannot completely eliminate factors that interfere with the data due to limitations in processing algorithms. These interfering factors include the covering effect of near-surface low-resistivity anomalies and the volume effect of methodological techniques, among others. These effects are superimposed on the maps obtained from geophysical exploration, thus affecting the identification of geological anomaly boundaries. Therefore, it is necessary to provide a method that can accurately extract the boundaries of geological anomalies. Summary of the Invention
[0004] To overcome the shortcomings of existing technologies, this invention provides a method for extracting anomalous boundaries in interpretation maps based on geophysical exploration. This method addresses the technical problem that existing technologies cannot accurately identify the boundaries of geological anomalies, thereby improving the ability to identify geological anomalies and providing a basis for boundary delineation.
[0005] To solve the above problems, the technical solution adopted by the present invention is as follows:
[0006] A method for extracting anomaly boundaries in interpretation maps based on geophysical exploration includes the following steps:
[0007] Obtain resistivity data of the survey area, and perform inversion on the resistivity data to obtain inverted resistivity data;
[0008] A resistivity space grid file is obtained based on the inverted resistivity data;
[0009] Obtain the median value of the resistivity of any two adjacent grid nodes in the resistivity space grid file;
[0010] Obtain the normalized derivative operator of the resistivity of any two adjacent grid nodes in the resistivity space grid file;
[0011] After redefining the resistivity of each grid node in the resistivity space grid file based on the median value and the normalized derivative operator, the normalized resistivity of each grid node is obtained.
[0012] A normalized two-dimensional graph is plotted based on the spatial relationship of the normalized resistivity, and the abnormal boundaries are clearly displayed through the normalized two-dimensional graph.
[0013] In a preferred embodiment of the present invention, obtaining a resistivity spatial grid file based on the inverted resistivity data includes:
[0014] Different detection points are set according to the size of the survey area, and the measured resistivity and measured potential values of the different detection points are collected to establish an initial resistivity model.
[0015] The model resistivity of the survey area is obtained based on the initial resistivity model;
[0016] The inverted resistivity data is obtained by inverting the measured resistivity, the model resistivity, and the initial resistivity model.
[0017] The resistivity space grid file is obtained based on the inverted resistivity data.
[0018] In a preferred embodiment of the present invention, obtaining the inverted resistivity data includes:
[0019] An inversion equation is established based on the measured resistivity, the model resistivity, and the initial resistivity model, as shown in Formula 1:
[0020] (1);
[0021] In the formula, S is the sensitivity matrix, c is the measured resistivity, ρ is the model resistivity, G is the forward modeling operator, and R... dd Let R be the covariance matrix of the measured resistivity. mm Let ρ0 be the covariance matrix of the model resistivity, and ρ0 be the initial resistivity model. kThe resistivity of the model after the k-th iteration;
[0022] The inversion resistivity data are obtained based on the inversion equation.
[0023] In a preferred embodiment of the present invention, when setting different detection points according to the size of the survey area, the method includes:
[0024] The first and second setting points are determined based on the size of the survey area.
[0025] Install the electrode device at the first setting point and the second setting point;
[0026] Between the first setting point and the second setting point, multiple detection lines are set according to a preset detection line distance;
[0027] On each testing line, multiple testing points are set according to the preset distance between testing points.
[0028] In a preferred embodiment of the present invention, the initial resistivity model is established by including:
[0029] The first detection electrode of the detection device is set at the first detection point of each detection line;
[0030] The second detection electrode of the detection device is set at the second detection point of each detection line;
[0031] Obtain the measured resistivity of the first detection point and the second detection point, as well as the measured potential value between the detection points;
[0032] The resistivity between the first detection point and the second detection point is obtained based on the measured resistivity and the measured potential value.
[0033] The initial resistivity model is established based on the resistivity between the first detection point and the second detection point;
[0034] The second detection point is the remaining detection point on each detection line other than the first detection point, and the second detection electrode is sequentially set on the remaining detection points far away from the first detection point.
[0035] In a preferred embodiment of the present invention, when obtaining the resistivity space grid file based on the inverted resistivity data, the method further includes:
[0036] Based on the distribution of the inverted resistivity data, the region and grid size to be interpolated are determined, the region is gridded to simulate a preliminary two-dimensional grid surface, and each grid node in the preliminary two-dimensional grid surface is obtained.
[0037] The data of each grid node is checked and analyzed to obtain the distance value between each grid node. If there is a distance value less than a preset distance value, the grid node is removed.
[0038] A mutation function is selected, and for each grid node, the Kriging estimate corresponding to each grid node is obtained based on the mutation function and the surrounding points of each grid node;
[0039] Obtain the squared error between the Kriging estimate and the corresponding grid node data, and determine whether the squared error is less than a preset error value;
[0040] If so, it is assumed that there are constraints between the grid node data, and the weighting coefficients are obtained through the Kriging equations.
[0041] The estimated point value corresponding to each grid node is obtained based on the weighting coefficients and the data of each grid node;
[0042] Each estimated point value is interpolated onto the corresponding grid node to obtain the resistivity space grid file;
[0043] The data for each grid node is measured data.
[0044] In a preferred embodiment of the present invention, obtaining the weighting coefficients through the Kriging equations includes:
[0045] Assuming that each estimated point value follows an intrinsic assumption, the Kriging equations are as shown in Equation 2:
[0046] (2)
[0047] In the formula, γ(x) i ,x j ) is the grid node x i and x j The variogram values between them, where μ is the Lagrange daily number;
[0048] The weighting coefficients are obtained according to Formula 2.
[0049] In a preferred embodiment of the present invention, obtaining the estimated point value corresponding to each grid node includes:
[0050] For the regional variation Q(x), let it be at a series of sampling points x i The observed value Q(x) at (i=1,2,.......,n) i (i=1,2,.......,n) Then the estimated value Q(x0) at a certain grid node Xn in the region can be estimated by a linear equation, as shown in Formula 3:
[0051] (3)
[0052] In the formula, λ i These are the weighting coefficients.
[0053] In a preferred embodiment of the present invention, obtaining the median resistivity of any two adjacent grid nodes includes:
[0054] The median resistivity of any two adjacent grid nodes is obtained using the median value acquisition function, which is specifically shown in Formula 4 or Formula 5:
[0055] (4)
[0056] (5)
[0057] In the formula, ρ med The median resistivity of any two adjacent grid nodes, ρ i,j The resistivity of any grid node in the resistivity space grid file;
[0058] When obtaining the normalized derivative operator of the resistivity of any two adjacent grid nodes and the normalized resistivity of each grid node, the following is included:
[0059] The normalized derivative operator is obtained by using the horizontal and vertical cell scales of the grid nodes in the resistivity space grid file, as shown in formula (5):
[0060] (6)
[0061] In the formula, δ p For the normalized derivative operator, Δx is the difference in the horizontal element scale between any two adjacent grid nodes, and Δy is the difference in the vertical element scale between any two adjacent grid nodes;
[0062] The resistivity of each grid node is redefined based on the normalized derivative operator obtained by formula (6), and the normalized resistivity of each grid node is obtained as shown in formula (7):
[0063] (7)
[0064] In the formula, ρ T This is the normalized resistivity.
[0065] Compared with the prior art, the beneficial effects of the present invention are as follows:
[0066] (1) The extraction method provided by the present invention can clearly display the abnormal boundaries in the resistivity cross-sectional diagram and resistivity plane diagram in the results of DC electrical method and electromagnetic method exploration;
[0067] (2) The extraction method provided by the present invention improves the ability to identify geological anomalies and serves as the basis for boundary delineation, thereby providing a basis for the accurate positioning of underground geological bodies.
[0068] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments. Attached Figure Description
[0069] Figure 1 - This is a diagram showing the spatial relationship between the median resistivity and resistivity of two adjacent grid nodes in an embodiment of the present invention.
[0070] Figure 2 - is the original resistivity plane contour map of an embodiment of the present invention;
[0071] Figure 3 - is a normalized resistivity plane contour map of an embodiment of the present invention;
[0072] Figure 4 - This is a step diagram illustrating the method for extracting anomalous boundaries in interpretation maps based on geophysical exploration, according to an embodiment of the present invention. Detailed Implementation
[0073] The method for extracting anomaly boundaries in interpretation maps based on geophysical exploration provided by this invention, such as... Figure 4 As shown, it includes the following steps:
[0074] Step S1: Obtain resistivity data of the survey area, perform inversion on the resistivity data, and obtain inverted resistivity data;
[0075] Step S2: Obtain the resistivity space grid file based on the inverted resistivity data;
[0076] Step S3: Obtain the median resistivity of any two adjacent grid nodes in the resistivity space grid file;
[0077] Step S4: Obtain the normalized derivative operator of the resistivity of any two adjacent grid nodes in the resistivity space grid file;
[0078] Step S5: After redefining the resistivity of each grid node in the resistivity space grid file according to the median value and the normalized derivative operator, the normalized resistivity of each grid node is obtained.
[0079] Step S6: Draw a normalized two-dimensional graph based on the spatial relationship of the normalized resistivity, and clearly display the abnormal boundaries through the normalized two-dimensional graph.
[0080] Specifically, Figure 2 It is a raw resistivity plane contour map obtained using existing geophysical exploration methods. Figure 3 The normalized resistivity plane contour map is obtained using the extraction method provided by this invention. Figure 2 and Figure 3 The comparison shows that the extraction method provided by this invention can clearly display abnormal boundaries.
[0081] In step S2 above, when obtaining the resistivity space grid file based on the inverted resistivity data, the following is included:
[0082] Different detection points are set according to the size of the survey area, and the measured resistivity and measured potential values of different detection points are collected to establish an initial resistivity model.
[0083] The model resistivity within the survey area is obtained based on the initial resistivity model;
[0084] Inversion resistivity data are obtained by inverting the measured resistivity, model resistivity, and initial resistivity model.
[0085] Obtain a resistivity space grid file based on the inverted resistivity data.
[0086] Specifically, the exploration results can be greatly improved by inverting resistivity data.
[0087] Furthermore, when obtaining the inverted resistivity data, the following are included:
[0088] An inversion equation is established based on the measured resistivity, model resistivity, and initial resistivity model, as shown in Formula 1:
[0089] (1);
[0090] In the formula, S is the sensitivity matrix, c is the measured resistivity, ρ is the model resistivity, G is the forward modeling operator, and R... dd Let R be the covariance matrix of the measured resistivity. mm Let ρ0 be the covariance matrix of the model resistivity, and ρ0 be the initial resistivity model. k The resistivity of the model after the k-th iteration;
[0091] The inversion resistivity data are obtained based on the inversion equation.
[0092] Furthermore, when setting different detection points according to the size of the survey area, this includes:
[0093] The first and second setting points are determined based on the size of the surveyed area.
[0094] Install the electrode device at the first setting point and the second setting point;
[0095] Multiple detection lines are set between the first and second setting points according to the preset detection line distance;
[0096] On each testing line, multiple testing points are set according to the preset distance between testing points.
[0097] Furthermore, when establishing the initial resistivity model, the following are included:
[0098] The first detection electrode of the detection device is set at the first detection point of each detection line;
[0099] The second detection electrode of the detection device is set at the second detection point of each detection line;
[0100] Obtain the measured resistivity of the first and second detection points, as well as the measured potential value between the detection points;
[0101] The resistivity between the first and second detection points is obtained based on the measured resistivity and measured potential values.
[0102] An initial resistivity model is established based on the resistivity between the first and second detection points.
[0103] The second detection point is the remaining detection point on each detection line excluding the first detection point, and the second detection electrode is sequentially set on the remaining detection points far away from the first detection point.
[0104] Specifically, the above detection method can accurately obtain the resistivity between the first and second detection points, thereby establishing an effective initial resistivity model.
[0105] In step S2 above, when obtaining the resistivity space grid file based on the inverted resistivity data, the following is also included:
[0106] Based on the distribution of the inverted resistivity data, the region and grid size to be interpolated are determined, the region is gridded, a preliminary two-dimensional grid surface is simulated, and each grid node in the preliminary two-dimensional grid surface is obtained.
[0107] The data of each grid node is checked and analyzed to obtain the distance value between each grid node. If there is a distance value less than the preset distance value, the grid node is removed.
[0108] Select a mutation function, and for each grid node, obtain the Kriging estimate corresponding to each grid node based on the mutation function and the surrounding points of each grid node;
[0109] Obtain the squared error between the Kriging estimate and the corresponding grid node data, and determine whether the squared error is less than the preset error value;
[0110] If so, it is assumed that there are constraints between the grid node data, and the weighting coefficients are obtained through the Kriging equations.
[0111] The estimated point value for each grid node is obtained based on the weighting coefficients and the data of each grid node.
[0112] Each estimated point value is interpolated onto the corresponding grid node to obtain a resistivity space grid file;
[0113] The data for each grid node is actual measured data.
[0114] Specifically, this invention compares the distance values between each grid node with preset distance values, eliminating invalid grid nodes to obtain a more effective preliminary two-dimensional grid surface. Furthermore, by selecting an appropriate variogram, the kriging estimate corresponding to each grid node is accurately obtained. The validity of the obtained kriging estimate is judged based on the squared error between the kriging estimate and the corresponding grid node data. Weighting coefficients are obtained from the valid kriging estimates, and the estimated point value corresponding to each grid node is obtained based on the weighting coefficients. Interpolation is then performed on each grid node to obtain an accurate and effective resistivity space grid file.
[0115] Furthermore, when obtaining the weighting coefficients through the Kriging equations, the following steps are included:
[0116] Assuming that each estimated point value follows an intrinsic assumption, the Kriging equations are as shown in Equation 2:
[0117] (2)
[0118] In the formula, γ(x) i ,x j ) is the grid node x i and x j The variogram values between them, where μ is the Lagrange daily number;
[0119] The weighting coefficients are obtained according to Formula 2.
[0120] Furthermore, when obtaining the estimated point value corresponding to each grid node, the process includes:
[0121] For the regional variation Q(x), let it be at a series of sampling points x i The observed value Q(x) at (i=1,2,.......,n) i (i=1,2,.......,n) Then the estimated value Q(x0) at a certain grid node Xn in the region can be estimated by a linear equation, as shown in Formula 3:
[0122] (3)
[0123] In the formula, λ i These are the weighting coefficients.
[0124] In step S3 above, when obtaining the median resistivity of any two adjacent grid nodes, the following steps are included:
[0125] The median resistivity of any two adjacent grid nodes can be obtained using the median value acquisition function, which is shown in Equation 4 or Equation 5.
[0126] (4)
[0127] (5)
[0128] In the formula, ρ med The median resistivity of any two adjacent grid nodes, ρ i,j The resistivity is the resistivity of any grid node in the resistivity space grid file.
[0129] Specifically, the spatial relationship between the median resistivity of any two adjacent grid nodes and the resistivity is as follows: Figure 1 As shown.
[0130] In steps S4 and S5 above, obtaining the normalized derivative operator of the resistivity of any two adjacent grid nodes and the normalized resistivity of each grid node includes:
[0131] The normalized derivative operator is obtained by using the horizontal and vertical element scales of the grid nodes in the resistivity space grid file, as shown in Equation 6:
[0132] (6)
[0133] In the formula, δ p For the normalized derivative operator, Δx is the difference in the horizontal element scale between any two adjacent grid nodes, and Δy is the difference in the vertical element scale between any two adjacent grid nodes;
[0134] The resistivity of each grid node is redefined based on the normalized derivative operator obtained from Formula 6, resulting in the normalized resistivity of each grid node, as shown in Formula 7:
[0135] (7)
[0136] In the formula, ρ T This is the normalized resistivity.
[0137] Specifically, in geophysical exploration, the results of direct current resistivity and electromagnetic methods can be described using spatial information such as point numbers, line numbers, depth (elevation), and resistivity. Based on their methodological principles and construction methods, maps such as depth-resistivity curves, distance-depth-resistivity cross-sections, and planar position-resistivity slices can be generated. However, these maps cannot accurately display the boundaries of geological anomalies. This invention obtains normalized resistivity using the aforementioned formula 7 and then uses this normalized resistivity to regenerate distance-depth-resistivity cross-section maps or planar position-resistivity slice maps, highlighting the anomaly boundaries.
[0138] Compared with the prior art, the beneficial effects of the present invention are as follows:
[0139] (1) The extraction method provided by the present invention can clearly display the abnormal boundaries in the resistivity cross-sectional diagram and resistivity plane diagram in the results of DC electrical method and electromagnetic method exploration;
[0140] (2) The extraction method provided by the present invention improves the ability to identify geological anomalies and serves as the basis for boundary delineation, thereby providing a basis for the accurate positioning of underground geological bodies.
[0141] The above embodiments are merely preferred embodiments of the present invention and should not be construed as limiting the scope of protection of the present invention. Any non-substantial changes and substitutions made by those skilled in the art based on the present invention shall fall within the scope of protection claimed by the present invention.
Claims
1. A method for extracting anomaly boundaries in interpretation maps based on geophysical exploration, characterized in that, Includes the following steps: Obtain resistivity data of the survey area, and perform inversion on the resistivity data to obtain inverted resistivity data; A resistivity space grid file is obtained based on the inverted resistivity data; Obtain the median value of the resistivity of any two adjacent grid nodes in the resistivity space grid file; Obtain the normalized derivative operator of the resistivity of any two adjacent grid nodes in the resistivity space grid file; After redefining the resistivity of each grid node in the resistivity space grid file based on the median value and the normalized derivative operator, the normalized resistivity of each grid node is obtained. A normalized two-dimensional plot is drawn based on the spatial relationship of the normalized resistivity, and the abnormal boundaries are clearly displayed through the normalized two-dimensional plot; Among them, obtaining the median resistivity of any two adjacent grid nodes includes: The median resistivity of any two adjacent grid nodes is obtained using the median value acquisition function, which is specifically shown in formula (4) or formula (5): (4) (5) In the formula, ρ med The median resistivity of any two adjacent grid nodes, ρ i,j The resistivity of any grid node in the resistivity space grid file; When obtaining the normalized derivative operator of the resistivity of any two adjacent grid nodes and the normalized resistivity of each grid node, the following is included: The normalized derivative operator is obtained by using the horizontal and vertical cell scales of the grid nodes in the resistivity space grid file, as shown in formula (5): (6) In the formula, δ p For the normalized derivative operator, Δx is the difference in the horizontal element scale between any two adjacent grid nodes, and Δy is the difference in the vertical element scale between any two adjacent grid nodes; The resistivity of each grid node is redefined based on the normalized derivative operator obtained by formula (6), and the normalized resistivity of each grid node is obtained as shown in formula (7): (7) In the formula, ρ T Normalized resistivity; When obtaining a resistivity space grid file based on the inverted resistivity data, the following steps are included: Different detection points are set according to the size of the survey area, and the measured resistivity and measured potential values of the different detection points are collected to establish an initial resistivity model. The model resistivity of the survey area is obtained based on the initial resistivity model; The inverted resistivity data is obtained by inverting the measured resistivity, the model resistivity, and the initial resistivity model. The resistivity space grid file is obtained based on the inverted resistivity data; Obtaining the inverted resistivity data includes: An inversion equation is established based on the measured resistivity, the model resistivity, and the initial resistivity model, as shown in formula (1): (1); In the formula, S is the sensitivity matrix, c is the measured resistivity, ρ is the model resistivity, G is the forward modeling operator, and R... dd Let R be the covariance matrix of the measured resistivity. mm Let ρ0 be the covariance matrix of the model resistivity, and ρ0 be the initial resistivity model. k The resistivity of the model at the k-th iteration; The inversion resistivity data are obtained based on the inversion equation.
2. The method for extracting anomaly boundaries in interpretation maps based on geophysical exploration according to claim 1, characterized in that, When setting different detection points according to the size of the survey area, the following are included: The first and second setting points are determined based on the size of the survey area. Install the electrode device at the first setting point and the second setting point; Between the first setting point and the second setting point, multiple detection lines are set according to a preset detection line distance; On each testing line, multiple testing points are set according to the preset distance between testing points.
3. The method for extracting anomaly boundaries in interpretation maps based on geophysical exploration according to claim 2, characterized in that, When establishing the initial resistivity model, the following are included: The first detection electrode of the detection device is set at the first detection point of each detection line; The second detection electrode of the detection device is set at the second detection point of each detection line; Obtain the measured resistivity of the first detection point and the second detection point, as well as the measured potential value between the detection points; The resistivity between the first detection point and the second detection point is obtained based on the measured resistivity and the measured potential value. The initial resistivity model is established based on the resistivity between the first detection point and the second detection point; The second detection point is the remaining detection point on each detection line other than the first detection point, and the second detection electrode is sequentially set on the remaining detection points far away from the first detection point.
4. The method for extracting anomaly boundaries in interpretation maps based on geophysical exploration according to claim 1, characterized in that, When obtaining the resistivity space grid file based on the inverted resistivity data, the process also includes: Based on the distribution of the inverted resistivity data, the region and grid size to be interpolated are determined, the region is gridded to simulate a preliminary two-dimensional grid surface, and each grid node in the preliminary two-dimensional grid surface is obtained. The data of each grid node is checked and analyzed to obtain the distance value between each grid node. If there is a distance value less than a preset distance value, the grid node is removed. A mutation function is selected, and for each grid node, the Kriging estimate corresponding to each grid node is obtained based on the mutation function and the surrounding points of each grid node; Obtain the squared error between the Kriging estimate and the corresponding grid node data, and determine whether the squared error is less than a preset error value; If so, it is assumed that there are constraints between the grid node data, and the weighting coefficients are obtained through the Kriging equations. The estimated point value corresponding to each grid node is obtained based on the weighting coefficients and the data of each grid node; Each estimated point value is interpolated onto the corresponding grid node to obtain the resistivity space grid file; The data for each grid node is measured data.
5. The method for extracting anomaly boundaries in interpretation maps based on geophysical exploration according to claim 4, characterized in that, When obtaining the weighting coefficients through the Kriging equations, the following is included: Assuming that each estimated point value follows an intrinsic hypothesis, the Kriging equations are as shown in equation (2): (2) In the formula, γ(x) i ,x j ) is the grid node x i and x j The variogram values between them, where μ is the Lagrange daily number; The weighting coefficients are obtained according to the formula (2).
6. The method for extracting anomaly boundaries in interpretation maps based on geophysical exploration according to claim 5, characterized in that, When obtaining the estimated point value corresponding to each grid node, the following steps are included: For the regional variation Q(x), let it be at a series of sampling points x i The observed values Q(x) at i=1,2,.......,n i Given that i = 1, 2, ..., n, then a certain grid node X in the region... n The estimated value Q(x0) at point X can be estimated by a linear equation, as shown in formula (3): (3) In the formula, λ i These are the weighting coefficients.