A GNSS imaging method and device for detecting short-wavelength deformation

Abstract

The present invention provides a GNSS imaging method and device for detecting short-wavelength deformation. For each point to be estimated, a weight is first determined based on the spatial structure, and then an estimated vertical velocity is extracted to generate a GNSS image. Determining the weight based on the spatial structure includes first determining known stations connected to the point to be estimated based on the Delaunay triangle, then using Gaussian process regression to obtain preliminary weights, and then re-weighting the connected known stations based on velocity uncertainty to obtain the final weights. Extracting the estimated vertical velocity includes extracting the estimated vertical velocity at the point to be estimated using a generalized weighted median filter based on the final weights of the known stations connected to the point to be estimated. The present invention proposes a Gaussian process regression (GPR-VU) scheme that accounts for velocity uncertainty. This scheme uses a covariance matrix to construct the spatial structure between known and unknown points, generating GNSS images of crustal deformation within the study area, and enabling the detection of short-wavelength deformation in areas with complex geology and intense tectonic activity.

Description

Technical Field

[0001] The present invention belongs to the technical field of GNSS data precision processing, and in particular relates to a GNSS imaging technology solution for detecting short-wavelength deformation. Background Art

[0002] Vertical land motion is caused by long-term processes such as glacial isostatic adjustment (GIA), hydrological loading, and tectonic movement. A direct analysis of vertical land motion can directly explain geophysical processes and the underlying causes of crustal motion, Earth's internal structure, and global climate change. Over the past few decades, numerous researchers have conducted extensive observations of vertical motion in the solid Earth using various measurement techniques, such as leveling, synthetic aperture radar interferometry (InSAR), and tide gauges. However, leveling is difficult to automate and data acquisition is challenging; InSAR exhibits errors in long-wavelength deformation and is only sensitive to surface motion in the satellite-to-ground line of sight, making it difficult to accurately assess the magnitude and direction of surface deformation; and tide gauges can measure sea level changes but cannot observe crustal motion in inland areas. However, GNSS can automatically provide continuous, high-precision coordinates. The long-term trends (velocities) of GNSS station coordinate time series reflect phenomena such as tectonic movement, fault stress accumulation, and glacial isostatic adjustment, providing fundamental data for geophysical research.

[0003] With the increasing number of GNSS stations and the expansion of observation ranges, it has become possible to construct GNSS imagery reflecting spatially continuous deformation using time series of multiple GNSS station coordinates within a target study area. Hammond et al. (2016) first proposed a GNSS imaging method (GIM method). This method first defines a spatial structure function (SSF), then quantifies the "structure" between station pairs. Velocities are then interpolated based on the SSF and velocity uncertainty to generate GNSS imagery reflecting spatially continuous VLM. This method has been successfully used to monitor deformation in California and Nevada, USA. However, GNSS imagery derived from this GIM method can exhibit patchy patterns and suffer from oversmoothing of velocity peaks, making it difficult to detect short-wavelength deformation and potentially leading to inappropriate geophysical interpretation. Husson et al. (2018) used transdimensional regression (TR) to generate GNSS imagery for studying deformation caused by long-wavelength processes, such as GIA. However, this method is insensitive to local variations and cannot explore short-wavelength deformation. Zhou et al. (2020) used the GIM method with a modified spatial structure function (SSF) to generate detailed GNSS images for the western United States and China, verifying that GNSS images can accurately reflect the spatiotemporal distribution characteristics of VLM. However, the GIM method is highly dependent on the spatial density of the GNSS network to construct an SSF that can correctly reflect the spatial structure of the site. In addition, the SSF is a piecewise function, which assumes that the velocity changes within a certain range in a specific direction are smooth. This assumption is not sufficient to describe the actual situation when applying the GIM method, and may lead to significant differences in velocity estimates of adjacent areas in the results, forming fragmentation and possibly over-smoothing velocity peaks. Therefore, the SSF-based method has a certain degree of empiricism, and the detection effect of short-wavelength deformation is unstable.

[0004] Therefore, this field urgently needs to propose a short-wavelength deformation measurement scheme that is simpler to implement and more applicable. Summary of the Invention

[0005] In view of the shortcomings of the existing technology, the purpose of the present invention is to propose a GNSS imaging technology that uses multiple GNSS station coordinate time series to construct a regional crustal vertical deformation map to obtain small-scale, short-wavelength deformation.

[0006] To achieve the above objectives, the technical solution of the present invention provides a GNSS imaging method for detecting short-wavelength deformation. For each point to be estimated, a weight is first determined based on the spatial structure, and then an estimated value of the vertical velocity is extracted to generate a GNSS image.

[0007] Determining the weights based on the spatial structure includes first obtaining known stations connected to the point to be estimated based on the Delaunay triangle, then using Gaussian process regression to obtain preliminary weights, and then re-weighting based on velocity uncertainty to obtain the final weights of the connected known stations;

[0008] The extracting of the estimated value of the vertical velocity includes extracting the estimated value of the vertical velocity of the point to be estimated by using a generalized weighted median filter based on the final weights of the known measuring stations connected to the point to be estimated.

[0009] Moreover, the study area is divided into regular grids, and each grid point is used as a point to be estimated, thus obtaining a set of points to be estimated. After extracting the vertical velocity estimate of each grid point, a GNSS image of the crustal deformation in the study area is generated.

[0010] Moreover, the vertical velocity estimation value of one to-be-estimated point in the set of to-be-estimated points is extracted each time.

[0011] Furthermore, preliminary weights are obtained by using Gaussian process regression to describe the spatial relationships between station pairs within the GNSS network.

[0012] Furthermore, the re-weighting based on velocity uncertainty is implemented by re-weighting the vertices based on the velocity uncertainty of the measuring station at each Delaunay triangle vertex.

[0013] Moreover, the discrete velocity field is converted into a continuous velocity field by re-weighting the known stations based on the velocity uncertainty.

[0014] Furthermore, the generalized weighted median filter is implemented by using the vertical velocity values ​​of known measuring stations as samples, transferring the positive and negative values ​​of the weights to the samples. The resulting samples are called signed samples. The signed samples and the corresponding absolute value weights are input into the generalized weighted median filter, and the vertical velocity estimate of the point to be estimated is output.

[0015] On the other hand, the present invention also provides an electronic device, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein when the processor executes the program, the GNSS imaging method for measuring short-wavelength deformation as described above is implemented.

[0016] On the other hand, the present invention further provides a non-transitory computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, implements the GNSS imaging method for measuring short-wavelength deformation as described above.

[0017] On the other hand, the present invention further provides a computer program product, comprising a computer program, which, when executed by a processor, implements the GNSS imaging method for measuring short-wavelength deformation as described above.

[0018] To address the problems of GNSS image fragmentation and over-smoothing of velocity peaks caused by sparse GNSS networks, this paper proposes a Gaussian process regression (GPR-VU) method that considers velocity uncertainty. The method uses the covariance matrix to construct the spatial structure between known and unknown points, breaking away from the empirical SSF construction in the traditional GIM method and inferring a posterior velocity estimate with a probabilistic advantage. Based on this, a GNSS image of crustal deformation in the study area is generated, enabling the detection of short-wavelength deformation in areas with complex geology and intense tectonic activity.

[0019] The solution of the present invention is simple and convenient to implement and has strong practicality. It solves the problems of low practicality and inconvenience in actual application existing in related technologies, can improve user experience, and has important market value. DETAILED DESCRIPTION

[0020] The following will further illustrate the concept, specific structure and technical effects of the present invention in conjunction with embodiments, so as to fully understand the purpose, characteristics and effects of the present invention.

[0021] The embodiment of the present invention proposes a GNSS imaging method for detecting short-wavelength deformation, which is implemented based on the following parts:

[0022] (1) Extracting vertical velocities from known stations: For implementation, vertical GNSS coordinate time series of multiple stations can be obtained from public datasets or calculated independently. MIDAS (a nonparametric method based on medians and accounting for interannual variability) is then used to calculate vertical velocities at these stations, serving as the basis for the known point data. These time series must meet the following criteria: 1) be longer than five years; 2) have a data availability greater than 80%; and 3) have a velocity uncertainty less than 5 mm / yr (millimeter per year).

[0023] (2) Constructing a spatial structure description between pairs of stations within the GNSS network, specifically implemented as follows:

[0024] 2.1) Use the covariance matrix of Gaussian process regression (GPR) to construct the spatial structure between known and unknown points.

[0025] Furthermore, to facilitate understanding of the technical solution of the present invention, a construction principle description is provided, which includes the following sub-steps:

[0026] 2.1.1) Under the kernel function metric, define the corresponding coordinate points x of any two stations i ,x j The covariance function k(x i ,x j ):

[0027]

[0028] Among them, ‖x i -x j ‖ is any two points x i ,x j The distance between them is in units of spherical angular distance on a great circle, i and j are station numbers; σ f 2 and l are the parameters to be estimated, σ f 2 is the scale factor, l is the distance scaling factor; exp() is the natural exponent.

[0029] 2.1.2) Assume that there are n known GNSS station coordinates, denoted as vector X = (x1, x2, …, x n ), defines the covariance function between the coordinates of n GNSS stations, which is an n×n matrix:

[0030]

[0031] Among them, k(x i ,x j ) is the Gaussian kernel function in formula (1).

[0032] 2.1.3) Under the Gaussian process regression assumption, the vertical velocities corresponding to n known stations are expressed as a vector Y, written as:

[0033] Y=f(X)+ε (3)

[0034] Y~N(0,K(X,X)+σ n 2 I n ) (4)

[0035] Among them, ε is the noise term, N(,) represents the normal distribution, σ n 2 is the prior variance of the noise term, i.e. ε~N(0,σ n 2 I n ) means that ε obeys the expectation of 0 and the prior variance is σ n 2 I n Normal distribution of I n represents the n×n identity matrix; f(X) represents the function value corresponding to X without considering noise, Y~N(,) represents that Y obeys the expectation of 0, and the variance is (K(X,X)+σ n 2 I n ) is normally distributed.

[0036] Therefore, the data set of known points is D = {X, Y}, and the posterior distribution p(f|D) can be obtained.

[0037] 2.2) For new input points (points to be estimated), the conditional distribution is calculated based on the posterior distribution to estimate the function value.

[0038] Furthermore, to facilitate understanding of the technical solution of the present invention, a corresponding implementation principle description is provided, which includes the following sub-steps:

[0039] 2.2.1) Introduce the vector X of the set of points to be estimated * You can get the corresponding function value Y * , that is, the vertical velocity value to be estimated, assuming that the observed value Y and the unknown function value Y * They satisfy the joint multivariate normal distribution and are recorded as:

[0040]

[0041] Among them, K(X,X) is the covariance matrix between the corresponding vectors X of all known observation points, K(X,X * ) is the covariance matrix between the newly introduced estimated point and the corresponding vector of the known point, K(X * ,X * ) is the covariance matrix between the corresponding vectors of the newly introduced points to be estimated; Represents the observed value Y and the unknown function value Y * The joint random variable of .

[0042] 2.2.2) Taking the maximization of the maximum likelihood function as the goal, the maximum likelihood estimation method is used to solve the optimal parameters to be estimated The expression is as follows:

[0043]

[0044] Where Ln() is the logarithmic symbol, and p(Y|X,θ) is the probability distribution of Y under the conditions of X and θ.

[0045] 2.2.3) According to the standard process, calculate the conditional probability distribution p(Y * |X * ,X,Y):

[0046] E(Y * |Y)=K(X * ,X)[K(X,X)+σ n 2 I n ] -1 Y (7)

[0047] D(Y * |Y)=K(X* ,X * )-K(X * ,X)[K(X,X)+σ n 2 I n ] -1 K(X,X * ) (8)

[0048] Among them, E(Y * |Y) is Y when Y is known * The conditional expectation of the matrix, the superscript -1 indicates the inverse of the matrix, D(Y * |Y) is Y when Y is known * conditional variance.

[0049] (3) Perform Delaunay triangulation on the site locations on the sphere and construct a weight set.

[0050] The spherical surface is the surface of the earth. In order to speed up the calculation efficiency and find the optimal parameters during the GPR process, the present invention further proposes a preferred solution of estimating only the set of points to be estimated X at a time. * Therefore, the conditional mean can be regarded as the output value y * A point x to be estimated * The output y at * , can be obtained using the following simplified formula:

[0051] y * =K(x * ,X)[K(X,X)+σ n 2 I n ] -1 Y=PY (9)

[0052] Where P = K(x * ,X)[K(X,X)+σ n 2 I n ] -1 , corresponds to the weight of the station connected to the point to be estimated, which is a 1×n matrix, where n represents the number of stations connected to the point to be estimated in Delaunay; here x * is the coordinate of a point to be estimated; X and Y are the coordinates of a known point and the corresponding vertical velocity vector respectively.

[0053] (4) Recalculate the weights, specifically:

[0054] Calculate the weight of each known station based on the GPR and velocity uncertainty. Within each Delaunay triangle, recalculate the weight of the known station taking into account the velocity uncertainty.

[0055] In this embodiment, the velocity uncertainty of each station on the vertex of the Delaunay triangle is considered, and the vertex is re-weighted. The j-th known point is re-weighted to weight w j , calculated by the following formula:

[0056]

[0057] Where, j = 1, 2, ..., n, P j is the weight corresponding to the jth known point (station) in GPR, σ j is the velocity uncertainty of the j-th known point, which can be estimated by the MIDAS method in specific implementation. The MIDAS method is a prior art and will not be described in detail in this invention.

[0058] (5) Since the use of a general weighted median filter may include negative weights, the present invention proposes to use a generalized weighted median filter to process the sites; when the weights are non-negative, the filter degenerates into a general weighted median filter; median filtering is performed on known sites to reduce the impact of outliers and enhance regional common features.

[0059] The embodiment uses the vertical velocity value of the known measuring station as a sample, transfers the positive and negative properties of the weight to the sample, and the result is called a symbol sample. The specific implementation is: the symbol sample and its corresponding absolute value weight are represented as S i and |w i |, where S i =sign(w i )·Y i .

[0060] Then, the sorted symbol samples and their absolute value weights are represented as S (i) and |w (i) |, that is, S (1) ≤S (2) ≤…≤S (n) The output v of the generalized weighted median filter is the estimated vertical velocity of the point to be estimated, which can be written as:

[0061]

[0062] in is the threshold, S (k) refers to the kth symbol sample value after sorting, It refers to the smallest symbol sample value obtained by adding the absolute value weights of the symbol samples in ascending order to a value greater than the threshold T0.

[0063] Based on the above theoretical foundation, the present invention proposes a GNSS imaging method for detecting short-wavelength deformation to construct a GPR-based GNSS image. The specific implementation is as follows:

[0064] 1) Divide the study area into regular grids.

[0065] After division, each grid point is a point to be estimated, forming a set of points to be estimated.

[0066] 2) extracting the vertical velocity estimate for each point to be estimated, including processing each point to be estimated separately, first determining the weight based on the spatial structure, and then extracting the vertical velocity estimate;

[0067] The implementation method of determining the weight based on the spatial structure is as follows:

[0068] First, the Delaunay triangle is used to obtain the known stations connected to the point to be estimated. Then, the Gaussian process regression (Formula (9)) is used to obtain the preliminary weights. Finally, Formula (10) is used to obtain the final weights of the connected stations after the velocity uncertainty is re-weighted.

[0069] The method for extracting the estimated value of the vertical velocity is to obtain the estimated value of the vertical velocity of the point to be estimated by using the generalized weighted median method (i.e., using the generalized weighted median filter shown in formula (11)) based on the final weights of the connected measuring stations.

[0070] In specific implementation, depending on computing resources, the above processing can be performed on multiple points to be estimated in parallel, or the above processing can be performed on each point to be estimated in sequence.

[0071] 3) Use the weighted median processing results to construct the GNSS image.

[0072] By combining the probabilistic advantages of Gaussian process regression and the corresponding velocity uncertainty, we determine the final weights for the known stations. Finally, we use the generalized weighted median to estimate the vertical velocity of the target point. Based on the vertical velocity estimate for each grid point, we can then map the estimated results for all grid points in the study area and construct a GNSS image.

[0073] To address the problem of over-smoothing in the traditional GNSS imaging method (GIM) when processing velocity peaks, which makes it difficult to detect short-wavelength deformations, the present invention innovatively utilizes Gaussian process regression (GPR) to describe the spatial relationship between adjacent station pairs and uses their uncertainties to reweight the velocities of known stations, converting the discrete velocity field into a continuous velocity field. To illustrate the technical effects of the present invention, the statistical results of checkerboard simulation experiments using the GIM method and the embodiment method of the present invention (GPR-VU) in the southwestern United States and the eastern Qinghai-Tibet Plateau are provided as shown in Tables 1 and 2:

[0074] Table 1. Experimental results of 1°×1° to 4°×4° checkerboard grids in the southwestern United States.

[0075]

[0076] Table 2 Experimental results of the 3°×3° to 7°×7° checkerboard grid in the eastern Qinghai-Tibet Plateau

[0077]

[0078] The method has been demonstrated in simulation experiments and experiments with field data from the southwestern United States and the eastern Tibetan Plateau of China. GPR-VU performs better than the GIM method in detecting short-wavelength deformations.

[0079] In specific implementation, the method proposed in the technical solution of the present invention can be automatically run by those skilled in the art using computer software technology. System devices that implement the method, such as computer-readable storage media that store the corresponding computer program of the technical solution of the present invention and computer equipment that runs the corresponding computer program, should also be within the scope of protection of the present invention.

[0080] The following describes the GNSS imaging electronic device for measuring short-wavelength deformation provided by the present invention. The GNSS imaging electronic device for measuring short-wavelength deformation described below and the GNSS imaging method for measuring short-wavelength deformation described above can be referenced to each other.

[0081] The electronic device may include a processor, a communications interface, memory, and a communications bus. The processor, communications interface, and memory communicate with each other via the communications bus. The processor may invoke logic instructions in the memory to execute the GNSS imaging method for measuring short-wavelength deformation, primarily including the software processing portion described above.

[0082] In addition, the logical instructions in the above-mentioned memory can be implemented in the form of a software functional unit and can be stored in a computer-readable storage medium when sold or used as an independent product. Based on this understanding, the technical solution of the present invention, or the part that contributes to the prior art, or the part of the technical solution, can be embodied in the form of a software product. The computer software product is stored in a storage medium and includes several instructions for enabling a computer device (which can be a personal computer, a server, or a network device, etc.) to execute all or part of the steps of the method described in each embodiment of the present invention. The aforementioned storage medium includes: various media that can store program codes, such as a USB flash drive, a mobile hard disk, a read-only memory (ROM), a random access memory (RAM), a magnetic disk or an optical disk.

[0083] In another embodiment, the present invention also provides a computer program product, which includes a computer program. The computer program can be stored on a non-transitory computer-readable storage medium. When the computer program is executed by a processor, the computer can execute the software processing part of the GNSS imaging method for measuring short-wavelength deformation provided by the above-mentioned methods.

[0084] In another embodiment, the present invention further provides a non-transitory computer-readable storage medium having a computer program stored thereon, which, when executed by a processor, is implemented to execute the software processing portion of the GNSS imaging method for measuring short-wavelength deformation provided by the above-mentioned methods.

[0085] The device embodiments described above are merely illustrative. The units described as separate components may or may not be physically separate, and the components shown as units may or may not be physical units, i.e., they may be located in one location or distributed across multiple network units. Some or all of the modules may be selected based on actual needs to achieve the objectives of the present embodiment. Persons of ordinary skill in the art will be able to understand and implement the present invention without inventive effort.

[0086] Through the description of the above embodiments, those skilled in the art can clearly understand that each embodiment can be implemented by means of software plus a necessary general hardware platform, or of course, by hardware. Based on this understanding, the essence of the above technical solution or the part that contributes to the existing technology can be embodied in the form of a software product. The computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, a magnetic disk, an optical disk, etc., and includes a number of instructions for enabling a computer device (which can be a personal computer, a server, or a network device, etc.) to execute the methods described in each embodiment or certain parts of the embodiments.

[0087] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, rather than to limit it. Although the present invention has been described in detail with reference to the aforementioned embodiments, those skilled in the art should understand that they can still modify the technical solutions described in the aforementioned embodiments, or make equivalent replacements for some of the technical features therein. However, these modifications or replacements do not deviate the essence of the corresponding technical solutions from the spirit and scope of the technical solutions of the various embodiments of the present invention.

Claims

1. A GNSS imaging method for detecting short-wavelength deformation, characterized by: For each point to be estimated, a weight is first determined based on the spatial structure, and then the estimated value of the vertical velocity is extracted to generate a GNSS image. Determining the weights based on the spatial structure includes first obtaining known stations connected to the point to be estimated based on the Delaunay triangle, then using Gaussian process regression to obtain preliminary weights, and then re-weighting based on velocity uncertainty to obtain the final weights of the connected known stations; Extracting the estimated vertical velocity includes extracting the estimated vertical velocity of the point to be estimated using a generalized weighted median filter based on the final weights of known measuring stations connected to the point to be estimated. The generalized weighted median filter is implemented by using the vertical velocity values ​​of the known measuring stations as samples, transferring the positive and negative values ​​of the weights to the samples, and obtaining a result called a signed sample. The signed sample and the corresponding absolute value weight are input into the generalized weighted median filter to output the estimated vertical velocity of the point to be estimated.

2. The GNSS imaging method for detecting short-wavelength deformation according to claim 1, characterized in that: The study area is divided into a regular grid, and each grid point is used as a point to be estimated, thus obtaining a set of points to be estimated. The vertical velocity estimate of each grid point is extracted and then a GNSS image of the crustal deformation in the study area is generated.

3. The GNSS imaging method for detecting short-wavelength deformation according to claim 2, characterized in that: Each time, the vertical velocity estimation value of one point to be estimated in the set of points to be estimated is extracted.

4. The GNSS imaging method for detecting short-wavelength deformation according to claim 1, characterized in that: The preliminary weights are obtained by using Gaussian process regression to describe the spatial relationship between station pairs in the GNSS network.

5. The GNSS imaging method for detecting short-wavelength deformation according to claim 1, characterized in that: The re-weighting based on velocity uncertainty is implemented by re-weighting the vertices based on the velocity uncertainty of the measuring station at each Delaunay triangle vertex.

6. The GNSS imaging method for detecting short-wavelength deformation according to claim 1, characterized in that: The discrete velocity field is converted into a continuous velocity field by re-weighting the known stations based on the velocity uncertainty.

7. An electronic device comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein: When the processor executes the program, the GNSS imaging method for detecting short-wavelength deformation according to any one of claims 1 to 6 is implemented.

8. A non-transitory computer-readable storage medium having a computer program stored thereon, characterized in that: When the computer program is executed by a processor, the GNSS imaging method for detecting short-wavelength deformation according to any one of claims 1 to 6 is implemented.

9. A computer program product comprising a computer program, characterized in that: When the computer program is executed by a processor, the GNSS imaging method for detecting short-wavelength deformation according to any one of claims 1 to 6 is implemented.

Images