Spatial Filtering Method Based on Principal Component Analysis of Delaunay Triangulation Network

Abstract

The present invention relates to a common mode error extraction technology. The present invention aims to solve the problem that the method for extracting CME is only effective when the spatial response of CME is close to uniform, and provides a spatial filtering method based on the principal component analysis method of Delaunay triangulation. The technical solution can be summarized as follows: First, perform Delaunay triangulation on each site in the region to obtain multiple small triangulations, and then use the principal component analysis method to extract the common mode error in each divided small triangulation. The beneficial effect of the present invention is that it takes into account the irregularity of the common mode error in a large region, making the calculation of the common mode error in the large region more accurate and reasonable, and is applicable to the extraction of the common mode error.

Description

Technical Field

[0001] The present invention relates to CME (Common Mode Error) extraction technology, and particularly to a spatial filtering method based on the principal component analysis method of Delaunay triangulation network. Background Art

[0002] The errors in the regional GPS network time series are caused by various biases and noises, some of which are site-specific and others are region-related. However, the region-related errors are not negligible and are usually referred to as common mode errors (CME). Generally, CME will affect the reliability of GPS applications and even lead to misinterpretations of certain geophysical phenomena, such as transient deformation, subtle tectonics and displacements. Although CME is well-known in the regional GPS network, its origin and spatial characteristics are still under investigation. Currently, most studies focus on uncorrected environmental loading effects, GPS satellite orbit deviations, inaccurate atmospheric modeling, systematic errors and any combination thereof for CME.

[0003] As is well-known, there are several traditional global spatial filtering methods that have been used to remove CME. Wdowinski et al. first introduced the stacking filtering technique to accurately detect the coseismic and post-seismic deformations of the 1992 Mw 7.3 Landers earthquake in Southern California. Then, a weighted stacking method was proposed to calculate CME. Márquez Asua and De Mets also used the baseline length and time series as weights when studying the large geographic aperture in Mexico. Tian et al. developed the weighted spatial filtering (CWSF) method to extract the common mode components from a dense regional GPS network. The correlation between sites as weights focuses on filtering a larger regional network. In addition, the principal component analysis method (PCA) was proposed, which can extract CME with spatial response and is applicable to large networks. Based on PCA, Shen et al. proposed an improved PCA to extract the common mode error from incomplete position time series, and Li proposed a weighted spatio-temporal filtering (WSF) method based on PCA, which takes into account the positioning form error of the daily solution and the data gaps in the time series. In addition, some scientists regard CME extraction as a blind signal separation problem and use independent component analysis (ICA) to decompose the residual time series into components that are as statistically independent of each other as possible. The above methods are only effective when the spatial response of CME is close to uniform. However, recent studies have shown that CME has non-uniform behavior because in a large area, the common spatial correlation may not be strictly uniform, so a new filtering method is needed to extract CME that is considered to have non-uniform behavior in the regional GPS network. Summary of the Invention

[0004] The object of the present invention is to solve the problem that the current method for extracting CME is only effective when the spatial response of CME is close to uniform, and a spatial filtering method based on the principal component analysis method of Delaunay triangulation is provided.

[0005] The technical solution adopted by the present invention to solve the above technical problems is a spatial filtering method based on the principal component analysis method of Delaunay triangulation, including the following steps: Step 1, perform Delaunay triangulation on each station in the area to obtain a plurality of small triangulations;

[0006] Step 2, adopt the principal component analysis method to extract the common mode error in each divided small triangulation.

[0007] Specifically, to explain how to perform Delaunay triangulation on each station in the area, in Step 1, the method for performing Delaunay triangulation on each station in the area is as follows: Given the three-dimensional coordinates of each station in the area, assume that a set of position points P = {p 1 , p 2 , p 3 ,..., p n} is given in the Euclidean space, and n > 3, where the set P is non-linear. If x ∈ P, then there is:

[0008]

[0009] Among them, R k , as the cell element of the Thiessen polygon, is a polygon formed by connecting the perpendicular bisectors of the side lengths of each small triangulation formed, and d(x, p i ) is the Euclidean distance between x and p i , and x represents the element in P;

[0010] For any station, obtain the stations corresponding to each R k adjacent to this R k as the adjacent stations of this station, connect this station and the two adjacent stations closest to it to form a triangulation, and obtain a plurality of small triangulations.

[0011] Further, to explain how to adopt the principal component analysis method to extract the common mode error in each divided small triangulation, in Step 2, the method for adopting the principal component analysis method to extract the common mode error in each divided small triangulation is as follows:

[0012] Step 201, construct an (m × n) matrix X of the cleaned residual sequence for n stations and m epochs in the Delaunay triangulation, where m ≥ n, and the epoch refers to the time;

[0013] Step 202. For any small triangular network, there are three matrices X, each X having n columns and m rows, where n = 3. Each matrix X represents the north-direction coordinate component, east-direction component, and vertical component of m epochs of n stations. After performing the orthogonal transformation in the principal component analysis method, the principal component analysis (PCs) is decomposed into:

[0014]

[0015] where α k is the k-th principal component of X, j represents the j-th station, i represents the i-th epoch. Perform eigen-decomposition on B to calculate the eigenvalues and corresponding eigenvectors. v k is the k-th eigenvector of the positive definite symmetric matrix B, and the definition of B is: , t k , t i respectively refer to the k-th row and the i-th row of X;

[0016] Arrange the eigenvalues in descending order to obtain the number of the principal component energy spectrum, denoted by p;

[0017] Calculate the common mode error of this small triangular network. The calculation formula is:

[0018]

[0019] Specifically, since a station may be located in two or more small triangular networks, there may be different CMEs at the stations of the overlapping points. Therefore, to calculate the CME of the overlapping points, step 2 also includes the following steps:

[0020] Step 203. Calculate the distance from the overlapping point to the centroid of the small triangular network containing this overlapping point, denoted as: , where (x i , y i , z i ) represents the coordinates of the overlapping point, and (x j , y j , z j ) represents the coordinates of the centroid of this small triangular network;

[0021] Step 204. Calculate the common mode error CME of the overlapping point through the formula. The calculation formula is:

[0022]

[0023] where C j is the common mode error of the j-th small triangular network, and n represents the number of small triangular networks to which this overlapping point belongs.

[0024] The beneficial effects of the present invention are as follows. In the solution of the present invention, by dividing a large area according to the Delaunay triangulation rule, then calculating the common-mode error in each triangulation respectively, and finally distributing the common-mode error in the triangulation to each site, such an approach takes into account the irregularity of the common-mode error in the large area, making the calculation of the common-mode error in the large area more accurate and reasonable. Description of the Drawings

[0025] Figure 1 is a flowchart of the spatial filtering method of the principal component analysis method based on the Delaunay triangulation of the present invention.

[0026] Figure 2 is a schematic diagram of the Thiessen polygon in the embodiment of the present invention.

[0027] Figure 3 is a schematic diagram of the triangulation network divided in the embodiment of the present invention. Detailed Embodiment

[0028] The technical solution of the present invention will be described in detail below in conjunction with the embodiments and the drawings.

[0029] The spatial filtering method of the principal component analysis method based on the Delaunay triangulation of the present invention includes the following steps:

[0030] Step 1: Perform Delaunay triangulation on each site in the area to obtain a plurality of small triangulations.

[0031] To explain how to perform Delaunay triangulation on each site in the area, in this step, the method for performing Delaunay triangulation on each site in the area is preferably: Given the three-dimensional coordinates of each site in the area, assuming that a set of position point sets P = {p 1 , p 2 , p 3 ,..., p n} is given in the Euclidean space, and n > 3, where the set P is non-linear. If x ∈ P, then there is:

[0032]

[0033] Among them, R k , as the cell element of the Thiessen polygon, is a polygon formed by connecting the perpendicular bisectors of the side lengths of each small triangulation formed. For any site, it is intended to connect the site with other sites respectively, obtain each connecting line segment, and take the perpendicular bisectors of each connecting line segment respectively. The intersection of the perpendicular bisectors obtains a polygon that only contains the site, which is R k , and d(x, p i ) is the Euclidean distance between x and p i , and x represents an element in P;

[0034] For any site, obtain the R corresponding to that site k Each adjacent R k The corresponding sites are used as the adjacent sites of that site. Connect the site and the two adjacent sites closest to it to form a triangular network, that is, obtain the divided triangular network, and get multiple small triangular networks, as Figure 3 shown

[0035] such as Figure 2 shown, the formation process of R k takes the site 1 in Figure 2 as an example. First, connect site 1 to site 2, site 3, site 4, site 5, and site 6 respectively to get 6 connecting lines. Then, take the perpendicular bisectors of the 6 connecting lines respectively, extend each perpendicular bisector to make them intersect, and 6 intersection points can be obtained. Connect the 6 intersection points to get a polygon, and there is only one site 1 in this polygon. Therefore, this polygon is the R k of site 1. It can be seen that the line segments connecting the 6 intersection points are the corresponding partial line segments of the above perpendicular bisectors. From this, it can be known that site 2, site 3, site 4, site 5, and site 6 are the adjacent sites of site 1

[0036] Delaunay triangulation (DT) is used to obtain tetrahedrons from dense positions in the WNAM network for two main reasons: 1) According to the three-dimensional coordinates of each site in the network, divide the large coverage area into different spatial scales; 2) Small tetrahedrons are generated by dense site stations, that is, the external environment of the site in any small tetrahedron is almost unchanged. In addition, there are no any generating points in any small tetrahedron, and DT is unique in the Euclidean plane, which means that DT is very suitable for dividing spatial patterns

[0037] Step 2: In each divided small triangular network, use the principal component analysis method to extract the common mode error

[0038] The principal component analysis method (PCA) is a tool that uses orthogonal transformation to convert the residual time series into principal components (PCs) and reconstruct the residual time series

[0039] To explain how to use the principal component analysis method to extract the common mode error in each divided small triangular network, in this step, the method of using the principal component analysis method to extract the common mode error in each divided small triangular network is preferably

[0040] Step 201: Construct an (m×n) matrix X of the cleaned residual series for n stations and m epochs in the Delaunay triangular network, where m≥n, and the epoch refers to the time

[0041] Step 202. For any small triangular network, there are three matrices X, each with n columns and m rows, where n = 3. Each matrix X represents the north coordinate component, east coordinate component, and vertical component of m epochs of n stations. After performing the orthogonal transformation in the principal component analysis method, the principal component analysis (PCs) is decomposed into:

[0042]

[0043] where α k is the k-th principal component of X, j represents the j-th station, i represents the i-th epoch. Perform eigenvalue decomposition on B to calculate the eigenvalues and corresponding eigenvectors. v k is the k-th eigenvector of the positive definite symmetric matrix B, and the definition of B is: , t k , t i respectively refer to the k-th row and the i-th row of X;

[0044] Arrange the eigenvalues in descending order to obtain the number of the principal component energy spectrum, denoted by p;

[0045] Calculate the common mode error of this small triangular network. The calculation formula is:

[0046] .

[0047] Since a station may be located in two or more small triangular networks, there may be different CMEs at the stations of the overlapping points. Therefore, it is necessary to calculate the CME of the overlapping points. Also, because the CME is related to spatial correlation, and the correlation is approximately inversely proportional to the distance from the overlapping point to the centroid of the small tetrahedron containing the overlapping point. So when calculating the CME of the overlapping points, the distance is regarded as a weighting factor. Then in this step 2, the following steps may also be included:

[0048] Step 203. Calculate the distance from the overlapping point to the centroid of the small triangular network containing this overlapping point, denoted as; , where (x i , y i , z i ) represents the coordinates of the overlapping point, and (x j , y j , z j ) represents the coordinates of the centroid of this small triangular network;

[0049] Step 204. Calculate the common mode error CME of the overlapping point through the formula. The calculation formula is:

[0050]

[0051] where C jis the common mode error of the j-th small triangular network, and n represents the number of small triangular networks to which the overlapping point belongs.

[0052] The following is the experimental analysis of the spatial filtering method using the above-mentioned principal component analysis method based on the Delaunay triangular network:

[0053] Two regional GPS network time series sets were used, including data from 881 stations between 2010 and 2016: the cleaned residual series and the filtered residual series. They are from SOPAC [(http: / / garner.ucsd.edu / pub / timeseries / measures / ats / )] at UCSD, California, USA. The tectonic signals, non-tectonic signals, instrument jumps, and post-seismic displacements of some strong earthquakes have been removed from the cleaned residual series of the WNAM network. Since the missing rate of the residual time series reaches 2.96%, the third-order spline interpolation method is used to fit the missing values. The common mode error (CME) is extracted from the cleaned residual series using PCA, and then compared with the filtered residual series.

[0054] In PCA, low-order PCs may not be able to accurately represent the CME information of the original residual sequence, while high-order PCs may over-filter the residual time series, even leading to incorrect CMEs. To extract appropriate CMEs, refer to Dong's criterion (see "Dong D, Fang P, Bock Y, et al. Spatiotemporal filtering using principal component analysis and Karhunen-Loeve expansion approaches for regional GPS network analysis[J]. Journal of Geophysical Research, 2006, 111(B3):B03405"), and take the first PC as the common-mode component, dividing each eigenvector by the maximum element in the tetrahedron (the element with the largest absolute value). The results show that the average values of the first PCs in all tetrahedrons are 70.0%, 62.4%, and 72.1% of the total variance of the NEU components respectively. Using PCA filtering, the first PC (the first PC refers to the first principal component obtained by the principal component analysis method, i.e., the principal component energy spectrum p = 1. In fact, in this invention, multiple principal components need to be calculated together, i.e., the principal component energy spectrum p > 1) accounts for 43.5%, 28.2%, and 41.5% respectively. Obviously, the first PC using PCA filtering may not be able to accurately represent the common-mode component of the original residual sequence. In addition, the maximum spatial responses appear in the NEU components of 32.8%, 20.3%, and 46.9% of the tetrahedrons respectively. Therefore, it can be concluded that the spatial pattern is non-uniform among the PCs of all tetrahedrons.

[0055] For comparison, the residual sequences of DT-PCA filtering (i.e., the spatial filtering method of the principal component analysis method based on Delaunay triangulation of this invention) and PCA filtering related to the TBLP station are selected. By comparison, it can be seen that the residual sequence using the PCA filtering method is greater than that of the DT-PCA filtering method. DT-PCA filtering significantly weakens the seasonal variation, especially in the vertical direction. In PCA and DT-PCA filtering, the root mean square (RMS) of the vertical residual sequence is 4.16 mm and 3.15 mm respectively. In other words, compared with the original residual time series, the dispersion of the vertical residual sequence is reduced by 43.6% and 57.3% respectively. This means that the DT-PCA filtering method can more effectively reduce the regional effect. The number of stations where DT-PCA filtering has better performance than PCA filtering is 70.6%, 80.5%, and 77.6% respectively. This means that DT-PCA filtering has better performance across the entire network.

[0056] According to the above analysis, it can be seen that CME is a significant differential behavior in terms of spatial scale. This DT-PCA filtering method can better separate CME, more effectively suppress regional effects; at the same time, it can provide high-precision time series. The DT-PCA filtering method may have better performance across the entire network.

[0057] To further statistically prove the performance of the DT-PCA filtering method, the -RMS and correlation of the entire network are shown. In terms of RMS, the average RMS of 881 stations is reduced by 49.5% (north), 42.3% (east), and 51.4% (vertical), respectively. Since 46.9% of the tetrahedrons with the largest spatial response appear in the vertical component, the RMS is concentrated on the vertical component. Obviously, through DT-PCA filtering, the root mean square values of most stations will be greatly reduced. It can be seen that the root mean square values of PCA filtering are concentrated in the range of 8 - 12 mm, while those of DT-PCA filtering are concentrated in the range of 0 - 6 mm. In addition, compared with the RMS (vertical) map corresponding to PCA filtering, there are no abnormal parts in the RMS (vertical) map corresponding to DT-PCA filtering.

[0058] In addition, the Pearson correlation coefficients of the residual sequences between the TBLP of the three components and other stations using PCA and DT-PCA filtering are obtained and linearly fitted. Compared with other methods, linear fitting is less sensitive to data outliers. As for PCA filtering, the inter-station correlation after its linear fitting shows an approximate inverse proportion to the distance. Compared with PCA filtering, the inter-station correlations using DT-PCA filtering in the three directions are significantly smaller and close to zero, which means that the common mode error has been accurately eliminated. Obviously, DT-PCA filtering is more suitable for eliminating the correlation of residual sequences across the entire network.

[0059] From the analysis of the root mean square value and correlation, it can be seen that most of the filtered residual sequences have more reliable accuracy, and the inter-station correlation is more effectively eliminated when using DT-PCA filtering. Obviously, when no strictly homogeneous spatial pattern is found in the regional network, DT-PCA filtering is superior to PCA filtering.

[0060] To sum up, in this example, it is first assumed that CME has non-uniform behavior in terms of spatial scale, and the DT-PCA filtering method is proposed to extract CME. The performance of DT-PCA filtering is first manifested in the filtered residual sequences and power spectra related to TBLP. Then, the 881 GPS residual time series within 6 years of WNAM are analyzed using PCA filtering and DT-PCA filtering methods, and the RMS and inter-station correlation are compared. The conclusions are as follows:

[0061] (1)The results show that CME has significant non-uniformity on the spatial scale. Through DT-PCA filtering, the periodic error of the residual time series can be more effectively weakened. It is the best choice for high-precision GPS applications.

[0062] (2)Compared with PCA filtering, DT-PCA filtering can remove CME in the regional GPS network more reliably and accurately. When the common spatial correlation of the entire region may not be strictly consistent, it is very suitable for fitting CME.

Claims

1. A spatial filtering method based on the principal component analysis method of Delaunay triangulation, characterized in that, it includes the following steps: Step 1: Divide Delaunay triangulation for each station in the area to obtain multiple small triangulations; Step 2: Use the principal component analysis method to extract the common mode error in each divided small triangulation; In the said Step 1, the method for dividing Delaunay triangulation for each station in the area is: The three-dimensional coordinates of each site within a known area. Assume that a set of position point sets P = {p 1 , p 2 , p 3 ,..., p n} is given in a Euclidean space, and n > 3. The set P is non-linear. If x ∈ P, then there is: Among them, R k As the cell element of the Thiessen polygon, it is a polygon formed by connecting the perpendicular bisectors of the side lengths of each small triangular network. d(x, p i ) is the Euclidean distance between x and p i , where x represents an element in P; For any site, obtain the R corresponding to that site k Adjacent Rs k The corresponding sites are used as the adjacent sites of that site. Connect the site with its two closest adjacent sites to form a triangular network, obtaining multiple small triangular networks; In the said Step 2, the method for using the principal component analysis method to extract the common mode error in each divided small triangulation is: Step 201: Construct an m×n matrix X of the cleaned residual sequence for n stations and m epochs in the Delaunay triangulation, where m≥n, and the epoch refers to the time; Step 202: For any one small triangulation, there are three matrices X, each X has n columns and m rows, n = 3, and each matrix X represents the north direction coordinate component, east direction component and vertical component of m epochs of n stations respectively. After performing the orthogonal transformation in the principal component analysis method, the principal component analysis is decomposed into: where α k is the k-th principal component of X, j represents the j-th site, i represents the i-th epoch, perform eigen decomposition on B to calculate the eigenvalues and corresponding eigenvectors, v k is the k-th eigenvector of the positive definite symmetric matrix B, and the definition of B is: , t k and t i refer to the k-th row and the i-th row of X respectively; Arrange the eigenvalues in descending order to obtain the number of the principal component energy spectrum, denoted by p; Calculate the common mode error of this small triangulation, and the calculation formula is: 。 2. The spatial filtering method based on the principal component analysis method of Delaunay triangulation according to claim 1, characterized in that, in Step 2, it further includes the following steps: Step 203: Calculate the distance from the overlapping point to the centroid of the small triangular mesh containing the overlapping point, denoted as: , where (x i , y i , z i ) represents the coordinates of the overlapping point, and (x j , y j , z j ) represents the coordinates of the centroid of the small triangular mesh; Step 204: Calculate the common mode error CME of the overlapping points through the formula, and the calculation formula is: Among them, C j is the common-mode error of the j-th small triangular network, and n represents the number of small triangular networks to which the overlapping point belongs.

Images