GNSS imaging method and device for detecting short wavelength deformation

Abstract

The invention provides a GNSS imaging method and device for detecting short-wavelength deformation, and the method comprises the steps: firstly determining the weight of each to-be-estimated point based on a space structure, then extracting an estimated value of a vertical speed, and generating a GNSS image; the step of determining the weight based on the space structure comprises the steps of firstly obtaining a known observation station connected with a to-be-estimated point based on a Delaunay triangle, then obtaining a preliminary weight by using Gaussian process regression, and obtaining a final weight of the connected known observation station after re-weighting based on velocity uncertainty; the step of extracting the estimated value of the vertical velocity comprises the step of extracting the estimated value of the vertical velocity of the to-be-estimated point by adopting generalized weighted median filtering processing based on the final weight of a known observation station connected with the to-be-estimated point. The invention provides a Gaussian process regression GPR-VU scheme considering speed uncertainty, a space structure between a known point and an unknown point is constructed by using a covariance matrix, an earth crust deformation GNSS image in a research area is generated, and detection of short wavelength deformation in an area with complex geology and strong structural activity is realized.

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 loads, and tectonic movements. Intuitive analysis of vertical land motion can directly explain geophysical processes and potential causes of crustal movement, Earth's internal structure, and global climate change. In the past few decades, many scholars have made extensive observations of vertical motion of 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 difficult; InSAR has errors in long-wavelength deformation and is only sensitive to surface motion in the direction of the satellite-to-ground line of sight, making it difficult to accurately assess the magnitude and direction of surface deformation; tide gauges can measure sea level changes but cannot observe crustal movement in inland areas. However, GNSS can automatically provide continuous, high-precision coordinates. The long-term trend (velocity) of the GNSS station coordinate time series reflects phenomena such as tectonic movement, fault stress accumulation, and glacial isostatic adjustment, providing basic data for geophysical research.

[0003] With the increase in the number of GNSS stations and the expansion of the observation range, it has become a reality to construct GNSS images reflecting spatially continuous deformation using multiple GNSS station coordinate time series in the target study area. Hammond et al. (2016) first proposed a GNSS image method (GIM method), which first defined a spatial structure function (SSF), then quantified the "structure" between station pairs, and then interpolated the velocity based on the SSF and velocity uncertainty to generate GNSS images reflecting spatially continuous VLM. It was successfully used to monitor deformation in California and Nevada, USA. The GNSS images derived from the GIM method may have fragmented patches and the problem of excessive smoothing of velocity peaks, which makes it difficult to detect short-wavelength deformation and improper geophysical interpretation. Husson et al. (2018) used cross-dimensional regression (TR) to generate GNSS images for studying long-wavelength processes, such as deformation caused by GIA. However, it is not sensitive to local regional changes 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 change within a certain range in a specific direction is 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 in adjacent areas in the results, forming fragmentation and possibly over-smoothing velocity peaks. Therefore, the SSF-based method is somewhat empirical and the detection effect of short-wavelength deformation is unstable.

[0004] Therefore, there is an urgent need in this field 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 prior art, the purpose of the present invention is to propose a GNSS imaging technology that utilizes 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-mentioned purpose, the technical solution of the present invention provides a GNSS imaging method for detecting short-wavelength deformation, wherein a weight is first determined for each point to be estimated based on the spatial structure, and then an estimated value of the vertical velocity is extracted to generate a GNSS image;

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

[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 process 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 to-be-estimated point set is extracted each time.

[0011] Furthermore, preliminary weights are obtained by using Gaussian process regression to describe the spatial relationship 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 on each Delaunay triangle vertex.

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

[0014] Moreover, the implementation method of the generalized weighted median filtering process is to use the vertical velocity value of the known measuring station as a sample, transfer the positive and negative properties of the weight to the sample, and the result is called a sign sample. The sign sample and the corresponding absolute value weight 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 also 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] Aiming at the problems of fragmentation and over-smoothing of velocity peaks caused by sparse GNSS network in GNSS images, this paper proposes a Gaussian process regression (GPR-VU) method considering velocity uncertainty, which uses the covariance matrix to construct the spatial structure between known points and unknown points, gets rid of the empiricism of SSF construction in the traditional GIM method, and infers the posterior velocity estimate with probability advantage. Based on this, the GNSS image of crustal deformation in the study area is generated, realizing 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, has strong practicality, 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 concept, specific structure and technical effects of the present invention will be further described below in conjunction with embodiments, so as to fully understand the purpose, characteristics and effects of the present invention.

[0021] A GNSS imaging method for detecting short-wavelength deformation proposed in an embodiment of the present invention is implemented based on the following parts:

[0022] (I) Extracting vertical velocities of known stations: In specific implementation, the vertical GNSS coordinate time series of multiple stations can be obtained from public data sets or calculated by oneself, and the vertical velocities of these stations can be calculated using the MIDAS method (a non-parametric method based on the median and taking into account interannual differences) as the basis for known point data. These time series must meet the following screening conditions: 1) the time series is more than 5 years; 2) the data availability is greater than 80%; 3) the velocity uncertainty is less than 5mm / yr (millimeter per year).

[0023] (II) Constructing the spatial structure description between the station pairs in the GNSS network, which is specifically implemented as follows:

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

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

[0026] 2.1.1) Under the kernel function metric, define the corresponding coordinates 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 measured in 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 involved, denoted by 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 velocity corresponding to n known stations is 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, that is, ε~N(0,σ n 2 I n ) means that ε has an expectation of 0 and a prior variance of σ n 2 I n Normal distribution of I n represents the n×n unit 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 a normal distribution.

[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, in order to facilitate understanding of the technical solution of the present invention, a corresponding implementation principle description is provided, including the following sub-steps:

[0039] 2.2.1) Introduce the vector X of the set of points to be estimated * The corresponding function value Y can be obtained * , 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 * is a 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, 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 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 each time. * Therefore, the conditional mean can be regarded as the output value y * . A certain 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, which is specifically implemented as follows:

[0054] The weight of each known station is calculated based on the GPR and velocity uncertainty. Within each Delaunay triangle, the weights of the known stations are recalculated taking into account the velocity uncertainty.

[0055] In the embodiment, the velocity uncertainty of the measuring station on each Delaunay triangle vertex is considered, and the vertices are re-weighted. The jth 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 the present 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 means adding the absolute value weights of the symbol samples in ascending order to the smallest symbol sample value greater than the threshold T0.

[0063] Based on the above theoretical foundation, the embodiment of the present invention proposes a GNSS imaging method for detecting short-wavelength deformation to construct a GNSS image based on GPR, which is specifically implemented 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 estimated value of the vertical velocity of 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 estimated value of the vertical velocity;

[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, i.e., 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, according to computing resources, the above processing can be performed on multiple points to be estimated in parallel, or one point to be estimated can be sequentially processed.

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

[0072] The final weights of the known stations are obtained by combining the probability advantage of Gaussian process regression and the corresponding velocity uncertainty, and finally the vertical velocity estimate of the point to be estimated is obtained using the generalized weighted median. Based on the estimated vertical velocity of each grid point, the estimated results of all grid points in the study area can be plotted to construct a GNSS image.

[0073] In view of the problem that the traditional GNSS imaging method (GIM) has excessive smoothing when processing velocity peaks, which makes it difficult to detect short-wavelength deformation, the present invention innovatively uses Gaussian process regression (GPR) to describe the spatial relationship between adjacent station pairs, and uses their uncertainties to reweight the velocities of known stations to convert discrete velocity fields into continuous velocity fields. To facilitate the explanation of the technical effects of the present invention, the results of the chessboard simulation experiment of 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 checkerboards from 1°×1° to 4°×4° in the southwestern United States

[0075]

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

[0077]

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

[0079] In specific implementation, the method proposed in the technical solution of the present invention can be implemented by technical personnel in this field using computer software technology to realize automatic operation process. System devices for implementing the method, such as computer-readable storage media storing the corresponding computer program of the technical solution of the present invention and computer equipment running the corresponding computer program, should also be within the scope of protection of the present invention.

[0080] The GNSS imaging electronic device for measuring short wavelength deformation provided by the present invention is described below. 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 communication interface, a memory, and a communication bus, wherein the processor, the communication interface, and the memory communicate with each other through the communication bus. The processor may call logic instructions in the memory to execute the GNSS imaging method for measuring short-wavelength deformation, which mainly includes the software processing part in the above steps.

[0082] In addition, the logic instructions in the above-mentioned memory can be implemented in the form of software functional units and can be stored in a computer-readable storage medium when sold or used as an independent product. Based on such an understanding, the technical solution of the present invention, in essence, 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, which is stored in a storage medium and includes several instructions for a computer device (which can be a personal computer, a server, or a network device, etc.) to perform 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 further 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, wherein the units described as separate components may or may not be physically separated, and the components displayed as units may or may not be physical units, that is, they may be located in one place, or they may be distributed on multiple network units. Some or all of the modules may be selected according to actual needs to achieve the purpose of the scheme of this embodiment. Ordinary technicians in this field can understand and implement it without paying creative labor.

[0086] Through the description of the above implementation methods, those skilled in the art can clearly understand that each implementation method can be implemented by means of software plus a necessary general hardware platform, and of course, it can also be implemented by hardware. Based on this understanding, the above technical solution is essentially or the part that contributes to the prior art can be embodied in the form of a software product, and the computer software product can be stored in a computer-readable storage medium, such as ROM / RAM, a disk, an optical disk, etc., including a number of instructions for 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 some 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 embodiments of the present invention.

Claims

1. A GNSS imaging method for detecting short-wavelength deformation, characterized in that: For each point to be estimated, the 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 weight based on the spatial structure includes first obtaining known stations connected to the point to be estimated based on the Delaunay triangle, then obtaining preliminary weights using Gaussian process regression, and obtaining final weights of the connected known stations after re-weighting based on velocity uncertainty; 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 process based on the final weights of the known measuring stations connected to 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 regular grids, and each grid point is used as a point to be estimated, thus obtaining a set of points to be estimated. After the vertical velocity estimate of each grid point is extracted, 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 on 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. The GNSS imaging method for detecting short-wavelength deformation according to claim 1, characterized in that: The implementation method of the generalized weighted median filter processing is to use the vertical velocity value of the known measuring station as a sample, transfer the positive and negative properties of the weight to the sample, and the result is called a sign sample. The sign sample and the corresponding absolute value weight are input into the generalized weighted median filter to output the vertical velocity estimate of the point to be estimated.

8. 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 as described in any one of claims 1 to 7 is implemented.

9. 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 as claimed in any one of claims 1 to 7 is implemented.

10. 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 as claimed in any one of claims 1 to 7 is implemented.

Images