A gravity background field construction method based on iterative optimization kriging

Abstract

The present application relates to a kind of gravity background field construction method based on iterative optimization kriging, comprising the following steps: S1. obtain local gravity measurement data as known point set;S2. based on the known point set, preliminary interpolation is carried out to the grid point of grid background map using kriging method, and the kriging estimated value of each grid point is obtained;S3. by introducing local gravity anomaly covariance and kriging variation function, the random uncertainty of kriging estimated value in interpolation process is evaluated, and each grid point is divided into reference point and to be modified point based on random uncertainty;S4. reference point is combined with known point set to construct continuous distribution constraint surface, and the continuous distribution constraint surface is locally continuously modified based on the obtained reference point, and the continuous distribution constraint surface is iteratively extended based on the locally modified constraint surface until the continuous distribution constraint surface is continuously modified, to realize the construction of gravity background field.

Description

Technical Field

[0001] This invention relates to the field of navigation technology, and in particular to a method for constructing a gravity background field based on iterative optimization kriging. Background Technology

[0002] Gravity matching and gravity-assisted navigation technologies have become important development directions in the field of high-precision navigation due to their significant advantages such as passivity and anti-interference. High-precision, high-resolution gravity background fields are a prerequisite for the implementation and effectiveness of these two technologies. Gravity background fields are usually provided in the form of gravity anomaly fields or gravity gradient fields. This paper selects gravity anomaly fields as the background field, and its construction mainly relies on two types of data sources: one is obtained through satellite measurements, and the other is obtained through on-site measurements such as aerial and shipborne measurements. Currently, satellite measurements have achieved background field coverage in approximately 70% of the global area, with a nominal grid resolution of 1′×1′ (approximately 1.8km). However, due to limitations in data source characteristics and processing methods, the actual effective grid resolution in some areas is only 7-12km, with an accuracy of approximately 1-3mGal in open sea areas and decreasing to 5-8mGal in complex areas such as shallow waters and island / reef areas. Related research indicates that gravity-matched navigation requires a background field containing high-frequency information with wavelengths less than 10 km, with an upper limit to its positioning accuracy approximately √6 / 6 times the grid size. Gravity-assisted navigation requires a background field accuracy of at least 2 mGal to meet the deviation correction requirements of inertial navigation systems. Clearly, the background field constructed by satellite measurements is insufficient in terms of resolution and accuracy to directly support the large-scale application of both navigation technologies. Currently, the accuracy of gravity anomalies measured in-situ can reach 0.1-2 mGal, and the resolution can reach 0.5-2 km, which can meet application requirements. However, data acquisition is constrained by factors such as cost, time, and weather, and sampling points exhibit non-grid distribution, especially in environments such as plateaus and open seas, where data gaps are more pronounced. Therefore, how to convert in-situ measured gravity anomaly data into a regular gravity background field is a common problem that urgently needs to be solved to achieve gravity matching and gravity-assisted navigation.

[0003] Currently, mainstream modeling methods can be divided into two main categories: deterministic interpolation and geostatistical interpolation. Spline function interpolation is a typical example of deterministic interpolation, achieving high-order continuity between data points through piecewise polynomials; common types include cubic splines, B-splines, and natural splines. In the geosciences field, the commonly used GMT software employs a modified tension spline method, introducing a tension parameter to achieve a balance between interpolation smoothness and fidelity to the original data, thereby effectively suppressing boundary oscillations and overshoot. A typical example of geostatistical interpolation is Kriging interpolation, which is widely integrated and applied in mainstream software such as ArcGIS and QGIS. Scholars have made a series of improvements to the classic kriging method for different application scenarios: Universal kriging introduces a deterministic drift term in addition to the random term, usually using a low-order polynomial to characterize the global trend of the study area, so that the model can still give an unbiased estimate when there is a system gradient; Regression kriging combines multiple regression with residual kriging, first using the regression model to explain the influence of large-scale environmental covariates, and then performing local spatial interpolation on the residuals, thereby capturing both complex trends and local variations; Co-kriging is geared towards multivariate scenarios, using cross-variogram functions to quantitatively describe the spatial synergistic relationship between main and auxiliary variables, and even if the main variable sampling is sparse, it can significantly improve the interpolation accuracy and spatial resolution by using high-frequency sampled auxiliary variables. Spline function method generates highly smooth continuous surfaces through simple algorithms, which are sensitive to subtle changes in data, but are fragile to outliers and are prone to oversmoothing and loss of high-frequency information, and cannot provide quantitative accuracy indicators. The advantage of the Kriging method lies in its high accuracy in local spatial modeling, which can provide quantitative accuracy indicators. However, it has limitations such as high computational cost, cumbersome selection of variogram parameters, poor adaptability to non-stationary data, and does not consider the continuity between unknown points. Summary of the Invention

[0004] The technical problem to be solved by the present invention is to provide a method for constructing a gravity background field based on iterative optimization kriging.

[0005] To achieve the above-mentioned objectives, this invention provides a method for constructing a gravity background field based on iterative optimization kriging, comprising the following steps:

[0006] S1. Obtain on-site gravity measurement data as a set of known points;

[0007] S2. Based on the known set of points, the kriging method is used to perform preliminary interpolation on the grid points of the gridded background map to obtain the kriging estimate of each grid point;

[0008] S3. By introducing the local gravity anomaly covariance and the Kriging variogram, the random uncertainty of the Kriging estimate during the interpolation process is evaluated, and the reference point and the point to be corrected are divided into grid points based on the random uncertainty.

[0009] S4. Merge the set of reference points and known points to construct a continuous distributed constraint surface, and perform local continuity correction on the continuous distributed constraint surface based on the obtained reference points. Iterate and expand the corrected local constraint surface until the entire continuous distributed constraint surface is completed and the gravity background field is constructed.

[0010] According to one aspect of the present invention, step S2, which involves performing preliminary interpolation on the grid points of the gridded background image using the kriging method to obtain the kriging estimate of each grid point, includes:

[0011] S21. Define the Kriging estimation model for grid points;

[0012] S22. Construct an objective function for solving the Kriging estimator model based on the Lagrange multiplier method;

[0013] S23. Under the condition of uniform spatial properties, establish a semivariance model;

[0014] S24. The first law of geosciences transforms the semivariogram into a function of distance, which is then used to fit and obtain the corresponding variogram.

[0015] S25. Introduce the fitted variogram into the objective function to solve for the Kriging estimate of the model and complete the initial interpolation of the grid points.

[0016] According to one aspect of the present invention, in step S21, the step of setting the kriging estimation model for the grid points is expressed as follows:

[0017] ;

[0018] in, Indicates the estimated point Kriging's estimate, Represents a known point The value, Indicates the estimated point Given points value The weight it accounts for This represents the number of grid points to be interpolated in the gridded background map, i.e., the estimated points. Quantity;

[0019] In step S22, the objective function for solving the Kriging estimation model is constructed based on the Lagrange multiplier method. The objective function is expressed as:

[0020] ;

[0021] in, Indicates the error variance. express The truth value of , Represents the Lagrange multipliers. Indicates the estimated point Given points value The weight it accounts for Indicates the number of known points;

[0022] In step S23, under the condition of uniform spatial attributes, the step of establishing a semivariance model is expressed as follows:

[0023] ;

[0024] in, Represents a known point With known points The semivariance between them, and , Let (x, y) represent the value of the point (x, y), and , indicating a known point With known points Covariance of values;

[0025] In step S24, the semivariogram is converted into a function of distance using the first law of geosciences, and the corresponding variogram is obtained by fitting the semivariogram represented by distance. The discrete variogram values ​​of the known point pairs formed by the known points are obtained based on the semivariogram represented by distance. The variogram is then fitted based on the obtained discrete variogram values ​​and the selected fitting function model.

[0026] In step S25, the fitted variogram is introduced into the objective function to solve for the Kriging estimate of the Kriging estimation model. This step, which completes the initial interpolation of the grid points, includes:

[0027] Based on the fitted variogram, an objective function is introduced to construct an extended matrix equation, and the weights in the Kriging estimation model are solved based on the extended matrix equation. The extended matrix equation is expressed as:

[0028] ;

[0029] in, The distance between the first known point and the second known point is used to calculate the semivariogram value. The distance between the 1st known point and the nth known point is used to calculate the semivariogram value using the semivariogram function. The distance between the nth known point and the 1st known point is used to calculate the semivariogram value using the semivariogram function. The distance between the nth known point and the nth known point is used to calculate the semivariogram value using the semivariogram function. This represents the weight of the first known point when estimating the k-th grid point. This represents the weight of the nth known point when estimating the kth grid point. This represents the semivariogram value calculated by substituting the distance between the first known point and the kth grid point into the semivariogram function. This represents the semi-variogram value calculated by substituting the distance between the nth known point and the kth grid point into the semi-variogram function;

[0030] Based on the obtained weights By substituting the values ​​into the Kriging estimation model, the Kriging estimate of the Kriging estimation model can be obtained, thus completing the initial interpolation of the grid points.

[0031] According to one aspect of the present invention, in step S3, the random uncertainty index used to evaluate the random uncertainty of the Kriging estimate during the interpolation process by introducing the local gravity anomaly covariance and the Kriging variogram is expressed as follows:

[0032] ;

[0033] in, Indicators representing random uncertainty Represents a known point With the estimated point The distance between them This represents the range in the variogram.

[0034] According to one aspect of the present invention, in step S3, in the step of dividing each grid point into a reference point and a point to be corrected based on random uncertainty, if the random uncertainty index of the grid point is... Less than or equal to the first threshold If the grid point's random uncertainty index is used as a benchmark, then it is classified as a benchmark point. Less than the second threshold If the value is 0, it is classified as a point to be corrected; among them, , , express The mean, subscript Indicates the number of iterations.

[0035] According to one aspect of the present invention, in step S4, the step of merging the reference points and the set of known points to construct a continuous distributed constraint surface, and performing local continuity correction on the continuous distributed constraint surface based on the obtained reference points, involves calculating the correction amount for the adjacent points to be corrected around the obtained reference points using an iterative calculation method, and updating the estimated value of the points to be corrected based on the obtained correction amount to achieve local continuity correction of the continuous distributed constraint surface.

[0036] According to one aspect of the invention, in the step of calculating the correction amount for the neighboring points to be corrected based on the obtained reference point using an iterative calculation method, the correction amount is expressed as:

[0037] ;

[0038] ;

[0039] ;

[0040] in, Indicates the estimated point The correction amount, Indicates the estimated point x-coordinate, Represents a known point x-coordinate, Indicates the estimated point The y-coordinate, Represents a known point The y-coordinate, Represents the Hermite basis matrix. Represents the derivative information matrix. , , , These represent intermediate variables.

[0041] According to one aspect of the present invention, in the step of updating the estimated value of the point to be corrected based on the obtained correction amount to achieve local continuity correction of the continuously distributed constrained surface, the result of updating the estimated value of the point to be corrected based on the correction amount is expressed as follows:

[0042] ;

[0043] ;

[0044] in, This represents the updated estimated value of the point to be corrected. Indicates the weight of the correction value. express The set that constitutes the composition.

[0045] According to one aspect of the present invention, in step S4, the step of merging the reference points and the set of known points to construct a continuous distributed constraint surface, and performing local continuity correction on the continuous distributed constraint surface based on the obtained reference points, and iteratively expanding based on the corrected local constraint surface until the entire continuous distributed constraint surface completes the continuity correction, thereby realizing the construction of the gravity background field, includes:

[0046] S41. Let the known set of points be... Let the set of grid points of the gridded background map interpolated using the Kriging method be... And, constructing a set of estimated points to be partitioned for dynamic updating. and the reference point set used to store reference points ;

[0047] S42. Initialize the point set, where, Initialize to , Initialize to ;

[0048] S43. Solving for the estimated point set to be partitioned using the Kriging method. Middle Estimation Point The Kriging estimate is obtained, and the stochastic uncertainty of the Kriging estimate is assessed.

[0049] S44. Based on random uncertainty, divide each grid point into reference points and points to be corrected, and then perform a set of reference points based on the reference points. Update and determine the reference point set. Does it cover the set of grid points? If not, then based on the known point set and reference point set Construct the first set of fusion points And it is represented as: ;

[0050] S45. Based on the first fusion point set Extract the reference point set separately and the set of points to be corrected ;

[0051] S46. Based on the first fusion point set Extract the reference points and use iterative calculation to calculate the set of points to be corrected. The correction amount is calculated for the adjacent points to be corrected, and the estimated value of the points to be corrected is updated based on the obtained correction amount;

[0052] S47. Let the set of points to be corrected be completed. For the reference point set and estimated point set Update them separately;

[0053] S48. Repeat steps S43 to S47 to progressively expand the local constraint surface until the set of reference points is reached. Covering grid point set Output complete background field modeling results.

[0054] According to one aspect of the invention, in step S45, based on the first fusion point set... Extract the reference point set separately and the set of points to be corrected In the steps, the reference point set is adopted. Represented as:

[0055] ;

[0056] Point set to be corrected Represented as:

[0057] ;

[0058] In step S47, the set of points to be corrected is completed. For the reference point set and estimated point set In the separate update steps,

[0059] Updated reference point set Represented as:

[0060] ;

[0061] Updated estimated point set Represented as:

[0062] .

[0063] According to one aspect of the present invention, the method first combines the advantages of Kriging in local spatial correlation modeling and quantitative error assessment to construct an initial local constraint surface; then, a spline function is introduced to continuously correct the estimated points, improving the smoothness and physical consistency of the surface in the transition region; further, the correction results are used to optimize the variogram parameters, realizing the spatial adaptive update of the local structure, and gradually expanding the constraint surface to the entire modeling region, ultimately completing the high-precision construction of the background field.

[0064] According to one embodiment of the present invention, the proposed method has been verified by the EIGEN-6C4 model and airborne gravity measurement data. The RMSE of the background field established by this method is better than that of the single Kriging method and the spline function method. In particular, it shows stronger stability and adaptability in background field modeling with many local extrema, providing a new approach to improve the high-precision construction of regional gravity background fields.

[0065] According to one aspect of the present invention, this approach outperforms traditional single methods in terms of modeling accuracy and stability, and effectively avoids the problems of spurious oscillations and missing results in the spline function method. Attached Figure Description

[0066] Figure 1 This is a step diagram of the gravity background field construction method based on iterative optimization kriging of the present invention;

[0067] Figure 2 The flowchart shows the gravity background field construction method based on iterative optimization kriging of the present invention.

[0068] Figure 3 This is a gravity background image of different regions in Embodiment 1 of the present invention, wherein, Figure 3 (a) shows the gravity background map of region A1. Figure 3 (b) shows the gravity background map of region A2. Figure 3 (c) shows the gravity background map of region A3. Figure 3 (d) shows the gravity background map of region A4. Figure 3 (e) shows the gravity background map of region A5;

[0069] Figure 4 This is a distribution diagram of the survey lines and intersections in Embodiment 2 of the present invention. Detailed Implementation

[0070] The present invention will now be described in detail with reference to the accompanying drawings and specific embodiments. The embodiments cannot be described in detail here, but the embodiments of the present invention are not limited to the following embodiments.

[0071] Combination Figure 1 and Figure 2 As shown, according to one embodiment of the present invention, a method for constructing a gravity background field based on iterative optimization kriging includes the following steps:

[0072] S1. Obtain on-site gravity measurement data as a set of known points;

[0073] S2. Based on the known set of points, the kriging method is used to perform preliminary interpolation on the grid points of the gridded background map to obtain the kriging estimate of each grid point;

[0074] S3. By introducing the local gravity anomaly covariance and the Kriging variogram, the random uncertainty of the Kriging estimate during the interpolation process is evaluated, and the reference point and the point to be corrected are divided into grid points based on the random uncertainty.

[0075] S4. Merge the set of reference points and known points to construct a continuous distributed constraint surface, and perform local continuity correction on the continuous distributed constraint surface based on the obtained reference points. Iterate and expand the corrected local constraint surface until the entire continuous distributed constraint surface is completed and the gravity background field is constructed.

[0076] According to one embodiment of the present invention, in step S1, the step of acquiring local gravity measurement data as a set of known points, the set of known points... Represented as:

[0077] ;

[0078] in, This represents the known points of the local gravity measurement data. , Represents a known point coordinates The number of known points; in this embodiment, The corresponding value is known, and it is used express.

[0079] Combination Figure 1 and Figure 2 As shown, according to one embodiment of the present invention, step S2, which involves using the kriging method to perform preliminary interpolation on the grid points of the gridded background map to obtain the kriging estimate of each grid point, includes:

[0080] S21. Define the Kriging estimation model for grid points;

[0081] S22. Construct an objective function for solving the Kriging estimator model based on the Lagrange multiplier method;

[0082] S23. Under the condition of uniform spatial properties, establish a semivariance model;

[0083] S24. The first law of geosciences transforms the semivariogram into a function of distance, which is then used to fit and obtain the corresponding variogram.

[0084] S25. Introduce the fitted variogram into the objective function to solve for the Kriging estimate of the model and complete the initial interpolation of the grid points.

[0085] According to one embodiment of the present invention, in step S21, the step of setting the Kriging estimation model for grid points is expressed as follows:

[0086] ;

[0087] in, Indicates the estimated point Kriging's estimate, Represents a known point The value, Indicates the estimated point Given points value The weight it accounts for This represents the number of grid points to be interpolated in the gridded background map, i.e., the estimated points. Quantity;

[0088] Furthermore, in step S22, the step of constructing the objective function for solving the Kriging estimation model based on the Lagrange multiplier method, uses the constructed objective function to obtain information about the estimation point. The set of weights whose estimated variance is minimized; where the objective function is expressed as:

[0089] ;

[0090] in, Indicates the error variance. express The truth value of , Represents the Lagrange multipliers. Indicates the estimated point Given points value The weight it accounts for Indicates the number of known points;

[0091] Furthermore, in step S23, under the condition of uniform spatial attributes, in the step of establishing the semivariance model, it is assumed that the spatial attributes are uniform. , Therefore, the semivariance model can be defined as follows:

[0092] ;

[0093] in, Expressing expectations, , Represents a known point in space With known points The spatial properties exhibit random fluctuations.

[0094] Therefore, the semivariance model can be expressed as:

[0095] ;

[0096] in, Represents a known point With known points The semivariance between them, and , Let (x, y) represent the value of the point (x, y), and , indicating a known point With known points Covariance of values;

[0097] Furthermore, in step S24, where the semivariogram is converted into a function of distance using the first law of geosciences, and used to fit and obtain the corresponding variogram, the discrete variogram values ​​of known point pairs formed by known points are obtained based on the semivariogram expressed as distance. The variogram is then fitted using the obtained discrete variogram values ​​and the selected fitting function model. Specifically, according to the first law of geosciences, the semivariogram can be expressed as the distance between known points. With known points Distance between The function, therefore, for the mutated function, can be obtained from The fitting yields the result, where the variogram function can be defined as follows: .

[0098] Furthermore, since the variogram obtained through fitting is a discrete value, it is necessary to further fit the discrete results to construct the variogram; specifically, the variogram is... The variogram is composed of three main parameters: sill, nugget, and range. The fitting function model used can be a spherical function model. Therefore, the general expression of the variogram is:

[0099] ;

[0100] in, Indicates the value of the nugget. Indicates the base value. Represents distance, used for substituting known points. With known points Distance between , This represents the range in the variogram.

[0101] Furthermore, in step S25, the fitted variogram is introduced into the objective function to solve for the Kriging estimate of the Kriging estimation model. This step, which completes the initial interpolation of the grid points, includes:

[0102] Based on the fitted variogram, an objective function is introduced to construct an extended matrix equation, and the weights in the Kriging estimation model are solved based on the extended matrix equation. The extended matrix equation is expressed as:

[0103] ;

[0104] in, The distance between the first known point and the second known point is used to calculate the semivariogram value. The distance between the 1st known point and the nth known point is used to calculate the semivariogram value using the semivariogram function. The distance between the nth known point and the 1st known point is used to calculate the semivariogram value using the semivariogram function. The distance between the nth known point and the nth known point is used to calculate the semivariogram value using the semivariogram function. This represents the weight of the first known point when estimating the k-th grid point. This represents the weight of the nth known point when estimating the kth grid point. This represents the semivariogram value calculated by substituting the distance between the first known point and the kth grid point into the semivariogram function. This represents the semi-variogram value calculated by substituting the distance between the nth known point and the kth grid point into the semi-variogram function;

[0105] Based on the obtained weights By substituting the values ​​into the Kriging estimation model, the Kriging estimate of the Kriging estimation model can be obtained, thus completing the initial interpolation of the grid points.

[0106] According to one embodiment of the present invention, in step S3, the random uncertainty index used to evaluate the random uncertainty of the Kriging estimate during the interpolation process by introducing the local gravity anomaly covariance and the Kriging variogram is expressed as follows:

[0107] ;

[0108] in, Indicators representing random uncertainty Represents a known point With the estimated point The distance between them This represents the range in the variogram.

[0109] In this embodiment, the process of constructing the stochastic uncertainty index is as follows:

[0110] The random uncertainty of the Kriging interpolation result can be assessed by the Kriging variance, and can be expressed as:

[0111] ;

[0112] in, Indicates the estimated point The estimated value of variance. Indicates distance The value of the mutation function.

[0113] Under the simplified model, It can be approximated as the distance Given the relevant first-order piecewise linear function, the Kriging variance can be viewed as... The relevant first-order linear combination. Since the error propagation relationship of gravity anomalies with distance in space exhibits a significant nonlinearity, the Kriging variance derived purely mathematically cannot accurately assess random uncertainty.

[0114] Furthermore, by using an isotropic empirical spatial covariance function to characterize the spatial correlation structure of regionalized variables, the empirical spatial covariance function is expressed as:

[0115] ;

[0116] in, This represents the prior variance within the study area and can be considered a constant. The model parameter used to control the correlation decay rate is its reciprocal. Indicates practical range variation. The larger the spatial correlation, the slower the decay and the slower the spread of random uncertainty. When the distance between two points approaches zero, the positive correlation is strongest, and the random uncertainty approaches zero; as the distance increases, the positive correlation weakens, and the random uncertainty increases accordingly.

[0117] Based on this, As a measure of stochastic uncertainty, the prior variance is treated as a constant. Combined with the weighting relationships obtained through the Kriging method, and based on the form of the empirical space covariance function, a stochastic uncertainty index is further derived. And it is expressed as:

[0118] .

[0119] According to one embodiment of the present invention, in step S3, in the step of dividing each grid point into a reference point and a point to be corrected based on random uncertainty, if the random uncertainty index of the grid point... Less than or equal to the first threshold If the grid point's random uncertainty index is used as a benchmark, then it is classified as a benchmark point. Less than the second threshold If the value is 0, it is classified as a point to be corrected; among them, , , express The mean, subscript Indicates the number of iterations.

[0120] According to one embodiment of the present invention, in step S4, the step of merging the reference points and the set of known points to construct a continuous distributed constraint surface, and performing local continuity correction on the continuous distributed constraint surface based on the obtained reference points, is to calculate the correction amount of the adjacent points to be corrected around the obtained reference points using an iterative calculation method, and update the estimated value of the points to be corrected based on the obtained correction amount to realize the local continuity correction of the continuous distributed constraint surface.

[0121] In this embodiment, in the step of calculating the correction amount for the neighboring points to be corrected based on the obtained reference point using an iterative calculation method, the correction amount is expressed as:

[0122] ;

[0123] ;

[0124] ;

[0125] in, Indicates the estimated point The correction amount, Indicates the estimated point x-coordinate, Represents a known point x-coordinate, Indicates the estimated point The y-coordinate, Represents a known point The y-coordinate, Represents the Hermite basis matrix. Represents the derivative information matrix. , , , These represent intermediate variables.

[0126] In this embodiment, the derivative information matrix The partial differential value is from The derivative information matrix is ​​obtained by difference calculus. Represented as:

[0127] ;

[0128] ;

[0129] ;

[0130] ;

[0131] .

[0132] According to one embodiment of the present invention, in the step of updating the estimated value of the point to be corrected based on the obtained correction amount to achieve local continuity correction of the continuously distributed constrained surface, the result of updating the estimated value of the point to be corrected based on the correction amount is expressed as follows:

[0133] ;

[0134] ;

[0135] in, This represents the updated estimated value of the point to be corrected. Indicates the weight of the correction value. express The set that constitutes the composition.

[0136] According to one embodiment of the present invention, in step S4, the steps of merging the reference points and the set of known points to construct a continuous distributed constraint surface, and performing local continuity correction on the continuous distributed constraint surface based on the obtained reference points, and iteratively expanding based on the corrected local constraint surface until the entire continuous distributed constraint surface completes the continuity correction to realize the construction of the gravity background field, include:

[0137] S41. Let the known set of points be... Let the set of grid points of the gridded background map interpolated using the Kriging method be... And, constructing a set of estimated points to be partitioned for dynamic updating. and the reference point set used to store reference points In this embodiment, the set of grid points is: Represented as:

[0138] ;

[0139] in, The number of grid points; in this embodiment, the grid point set The values ​​of the grid points can be unknown, or they can be reference values ​​used for verification. If reference values ​​exist, they are applicable to the application. express.

[0140] S42. Initialize the point set, where, Initialize to , Initialize to ;

[0141] S43. Solving for the estimated point set to be partitioned using the Kriging method. Middle Estimation Point The Kriging estimate is obtained, and the stochastic uncertainty of the Kriging estimate is assessed.

[0142] S44. Based on random uncertainty, divide each grid point into reference points and points to be corrected, and then perform a set of reference points based on the reference points. Update and determine the reference point set. Does it cover the set of grid points? If not, then based on the known point set and reference point set Construct the first set of fusion points And it is represented as: ;

[0143] S45. Based on the first fusion point set Extract the reference point set separately and the set of points to be corrected Among them, the set of reference points. Represented as:

[0144] ;

[0145] Point set to be corrected Represented as:

[0146] ;

[0147] S46. Based on the first fusion point set Extract the reference points and use iterative calculation to calculate the set of points to be corrected. The correction amount is calculated for the adjacent points to be corrected, and the estimated value of the points to be corrected is updated based on the obtained correction amount;

[0148] S47. Let the set of points to be corrected be completed. For the reference point set and estimated point set Update them separately; among them, the updated set of reference points Represented as:

[0149] ;

[0150] Updated estimated point set Represented as:

[0151] .

[0152] S48. Repeat steps S43 to S47 to progressively expand the local constraint surface until the set of reference points is reached. Covering grid point set Output complete background field modeling results.

[0153] According to one embodiment of the present invention, the gravity background field construction method based on iterative optimization kriging can be implemented using a gravity background field construction device. Specific limitations on the gravity background field construction device can be found in the above-described limitations of the gravity background field construction method based on iterative optimization kriging, and will not be repeated here. Furthermore, the various modules involved in the gravity background field construction device based on iterative optimization kriging can be implemented entirely or partially through software, hardware, or a combination thereof. These modules can be embedded in or independent of the processor in a computer device in hardware form, or stored in the memory of a computer device in software form, so that the processor can call and execute the operations corresponding to each module.

[0154] In this embodiment, the memory may be, but is not limited to, Random Access Memory (RAM), Read Only Memory (ROM), Programmable Read-Only Memory (PROM), Erasable Programmable Read-Only Memory (EPROM), Electrically Erasable Programmable Read-Only Memory (EEPROM), etc.

[0155] In this embodiment, the processor can be an integrated circuit chip with signal processing capabilities. The processor can be a general-purpose processor, including a central processing unit (CPU), a network processor (NP), etc.; it can also be a digital signal processor (DSP), an application-specific integrated circuit (ASIC), a field-programmable gate array (FPGA), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components.

[0156] To further illustrate this plan, further examples will be provided.

[0157] Example 1: EIGEN-6C4 model test

[0158] Data source

[0159] To systematically verify the adaptability of the proposed method under different terrain conditions, five regions with different characteristics were selected: desert (A1), plateau (A2), karst mountains (A3), sea area (A4), and plain (A5). A 2°×2° rectangular area was selected as the study area for each region.

[0160] Table 1 shows the statistical information on gravity anomalies in the study area. From the statistical characteristics of gravity anomalies, the extreme values ​​of gravity anomalies in regions A2 and A3 have a large range, reflecting the significant disturbance effect of the strong crustal uplift in region A2 and the complex lithological distribution in region A3 on the gravity field; regions A1, A4 and A5, on the other hand, show a weak and stable fluctuation characteristic.

[0161] Table 1. Statistical information on gravity anomalies in the study area

[0162]

[0163] The experimental data were obtained from the calculation service of the International Centre for Global Gravity Models (ICGEM) website, using the EIGEN-6C4 gravity field model. The data grid resolution was set to 0.5′×0.5′, and 58,081 (241×241) gravity anomaly reference points were obtained for each 2°×2° study area, forming the ground truth dataset of the gravity anomaly background field.

[0164] The distribution of gravity anomalies in the study area is shown in [reference]. Figure 3 From the spatial characteristics of gravity anomalies, region A1 shows a large-scale, gradually changing regional anomaly, regions A2 and A3 show complex local disturbances, and regions A4 and A5 are a superposition of gradually changing anomalies and local disturbances.

[0165] To simulate the actual gravity measurement scenario, the sampling strategy followed the network layout requirements in the "Technical Specification for Airborne Gravity Measurement (DZ / T 0381-2021)". The spacing between survey lines was 500m to 2000m, and the spacing between control lines was 10 times the spacing between survey lines. Calculations showed that within a 2°×2° study area, the number of intersections was approximately 20,000 when the survey line spacing was 500m and approximately 1,200 when the spacing was 2000m. Based on this, 1200, 2500, 5000, 8000, 10000, 15000, and 20000 points were randomly sampled from the ground truth points in each study area to construct a simulated airborne gravity measurement dataset. This dataset served as known points for constructing the background field and was used for method validation.

[0166] Results Comparison

[0167] Table 2 presents the RMSE statistics of this model experiment. The data shows that when the number of sampling points is greater than 5000 (corresponding to a survey line spacing of approximately 1 km), this scheme demonstrates a stable error reduction effect in all five study areas. Regardless of whether facing large-scale gradual changes, local strong disturbances, or background superimposed anomalies, this scheme significantly outperforms single methods in terms of modeling accuracy and spatial feature fidelity. Especially in scenarios with complex anomalies and variable sampling densities, it can more accurately reconstruct the spatial distribution pattern of gravity anomalies, indicating that this method can effectively improve the accuracy of constructing the gravity background field from on-site measurement data.

[0168] Table 2 RMSE statistics of model experiment

[0169]

[0170] Example 2: Application of Airborne Gravity Measurement Data

[0171] Data source

[0172] The measured data comes from an airborne gravity survey conducted in a certain location, with a longitude range of approximately 0.4° and a latitude range of approximately 0.35°. Data was collected using a SGAWZ03 strapdown airborne gravimeter, yielding a total of 71 survey lines: 63 east-west lines with a spacing of approximately 617.7m, and 8 north-south lines with a spacing of approximately 5559.7m. A total of 504 intersection points were obtained, with an elevation difference of no more than 22m between intersection points. The layout of the survey lines and intersection points conforms to the specifications, and the intersection point accuracy is better than 1.06mGal. The distribution of survey lines and points is shown in [reference needed]. Figure 4 .

[0173] The experiment used a 10-fold cross-validation method, randomly dividing the data into 10 groups. Nine groups were used as known points to construct the background field, and one group was used as the unknown points to verify the calculation of RMSE.

[0174] Results Analysis

[0175] Table 3 shows the RMSE statistics for 10-fold cross-validation. The data shows that this scheme has the lowest RMSE in cross-validation, with an average RMSE of only 0.6842 mGal, which is about 2.8% lower than the Kriging method and 44.2% lower than the spline function method, indicating the applicability of this scheme in constructing gravity background fields from measured data.

[0176] Table 3. RMSE statistics for 10-fold cross-validation

[0177]

[0178] Therefore, compared with the Kriging method and the spline method, this method reduces the RMSE by an average of 18.7% and 27.3% in the model test and by an average of 2.8% and 42.2% in the experimental data.

[0179] This scheme constructs a local spatial correlation constraint based on the variogram using the Kriging method, quantifies the confidence boundary of the interpolation results, and introduces the global smoothness property of spline interpolation to correct local surface deviations. Through iterative optimization, it achieves a deep integration of the advantages of the two methods, effectively reducing the error in the construction of the gravity background field.

[0180] The above description is merely an example of a specific solution of the present invention. For any devices and structures not described in detail herein, it should be understood that they are implemented using common devices and methods already available in the art.

[0181] The above description is merely one embodiment of the present invention and is not intended to limit the invention. Various modifications and variations can be made to the invention by those skilled in the art. Any modifications, equivalent substitutions, or improvements made within the spirit and principles of the present invention should be included within the scope of protection of the present invention.

Claims

1. A method for constructing a gravity background field based on iterative optimization kriging, characterized in that, Includes the following steps: S1. Obtain on-site gravity measurement data as a set of known points; S2. Based on the known set of points, the kriging method is used to perform preliminary interpolation on the grid points of the gridded background map to obtain the kriging estimate of each grid point; S3. By introducing the local gravity anomaly covariance and the Kriging variogram, the random uncertainty of the Kriging estimate during the interpolation process is evaluated, and the reference point and the point to be corrected are divided into grid points based on the random uncertainty; wherein, the random uncertainty index used to evaluate the random uncertainty of the Kriging estimate during the interpolation process is expressed as: in, Indicators representing random uncertainty Represents a known point With the estimated point The distance between them This represents the range in the variogram; If the random uncertainty index of grid points Less than or equal to the first threshold If the grid point's random uncertainty index is used as a benchmark, then it is classified as a benchmark point. Less than the second threshold If the value is 0, it is classified as a point to be corrected; among them, , , express The mean, subscript Indicates the number of iterations; S4. Merge the set of reference points and known points to construct a continuous distributed constraint surface, and perform local continuity correction on the continuous distributed constraint surface based on the obtained reference points. Iterate and expand the corrected local constraint surface until the entire continuous distributed constraint surface is completed and the gravity background field is constructed.

2. The method for constructing a gravity background field based on iterative optimization kriging according to claim 1, characterized in that, Step S2, which involves using the kriging method to perform preliminary interpolation on the grid points of the gridded background map to obtain the kriging estimate of each grid point, includes: S21. Define the Kriging estimation model for grid points; S22. Construct an objective function for solving the Kriging estimator model based on the Lagrange multiplier method; S23. Under the condition of uniform spatial properties, establish a semivariance model; S24. The first law of geosciences transforms the semivariogram into a function of distance, which is then used to fit and obtain the corresponding variogram. S25. Introduce the fitted variogram into the objective function to solve for the Kriging estimate of the model and complete the initial interpolation of the grid points.

3. The method for constructing a gravity background field based on iterative optimization kriging according to claim 2, characterized in that, In step S21, the step of setting the kriging estimation model for grid points is expressed as follows: in, Indicates the estimated point Kriging's estimate, Represents a known point The value, Indicates the estimated point Given points value The weight it accounts for This represents the number of grid points to be interpolated in the gridded background map, i.e., the estimated points. Quantity; In step S22, the objective function for solving the Kriging estimation model is constructed based on the Lagrange multiplier method. The objective function is expressed as: in, Indicates the error variance. express The truth value of , Represents the Lagrange multipliers. Indicates the estimated point Given points value The weight it accounts for Indicates the number of known points; In step S23, under the condition of uniform spatial attributes, the step of establishing a semivariance model is expressed as follows: in, Represents a known point With known points The semivariance between them, and , Let (x, y) represent the value of the point (x, y), and , indicating a known point With known points Covariance of values; In step S24, the semivariogram is converted into a function of distance using the first law of geosciences, and the corresponding variogram is obtained by fitting the semivariogram represented by distance. The discrete variogram values ​​of the known point pairs formed by the known points are obtained based on the semivariogram represented by distance. The variogram is then fitted based on the obtained discrete variogram values ​​and the selected fitting function model. In step S25, the fitted variogram is introduced into the objective function to solve for the Kriging estimate of the Kriging estimation model. This step, which completes the initial interpolation of the grid points, includes: Based on the fitted variogram, an objective function is introduced to construct an extended matrix equation, and the weights in the Kriging estimation model are solved based on the extended matrix equation. The extended matrix equation is expressed as: in, The distance between the first known point and the second known point is used to calculate the semivariogram value. The distance between the 1st known point and the nth known point is used to calculate the semivariogram value using the semivariogram function. The distance between the nth known point and the 1st known point is used to calculate the semivariogram value using the semivariogram function. The distance between the nth known point and the nth known point is used to calculate the semivariogram value using the semivariogram function. This represents the weight of the first known point when estimating the k-th grid point. This represents the weight of the nth known point when estimating the kth grid point. This represents the semivariogram value calculated by substituting the distance between the first known point and the kth grid point into the semivariogram function. This represents the semi-variogram value calculated by substituting the distance between the nth known point and the kth grid point into the semi-variogram function; Based on the obtained weights By substituting the values ​​into the Kriging estimation model, the Kriging estimate of the Kriging estimation model can be obtained, thus completing the initial interpolation of the grid points.

4. The method for constructing a gravity background field based on iterative optimization kriging according to any one of claims 1 to 3, characterized in that, In step S4, the process of merging the reference points and the set of known points to construct a continuous distributed constraint surface, and performing local continuity correction on the continuous distributed constraint surface based on the obtained reference points, involves calculating the correction amount for the adjacent points to be corrected around the obtained reference points using an iterative calculation method, and updating the estimated value of the points to be corrected based on the obtained correction amount to achieve local continuity correction of the continuous distributed constraint surface.

5. The method for constructing a gravity background field based on iterative optimization kriging according to claim 4, characterized in that, In the step of calculating the correction amount for the adjacent points to be corrected based on the obtained reference point using an iterative calculation method, the correction amount is expressed as: in, Indicates the estimated point The correction amount, Indicates the estimated point x-coordinate, Represents a known point x-coordinate, Indicates the estimated point The y-coordinate, Represents a known point The y-coordinate, Represents the Hermite basis matrix. Represents the derivative information matrix. , , , These represent intermediate variables.

6. The method for constructing a gravity background field based on iterative optimization kriging according to claim 5, characterized in that, In the step of updating the estimated value of the point to be corrected based on the obtained correction amount to achieve local continuity correction of the continuously distributed constrained surface, the result of updating the estimated value of the point to be corrected based on the correction amount is expressed as follows: in, This represents the updated estimated value of the point to be corrected. Indicates the weight of the correction value. express The set that constitutes the composition.

7. The method for constructing a gravity background field based on iterative optimization kriging according to claim 6, characterized in that, In step S4, the steps of merging the reference points and the set of known points to construct a continuous distributed constraint surface, and performing local continuity correction on the continuous distributed constraint surface based on the obtained reference points, and iteratively expanding the corrected local constraint surface until the entire continuous distributed constraint surface is continuously corrected to realize the construction of the gravity background field, include: S41. Let the known set of points be... Let the set of grid points of the gridded background map interpolated using the Kriging method be... And, constructing a set of estimated points to be partitioned for dynamic updating. and the reference point set used to store reference points ; S42. Initialize the point set, where, Initialize to , Initialize to ; S43. Solving for the estimated point set to be partitioned using the Kriging method. Middle Estimation Point The Kriging estimate is obtained, and the stochastic uncertainty of the Kriging estimate is assessed. S44. Based on random uncertainty, divide each grid point into reference points and points to be corrected, and then perform a set of reference points based on the reference points. Update and determine the reference point set. Does it cover the set of grid points? If not, then based on the known point set and reference point set Construct the first set of fusion points And it is represented as: ; S45. Based on the first fusion point set Extract the reference point set separately and the set of points to be corrected ; S46. Based on the first fusion point set Extract the reference points and use iterative calculation to calculate the set of points to be corrected. The correction amount is calculated for the adjacent points to be corrected, and the estimated value of the points to be corrected is updated based on the obtained correction amount; S47. Let the set of points to be corrected be completed. For the reference point set and estimated point set Update them separately; S48. Repeat steps S43 to S47 to progressively expand the local constraint surface until the set of reference points is reached. Covering grid point set Output complete background field modeling results.

8. The method for constructing a gravity background field based on iterative optimization kriging according to claim 7, characterized in that, In step S45, based on the first fusion point set Extract the reference point set separately and the set of points to be corrected In the steps, the reference point set is adopted. Represented as: Point set to be corrected Represented as: ； In step S47, the set of points to be corrected is completed. For the reference point set and estimated point set In the separate update steps, Updated reference point set Represented as: ； Updated estimated point set Represented as: 。

Images