Deformation Information Extraction from Multi-GNSS Coordinate Series Based on EWT-ICA-R

Abstract

Global navigation satellite system (GNSS) has been widely used in many deformation monitoring fields in recent years and can achieve centimeter-level or even sub-centimeter-level real-time monitoring accuracy through the carrier phase double-differenced technique. However, this technique cannot eliminate or weaken multipath errors, which become the main error source for GNSS deformation monitoring. Therefore, extracting deformation information from coordinate series mixed with multipath errors has become a key issue for further improving the accuracy of GNSS deformation monitoring. In this paper, we propose an approach to overcome this issue called empirical wavelet transform-independent component analysis with reference (EWT-ICA-R). The specific process is as follows. First, EWT is employed to model the multipath errors from a priori GNSS coordinate series, and the model is input to ICA-R as a reference signal. Then, the GNSS deformation monitoring series mixed with multipath errors and deformation information is decomposed into sub-series of different scales using EWT, and these sub-series are input to ICA-R as multi-channel signals. Finally, ICA-R is used to calculate the input signals together to obtain the multipath errors in the GNSS deformation monitoring series and then subtract the multipath errors from the GNSS deformation monitoring series to obtain accurate deformation information. Experiments show the following: (1) For the vibration deformation experiments, the correlation coefficients between the deformation information extracted by the proposed method and the real values reached 0.981, 0.981, and 0.885 in the E, N, and U directions, respectively, and the corresponding root mean square errors decrease to 0.694 mm, 0.694 mm, 1.852 mm, respectively. (2) For the slow-deformation experiment, the correlation coefficients in the three directions were all higher than 0.98, and the corresponding root mean square errors decrease to 1.345 mm, 1.546 mm, and 3.866 mm, respectively. The experiments verified the feasibility of the proposed method to accurately extract deformation information, which makes it possible to obtain sub-millimeter GNSS deformation information and provide effective technical support for deformation monitoring in related fields.

1. Introduction

Deformation disasters occur widely in many fields, such as large buildings or structures, mining surfaces, and landslides. When the deformation reaches the bearing capacity of the deformable body, the deformation disaster seriously threatens the safety of people’s lives and properties. Therefore, it is necessary to use technical means to monitor deformation, obtain deformation information, and use a certain method to mine the deformation evolution mechanism to guide the prevention and control of deformation disasters [1,2,3].

As a high-precision positioning technology, the global navigation satellite system (GNSS) can provide continuous real-time positioning information, and has been widely used in related deformation monitoring fields [4,5]. Lovse et al. [6] pioneered the use of the global positioning system (GPS) to measure the structural vibrations of the Calgary Tower in Calgary, AB, Canada, and experimentally verified the feasibility of GPS as a structural monitoring technology. Guo and Ge [7] used GPS to successfully obtain the displacement and vibration frequencies of the Di Wang Tower in Shenzhen, Guangzhou, China. Gili et al. [8] discussed the applicability of GPS for monitoring the displacement of landslide surfaces, which has larger coverage and higher work efficiency than classical measurement methods. Lü et al. [9] employed GPS for land subsidence monitoring and obtained the GPS control network of the observation stations with sufficient accuracy. At the same time, GPS provides reliable data for studying the laws of surface deformation caused by underground coal mining.

Usually, GNSS real-time deformation monitoring is performed using the carrier phase double-differenced technique. This technique can eliminate the satellite clock error and the receiver clock error and greatly weaken the tropospheric error and ionospheric error (especially for short baseline solutions), but cannot mitigate the influence of multipath errors on positioning accuracy [10]. Therefore, multipath errors become the main error in the real-time deformation monitoring of GNSS. Owing to the spatiotemporal repeatability of GPS (or GLONASS, Galileo) operations, the sidereal filtering (SF) technique and its modified methods have been proposed for the real-time mitigation of multipath errors [11,12,13]. However, the constellation structures of different navigation satellite systems and their operating periods are different. The spatiotemporal repeatability of multi-GNSS satellites is difficult to obtain, and the standard SF technique is no longer strictly applicable [14]. Ye et al. [15] proposed an improved SF technique based on the “zero-average” hypothesis, analyzed the differences in the return period and multipath signals of medium earth orbit (MEO), geosynchronous orbit (GEO), and inclined geosynchronous orbit (IGSO) satellites, and solved the constellation difference of the BeiDou navigation satellite system (BDS). This study makes it possible to reduce the multipath errors of multi-GNSS using the SF technique. Dong et al. [16] proposed a multipath hemispherical map (MHM) to mitigate multipath errors based on the single-receiver dual-antenna technique and compared and analyzed it using the improved SF technique. The experimental results demonstrate that the MHM model can mitigate the influence of multipath errors more effectively. However, this model has higher requirements for GNSS hardware equipment and its practicability needs to be improved. Liu et al. [17] employed the MHM model and nearest-neighbor search to mitigate the multipath errors of BDS-2/BDS-3, and achieved a good effect for multiple consecutive days. Tao et al. [18] considered the highly nonlinear characteristics of multipath errors in the GNSS coordinate series and established a multipath error model based on the convolutional neural network (CNN)-long short-term memory (LSTM) method, which realizes the real-time mitigation of multipath errors.

Multipath errors are mainly affected by the geometric structure between the satellite, receiver, and reflector positions. In the application of deformation monitoring, owing to the continuous deformation of the monitored object, the original geometric structure is destroyed, and the multipath error model established by different methods is not strictly applicable [19,20]. Zheng et al. [21] used a Vondrak filter and cross-validation to mitigate the effect of multipath errors in GPS coordinate series mixed with vibration deformation and improved the final deformation accuracy by 20–40%. Dai et al. [22] performed multi-scale decomposition on the GNSS coordinate series mixed with multipath errors and deformation information based on empirical mode decomposition (EMD); that is, multiple virtual channel signals were obtained, and then the virtual channel signals and the prior GNSS coordinate series (regarded as a reference signal) were used for the independent component analysis with reference (ICA-R) calculation to mitigate the multipath errors. Experimental results verified the effectiveness of the proposed strategy. Although EMD can perform multi-scale decomposition of GNSS coordinate series and has been widely used in many fields, it still has some insurmountable defects, such as frequency aliasing, endpoint effects, and uncertainty in the number of decompositions, especially the latter, which seriously restricts the use of EMD for automatic or batch processing [23,24,25]. Gilles [26] proposed the empirical wavelet transform (EWT), a method capable of constructing adaptive wavelets to segment different modes of a signal by designing appropriate wavelet filters with a clear mathematical theoretical basis, which has the characteristics of high computational efficiency, good adaptivity, and theoretical adequacy compared with EMD. The EWT has been widely applied to the fault detection of mechanical bearings [27], seismic data mining [28], medical image monitoring [29], and other fields [30].

In view of the excellent multi-scale decomposition ability of EWT, we propose the use of EWT to perform multi-scale decomposition of the GNSS coordinate series mixed with deformation information and then use the ICA-R method to separate deformation information and multipath errors to achieve higher-precision deformation information acquisition.

GNSS deformation monitoring series usually contain information at different frequencies, which mainly consist of specific deformation information and GNSS system errors. However, systematic errors in GNSS short-range relative positioning usually include multipath errors and other non-modeled systematic errors, which usually exhibit significant frequency aliasing. Therefore, according to the characteristics of structural deformation and land subsidence, this paper simulates vibration deformation and slow deformation respectively, hoping to verify the separation effect of the proposed method on deformation information under different deformation scenarios and different frequency aliasing conditions. It is precisely because of this that the deformation information in two different scenarios is not superimposed.

2. Methods

2.1. Multipath Errors

In the process of signal propagation, the signal wave encounters a multipath error effect, such as the reflection of light and electromagnetic waves. In satellite positioning, the satellite signal propagates in the form of electromagnetic waves. During the propagation of GNSS signals, the vast majority of the signals enter the receiver in the form of direct transmission, but there is still an inevitable small part of the signals that enter the receiver after being reflected by the reflectors around the station, as shown in Figure 1. The direct and reflected signals enter the receiver simultaneously, and the two are superimposed to form an interference signal. There is a certain deviation between the interference and direct signals. This phenomenon is called the multipath error, which reduces the positioning accuracy of the GNSS to a certain extent.

The path deviation $S_{\mathrm{m}}$ of the satellite signal caused by the multipath effect can be expressed as [31]:

$$S_{\mathrm{m}}=\frac{\Delta {\phi }_{\mathrm{m}}}{2\pi }\lambda =\frac{\lambda }{2\pi }\arctan \left(\frac{\alpha \sin \left(\frac{4\pi d}{\lambda }\sin \theta \right)}{1+\alpha \cos \left(\frac{4\pi d}{\lambda }\sin \theta \right)}\right)$$

(1)

where $\Delta {\phi }_{\mathrm{m}}$ is the multipath error in the carrier phase measurement; $\lambda$ is the wavelength of the reflected signal; $\alpha$ is the reflection coefficient; d is the distance between the receiver antenna and the surrounding reflector; $\theta$ is the incident angle of the reflected signal.

From Equation (1), we find that multipath is mainly affected by d, $\theta$, and $\alpha$. If the receiver position does not change significantly and the surrounding environment of the receiver remains unchanged, the multipath error is mainly affected by $\theta$. In other words, multipath errors are strongly related to the satellite position. Because GNSS satellites (except GEO satellites of BDS) run periodically on the set orbits, multipath errors usually exhibit periodic variation characteristics. The operating periods of different satellite navigation systems are different, and the frequency bands of signals broadcast by different satellite navigation systems are also different, which causes multipath errors in multi-GNSS system positioning to be no longer periodic.

2.2. EWT

Gilles [26] proposed EWT, which essentially utilizes the adaptive characteristics of EMD and the controllable number of scales of the wavelet transform. In contrast to wavelet transform, EWT adaptively constructs basis functions based on the basis pursuit method and then uses the framework of EMD to perform wavelet transform. Essentially, the EWT method is an adaptive segmentation and fully adaptive multi-scale decomposition method.

The main steps of the multiscale decomposition of the GNSS coordinate series mixed with deformation information using EWT are as follows.

(1)

Spectral analysis of the coordinate series $x(t)$ mixed with deformation information to obtain Fourier spectrum $F\left(\omega \right)$.

(2)

Divide $F\left(\omega \right)$ into K continuous segments and determine its boundary $\mathsf{\Omega }={\left\{{\omega }^k\right\}}_{k=0,1,\cdots ,K}$, where $\omega$ is the frequency, and ${\omega }^k$ is the kth boundary frequency (where ${\omega }^0=0$ and ${\omega }^K=\pi$).

(3)

Constructing low-pass and band-pass filters with the empirical scale function ${\phi }_k\left(\omega \right)$ and empirical wavelet function ${\psi }_k\left(\omega \right)$, as shown in Equations (2) and (3).

$${\phi }_k\left(\omega \right)=\left\{\begin{array}{cc}1,& \left|\omega \right|\le \left(1-\gamma \right){\omega }^k\\ \cos \left[\frac{\pi }{2}R\left(\gamma ,{\omega }^k\right)\right],& \left(1-\gamma \right){\omega }^k\le \left|\omega \right|\le \left(1-\gamma \right)\\ 0,& \mathrm{others}\end{array}\right.$$

(2)

$${\psi }_k\left(\omega \right)=\left\{\begin{array}{cc}1,& \left(1+\gamma \right){\omega }^k\le \left|\omega \right|\le \left(1-\gamma \right){\omega }^{k+1}\\ \cos \left[\frac{\pi }{2}R\left(\gamma ,{\omega }^{k+1}\right)\right],& \left(1-\gamma \right){\omega }^{k+1}\le \left|\omega \right|\le \left(1-\gamma \right){\omega }^{k+1}\\ \sin \left[\frac{\pi }{2}R\left(\gamma ,{\omega }^k\right)\right],& \left(1-\gamma \right){\omega }^k\le \left|\omega \right|\le \left(1+\gamma \right){\omega }^k\\ 0,& \mathrm{others}\end{array}\right.$$

(3)

where $\gamma$ is a parameter that ensures that the two continuous transforms do not overlap, and $R\left(\gamma ,{\omega }^k\right)$ is a function that can be expressed as

$$R\left(\gamma ,{\omega }^k\right)=\beta \left(\frac{1}{2\gamma {\omega }^k}\left(\left|\omega \right|-\left(1-\gamma \right){\omega }^k\right)\right)$$

(4)

where

$$\beta \left(q\right)=\left\{\begin{array}{cc}0,& q\le 0\\ 1,& q\ge 1\\ q^4\left(35-84q+70q^2-20q^3\right),& q\in \left(0,1\right)\end{array}\right.$$

(5)

(4)

The inner product of the empirical scale function (Equation (2)) and the GNSS coordinate series is the scale coefficient $W\left(0,t\right)$, which can be expressed as

$$W\left(0,t\right)=〈f\left(t\right),{\phi }_1\left(t\right)〉={\displaystyle \int f\left(\tau \right)}·\overline{{\phi }_1\left(\tau -t\right)}d\tau =F^{-1}\left[f\left(\omega \right)\widehat{{\phi }_1}\left(\omega \right)\right]$$

(6)

where $\tau$ is the time node, $〈·〉$ is the inner product operation, $\widehat{{\phi }_1}\left(\omega \right)$ is the Fourier transform of ${\phi }_1\left(t\right)$, and $\overline{{\phi }_1\left(\tau -t\right)}$ is the complex conjugate of ${\phi }_1\left(\tau -t\right)$.

In addition, the inner product of the empirical wavelet function (Equation (3)) and the GNSS coordinate series is the empirical wavelet coefficient $W\left(k,t\right)$, which can be expressed as

$$W\left(k,t\right)=〈f\left(t\right),{\psi }_k\left(t\right)〉={\displaystyle \int f}\left(\tau \right)·\overline{{\psi }_k\left(\tau -t\right)}d\tau =F^{-1}\left(\omega \right)\left[f\left(\omega \right)\widehat{{\psi }_k}\left(\omega \right)\right]$$

(7)

where $\widehat{{\psi }_k}\left(\omega \right)$ is the Fourier transforms of ${\psi }_k\left(t\right)$; $\overline{{\psi }_k\left(\tau -t\right)}$ is the complex conjugate of ${\psi }_k\left(\tau -t\right)$.

(5)

Equation (8) is used to reconstruct the data.

$$f\left(t\right)=W\left(0,t\right)\ast {\phi }_1\left(t\right)+{\displaystyle \sum _{k=1}^{n}W\left(k,t\right)}\ast {\psi }_k\left(t\right)$$

(8)

where $\ast$ is the convolution operator symbol.

2.3. ICA-R

The Independent Component Analysis (ICA) is a computational method for separating multivariate signals into additive sub-components. The research object of ICA is independent non-Gaussian data, and the multi-channel mixed data are separated one by one, and each time the signal is separated into an independent component until the iterative component changes do not change or change very little [32]. The input signal to be processed by ICA must be multiple and greater than or equal to the number of output signals. The purpose of ICA is to determine the unmixing matrix w and then linearly transform the input signal to obtain the output signal. The main calculation steps for w are as follows.

(1)

Centering, white-noising, and decorrelating the input signal to obtain a modified input signal with zero mean.

(2)

Choose an initialized w, which should have a unit norm, and update w according to Equation (9).

$$w_{k+1}=E\left(xg\left(w_k{}^{\mathrm{T}}x\right)\right)-E\left(g^{\prime }\left(w_k{}^{\mathrm{T}}x\right)\right)w_k$$

(9)

where x is the observation signal, and $g\left(·\right)$ is the derivative of a non-quadratic function $G\left(·\right)$.

(3)

normalize $w_{k+1}$

$$w_{k+1}=\frac{w_{k+1}}{‖w_{k+1}‖}$$

(10)

(4)

Check whether the results converge. Otherwise, return to Step (3) to continue the update iteration.

ICA can be used to recover the independent components of the source signal using only the mixed signal of the source signal without knowing the statistical distribution type and mixing parameters of the source signal. However, traditional ICA usually has the defects of uncertain output-signal ordering and low estimation efficiency. Some scholars have introduced prior knowledge into blind source separation and obtained the result of improved performance, which is the so-called semi-blind ICA. Because of the diversity of prior information and its utilization, there are many types of semi-blind algorithms. As a typical semi-blind ICA, ICA-R has the advantage of determining the output order and has high estimation efficiency [33]. The basic idea is to introduce the reference signal into the traditional ICA algorithm to transform the blind source separation problem into a constrained optimal solution problem, in which the mixed signal is obtained through the iterative Newton approximation process. Finally, a certain certainty is extracted from the mixed signal. The method used in this study was a single ICA-R algorithm, and its algorithmic model was [34].

$$\mathrm{Objective} \mathrm{function}:J\left(w\right)\approx \rho {\left[E\left\{G\left(w^{\mathrm{T}}z\right)\right\}-E\left\{G\left(v\right)\right\}\right]}^2$$

(11)

$$\mathrm{Restrictions}: \left\{\begin{array}{l}s\left(w\right)=E\left\{y^2\right\}-1=0\\ g\left(w\right)=\epsilon \left(y,r\right)-\xi \le 0\end{array}\right.$$

(12)

where $y=w^{T}z$ is the output signal; $w$ is the initial weight; r is the reference signal; $\rho$ is the normal quantity; v is a Gaussian variable with zero mean and unit variance; $J\left(w\right)$ is the negative entropy of $y$; $G\left(·\right)$ is any non-quadratic equation; $g\left(·\right)$ is the derivative of $G\left(·\right)$; $s\left(w\right)$ is a constraint function that ensures that $J\left(w\right)$ and $w$ are bounded; $\xi$ is its threshold value to distinguish the target signal from other signals; and $\epsilon \left(y,r\right)$ is the generalized distance function between y and r.

To obtain the optimal solution of the objective function in Equation (12) under the restriction of Equation (11), a Lagrangian function is established:

$$H\left(w,\mu ,u,z\right)=J\left(y\right)+{\mu }^T\widehat{g}\left(y\right)+\frac{1}{2}\eta {‖\widehat{g}\left(y\right)‖}^2+u^{T}s\left(y\right)+\frac{1}{2}\eta {‖\widehat{s}\left(y\right)‖}^2$$

(13)

where $\mu$ and $u$ are the Lagrangian factors. The ICA-R cell is obtained by polarizing the Lagrangian function in Equation (13), and the algorithm is given in Equation (14).

The maximum value operation is performed using Equation (13), and the unit value of ICA-R can be obtained as follows:

$$\Delta w=-\eta \left(\mathsf{\Phi }\left(y,\mu ,u\right)+\mathsf{\Gamma }\left(y,\mu ,u\right)+\mathsf{\Psi }\left(y,\mu ,u\right)\right){\displaystyle \sum x{x}^{-1}}$$

(14)

where

$$\left\{\begin{array}{c}\mathsf{\Phi }\left(y,\mu ,u\right)=\frac{E\left\{G^{\prime }\left(y\right)x\right\}}{d\left(y,\mu ,u\right)}\\ \mathsf{\Gamma }\left(y,\mu ,u\right)=\frac{E\left\{G^{\prime }\left(y,\mu \right)x\right\}}{d\left(y,\mu ,u\right)}\\ \mathsf{\Psi }\left(y,\mu ,u\right)=4u\frac{\left(E\left\{y^2\right\}-1\right)E\left\{yx\right\}}{d\left(y,\mu ,u\right)}\end{array}\right.$$

(15)

where $d\left(y,\mu ,u\right)=E\left\{\mu g^{″}\left(y\right)\right\}+8u-E\left\{\rho G^{″}\left(y\right)\right\}$; $G^{\prime }\left(y\right)$, and $G^{″}\left(y\right)$ are the first- and second-order derivatives of $G\left(y\right)$, respectively, and $g^{″}\left(y\right)$ is the second-order derivative of $g\left(y\right)$.

2.4. Deformation Information Extraction Based on EWT-ICA-R

To accurately extract the deformation information in the coordinate series of the multi-GNSS, a specific data processing flow is designed based on EWT-ICA-R, as shown in Figure 2. First, multipath modeling is performed on the a priori GNSS coordinate series based on the EWT, and the obtained multipath model is input to ICA-R as the reference signal. The method in this paper separates and models the multipath errors from the a priori coordinate data, where the simulated vibrational and slow deformations are not added to the a priori data and are added to the subsequent GNSS coordinate sequence of DOY135-147. Then, the GNSS deformation monitoring series mixed with deformation information is decomposed into sub-series of different scales using EWT, and these sub-series are input to ICA-R as multi-channel signals. Finally, ICA-R is used to calculate the input signals together to obtain the multipath errors in the GNSS deformation monitoring series and then subtract the multipath errors from the GNSS deformation monitoring series to extract accurate deformation information. The multipath error model established in the first step as the input of the subsequent deformation information separation will directly affect the accuracy of its separation. Because of this, we need to use a method with a good noise reduction effect (e.g., EWT) to model the multipath error in the first step. As evidenced by the many kinds of literature, EWT has relatively excellent noise reduction ability, so we adopt this method to model multipath error in this paper.

With the excellent multi-scale decomposition capability of EWT, the GNSS coordinate series mixed with deformation information can be adaptively decomposed into multiple IMFs, and the segmentation mechanism of the Fourier spectrum is the key issue for multi-scale decomposition. The segmentation strategies include Locmax, Locamaxmin, Adaptive, and Scalespace. Considering the excellent frequency division capability of Scalespace, it was adopted in this study for subsequent data processing.

When using EWT-ICA-R for multipath error reconstruction, there is a lack of uncertainty in the amplitude of the multipath errors [31]. Therefore, a scale coefficient $a=\mathrm{std}\left(r\right)/\mathrm{std}\left(y\right)$ must be set to restore the correct amplitude, where $\mathrm{std}\left(·\right)$ is the standard deviation. After the multipath errors s are recovered (by $s=a\ast y$), they can be used as a subsequent reference signal to weaken them and accurately extract the deformation information.

3. Experiments and Results

3.1. Data Collection

The experiment was conducted on the roof of the building of the School of Spatial Information and Geomatics Engineering, Anhui University of Science and Technology, China. Two GNSS receivers of the same type were used for continuous long-term data collection. The distance between the two receivers was approximately 13 m, the height of the receiver antennas was approximately 2.3 m, and the height of the walls around the receivers was approximately 1.5 m (Figure 3). The receiver parameters are the Trimble BD980 board and AT300 antenna, which can receive GPS triple-frequency signals, GLONASS, and BDS dual-frequency signals. The base station and rover station can receive data from three systems, GPS, GLONASS, and BDS, simultaneously, with a sampling interval of 1 Hz and a satellite elevation cutoff of 15°. The original observation data of two GNSS receivers from DOY134 to DOY147 in 2018 were collected, of which DOY142 was unavailable owing to a power outage, with a total of 13 days of observation data.

The RTKlib software (http://www.rtklib.com/, accessed on 14 May 2018) was used as the data solution to obtain the GNSS coordinate series. Because the baseline distance is short, the height difference is extremely small, and the receiver types are the same, the tropospheric delay, ionospheric delay, and antenna phase deviation can be eliminated using the carrier phase double-differenced technique. Therefore, the multi-GNSS coordinate series is mainly composed of noise and multipath errors. Meanwhile, considering the periodic characteristics of multipath errors, we take the average value of the long-term coordinate series as the true value for subsequent data validation.

3.2. Multipath Errors Modeling

In this section, EWT is used to denoise the GNSS coordinate series of DOY134, and the multipath error model of DOY134 was established. The model was used as a reference signal for the identification of multipath errors and the extraction of deformation information in the following days.

The established multipath error model of DOY134 and its power spectral density are shown in Figure 4. The multipath error model is consistent with the original GNSS coordinate series trend, which intuitively shows that the EWT can effectively eliminate the influence of noise in the GNSS coordinate series. In terms of power spectral density, the characteristics of the multipath error model and the low-frequency region of the original GNSS coordinate series are almost the same, which is consistent with the characteristics that the multipath error is mainly low-frequency.

3.3. Vibration Deformation Experiments and Results

To test the performance of the proposed method for the accurate extraction of deformation information, we designed two typical deformation scenarios: a vibration deformation experiment and a slow deformation experiment. In this section, we test how to efficiently extract the vibration deformation in a GNSS coordinate series mixed with the vibration deformation, multipath errors, and noise.

To accurately verify the extraction effect of deformation information, a harmonic signal (simulating vibration deformation) that can produce frequency aliasing with multipath errors is added to the GNSS coordinate series from DOY135 to DOY147, with an amplitude of 5 mm and a frequency of 0.015 Hz. Owing to space limitations, the E direction of DOY135 is taken as an example to introduce the relevant processing process, as shown in Figure 5. Figure 5 shows the GNSS deformation monitoring sequence of DOY135 after adding the simulated vibration deformation, and the frequency aliasing with multipath errors can be seen after zooming in locally.

The key issue in accurately extracting deformation information is the separation of the multipath errors. Therefore, we separate the multipath errors in the GNSS coordinate sequences mixed with the vibration deformation from DOY135 to DOY147.

To facilitate the observation of the accuracy of the separated multipath errors, they are compared with the corresponding original GNSS coordinate series. Figure 6 presents the original GNSS coordinate series and the separated multipath errors from DOY135 to DOY147 in the E direction; the white line in the figure is the separated multipath errors, and a certain value is added to the pair-by-day to facilitate the observation. As shown in Figure 6, the separated multipath errors from DOY135 to DOY147 are generally consistent with the variation law of the corresponding GNSS original coordinate series. However, they become smoother daily in the extracted multipath errors, especially in the series after DOY142, the same phenomenon occurs in the N and U directions. There are two primary reasons for this phenomenon. First, the separation accuracy of multipath errors continues to decrease with time, which is a common phenomenon, and owing to the lack of DOY142 GNSS data, DOY141 is required to separate the multipath errors of DOY143. Therefore, the multipath error separation accuracy of DOY143 and the following days are reduced to a certain extent.

After the separation of GNSS multipath errors, a vibration deformation series mixed with noise was obtained. To analyze and compare the obtained vibration deformation with its reference value (harmonic signal), the obtained vibration deformation series mixed with noises was denoised by the EWT, and the denoised series was used as the extracted vibration deformation. Figure 7 shows the extracted vibration deformation in the E direction of the first and last two days (DOY135 and DOY147) and the difference between it and the reference value. In general, the difference between the extracted vibration deformation information and the reference value is within 0.2–0.3 mm, reaching a sub-mm-level deformation information extraction. In addition, the data processing accuracies for the other time periods and directions were comparable.

Table 1. Root mean square errors of the extracted vibration deformation and their correlation coefficients with the corresponding harmonic signal.

DOY	E Direction		N Direction		U Direction
	R	RMS/mm	R	RMS/mm	R	RMS/mm
135	0.982	0.686	0.980	0.710	0.875	1.937
136	0.980	0.709	0.981	0.705	0.882	1.875
137	0.982	0.686	0.981	0.709	0.874	1.964
138	0.983	0.669	0.981	0.699	0.889	1.818
139	0.980	0.718	0.981	0.693	0.888	1.843
140	0.982	0.680	0.983	0.659	0.892	1.778
141	0.982	0.682	0.982	0.685	0.889	1.811
143	0.980	0.713	0.981	0.701	0.880	1.880
144	0.981	0.699	0.982	0.680	0.889	1.814
145	0.981	0.705	0.981	0.700	0.883	1.861
146	0.982	0.679	0.981	0.690	0.878	1.907
147	0.980	0.703	0.981	0.692	0.897	1.740
Mean	0.981	0.694	0.981	0.694	0.885	1.852

presents the statistical results for the accuracy and correlation coefficients of the extracted vibration deformations. It can be seen that the correlation coefficient between the extracted vibration deformation of DOY135 to DOY147 in the E and N directions and the corresponding reference values are maintained above 0.98 and that the U direction is maintained at approximately 0.88. The root mean square errors of the E and N directions are approximately 0.70 mm, whereas those in the U direction reach approximately 1.8 mm.

The correlation coefficients and root mean square errors in the U-direction were significantly worse than those in the other two directions. This is mainly because the amplitude of the multipath errors in the original GNSS coordinate series in the U-direction is significantly larger than that in the other two directions, and the power spectral density of the coordinate series in the U-direction is significantly larger than that in the other two directions in the low-frequency region (Figure 4). This increases the degree of aliasing between multipath errors and vibration deformations in the coordinate series, which is not conducive to the accurate separation of multipath errors and the accurate extraction of vibration deformations.

3.4. Slow Deformation Experiments and Results

To verify the deformation information extraction performance of EWT-ICA-R in the slow deformation application scenario, we added slow deformation information formed by the superposition of linear functions and trigonometric functions to the original GNSS coordinate series. The equation used is as follows:

$$y_{def}=\frac{-0.02\left(100+t\right)}{0.35n}+0.04\left(\sin \left(\frac{t}{2000}\right)+\cos \left(2+\frac{t}{3000}\right)\right)$$

(16)

where $y_{def}$ is the simulated slow deformation and $t\in [1,86400]$; n = 80. The simulated slow deformation was approximately 60 mm/day.

The simulated slow deformation information was added to the original GNSS coordinate series of DOY135-DOY147 respectively to form a simulated deformation monitoring series. Figure 8 shows the original GNSS coordinate series and simulated deformation monitoring series in the E direction of DOY135. Meanwhile, the multipath errors of DOY134 were used as a reference signal for multipath error separation and deformation information extraction in the subsequent simulation deformation series.

Based on the EWT-ICA-R method, multipath errors were separated from the simulated deformation monitoring series, and a deformation monitoring series mixed with noise was obtained. To verify the extraction effect of deformation information, we performed noise reduction processing on the obtained deformation monitoring series mixed with noise and compared and analyzed it with the simulated deformation series. Figure 9 shows the extracted deformation series of the first and last days (DOY135 and DOY147), the simulated deformation series, and the difference series between them. In general, the change trends of the extracted deformation series and the simulated deformation series were similar, and the differences were maintained within 5 mm. In addition, the difference between the extracted deformation series and the simulated deformation series of the last day (DOY147) was significantly greater than that of the first day (DOY135).

Table 2. Root mean square errors of the extracted deformation series and their correlation coefficients with the corresponding simulated deformation series.

DOY	E		N		U
	R	RMS/mm	R	RMS/mm	R	RMS/mm
135	0.998	0.992	0.998	1.052	0.988	2.675
136	0.998	1.342	0.997	1.336	0.980	3.370
137	0.997	1.246	0.997	1.396	0.984	3.190
138	0.997	1.363	0.996	1.544	0.988	2.906
139	0.997	1.474	0.995	1.597	0.987	2.995
140	0.997	1.488	0.994	1.871	0.986	3.117
141	0.998	1.232	0.997	1.389	0.976	4.256
143	0.997	1.392	0.997	1.812	0.980	4.973
144	0.998	1.089	0.995	1.644	0.980	4.475
145	0.995	1.612	0.997	1.372	0.978	5.082
146	0.998	1.379	0.994	1.928	0.987	4.537
147	0.996	1.525	0.997	1.605	0.986	4.815
Mean	0.997	1.345	0.996	1.546	0.983	3.866

presents the root mean square errors of the extracted deformation series and their correlation coefficients with the corresponding simulated deformation series in the three directions of E/N/U from DOY135 to DOY147. It can be seen that the correlation coefficients are all above 0.97, indicating that the extracted deformation series and the simulated deformation series maintain a high degree of consistency in shape, and the correlation coefficient does not change significantly over time. In addition, the RMS in the three directions of E/N/U increases with time, which indicates that the separation effect of the slow deformation information based on the EWT-ICA-R method decreases slowly with time.

4. Discussion

At present, the main constraint limiting GNSS millimeter-scale deformation monitoring is multipath errors. However, in the multi-GNSS coordinate domain, multipath errors no longer exhibit temporal repeatability; therefore, traditional SF is no longer applicable. In addition, when the position of the GNSS monitoring station changes continuously with the deformation of the monitored object, the geometric relationship between the original satellite, reflector, and GNSS receiver is destroyed, and the established multipath error model is no longer strictly applicable [14]. In this context, we comprehensively utilize the multi-scale decomposition and reconstruction characteristics of the EWT [26] and the blind source separation characteristics of the ICA-R [33] method and propose the use of the EWT-ICA-R method [22] to separate the multipath errors and deformation information in the deformation monitoring series.

The multipath errors were separated from the coordinate series by the EWT, and their power spectral density is shown in Figure 4. In the low-frequency region, the power spectral densities of the multipath errors are consistent with those of the original GNSS coordinate series, whereas, in the high-frequency region, the power spectral density of the multipath errors is significantly smaller than that of the original GNSS coordinate series, indicating that the EWT can accurately extract the multipath errors in the GNSS coordinate series.

Based on the EWT-ICA-R method, the multipath errors and vibration deformation were separated from the GNSS coordinate series mixed with the simulated vibration information, and the final accurate vibration deformation information was obtained. The correlation coefficient (approximately 0.98) between the vibration deformation extracted in the E and N directions and the simulated vibration information is significantly greater than that in the U direction (approximately 0.88), which is directly related to the fact that the geometric structure in the E and N directions is better than that in the U direction owing to the GNSS satellite constellation. In addition, owing to the difference in the geometric structure in different directions caused by the satellite constellation, the RMS in the U direction is significantly greater than that in the E and N directions. In addition, there was no significant difference in the correlation coefficient and RMS between DOY135 and DOY 147. The main reason for this is that the added vibration deformation was only 5 mm, and the multipath errors and noises dominated the GNSS coordinate series mixed with the simulated vibration information.

Slow deformation monitoring yields results similar to those of vibration deformation monitoring. The correlation coefficients between the extracted deformation information and the simulated deformation information are above 0.97, among which the E and N directions are above 0.99, which is slightly better than that of the U direction. In addition, with the continuous accumulation of deformation, the RMS of the three directions from DOY135 to DOY147 increases continuously, and the increase in the U direction was significantly larger than that in the other two directions. This shows that with the continuous accumulation of deformation, the accuracy of the multipath errors separated by EWT-ICA-R gradually decreases, resulting in a corresponding decrease in the accuracy of the extracted deformation information. For the same reason, the accuracy of the deformation information extracted in the slow deformation monitoring series is lower than that of the vibration deformation, and the decrease in accuracy increases with the continuous accumulation of deformation. To some extent, this phenomenon proves the validity of the method proposed in this paper, and it is similar to the process of processing other GNSS data. Limited by the length, we only chose the data with the longest time span to write this paper without adding other data.

As we all know, the main purpose of setting high-frequency deformation information is to simulate the deformation characteristics of structures or buildings and verify the effectiveness of the proposed method in such scenarios. Similarly, the purpose of setting the low-frequency deformation information is to simulate the deformation characteristics of the mining surface, urban surface, and landslide, and verify the effectiveness of the proposed method in such scenarios.

In addition, our research group has some GNSS and accelerometer data of bridge monitoring, and surface subsidence data from coal mining. For two reasons, this paper does not adopt such measured data: First, such data are highly professional, and the use of such data will add a certain length of background and relevant professional analysis, which will dilute the core of this paper to a certain extent. Second, the core of this paper is the proposed method and the biggest advantage of simulation deformation is that it has enough accurate verification data, which can more accurately verify the advantages of the new method. This convenience is not available in measured data.

5. Conclusions

Aiming at the separation of deformation information in GNSS deformation monitoring, the EWT-ICA-R algorithm was proposed to separate the vibration deformation and slow deformation from the coordinate series. Through theoretical analysis, data processing, and statistical analysis, the following conclusions were drawn.

For the vibration deformation with overlapping multipath frequencies, the difference between the extracted vibration deformation information and the reference value in the E direction on the first and last day is within 0.2–0.3 mm, which achieves a sub-millimeter level of deformation information extraction. The average correlation coefficients between the vibration deformation in the E, N, and U directions extracted for 12 consecutive days and the harmonic analog signal are 0.981, 0.981, and 0.885, respectively, and the average RMS values were 0.694 mm, 0.694 mm, and 1.852 mm, respectively. The deformation information extracted from the three directions maintains a high degree of similarity, and it is also believed that the method can accurately separate multipath errors from observations containing vibrational deformation.

For the slow deformation experiment, the differences between the extracted and simulated deformation sequences were kept within 5 mm, and the average correlation coefficients between the deformation information extracted for 12 consecutive days and the simulated deformation in the E, N, and U directions were 0.997, 0.996, and 0.983, respectively. In terms of extraction accuracy, the accurately separated deformation was different from the simulated data, and the RMS evaluation indicators in the E, N, and U directions were 1.345, 1.546, and 3.866 mm, respectively.

The results of this study indicate that the EWT-ICA-R method has excellent processing ability and wide application prospects for the deformation information separation problem of dynamic deformation monitoring, such as urban ground settlement deformation monitoring, dam deformation monitoring, bridge deformation monitoring, landslide monitoring, high-rise building deformation monitoring, mining deformation monitoring. The deformation information is obtained in real-time to guarantee the safety of human life and property. The method can separate multipath errors and extract vibration deformation and slow deformation information accurately in real-time. However, owing to the lack of DOY142 observation data in the middle, the accuracy of the vibration deformation extraction has certain fluctuations. Meanwhile, with the accumulation of deformation, the accuracy of multi-path error separation in the slow deformation experiment decreases, and the accuracy of the extracted slow deformation information decreases. Future work needs to solve the problem of deformation information separation under the condition of missing data and complex observation environments.