GNSS (Global Navigation Satellite System) earth crust strain field solving method based on gradient guided radial basis function

Abstract

The invention discloses a GNSS (Global Navigation Satellite System) earth crust strain field solving method based on a gradient-guided radial basis function, which belongs to the field of geophysics and space geodetic survey, and comprises the following steps: acquiring GNSS observation data of a target area and converting the GNSS observation data to a plane projection coordinate system to obtain the plane position of each observation station; estimating a local velocity gradient based on the plane position and the horizontal velocity component of each observation station; constructing an anisotropic constraint matrix based on the local velocity gradient to define a non-Euclidean distance metric; establishing a radial basis function interpolation model based on the non-Euclidean distance measurement and solving an interpolation coefficient to obtain a continuous velocity field; calculating a strain tensor component based on the continuous velocity field; and constructing a strain field based on the strain tensor component, and carrying out gridding and visualization processing on the strain field.

Description

Technical Field

[0001] This invention belongs to the field of geophysics and space geodesy, and particularly relates to a method for solving the crustal strain field in GNSS based on gradient-guided radial basis functions. Background Technology

[0002] In modern crustal deformation and seismic tectonic research, GNSS observations provide high-precision three-dimensional surface velocity information. However, due to factors such as uneven distribution of observation stations and sparse data sampling, the actual observed velocity field is discrete and irregularly distributed. To analyze regional crustal deformation patterns, it is necessary to convert these discrete point data into a continuous velocity field and further calculate the strain tensor field.

[0003] Currently, commonly used velocity-strain field inversion methods include least squares collocation, polynomial least squares fitting, locally linear or polynomial fitting, finite element method, and Bayesian or regularized inversion methods. While each method has its own advantages, they also have their limitations. Radial Basis Function (RBF) interpolation is an interpolation method introduced into the geosciences in recent years. This method constructs a kernel function using distance as a variable for smooth interpolation, exhibiting good continuity and differentiability. However, traditional RBF assumes spatial isotropy and does not consider the directional gradient information of the velocity field, making it unable to accurately characterize strain changes in regions with significant directionality, such as faults, leading to unstable strain estimation in anisotropic deformation regions. Therefore, it is necessary to propose an interpolation method that can maintain the smoothness and differentiability of RBF while reflecting the directional characteristics of the velocity field. Summary of the Invention

[0004] To address the aforementioned technical problems, this invention proposes a GNSS crustal strain field solution method based on gradient-guided radial basis functions, thereby resolving the issues present in the prior art.

[0005] To achieve the above objectives, this invention provides a method for solving the GNSS crustal strain field based on gradient-guided radial basis functions, comprising:

[0006] Acquire GNSS observation data of the target area and convert it to a plane projection coordinate system to obtain the plane position of each station;

[0007] Local velocity gradients are estimated based on the planar position and horizontal velocity components of each station.

[0008] An anisotropic constraint matrix is ​​constructed based on the local velocity gradient to define a non-Euclidean distance metric;

[0009] A radial basis function interpolation model is established based on the non-Euclidean distance metric, and the interpolation coefficients are solved to obtain a continuous velocity field.

[0010] Calculate the strain tensor components based on the continuous velocity field;

[0011] A strain field is constructed based on the strain tensor components, and the strain field is then meshed and visualized.

[0012] Optionally, the process of estimating the local velocity gradient includes: using a weighted linear fitting method, the local velocity gradient is obtained by solving the spatial derivative of the horizontal velocity component in the neighborhood of each observation point through least squares.

[0013] Optionally, the weights used in the weighted linear fitting are Gaussian weights, and the Euclidean distance between the observation points is processed by a Gaussian function to obtain weight values ​​that decay with distance.

[0014] The smoothing length parameter of the Gaussian function is taken as 0.5 to 1 times the average distance between measuring points in the target area.

[0015] Optionally, the process of constructing the anisotropic constraint matrix includes:

[0016] Calculation of symmetric strain rate tensor based on local velocity gradient;

[0017] Eigendecomposition of the symmetric strain rate tensor yields the principal variation direction and its corresponding eigenvalues;

[0018] Calculate the anisotropic scaling factor based on the eigenvalues;

[0019] Construct a rotation matrix based on the main direction of change;

[0020] Anisotropic constraint matrices are synthesized based on the rotation matrix and the anisotropic scaling factor.

[0021] Optionally, the anisotropic scaling factor is calculated by the ratio of the adjustment factor to the sum of the absolute values ​​of the eigenvalues, wherein the adjustment factor ranges from 0.2 to 1.0.

[0022] Optionally, the process of establishing a radial basis function interpolation model based on the non-Euclidean distance metric and solving for the interpolation coefficients to obtain the continuous velocity field includes:

[0023] The non-Euclidean distance metric is processed using a Gaussian kernel function to obtain the radial basis function relation matrix;

[0024] Based on the relationship matrix and the observed velocity values, a system of linear equations with polynomial constraints is constructed.

[0025] Solving the linear equations yields the radial basis function weighting coefficients and polynomial coefficients;

[0026] The expression for the continuous velocity field is obtained based on the weighting coefficients, polynomial coefficients, and radial basis functions.

[0027] Optionally, calculating the strain tensor components includes:

[0028] The spatial partial derivatives of the velocity field are analytically obtained based on the radial basis function interpolation model.

[0029] The components of the strain tensor, the maximum shear strain, the surface dilatation rate, and the principal strain rates are calculated using the spatial partial derivatives.

[0030] Compared with the prior art, the present invention has the following advantages and technical effects:

[0031] This invention constructs an anisotropic constraint matrix by introducing local velocity gradient information, thus improving the traditional isotropic radial basis function interpolation into anisotropic interpolation. This method can adaptively adjust the interpolation kernel shape according to the characteristics of velocity field changes, using a finer interpolation scale in directions of rapid velocity change and a more lenient scale in directions of gradual change. This gradient-guided mechanism enables the interpolation process to accurately capture the principal deformation directions of anisotropic deformation regions such as active fault zones and shear zones, significantly improving the stability and spatial resolution of strain field inversion. Compared with traditional methods, this invention eliminates the need to construct complex covariance models, reducing parameter dependence. While maintaining computational simplicity, it makes the strain field distribution more consistent with the actual crustal deformation pattern, providing more reliable technical support for crustal deformation analysis and tectonic activity research. Attached Figure Description

[0032] The accompanying drawings, which form part of this application, are used to provide a further understanding of this application. The illustrative embodiments and descriptions of this application are used to explain this application and do not constitute an undue limitation of this application. In the drawings:

[0033] Figure 1 This is a flowchart of an embodiment of the present invention. Detailed Implementation

[0034] It should be noted that, unless otherwise specified, the embodiments and features described in this application can be combined with each other. This application will now be described in detail with reference to the accompanying drawings and embodiments.

[0035] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0036] The main objective of this invention is to provide a method for solving the strain of GNSS velocity fields based on gradient-guided radial basis functions (RBFs). This method aims to address the problem that traditional interpolation methods based on radial basis functions (RBFs) do not adequately consider the local variation characteristics of the GNSS velocity field in strain calculations.

[0037] This invention innovatively incorporates local gradient information of the GNSS velocity field into the RBF interpolation kernel function, constructing a non-Euclidean anisotropic distance metric. This novel metric enables the interpolation process to automatically sense and adapt to the main direction of velocity field changes; that is, it employs a finer interpolation scale in regions with large velocity gradients (such as near fault-locked areas), while automatically widening the scale in regions with gentle velocity changes.

[0038] Compared to existing methods, this invention incorporates the physical constraints (gradient information) of the velocity field itself into the interpolation process, resulting in strain tensor estimations with stronger physical properties and a more precise characterization of deformation differences on both sides of the fault zone. Simultaneously, by utilizing the spatial gradient information of the velocity field, this invention significantly enhances the stability and reliability of strain inversion, especially in regions with uneven GNSS station distribution. Because the interpolation process automatically adapts to the spatial structure of the data, it effectively improves the spatial resolution and convergence speed of strain estimation, reduces uncertainties introduced by the interpolation method, and thus achieves high-precision and high-efficiency GNSS velocity field strain inversion, improving the strain model's analytical capability for local structural deformation and its overall physical rationality.

[0039] like Figure 1 As shown, this embodiment provides a method for solving the GNSS crustal strain field based on gradient-guided radial basis functions, including the following steps:

[0040] 1. Acquire GNSS observation data for the target area, extract the planar position and horizontal velocity components of each station, and convert the latitude and longitude coordinates into a unified planar projected coordinate system to ensure consistency in spatial calculations. Simultaneously, assign weights to different stations to reflect differences in observation accuracy.

[0041] 2. Local Velocity Gradient Estimation: Within the neighborhood of each observation point, weighted linear fitting or other smoothing estimation methods are used to obtain the changing trend of the local velocity field, thereby determining the main direction and intensity of velocity field variation in each direction. This process is equivalent to extracting the local structural features of the velocity field, providing a basis for subsequent directional constraints. Gradient estimation methods can use kernel regression or structural tensor methods instead of local least squares; simultaneously, strain calculation can be performed using either analytical differentiation or finite difference methods.

[0042] 3. Constructing gradient-guided anisotropic constraints: Based on the local velocity gradient results, a constraint matrix describing directional differences is established, making the interpolation influence range larger in directions with gentle velocity changes and smaller in directions with drastic changes. This forms a "non-Euclidean distance metric" with anisotropic characteristics, replacing the isotropic distance calculation method in traditional radial basis function interpolation.

[0043] 4. Radial Basis Function Interpolation Modeling: A spatial interpolation model is established around each observation point using a radial basis function kernel function with smoothing properties. The shape and influence range of the interpolation kernel are adjusted according to the aforementioned directional constraints. This step determines the interpolation weights by solving a system of linear equations, thereby obtaining a continuous, smooth, and directionally adaptive velocity field distribution model.

[0044] 5. Strain Field Calculation and Derivative Solving: Based on the interpolated continuous velocity field, the spatial derivative of the velocity field is calculated analytically or numerically to obtain strain parameters such as the principal strain components of the Earth's crust, shear strain, and surface dilatation rate. The resulting strain field exhibits good spatial continuity and physical consistency, effectively reflecting the principal directional characteristics of tectonic deformation.

[0045] 6. Results Analysis and Visualization: Spatial interpolation and gridding are performed on the obtained strain field to generate distribution maps of principal strain direction, shear strain rate and surface dilatation rate, which can assist in the quantitative analysis of crustal active zones, fault regions and strain concentration areas.

[0046] As a specific implementation method of this embodiment, taking GNSS observation data of a certain area as an example, the specific implementation steps are as follows:

[0047] Projection transformation converts latitude and longitude coordinates into planar coordinates. The planar coordinates after projection transformation are represented as follows: ;

[0048] Local velocity gradient estimation at each observation point Within the neighborhood of , the horizontal velocity component and It can be approximated as a linear function:

[0049] ;

[0050] ;

[0051] here: The velocity constant (local average velocity) represents the velocity at the observation point. This represents the local spatial derivative (gradient) of the velocity component. For another observation point in the neighborhood... , with horizontal velocity component For example, the objective function is:

[0052] ;

[0053] in, , ,calculate: .

[0054] in:

[0055] , ;

[0056] It is the weighted summation term in the locally weighted least squares solution, which originates from the weights of the neighborhood points. :

[0057] ;

[0058] exp(·) denotes the natural exponent, that is, the power of the natural constant ℮≈2.71828, i.e., exp(x) = .

[0059] These weights are calculated from the distance between observation points using a Gaussian or exponential function, and are used to make the local fit smoother and more physically reasonable. Gaussian weights mean that the weight decreases as the distance increases, used to smooth local fitting. The smoothing length controls the weight decay rate, typically taken as 0.5 to 1 times the average measurement point spacing, or adaptively determined based on local point density. According to the above formula, the solution can be obtained... Similarly, the horizontal velocity component can be derived. Unknown parameters These represent the local spatial derivatives (gradients) of the velocity components, reflecting the intensity of the velocity change in opposite directions.

[0060] The main innovation of this application is to construct gradient-guided anisotropic constraints, define a local principal direction matrix, and obtain the principal direction of velocity change through the local gradient results. The gradient vector is defined as follows:

[0061] Right now ;

[0062] Calculate the symmetric part using the defined gradient vector: .

[0063] Then to Perform eigenvalue decomposition: The eigenvectors obtained by solving , Indicates the principal direction of change of the velocity field. Eigenvalues , Indicates the intensity in a direction. Principal direction angle:

[0064] ;

[0065] Define an anisotropic scaling factor, where the scaling factor along the principal directions is a function of the inverse of the gradient, where... This is an adjustment factor, typically set to 0.2–1.0, or can be adaptively adjusted according to actual conditions. To prevent the use of tiny constants for division by zero, a value of 10 is generally chosen. -4 :

[0066] ;

[0067] Rotation matrix:

[0068] ;

[0069] Finally, an anisotropic matrix can be constructed:

[0070] ;

[0071] The principal direction can be determined by the eigenvalues ​​and eigenvectors of the gradient matrix, and is used to describe the local expansion and contraction directions of the velocity field. Traditional RBF uses Euclidean distance:

[0072] ;

[0073] Guided by gradients, we adjust the distance calculation based on the velocity gradient, causing the interpolation kernel to stretch or compress in the main direction of change. In this application, non-Euclidean distance is used.

[0074] ;

[0075] Radial basis function (RBF) interpolation modeling can use a kernel function other than the Gaussian kernel described in the application, such as a multiquadric or Matern kernel, to adjust the smoothness. In this embodiment, the RBF system uses a Gaussian kernel function. The non-Euclidean distances between observation points are used to construct the relationship matrix between all observation points:

[0076] ;

[0077] in It is a shape parameter that controls the smoothness of the kernel function. When the value is larger, the results are smoother; when the value is smaller, the interpolation is closer to the data points. Generally, the interpolation value is 1–3 times the average spacing of the measurement points. Construct the RBF system equations with polynomial constraints:

[0078] ;

[0079] Some parameters in the equation have a polynomial basis matrix:

[0080] ;

[0081] Constraints:

[0082] ;

[0083] Solving for interpolation coefficients , :

[0084] ;

[0085] Solved interpolation coefficients , Substituting into the radial basis interpolation function yields the following result. The expression for the interpolation result of the continuous velocity field of the velocity components on the plane. Similarly, the interpolation function is in the form of ( , (These are the radial basis weighting coefficients corresponding to the observation points).

[0086] ;

[0087] ;

[0088] Based on the results of the previous step, we have obtained a continuous velocity field, and can then calculate the strain components using the derivative of the velocity field. The interpolation results from the RBF model yield the following:

[0089] ;

[0090] ;

[0091] ;

[0092] ;

[0093] The strain tensor in the two-dimensional plane is:

[0094] ;

[0095] in:

[0096] ;

[0097] The maximum shear strain can then be calculated:

[0098] ;

[0099] Solve for the surface expansion rate:

[0100] ;

[0101] Principal strain rate:

[0102] ;

[0103] Visualize the calculation results. Further use GMT or Matplotlib for interpolation and gridding to plot the various strain field distributions calculated.

[0104] The core innovation of this invention lies in introducing local velocity gradient information into the radial basis function interpolation framework. By establishing a directional constraint mechanism based on the characteristics of velocity field changes, the interpolation kernel function can adaptively adjust its shape and range of action, realizing the transformation from isotropic interpolation to anisotropic interpolation. This design constructs a non-Euclidean distance metric related to the local gradient, allowing the interpolation kernel to expand in directions of gentle velocity change and contract in directions of drastic change. This effectively captures the main deformation direction in regions such as active fault zones and shear zones, maintaining the continuity and physical rationality of the velocity and strain fields.

[0105] This invention establishes a stable and deterministic strain inversion system through the aforementioned gradient guidance mechanism, achieving a unification of directional constraints and spatial smoothing without requiring assumptions about a covariance model or complex statistical parameter estimation. Theoretically, this method extends radial basis function interpolation from isotropic to anisotropic, and structurally forms a closed-loop strain inversion framework, which can be widely applied to high-precision strain estimation in GNSS velocity fields.

[0106] This application provides a GNSS crustal strain field solution method based on gradient-guided radial basis functions. By introducing local velocity gradient information into the interpolation kernel and establishing an anisotropic constraint model, the interpolation kernel function can adaptively conform to the main deformation direction of the crust, achieving accurate characterization of the spatial direction characteristics of the velocity field. This constraint effectively enhances the stability of the interpolation process in anisotropic deformation regions such as active fault zones and shear zones, avoiding common problems in traditional isotropic interpolation methods such as strain anomalies, gradient oscillations, and physical discontinuities.

[0107] Unlike traditional radial basis function or least squares configuration methods that rely solely on the isotropic assumption, this application uses the local velocity gradient as the core guiding factor, enabling the interpolation process to adaptively adjust the kernel function shape according to the actual directional variation characteristics of the observed data. Through this mechanism, the final strain field distribution is more consistent with the actual deformation mode and tectonic trend of the crust, significantly improving the spatial resolution and physical reliability of strain estimation.

[0108] Meanwhile, this method eliminates the need for complex covariance function models, reducing parameter dependencies and computational uncertainties. The algorithm boasts a simple structure and strong scalability. This method can be further extended to crustal deformation time-series analysis and strain evolution studies in tectonically active regions. By introducing gradient-guided constraints into the interpolation kernel, this application achieves enhanced direction recognition capabilities while maintaining smoothness, thereby improving strain inversion accuracy and computational stability.

[0109] The above are merely preferred embodiments of this application, but the scope of protection of this application is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application. Therefore, the scope of protection of this application should be determined by the scope of the claims.

Claims

1. A method for solving the GNSS crustal strain field based on gradient-guided radial basis functions, characterized in that, Includes the following steps: Acquire GNSS observation data of the target area and convert it to a plane projection coordinate system to obtain the plane position of each station; Local velocity gradients are estimated based on the planar position and horizontal velocity components of each station. An anisotropic constraint matrix is ​​constructed based on the local velocity gradient to define a non-Euclidean distance metric; A radial basis function interpolation model is established based on the non-Euclidean distance metric, and the interpolation coefficients are solved to obtain a continuous velocity field. Calculate the strain tensor components based on the continuous velocity field; A strain field is constructed based on the strain tensor components, and the strain field is then meshed and visualized.

2. The GNSS crustal strain field solution method based on gradient-guided radial basis functions according to claim 1, characterized in that, The process of estimating the local velocity gradient includes: using a weighted linear fitting method, the local velocity gradient is obtained by solving the spatial derivative of the horizontal velocity component in the neighborhood of each observation point through least squares.

3. The GNSS crustal strain field solution method based on gradient-guided radial basis functions according to claim 2, characterized in that, The weights used in the weighted linear fitting are Gaussian weights. A Gaussian function is used to process the Euclidean distance between observation points to obtain weight values ​​that decay with distance. The smoothing length parameter of the Gaussian function is taken as 0.5 to 1 times the average distance between measuring points in the target area.

4. The method for solving the GNSS crustal strain field based on gradient-guided radial basis functions according to claim 1, characterized in that, The process of constructing the anisotropic constraint matrix includes: Calculation of symmetric strain rate tensor based on local velocity gradient; Eigendecomposition of the symmetric strain rate tensor yields the principal variation direction and its corresponding eigenvalues; Calculate the anisotropic scaling factor based on the eigenvalues; Construct a rotation matrix based on the main direction of change; Anisotropic constraint matrices are synthesized based on the rotation matrix and the anisotropic scaling factor.

5. The GNSS crustal strain field solution method based on gradient-guided radial basis functions according to claim 4, characterized in that, The anisotropic scaling factor is calculated by the ratio of the adjustment factor to the sum of the absolute values ​​of the eigenvalues, wherein the adjustment factor ranges from 0.2 to 1.

0.

6. The method for solving the GNSS crustal strain field based on gradient-guided radial basis functions according to claim 1, characterized in that, The process of establishing a radial basis function interpolation model based on the aforementioned non-Euclidean distance metric and solving for the interpolation coefficients to obtain the continuous velocity field includes: The non-Euclidean distance metric is processed using a Gaussian kernel function to obtain the radial basis function relation matrix; Based on the relationship matrix and the observed velocity values, a system of linear equations with polynomial constraints is constructed. Solving the linear equations yields the radial basis function weighting coefficients and polynomial coefficients; The expression for the continuous velocity field is obtained based on the weighting coefficients, polynomial coefficients, and radial basis functions.

7. The GNSS crustal strain field solution method based on gradient-guided radial basis functions according to claim 1, characterized in that, The calculation of strain tensor components includes: The spatial partial derivatives of the velocity field are analytically obtained based on the radial basis function interpolation model. The components of the strain tensor, the maximum shear strain, the surface dilatation rate, and the principal strain rates are calculated using the spatial partial derivatives.

Images