Location via proxy:   [ UP ]  
[Report a bug]   [Manage cookies]                
Next Article in Journal
Soft Array Surface-Changing Compound Eye
Previous Article in Journal
Acceptance and Preferences of Using Ambient Sensor-Based Lifelogging Technologies in Home Environments
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

An Improved Adaptive IVMD-WPT-Based Noise Reduction Algorithm on GPS Height Time Series

1
School of Surveying and Mapping Engineering, East China University of Technology, 418 Guanglan Road, Nanchang 330013, China
2
Key Laboratory for Digital Land and Resources of Jiangxi Province, East China University of Technology, Nanchang 330013, China
3
Physikalisch-Meteorologisches Observatorium Davos/World Radiation Center (PMOD/WRC), Dorfstrasse 33, CH-7260 Davos, Switzerland
4
Space and Earth Geodetic Analysis Laboratory (SEGAL), University of Beira Interior, 6201-001 Covilhã, Portugal
5
School of Civil and Surveying & Mapping Engineering, Jiangxi University of Science and Technology, No. 86 Hongqi Ave., Ganzhou 341000, China
*
Author to whom correspondence should be addressed.
Sensors 2021, 21(24), 8295; https://doi.org/10.3390/s21248295
Submission received: 20 October 2021 / Revised: 6 December 2021 / Accepted: 8 December 2021 / Published: 11 December 2021
(This article belongs to the Section Remote Sensors)

Abstract

:
To improve the reliability of Global Positioning System (GPS) signal extraction, the traditional variational mode decomposition (VMD) method cannot determine the number of intrinsic modal functions or the value of the penalty factor in the process of noise reduction, which leads to inadequate or over-decomposition in time series analysis and will cause problems. Therefore, in this paper, a new approach using improved variational mode decomposition and wavelet packet transform (IVMD-WPT) was proposed, which takes the energy entropy mutual information as the objective function and uses the grasshopper optimisation algorithm to optimise the objective function to adaptively determine the number of modal decompositions and the value of the penalty factor to verify the validity of the IVMD-WPT algorithm. We performed a test experiment with two groups of simulation time series and three indicators: root mean square error (RMSE), correlation coefficient (CC) and signal-to-noise ratio (SNR). These indicators were used to evaluate the noise reduction effect. The simulation results showed that IVMD-WPT was better than the traditional empirical mode decomposition and improved variational mode decomposition (IVMD) methods and that the RMSE decreased by 0.084 and 0.0715 mm; CC and SNR increased by 0.0005 and 0.0004 dB, and 862.28 and 6.17 dB, respectively. The simulation experiments verify the effectiveness of the proposed algorithm. Finally, we performed an analysis with 100 real GPS height time series from the Crustal Movement Observation Network of China (CMONOC). The results showed that the RMSE decreased by 11.4648 and 6.7322 mm, and CC and SNR increased by 0.1458 and 0.0588 dB, and 32.6773 and 26.3918 dB, respectively. In summary, the IVMD-WPT algorithm can adaptively determine the number of decomposition modal functions of VMD and the optimal combination of penalty factors; it helps to further extract effective information for noise and can perfectly retain useful information in the original time series.

1. Introduction

With the rapid development of space observation technology, GPS has become an important observational approach in geodesy and geodynamics [1,2]. Globally distributed International GNSS Service (IGS, see Abbreviations) reference stations have accumulated nearly twenty years of coordinate time series, which provide valuable basic data for the study of the geodynamics and global tectonics of Earth’s lithosphere and mantle [3,4,5,6,7]. Moreover, GPS observable time series not only contain geophysical signals but also unmodelled errors and other nuisance parameters, making the GPS coordinate time series present a nonlinear variation that affects the performance of the estimation of site coordinates and velocity [8,9]. The study and analysis of the Global Navigation Satellite Systems (GNSS) time series are conducive to obtaining accurate positions and velocities of stations, reasonably understanding plate tectonic movements, and establishing and maintaining dynamic earth reference frames, and they contribute to the study of relevant geodynamic processes. Therefore, in the study of GNSS signal processing, how to effectively reduce the influence of various noises in the original timing signals has always been a hot research issue in GNSS timing analysis.
In terms of GNSS time series noise reduction, Ghaderpour and Pagiatakis developed a new method of spectral analysis, namely, the least-squares wavelet analysis (LSWA), which decomposes a time series into the time–frequency domain, allowing for the detection of short-duration signatures in the series [10,11,12]. In [13], they proposed a new method, the Hilbert–Huang Transform (HHT), which compared to the wavelet and Fourier analyses, offers much better temporal and frequency resolutions. In [14], the authors applied the Kalman filter to GPS data noise reduction, and the experimental results show that the Kalman filter has a good application effect on noise reduction of triple-difference observation data, but the accuracy of the system equation directly affects the filtering effect [14]. Mosavi et al. proposed the wavelet packet transform, which can decompose the low-frequency part and better process the high-frequency part of the signal. The wavelet packet transform improves the time–frequency resolution of the signal, but it cannot improve distortion phenomena such as blurring of the signal edge [15,16,17,18]. Huang et al. improved the empirical mode decomposition (EMD) algorithm and applied it to GPS time series noise reduction. The noise reduction effect was effectively improved, but certain endpoint effects and modal mixing phenomena occurred, affecting the noise reduction effect [19,20,21]. Through the optimisation of the EMD method, ensemble empirical mode decomposition (EEMD) [22,23] and the complementary ensemble empirical mode decomposition (CEEMD) method [24,25,26,27,28] are obtained. Although EEMD and CEEMD can effectively suppress the modal aliasing phenomenon, the calculation is complicated and large. Zhang et al. [29] improved the EEMD noise identification method based on the continuous mean square error criterion and verified that the method could correctly identify the boundary point between the signal and noise.
With the rapid development of time–frequency analysis methods, Dragomiretskiy et al. proposed a new signal multiscale time–frequency analysis and processing method, variational mode decomposition (VMD) [30]. This method is based on VMD to denoise mechanical signals, and the denoising effect is better than wavelet and EMD denoising methods to varying degrees. In view of the advantages of VMD in analysing complex nonlinear, multiscale and nonstationary data, its algorithm has good antinomies performance, but the number of modal functions and penalty factors in the VMD method needs to be set in advance, and the use of inappropriate parameter combinations will result in insufficient noise reduction, so it is not adaptive [31,32,33,34,35].
Saremi et al. proposed the grasshopper optimisation algorithm (GOA) in 2017 and compared it with a variety of optimisation algorithms. The results show that GOA has outstanding advantages in the optimisation of unimodal functions, multimodal functions and composite functions [36]. Since GOA considers a given optimisation problem as a black box and does not need any gradient information of the search space, this makes it a highly suitable optimisation technique for any properly formulated optimisation problem in different fields [37]. The GOA algorithm is not affected by the nonlinear or magnitude of a problem, where usually other global optimisation techniques show early convergence, it finds the best solution more efficient with faster convergence [38], and VMD is a non-recursive approach that can adaptively derive an ensemble of band-limited intrinsic mode functions (BLIMFs) from non-stationary and nonlinear signals simultaneously [39].
To solve the above problems, we propose an improved variational modal decomposition (VMD) algorithm combining wavelet packet and energy entropy mutual information as the objective function and combined it with data experiments and analysis to verify the effectiveness and universality of the proposed method.
The rest of the paper is organised as follow. Section 2 describes the background theory of VMD and WPD. The model of the IVMD and the flowchart of IVMD-WPT Algorithm are discussed in Section 3. The validity of the proposed method was verified by simulation time series and GPS height time series from the CMONOC in Section 4. The conclusions of the experiment are summarised in Section 5.

2. Principles and Methods

2.1. Basic Principles of the VMD Method

VMD is a new signal decomposition method that decomposes the input signal of f into K modal components with centre frequency ω k and reconstructs the input raw signal. Therefore, the process of VMD can be regarded as the construction and solution of the constrained variational problem described in Equation (1).
min μ k , ω k k = 1 K t δ t + j π t × μ k ( t ) e j ω k t 2 2 s . t . k μ k = f
where μ k ( t ) is the intrinsic modal function, ω k is the centre frequency of the modal function, δ t is the impulse function and e j ω k t is the estimated centre frequency of each analytic signal.
Here, a quadratic penalty factor α and Lagrange operator λ ( t ) are used to render the variational problem unconstrained. α can ensure the accuracy of reconstruction in the presence of Gaussian noise, and λ ( t ) can ensure the tightness of constraint conditions. Therefore, the extended Lagrange expression is:
L μ k ( t ) , ω k , λ ( t ) = α k = 1 K t δ t + j π t × μ k ( t ) e j ω k t 2 2 + f ( t ) k = 1 K μ k ( t ) 2 2 + λ ( t ) , f ( t ) k = 1 K μ k ( t )
Equation (2) is solved by using the alternating direction method of multipliers (ADMM) and μ k n + 1 ( t ) , ω k n + 1 and λ n + 1 ( t ) are updated to find the optimal solution of the original variational problem in Equation (2). Equations (3) and (4) give the iterative formulae of each modal function μ k ( t ) and corresponding central frequency ω k , respectively.
μ k n + 1 ( ω ) = f ^ ( ω ) i k μ ^ i ( ω ) + λ ^ ( ω ) 2 1 + 2 α ( ω ω k ) 2
ω k n + 1 = 0 ω μ ^ k ( ω ) 2 d ω 0 μ ^ k ( ω ) 2 d ω
where f ^ ( ω ) , μ ^ k ( ω ) , μ ^ k n + 1 ( ω ) and λ ^ ( ω ) represent the Fourier transform of f ^ ( t ) , μ k ( t ) , μ k n + 1 ( t ) , λ ( t ) and n is the number of iterations. Equation (5) is the iterative formula of the Lagrange operator.
λ ^ n + 1 ( ω ) = λ ^ n ( ω ) + τ f ^ ( ω ) k = 1 K μ ^ k n + 1 ( ω )
where τ is the iteration step, and Equation (6) is the convergence condition, ε is the convergence tolerance.
k = 1 K μ ^ k n + 1 μ ^ k n 2 2 μ ^ k n 2 2 < ε

2.2. Grasshopper Optimisation Algorithm

The grasshopper optimisation algorithm (GOA) is a nature-inspired algorithm with high search efficiency and fast convergence. It simulates the predation behaviour of natural grasshopper swarms. The process of searching for food sources can be divided into two steps: exploration and development. In the exploration process, the long distance of the swarming is conducive to the global search, while in the development process, it is local. The behaviour is mathematically modelled as follows:
X i = S i + G i + A i
where X i is the location of the ith grasshopper; S i is the interaction factor between the ith grasshopper and other grasshoppers; G i is the force of gravity on the first grasshopper; A i expresses wind advection. The calculation formula of S i is as follows:
S i = j = 1 , j i N s ( X j X i ) d i j
where N is the population size; X j X i is the distance between the ith and jth grasshoppers; d i j = X j X i / X j X i represents a unit vector from the ith to the jth grasshopper; s is the social forces between grasshoppers as shown below:
s ( r ) = f e r l e r
where f is the attraction intensity between grasshoppers and l is the ratio of the attraction length.
G i and A i can be calculated by Equations (10) and (11).
G i = g e g
A i = u e ω
where g is the gravity constant; e g is the unit vector pointing to the centre of the earth; u is the wind force constant; e ω is the unit vector of wind direction, so the X i expansion is as follows:
X i = j = 1 , j i N s ( X j X i ) X j X i d i j g e g + u e ω
In solving practical problems, gravity is usually not taken into account, the wind direction is always set to the optimisation target of T d and the parameter coordination global and local search process is introduced. Then, the position formula is as follows:
X i d = c j = 1 , j i N c U d L d 2 s ( X j X i ) X j X i d i j + T d
where U d and L d are the upper and lower limits of the function in dimensional space; T d is the optimal solution of the current grasshopper position; c is the attenuation coefficient as shown below:
c = c max l c max c min M
where c max and c min are the maximum and minimum values; l is the current iteration number; M is the maximum number of iterations.

2.3. Principle of the Wavelet Packet Algorithm

Wavelet decomposition decomposes the original time series into high frequency and low frequency through a set of high-pass and low-pass filters and then decomposes the low frequency part. By wavelet packet decomposition, the high-frequency part not involved in wavelet decomposition is further decomposed, and then the optimal wavelet basis function is selected. The time–frequency analysis effect is better than that of the wavelet function. The specific steps are as follows:
Step 1: Define φ ( x ) and ψ ( x ) as the orthogonal scaling function and its corresponding wavelet function. Set h ( k ) as the low-pass filter coefficient and g ( k ) as the high-pass filter coefficient.
φ ( x ) = 2 k z h ( k ) φ ( 2 x k ) ψ ( x ) = 2 k z g ( k ) φ ( 2 x k )
Let μ 0 = φ ( x ) and μ 1 = ψ ( x ) then:
μ 2 n ( x ) = 2 k z h ( k ) μ n ( 2 x k ) μ 2 n + 1 ( x ) = 2 k z g ( k ) μ n ( 2 x k )
Step 2: Assume subspace U j n as the closure space of function μ n ( x ) and subspace U j 2 n as the closure space of function μ 2 n ( x ) and g j n U j u . C l j , n is the coefficient of φ ( x ) in the subspace, so g j n ( x ) can be expressed as:
C l j , n = + φ ( x ) 2 j / 2 μ n ( 2 j t l ) d t
g j n = l C l j , n μ n ( 2 j x 1 )
The wavelet packet decomposition algorithm can be obtained as:
C l j , n = k h ( k 2 l ) C k j + 1 , n C l j , 2 n + 1 = k g ( k 2 l ) C k j + 1 , n
Step 3: Decompose the wavelet packet, perform the inverse operation and obtain the wavelet packet reconstruction expression as follows:
C l j + 1 , n = k h ( l 2 k ) C l j , 2 n + g ( l 2 k ) C l j , 2 n + 1

3. IVMD-WPT Algorithm

As discussed in the previous section, the key to VMD performing feature extraction on time series data is the determination of the decomposition modal number K and penalty factor α . If K is less than the number of useful components in the processed signal, it will cause insufficient data decomposition; conversely, it will cause an over-decomposition phenomenon. An improper value of the penalty factor α may lead to centre frequency overlap of the modal function; therefore, appropriate parameters must be selected.

3.1. Improved VMD Method

Since GPS time series data are polluted by stationary and nonstationary noise and a single indicator cannot be used to obtain signal features, the mixing of two or more single indicators provides stronger robustness [40,41]. Thus, to determine VMD parameters, the energy entropy mutual information (EEMI) index composed of energy entropy and mutual information is adopted in this paper. The sum of the EEMI of the previous two modal functions is taken as the objective function, and the VMD parameters are optimised by the GOA algorithm as shown in Equation (21).
f i t n e s s = min γ = ( K , α ) i = 1 2 E E M I K 2 , 8 α [ 1000 , 10000 ]
where f i t n e s s is the objective function, γ = ( K , α ) is the value range of the decomposition modal number K and penalty factor α and E E M I = E E * M I , E E and M I are the energy entropy and mutual information, respectively, which can be calculated by:
E i = + i m f i ( t ) 2 d t i , i = 1 , 2 , , K σ i = E i / E E E = i = 1 N ( σ i ) I n ( σ i )
M I ( X ; Y ) = y Y x X p ( x , y ) log p ( x , y ) p ( x ) p ( y )
where i m f i ( t ) i = 1 , , K are the modes of different frequency bands, E i = E 1 , E 2 , , E K is the energy distribution of the vibration signals in the frequency domain, E = E 1 + + E K , p ( x ) and p ( y ) are the edge probability distribution functions of X and Y , respectively, and p ( x , y ) denotes the joint probability distribution function of X and Y .

3.2. IVMD-WPT

Aiming at the problem that the VMD method cannot determine the number of modal functions K and the penalty factor α and that the removed noise contains any valid information, this paper adopted the energy entropy mutual information (EEMI) index as the parameter adaptive to VMD and wavelet packet of the objective function. The specific steps of the combined GNSS coordinate time series noise reduction algorithm are as follows, and the flowchart is depicted in Figure 1.
Step 1: The parameter range of the VMD algorithm was set and the parameters of the GOA algorithm were initialised in [34,35] the modal component number K 2 , 8 , K 2 , 10 and penalty factor α 1000 , 10000 ; however, in [40,41], K takes an integer in the interval of 2 , 8 , and this paper focused on VMD applying to the GPS; thus, the authors considered that the range of K and α was 2 , 8 , 1000 , 10000 , and the population number of the GOA algorithm N = 30 , and the maximum cycle number L = 10 [34,35].
Step 2: The method described in Section 3.1 was used to select the optimal parameters of the VMD method, and the VMD method with the optimal parameters was used to decompose the GNSS coordinate time series.
Step 3: The composite evaluation index T [42] was used to form the reconstructed time series by summing each modal component successively, and the composite evaluation index T value of each reconstructed time series was calculated. When the T value was the smallest, the corresponding reconstructed time series was a denoising time series, and the remaining IMF components were regarded as high-frequency noise.
Step 4: The decomposed noise was further denoised by wavelet packet transform (WPT). Finally, the effective signal component after denoising by wavelet packet transform and the denoising time series obtained in Step 3 were reconstructed into the final denoising time series.
To effectively evaluate the effectiveness of the text combination method, the signal-to-noise ratio, root mean square error (RMSE) and correlation coefficient (R) were selected as the evaluation indices for noise reduction, and the calculation formulae are shown as follows:
R s n = 10 × log n = 0 N 1 S n 2 n = 0 N 1 S n S ¯ n 2
R M S E = 1 n i = 1 n x i - x i ^ 2
r i = C o v x i , y i V a r X V a r Y , i = 1 , 2 , , n
where n and N denote the number of sampling points in the sequence; x i , S n and X denote the sequence of the signal; x i ^ , S ¯ n and Y denote the original sequence.

4. Experiment Analysis and Discussion

In this section, two groups of simulation data were used to verify the effectiveness of the IVMD-WPT method for denoising common periodic signals and time series signals, and GNSS measured elevation coordinate time series data were used to verify the universality of the IVMD-WPT denoising effect.

4.1. Simulation Experiment A

Simulation GNSS elevation time series 1 (Sim_1) consists of a trend term, seasonal term (periodic term) and noise term. First, Sim_1 TS containing three constant amplitude period terms and Gaussian white noise was generated by Equation (27) in which the sampling frequency was 1 Hz, the sampling number was 1024 and the signal-to-noise ratio was 6 dB. Figure 2 and Figure 3 are the simulated original time series diagram and simulated component waveform diagram, respectively.
y 1 = 5 sin 2 π t i 600 sin 2 π t i 350 y 2 = 7 sin 2 π t i 500 y 3 = 2 sin 2 π t i 50 ε = n o i s e f = y 1 + y 2 + y 3 + ε
The method in this paper was used to decompose simulated data A to select the optimal decomposition modal number K and penalty factor α of the VMD. The GOA parameters were set as follows: search agent n = 30 and maximum number of loops L = 10 . Figure 4 and Figure 5 are the historical values of K and α and the convergent variation diagram of the objective function. Figure 4 and Figure 5 show that the optimal parameter combination of VMD is K = 6 and α = 1009 . We prove the EEMI effectiveness by analysing the multiple indexes of the modes, including MI, EE and EEMI. As shown in Figure 6, the change in indexes MI was very small for the difference IMFs, while the indexes EE and EEMI obtained by GOA were larger than MI. It was indicated that EE and EEMI were sensitive to the two parameters of VMD, and the EEMI index was effective.
The VMD method with the optimal parameter combination was used to decompose simulated data A. Figure 6 shows the decomposition result graph and the corresponding spectrum graph. Analysis of the spectrum diagram in Figure 7 shows that the IVMD decomposed the simulated data I into six IMFs, among which the IMF1–IMF3 modal component extracted the main information in the simulated data, while IMF4–IMF6 contained noise and some useful information. Then, to better distinguish high-frequency noise from low-frequency useful signals, the composite evaluation index T was adopted, and each modal component (IMF1–IMF6) was successively accumulated to form reconstructed time series, and the composite evaluation index T value of each reconstructed time series was calculated. The corresponding reconstructed time series was a denoising signal when the T value was minimal, and the remaining IMF components were regarded as high-frequency noise. Table 1 shows the corresponding T value of each reconstructed time series. When the reconstructed time series is i = 1 2 I M F i , T = 0.1697 reached the minimum, so i = 1 2 I M F i was regarded as a denoising time series and I M F 3 I M F 6 was a high-frequency noise.
As the high-frequency noise I M F 3 I M F 6 also contained part of the original time series information, wavelet packet analysis was carried out for further noise reduction. The parameters of wavelet packet: thresholding rule = 4.2975, thresholding function was hard, decomposition level number = 4, wavelet function = db3. The autocorrelation coefficient threshold method combined with EMD and IVMD was used to denoise the analogue time series, and a comparative analysis was performed with the method in this paper. The effect of denoising was evaluated using three indicators: root mean square error, correlation coefficient, and signal-to-noise ratio. The statistical table of the analysis results is shown in Table 2.
Table 2 shows that this method is superior to the EMD and IVMD methods on the whole, and the root mean square error of the simulated noise reduction results was 0.1563 mm. Compared with the EMD and IVMD methods, the root mean square error decreased by 0.084 and 0.0715, respectively, and the correlation coefficient was 0.9997. Compared with EMD and IVMD, the improvement was 0.0005 and 0.0004, respectively, thus verifying the effectiveness of the method.

4.2. Simulation Experiment B

The method in this paper was used to conduct noise reduction analysis on elevation data of GNSS coordinate time series. A simulated time series containing a trend term, periodic term and noise term was constructed according to the GNSS single station and single component coordinate time series function model and random model, and the function expression is as follows:
y ( t i ) = a + b t i T + c sin ( 2 π t i T ) + d cos ( 2 π t i T ) + e sin ( 4 π t i T ) + f cos ( 4 π t i T ) + v i
where t i is the observation time in units of year; a is the starting position of the time series of the station; b is the linear speed of the station movement; c , d , e , f are the amplitudes of the annual and semi-annual movements of the station, respectively; v i is the noise in the GPS coordinate time series, which is more realistic and simulates GPS time series noise. In this paper, Equation (28) was used to generate white noise and coloured noise to realistically simulate GPS time series noise.
v ( t i ) = w ( t i ) , i 2 v ( t i ) = f 1 v ( t i 1 ) + f 2 v ( t i 2 ) + w ( t i ) , i 3
where w ( t i ) is white noise with zero mean and variance of 1; f 1 and f 2 are 0.2 and 0.2, respectively. According to Equations (28) and (29), this paper generated simulated time series data II with a length of 1826. Table 3 shows the parameters of the simulated data, and Figure 8 shows the simulated time series.
According to the steps in Section 2.1, the EMD, IVMD, and IVMD-WPT methods were used to analyse the noise reduction of analogue data II. Figure 9 and Figure 10 show the parameter changes determined by the IVMD method, and Table 4 shows the corresponding T values of reconstructed time series. The analysis of Figure 8 and Figure 9 shows that the optimal parameter combination of VMD of the simulated data II was K = 7 and α = 1000 . Combined with Table 4, the T index reached the minimum when the three IMF modal components were accumulated, so i = 1 3 I M F i will be used as the signal after denoising. Then, the root mean square error, correlation coefficient, and signal-to-noise ratio (SNR) after noise reduction using the EMD, IVMD and IVMD-WPT methods were calculated. Table 5 shows the statistical results of the three indicators.
By analysing the results in Table 5, it can be seen that the root mean square error decreased by 0.0087, and the correlation coefficient and the signal-to-noise ratio increased by 0.0007 and 0.7391, respectively. Therefore, for the GNSS coordinate time series data, the improved VMD method in this paper also had a better noise reduction effect than the EMD method. In the comparison between the IVMD-WPT method and IVMD method, the root mean square error decreased by 0.0003 and the signal-to-noise ratio increased by 0.0239, indicating that the IMF modal component that was removed also contained information. Therefore, the IVMD-WPT method proposed in this paper can extract effective information in the original time series more effectively.

4.3. Noise Reduction Analysis with Real GNSS Elevation Time Series

To further verify the reliability and applicability of the proposed method, this paper adopted the time series of the original elevation coordinates of 100 land state network reference stations in China to conduct noise reduction research (data came from the GNSS Data Product and Service Platform of the China Earthquake Administration). The observation epoch was a total of 5 years of elevation coordinate time series signals from 1 January 2010 to 1 January 2015, with a sampling interval of 1/365.25 a and a sampling frequency of 365.25 Hz. EMD, IVMD and the method in this paper were used to denoise the elevation data of 100 reference stations, and the relevant parameter settings were consistent with the simulation experiment. The BJFS station was taken as an example for detailed illustration. Figure 11 is the denoising effect diagram of the three methods for the BJFS station, and Table 6 is the statistical table of the denoising evaluation parameters of some stations.
As seen from Figure 10, compared with the EMD method, the IVMD-WPT method and IVMD method in this paper can better avoid the influence caused by the endpoint effect and have a better noise reduction effect at both ends of the original time series. The time series denoised by the IVMD-WPT method in this paper can better fit the original time series. It can effectively reflect local trend motion changes and small periodic oscillations. Table 3 summarises the root mean square error, signal-to-noise ratio and correlation coefficient of the results of the three algorithms. It can be seen from the table that the IVMD-WPT method in this paper is superior to the other two single methods on the whole. Compared with the other two methods, the IVMD-WPT method increased the root mean square error by 11.4648 and 6.7322 mm on average, respectively. The SNR index increased by 32.6773 and 26.3918 dB, and the correlation coefficient increased by 0.1458 and 0.0588, respectively, indicating that the noise reduction effect of the IVMD-WPT method presented in this paper was better and that the noise reduction results were more reliable.

5. Conclusions

Since the traditional VMD method cannot be sure of the number of decomposition mode functions (IMFs) and the value of the punishment factor in the process of noise reduction, resulting in inadequate decomposition or making the decomposition part a problem of denoising the time series effectively, a kind of energy entropy mutual information was proposed as the objective function to improve variational mode decomposition (VMD) combined with the wavelet packet denoising algorithm. The method was based on the traditional method of VMD, energy entropy mutual information as the objective function and the locusts optimised algorithm (GOA) to optimise the objective function; thus, the mode decomposition number and value of the punishment factor were determined adaptively, and the composite index T was used to determine the noise. Subsequently, the noise component was filtered using wavelet packet transform, and the filtered time series and denoising time series were reconstructed to obtain the final denoising time series. In this paper, two groups of analogue time series and 100 groups of measured time series of the land state network were used for research and analysis.
The main contributions of this approach are summarised below.
  • Compared with the traditional VMD method, this paper used the energy entropy mutual information as the objective function and used GOA to optimise the objective function to adaptively determine the number of decomposition mode functions (IMFs) and the value of the punishment factor and to improve the effect of noise reduction.
  • Compared with the single EMD method, the IVMD method can effectively weaken the influence of the endpoint effect, thus improving the noise reduction effect. The simulation results showed that the RMSE decreased by 0.0106 mm and the CC and SNR increased by 0.0004, and 428.42 dB, respectively.
  • Compared with the two single models of traditional EMD and IVMD proposed in this paper, the IVMD-WPT method was superior to the two single models in the three indicators of root mean square error, correlation coefficient and signal-to-noise ratio. The real results showed that the RMSE decreased by 11.4648 and 6.7322 mm and CC and SNR increased by 0.1458 and 0.0588 and 32.6773 and 26.3918 dB, respectively, thus verifying the effectiveness of the IVMD-WPT method in noise reduction. In addition, the local optimum problem of GOA (the appropriate parameters of WPT) need further exploration.

Author Contributions

H.X., T.L. and X.H. writing–original draft preparation and methodology; J.-P.M. reviewed and edited the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding

This work was sponsored by Key Laboratory for Digital Land and Resources of Jiangxi Province (DLLJ201909), National Natural Science Foundation of China (42061077, 42104023, 42064001), Jiangxi University of Science and Technology High-Level Talent Research Startup Project (205200100564) and Jiangxi Provincial Natural Science Foundation (20212BAB204030).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

GPSGlobal Positioning System
VMDVariational Mode Decomposition
IMFsIntrinsic Modal Functions
IVMD-WPTImproved Variational Mode Decomposition and Wavelet Packet Transform
EEMIEnergy Entropy Mutual Information
GOAGrasshopper Optimisation Algorithm
RMSERoot Mean Square Error
CCCorrelation Coefficient
SNRSignal-to-Noise Ratio
EMDEmpirical Mode Decomposition
IVMDImproved Variational Mode Decomposition
CMONOCCrustal Movement Observation Network of China
IGSInternational GNSS Service
GNSSGlobal Navigation Satellite System
LSWALeast-Squares Wavelet Analysis
HHTHilbert–Huang Transform
EEMDEnsemble Empirical Mode Decomposition
CEEMDComplementary Ensemble Empirical Mode Decomposition

References

  1. Bock, Y.; Melgar, D. Physical applications of GPS geodesy: A review. Rep. Prog. Phys. 2016, 79, 106801. [Google Scholar] [CrossRef] [PubMed]
  2. He, X.; Montillet, J.P.; Fernandes, R.; Bos, M.; Yu, K.; Hua, X.; Jiang, W. Review of current GPS methodologies for producing accurate time series and their error sources. J. Geodyn. 2017, 106, 12–29. [Google Scholar] [CrossRef]
  3. Serpelloni, E.; Faccenna, C.; Spada, G.; Dong, D.; Williams, S.D. Vertical GPS ground motion rates in the Euro-Mediterranean region: New evidence of velocity gradients at different spatial scales along the Nubia-Eurasia plate boundary. J. Geophys. Res. Solid Earth 2013, 118, 6003–6024. [Google Scholar] [CrossRef] [Green Version]
  4. Van Dam, T.M.; Blewitt, G.; Heflin, M.B. Atmospheric pressure loading effects on the global positioning system coordinate determinations. Geophys. Res. 1994, 99, 23939–23950. [Google Scholar] [CrossRef]
  5. Tregoning, P.; Watson, C. Atmospheric effects and spurious signals in GPS analyses. Geophys. Res. 2009, 114, B09403. [Google Scholar] [CrossRef] [Green Version]
  6. Soundy, A.W.; Panckhurst, B.J.; Brown, P.; Martin, A.; Molteno, T.C.; Schumayer, D. Comparison of Enhanced Noise Model Performance Based on Analysis of Civilian GPS Data. Sensors 2020, 20, 6050. [Google Scholar] [CrossRef] [PubMed]
  7. Shen, N.; Chen, L.; Liu, J.; Wang, L.; Tao, T.; Wu, D.; Chen, R. A Review of Global Navigation Satellite System (GNSS)-based Dynamic Monitoring Technologies for Structural Health Monitoring. Remote Sens. 2019, 11, 1001. [Google Scholar] [CrossRef] [Green Version]
  8. Han, M.; Liu, Y.; Xi, J.; Guo, W. Noise smoothing for nonlinear time series using wavelet soft threshold. IEEE Signal Process Let. 2007, 14, 62–65. [Google Scholar] [CrossRef]
  9. Altamimi, Z.; Rebischung, P.; Métivier, L.; Collilieux, X. ITRF2014: A new release of the International Terrestrial Reference Frame modelling nonlinear station motions. J. Geophys. Res. Solid Earth 2016, 121, 6109–6131. [Google Scholar] [CrossRef] [Green Version]
  10. Ghaderpour, E.; Pagiatakis, S.D. LSWAVE: A MATLAB software for the least-squares wavelet and cross-wavelet analyses. GPS Solut. 2019, 23, 1–8. [Google Scholar] [CrossRef]
  11. Kaczmarek, A.; Kontny, B. Identification of the noise model in the time series of GNSS stations coordinates using wavelet analysis. Remote Sens. 2018, 10, 1611. [Google Scholar] [CrossRef] [Green Version]
  12. Ghaderpour, E.; Pagiatakis, S.D. Least-squares wavelet analysis of unequally spaced and non-stationary time series and its applications. Math Geosci. 2017, 49, 819–844. [Google Scholar] [CrossRef]
  13. Huang, N.E.; Wu, M.L.; Qu, W.; Long, S.R.; Shen, S.S. Applications of Hilbert–Huang transform to non-stationary financial time series analysis. Appl. Stoch. Models Bus. 2003, 19, 245–268. [Google Scholar] [CrossRef]
  14. Baohong, F.; Yiqiang, G. Application of Kalman filter in GPS data preprocessing. Sci. Surv. Mapp. 2016, 41, 171–174. [Google Scholar]
  15. Mosavi, M.R.; Rezaei, M.J.; Pashaian, M.; Moghaddasi, M.S. A fast and accurate anti-jamming system based on wavelet packet transform for GPS receivers. Gps. Solut. 2017, 21, 415–426. [Google Scholar] [CrossRef]
  16. Sun, Z.; Chang, C.C. Structural damage assessment based on wavelet packet transform. J. Struct. Eng. 2002, 128, 1354–1361. [Google Scholar] [CrossRef]
  17. Cody, M.A. The wavelet packet transform: Extending the wavelet transform. DR Dobbs. J. 1994, 19, 44–46. [Google Scholar]
  18. Barros, J.; Diego, R.I. Analysis of harmonics in power systems using the wavelet-packet transform. IEEE T Instrum. Meas. 2007, 57, 63–69. [Google Scholar] [CrossRef]
  19. Huang, N.E.; Shen, Z.; Long, S.R.; Wu, M.C.; Shih, H.H.; Zheng, Q.; Yen, N.C.; Tung, C.C.; Liu, H.H. The empirical mode decomposition and the Hilbert spectrum for nonlinear and non-stationary time series analysis. Proc. R. Soc. Math. Phys. Eng. Sci. 1971, 454, 903–995. [Google Scholar] [CrossRef]
  20. Wang, T.; Zhang, M.; Yu, Q.; Zhang, H. Comparing the application of EMD and EEMD on time-frequency analysis of seimic signal. J. Appl. Geophys. 2012, 83, 29–34. [Google Scholar] [CrossRef]
  21. Liu, C.; Zhu, L.; Ni, C. The chatter identification in end milling based on combining EMD and WPD. Int. J. Adv. Manuf. Technol. 2017, 91, 3339–3348. [Google Scholar] [CrossRef]
  22. Zheng, J.D.; Cheng, J.S.; Yang, Y. Modified EEMD algorithm and its applications. J. Vib. 2013, 32, 21–26. [Google Scholar]
  23. Shen, W.; Peng, C. Detection of different-time-scale signals in the length of day variation based on EEMD analysis technique. Geod. Geodyn. 2016, 7, 180–186. [Google Scholar] [CrossRef] [Green Version]
  24. Wang, J.; He, X.; Ferreira, V.G. Ocean wave separation using CEEMD-Wavelet in GPS wave measurement. Sensors 2015, 15, 19416–19428. [Google Scholar] [CrossRef] [PubMed]
  25. Lu, J.; Chen, X.; Feng, S. A GPS Time Series Prediction Model Based on CEEMD. J. Adv. Comput. Netw. 2016, 4. [Google Scholar] [CrossRef] [Green Version]
  26. Zhang, J.; Tang, H.; Wen, T.; Ma, J.; Tan, Q.; Xia, D.; Zhang, Y. A hybrid landslide displacement prediction method based on CEEMD and DTW-ACO-SVR—Cases studied in the three gorges reservoir area. Sensors 2020, 20, 4287. [Google Scholar] [CrossRef] [PubMed]
  27. Zhang, J.; Tang, H.; Tannant, D.D.; Lin, C.; Xia, D.; Liu, X.; Ma, J. Combined forecasting model with CEEMD-LCSS reconstruction and the ABC-SVR method for landslide displacement prediction. J. Clean. Prod. 2021, 293, 126205. [Google Scholar] [CrossRef]
  28. He, X.; Yu, K.; Montillet, J.P.; Xiong, C.; Lu, T.; Zhou, S.; Ma, X.; Cui, H.; Ming, F. GNSS-TS-NRS: An Open-source MATLAB-Based GNSS time series noise reduction software. Remote Sens. 2020, 12, 3532. [Google Scholar] [CrossRef]
  29. Zhang, S.; Liu, H.; Hu, M.; Jiang, A.; Zhang, L.; Xu, F.; Hao, G. An adaptive CEEMDAN thresholding denoising method optimized by nonlocal means algorithm. IEEE Trans. Instrum. Meas. 2020, 69, 6891–6903. [Google Scholar] [CrossRef]
  30. Dragomiretskiy, K.; Zosso, D. Variational Mode Decomposition. IEEE Trans. Signal Process. 2014, 62, 531–544. [Google Scholar] [CrossRef]
  31. Shen, Y.; Zheng, W.; Yin, W.; Xu, A.; Zhu, H. Feature Extraction Algorithm Using a Correlation Coefficient Combined with the VMD and Its Application to the GPS and GRACE. IEEE Access 2021, 9, 17507–17519. [Google Scholar] [CrossRef]
  32. Zhang, X.; Miao, Q.; Zhang, H.; Wang, L. A parameter-adaptive VMD method based on grasshopper optimization algorithm to analyze vibration signals from rotating machinery. Mech. Syst. Signal PR 2018, 108, 58–72. [Google Scholar] [CrossRef]
  33. Ahmed, W.A.; Wu, F.; Marlia, D.; Zhao, Y. Mitigation of Ionospheric Scintillation Effects on GNSS Signals with VMD-MFDFA. Remote Sens. 2019, 11, 2867. [Google Scholar] [CrossRef] [Green Version]
  34. Zhou, C.; Ma, J.; Wu, J.; Yuan, X. An adaptive VMD method based on improved GOA to extract early fault feature of rolling bearings. Int. J. Innov. Comput. Inf. Control 2019, 15, 1485–1505. [Google Scholar]
  35. Li, C.; Liu, Y.; Liao, Y. An Improved Parameter-Adaptive Variational Mode Decomposition Method and Its Application in Fault Diagnosis of Rolling Bearings. Shock Vib. 2021, 2021, 2968488. [Google Scholar] [CrossRef]
  36. Saremi, S.; Mirjalili, S.; Lewis, A. Grasshopper optimisation algorithm: Theory and application. Adv. Eng. Softw. 2017, 105, 30–47. [Google Scholar] [CrossRef] [Green Version]
  37. Meng, H.; Han, Y.; Xun, Y.; Chen, J.; Du, N.; Cao, Y.; Zheng, Y. GPS/INS Integrated Navigation Based on Grasshopper Optimization Algorithm. IFAC-Papers Online 2019, 52, 29–34. [Google Scholar] [CrossRef]
  38. Mafarja, M.; Aljarah, I.; Heidari, A.A.; Hammouri, A.I.; Faris, H.; Ala’M, A.Z.; Mirjalili, S. Evolutionary population dynamics and grasshopper optimization approaches for feature selection problems. Knowl. Based Syst. 2018, 145, 25–45. [Google Scholar] [CrossRef] [Green Version]
  39. Jiang, X.; Shen, C.; Shi, J.; Zhu, Z. Initial center frequency-guided VMD for fault diagnosis of rotating machines. J. Sound Vib. 2018, 435, 36–55. [Google Scholar] [CrossRef]
  40. He, X.; Bos, M.S.; Montillet, J.P.; Fernandes, R.; Melbourne, T.; Jiang, W.; Li, W. Spatial Variations of Stochastic Noise Properties in GPS Time Series. Remote Sens. 2021, 13, 4534. [Google Scholar] [CrossRef]
  41. Silva, L.K.J.; Ramarakula, M. An efficient interference mitigation approach for NavIC receivers using improved variational mode decomposition and wavelet packet decomposition. Trans. Emerg. Telecommun. Technol. 2021, 32, e4242. [Google Scholar]
  42. Zhu, J.; Zhang, Z.; Kuang, C. A Reliable Evaluation Indicator of Wavelet Denoising. Geomat. Inf. Sci. Wuhan Univ. 2015, 40, 688–694. [Google Scholar]
Figure 1. IVMD-WPT algorithm flowchart.
Figure 1. IVMD-WPT algorithm flowchart.
Sensors 21 08295 g001
Figure 2. Simulation of the original time series.
Figure 2. Simulation of the original time series.
Sensors 21 08295 g002
Figure 3. Waveform diagram of each component of the analogue signal: (a) waveform diagram of y1, (b) waveform diagram of y2, (c) waveform diagram of y3, (d) waveform diagram of noise.
Figure 3. Waveform diagram of each component of the analogue signal: (a) waveform diagram of y1, (b) waveform diagram of y2, (c) waveform diagram of y3, (d) waveform diagram of noise.
Sensors 21 08295 g003
Figure 4. Historical values.
Figure 4. Historical values.
Sensors 21 08295 g004
Figure 5. Convergence of the objective function of data I.
Figure 5. Convergence of the objective function of data I.
Sensors 21 08295 g005
Figure 6. The three indexes for different IMFs.
Figure 6. The three indexes for different IMFs.
Sensors 21 08295 g006
Figure 7. VMD decomposition and spectrum diagram.
Figure 7. VMD decomposition and spectrum diagram.
Sensors 21 08295 g007
Figure 8. Simulated time series data.
Figure 8. Simulated time series data.
Sensors 21 08295 g008
Figure 9. Historical values of simulated data II.
Figure 9. Historical values of simulated data II.
Sensors 21 08295 g009
Figure 10. Convergence of the objective function of data Ⅱ.
Figure 10. Convergence of the objective function of data Ⅱ.
Sensors 21 08295 g010
Figure 11. Noise reduction effect of the three methods for the BJFS station.
Figure 11. Noise reduction effect of the three methods for the BJFS station.
Sensors 21 08295 g011
Table 1. T indicators corresponding to different IMF s of data I.
Table 1. T indicators corresponding to different IMF s of data I.
Reconstructed Time Series
Index i = 1 1 I M F i i = 1 2 I M F i i = 1 3 I M F i i = 1 4 I M F i i = 1 5 I M F i i = 1 6 I M F i
T 0.60150.16790.17310.20890.29610.3985
Table 2. Statistical table of the evaluation parameters for different noise reduction methods in data Ⅰ.
Table 2. Statistical table of the evaluation parameters for different noise reduction methods in data Ⅰ.
MethodsRMSE/mmCorrelation Coefficient (R)Signal-to-Noise Ratio
(Rsn/dB)
EMD0.24030.9992637.95
IVMD0.22780.99931494.06
IVMD-WPT0.15630.99971500.23
Table 3. Simulation data statistics.
Table 3. Simulation data statistics.
TimeIntercept (a)/mmLinear Velocity (b)/(mm/a)Annual Amplitude (c)/mmAnnual Amplitude (d)/mmHalf-Year Amplitude (e)/mmHalf-Year Amplitude (f)/mmPeriod (T)
2013–20171211.21.20.8200
Table 4. T indicators corresponding to different IMF s of data II.
Table 4. T indicators corresponding to different IMF s of data II.
Index Reconstructed Time Series
i = 1 1 I M F i i = 1 2 I M F i i = 1 3 I M F i i = 1 4 I M F i i = 1 5 I M F i i = 1 6 I M F i i = 1 7 I M F i
T 0.36160.31050.29310.31500.37730.49500.6384
Table 5. Statistical table of the evaluation parameters for different noise reduction methods in data II.
Table 5. Statistical table of the evaluation parameters for different noise reduction methods in data II.
MethodsRMSE/mmCorrelation
Coefficient (R)
Signal-to-Noise Ratio
(Rsn/dB)
EMD0.65460.978523.5371
IVMD0.64590.979224.2762
IVMD-WPT0.64560.979224.3001
Table 6. Statistical table of the noise reduction evaluation parameters.
Table 6. Statistical table of the noise reduction evaluation parameters.
SiteMethodsRMSE/mmCorrelation Coefficient (R)Signal-to-Noise Ratio
(Rsn/dB)
ARTUEMD4.87495.33180.9187
IVMD3.598410.34010.9568
IVMD-WPT2.584620.86240.9781
BJFSEMD8.809425.40940.9808
IVMD6.920441.31460.9882
IVMD-WPT3.2106194.66150.9975
CHANEMD5.877910.88990.9550
IVMD3.924023.99650.9800
IVMD-WPT2.538958.29240.9917
CHUNEMD5.34582.02000.8184
IVMD3.63525.04990.9215
IVMD-WPT2.95897.91850.9496
DLHAEMD142.72000.32290.0596
IVMD101.34570.26010.5831
IVMD-WPT36.82916.44180.9772
HRBNEMD5.792123.46290.9792
IVMD3.968450.09570.9903
IVMD-WPT2.8025101.14860.9952
KMINEMD7.52111.65700.7910
IVMD5.35863.77070.9017
IVMD-WPT3.183911.94090.9687
LUZHEMD4.08904.83020.9092
IVMD3.84655.26190.9203
IVMD-WPT2.476013.65320.9684
PIMOEMD4.528829.04760.9831
IVMD4.589627.67310.9827
IVMD-WPT4.118534.46990.9861
TAINEMD5.26845.00860.9108
IVMD3.86869.31590.9533
IVMD-WPT2.720119.52940.9776
WUSHEMD5.181510.86380.9562
IVMD4.315315.27260.9699
IVMD-WPT2.978832.76060.9860
XIAGEMD6.61841.38740.7279
IVMD4.46563.30620.8854
IVMD-WPT2.647510.67940.9629
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Xu, H.; Lu, T.; Montillet, J.-P.; He, X. An Improved Adaptive IVMD-WPT-Based Noise Reduction Algorithm on GPS Height Time Series. Sensors 2021, 21, 8295. https://doi.org/10.3390/s21248295

AMA Style

Xu H, Lu T, Montillet J-P, He X. An Improved Adaptive IVMD-WPT-Based Noise Reduction Algorithm on GPS Height Time Series. Sensors. 2021; 21(24):8295. https://doi.org/10.3390/s21248295

Chicago/Turabian Style

Xu, Huaqing, Tieding Lu, Jean-Philippe Montillet, and Xiaoxing He. 2021. "An Improved Adaptive IVMD-WPT-Based Noise Reduction Algorithm on GPS Height Time Series" Sensors 21, no. 24: 8295. https://doi.org/10.3390/s21248295

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop