COA-VMPE-WD: A Novel Dual-Denoising Method for GNSS Time Series Based on Permutation Entropy Constraint

1. Introduction

Due to factors such as the external environment of the monitoring station, unmodeled errors, and geophysical effects, Global Navigation Satellite System (GNSS) coordinate time series exhibit significant nonlinear variations, containing various signals [1] and colored noise [2,3]. The noise components in GNSS coordinate time series are complex, affecting the estimation of station velocity and its uncertainty, and even leading to misinterpretations of certain geophysical phenomena [4,5]. Therefore, in the analysis of GNSS coordinate time series, effectively reducing the impact of noise to obtain accurate station velocity and its uncertainty is of great significance for establishing high-precision velocity field models [6,7,8,9] and analyzing geophysical phenomena such as plate tectonic movements [10,11].Studies have shown that common methods for denoising GNSS coordinate time series [3,4,12,13,14] include wavelet analysis [15,16], singular spectrum analysis (SSA) [17,18], empirical mode decomposition (EMD), ensemble empirical mode decomposition(EEMD) and complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) et al. [19,20]. Kunpu and Yunzhong S. (2020) proposed a non-interpolated wavelet analysis algorithm that effectively extracts seasonal signals from GNSS coordinate time series [21]. However, in wavelet analysis, the selection of basis functions and decomposition levels lacks adaptability, and different basis functions and decomposition levels can significantly impact the denoising performance of the algorithm. Dong et al. (2002) employed wavelet decomposition (WD) and SSA to model the nonlinear variations in GNSS coordinate time series, reducing the likelihood of useful signals being mistakenly filtered out as noise [22]. Nevertheless, SSA analysis involves subjectivity in selecting the lag window, and different window lengths greatly affect the signal extraction results. Liu et al. (2023) used EMD to correct the periodic terms in GNSS continuous station time series, obtaining relatively reliable station velocities [23]. However, EMD is prone to endpoint effects and mode mixing. To mitigate this phenomenon, Yeh et al. (2010) proposed an optimized version of EMD called ensemble empirical mode decomposition (EEMD) [24], and Agnieszka and Dawid (2022) supplemented the ensemble EMD algorithm. However, both algorithms suffer from low computational efficiency [25], and their decomposition performance heavily depends on the number of ensemble trials and the amplitude of the added white noise.

Dragomiretskiy and Zosso (2013) proposed a non-recursive signal decomposition method called variational mode decomposition (VMD) [26]. Unlike EMD, which iteratively sifts and peels off signals [27], VMD essentially constructs and solves a variational problem, effectively avoiding the mode mixing and endpoint effects inherent in EMD. VMD possesses adaptive characteristics and circumvents the issues associated with EMD while also exhibiting advantages in noise robustness. As a result, VMD has been widely applied in signal denoising research [28,29,30]. However, VMD requires the predefinition of two critical parameters that influence its decomposition performance: the number of modes ($K$) and the penalty factor ($\alpha$). In practical applications, $K$ and $\alpha$ are often selected empirically, and inappropriate parameter choices can lead to over-decomposition or under-decomposition of the signal, adversely affecting denoising results [31]. Lu and Xie (2021) introduced a comprehensive index T to determine the optimal $K$value for VMD denoising [32]. This method fixes the $\alpha$ value and only considers the impact of $K$ on the denoising performance of the VMD algorithm, neglecting the interplay between $K$ and $\alpha$, and thus can only yield a relatively suboptimal parameter combination.

Based on the above, in this study, we propose a novel denoising method (COA-VMPE-WD) that combines Variational Mode Decomposition (VMD) optimized by the Crayfish Optimization Algorithm (COA) using envelope entropy as the fitness function, with multiscale permutation entropy (MPE) as a constraint condition (COA-VMPE), and wavelet decomposition (WD). This method is applied to denoise GNSS coordinate time series, and its effectiveness and reliability are validated through experiments on simulated signals and real GNSS observation data. The root means square error (RMSE), signal-to-noise ratio (SNR), mean absolute error (MAE), and cross-correlation coefficient (R) as the denoising evaluation indicators.

The organization of the work is as follows: Section 2 describes the GNSS data and mathematical methods. Section 3 shows the parameters of COA optimized VMD, optimal parameters selection, results of different denoising methods between simulated signals and real GNSS. Section 4 discusses the COA-VMPE-WD method and existing noise reduction methods for detecting the optimal noise model and analyzing the impact of GNSS station annual term and station velocity. Section 5 concludes this work.

2. Data and Methods

2.1. Data

In this study, to further verify the reliability and applicability of the COA-VMPE-WD method, we selected the original coordinate time series of 52 GNSS reference station networks in China for research and analysis. The observation time was GNSS coordinate time series data from 2010 to 2023, sourced from Crustal Movement Observation Network of China (CMONOC).The data we have selected meets the following conditions: 1) Data missing less than 10%; 2) The missing data is filled in by using the single spectrum analysis (SSA) method to fill in the differences; 3) We used the 3 interquartile range (3IQR) method for coarse error removal. The distribution of selected sites in the study is shown in Figure 1.

2.2. Methods

2.2.1. Crayfish Optimization Algorithm

Crayfish Optimization Algorithm (COA) is a metaheuristic algorithm inspired by the foraging and enemy avoidance behavior of crayfish COA has the ability to efficiently balance global exploration and local development [33]. It simulates the temperature convergence, foraging competition, and cave avoidance behavior of crayfish, and dynamically balances global exploration (moving away from the current optimal solution) and local development (finely searching nearby areas) through temperature factors (controlling individual movement range) and random disturbance strategies, avoiding premature convergence. As the iteration progresses, the algorithm automatically reduces the search range (similar to a cooling process), gradually shifting from extensive exploration to local optimization, and improving convergence accuracy.

The implementation of COA algorithm mainly consists of six steps [34]:

Step 1: Parameter definition and initialization of population

We define T is the number of iterations, N is population, dim is dimension, ub is upper bound, lb is lower bound, X is initializing population based on upper and lower bounds (formula1). X_i,j is the position of individual i in the j dimension, and X_i,j value is obtained from formula (1). and rand is a random number of Ezugwu et al. (2021) [35].

$$X=[X_1,X_2,\mathrm{...},X_N]=\left[\begin{array}{c}{X}_{1,1}\mathrm{....}{X}_{1,j}\mathrm{...}{X}_{1,\dim }\\ ⋮\\ X_{i,1}\mathrm{....}{X}_{i,j}\mathrm{...}{X}_{i,\dim }\\ ⋮\\ X_{N,1}\mathrm{....}{X}_{N,j}\mathrm{...}{X}_{N,\dim }\end{array}\right]$$

(1)

$$X_{i,j}=l{b}_j+(u{b}_j-l{b}_j)\times rand$$

(2)

Step 2: Define temperature

The ambient temperature of crayfish is defined according to formula (3) to make COA enter different stages. We determine the temperature(temp) and intake of crayfish p [36], as shown in formula (3) .

$$\begin{array}{l}temp=rand\times 15+20\\ p=C_1\times ({\displaystyle \frac{1}{\sqrt{2\times \pi }\times \sigma }}\times \exp (-{\displaystyle \frac{{(temp-\mu )}^2}{2{\sigma }^2}}))\end{array}$$

(3)

µ is the temperature most suitable for crayfish, σ and C₁ are used to control the intake of crayfish at different temperatures.

Step 3: Summer resort stage and competition stage

When the temperature is above 30 (temp>30°), the temperature is too high, crayfish will choose to spend their summer vacation in caves. The definition of X_shade cave is as follows [34]:

$$X_{shade}=(X_G+X_L)/2$$

(4)

X_G is the optimal position obtained so far by the number of iterations, and X_L is the optimal position of the current population. When rand<0.5, crawfish will directly enter the cave for summer vacation, and we have:

$$\begin{array}{l}{X}_{i,j}^{t+1}=X_{i,j}^t+C_2\times rand\times (X_{shade}-X_{i,j}^t)\\ C_2=2-(t/\mathrm{T})\end{array}$$

(5)

$$X_{i,j}^{t+1}=X_{i,j}^t-X_{z,j}^t+X_{shade}$$

(6)

C₂ is a decreasing curve, t is the current iteration number, and t+1 is the next generation iteration number. when rand<0.5, the updated formula(4) is formula (5). Rand ≥ 0.5, COA enters the competitive stage. the two crayfish will compete for the cave according to formula (6).

Step 4: Foraging stage

When temp ≤ 30°, the food intake p and the food size Q are defined according to formula (4) and formula (7).

If Q > (C₃+1)/ 2, shred the food according to formula(8). After that, obtain the new position through formula (9) and proceed to step 5.

$$Q=C_3\times rand\times (fitnes{s}_i/fitnes{s}_{food})$$

(7)

$$X_{food}=\exp (-{\displaystyle \frac{1}{Q}})\times X_{food}$$

(8)

$$X_{i,j}^{t+1}=X_{i,j}^t+X_{food}\times p\times (\cos (2\times \pi \times rand)-\sin (2\times \pi \times rand))$$

(9)

If Q ≤ (C₃+1)/ 2, obtain a new position through formula (10) and proceed to step 5.

$$X_{i,j}^{t+1}=(X_{i,j}^t-X_{food})\times p+p\times rand\times X_{i,j}^t$$

(10)

z is the random individual of crayfish, X_food is the food location, C₃ is the food factor, fitness_i is the fitness value of the ith crayfish, and fitness_food is the fitness value of the food location.

Step 5: Evaluation function

Evaluate the population and determine whether to exit the cycle. If not, return to step 2.

Step 6: Output the best parameter values value.

The core advantage of COA lies in its intelligent biomimetic behavior and dynamic balancing ability, especially in complex and multimodal optimization problems. The algorithm flowchart of COA is shown in Figure 2 [34].

2.2.2. COA optimization VMD

For the modes$K$and penalty factors$\alpha$that have the greatest impact on the VMD decomposition process, improper settings can have a serious impact on the decomposition results. Other parameters are generally set to default values. When set$K$too small, the signal will be under decomposed, and when set$K$too large, the signal will be over decomposed and mode$K$mixing will occur. Therefore, the key to the VMD algorithm is to find the optimal combination of decomposition layers and secondary penalty factors$\alpha$. Therefore, we use the crayfish optimization algorithm (COA) to optimize the parameters of VMD, which can quickly and accurately obtain the optimized parameters. The introduction of VMD principle can be found in reference [26].

When using COA algorithm to optimize VMD, it is particularly important to choose a suitable fitness function. In this paper, envelope entropy is selected as the fitness function optimized by COA(Figure 3). Envelope entropy can better reflect the sparsity and uncertainty of the original signal. When there is a lot of noise in the signal, the entropy value is larger, and conversely, the entropy value is smaller. The principle of envelope entropy [37]:

$$\left.\begin{array}{l}{{\displaystyle E}}_p=-{\displaystyle \sum _{j=1}^N{{\displaystyle p}}_j\lg {{\displaystyle p}}_j}\\ {{\displaystyle p}}_j=\alpha (j)/{\displaystyle \sum _{j-1}^{N}a(j)}\end{array}\right\}$$

(10)

$N$is the number of sampling points of the signal, $p_j$ is the normalized form of$\alpha (j)$, and$\alpha (j)$ is the envelope signal obtained by Hilbert demodulation of the signal$x(j)$.

2.2.3. Multi-scale Permutation Entropy Judgment of COA-VMD Reconstructed Signal

To achieve fast and accurate decomposition, we use the COA to optimize the two key parameters [$K$,$\alpha$]of VMD. For the IMF components decomposed by VMD, Multiscale Permutation Entropy (MPE) is used as the criterion for judging noise and signal. MPE is an improved method based on Permutation Entropy (PE), which calculates permutation entropy [38] at multiple time scales. It has better stability and stronger noise resistance than permutation entropy, and the specific calculation steps are referred to in reference [39]. After calculating the multi-scale permutation entropy of each IMF component, the low-frequency signal and high-frequency noise are determined by setting a threshold, so that, we propose the COA-VMPE algorithm (Figure 4). For the high-frequency noise part, the study does not directly remove it, but further extracts the signal using wavelet decomposition (WD). For a detailed introduction of the WD method, refer to Starck et al. (2017) and Soltani (2002) [40,41].

2.2.4. COA-VMPE-WD Algorithm

The wavelet analysis method overcomes the shortcomings of traditional Fourier transform and has good time-frequency local analysis and multi-resolution analysis performance. Therefore, it is widely used in the study of nonlinear and non-stationary signals [42]. Wavelet decomposition decomposes the original time series into low-frequency and high-frequency components through a set of high pass and low-pass filters, and then decomposes the low-frequency components again. The wavelet basis function is represented as [43]:

$$\left\{\begin{array}{l}{{\displaystyle W}}_f(a,b)=\left|{\displaystyle \frac{1}{\sqrt{a}}}\right|{\displaystyle {\int }_{-\infty }^{+\infty }f(t)\phi (\frac{t-b}{a})dt}\\ {{\displaystyle \phi }}_{a,b}(t)=\left|{\displaystyle \frac{1}{\sqrt{a}}}\right|\phi ({\displaystyle \frac{t-b}{a}})\end{array}\right.$$

(11)

$\phi (t)$is the basis wavelet or mother wavelet function, and the${\phi }_{a,b}(t)$transformed by scale factor$a$and translation factor$b$is collectively referred to as the wavelet.

The key to wavelet decomposition is to choose the appropriate wavelet basis function and determine the decomposition level. In this study, the db4 wavelet with good regularity was selected for decomposition, and the optimal decomposition level was determined by the method in reference [44,45]. We denoise the high-frequency noise obtained by decomposing COA-VMPE using WD to obtain COA-VMPE-WD.

The noise reduction process of COA-VMPE-WD is as follows:

Step1: Initialize the COA algorithm parameters, set the number of crayfish groups to 30, and the maximum number of iterations to 15. Based on considerations of computational efficiency and algorithm accuracy, we set$K$the value range of [3,15] and the $\alpha$value range of [100,3500].

Step2: According to the optimal parameter [$K$,$\alpha$] combination obtained in step 1), perform VMD decomposition on the original reference signal to obtain IMF modal components.

Step3: Calculate the multi-scale permutation entropy of each IMF component, set the threshold for MPE, determine the effective IMF components based on the threshold size and reconstruct them as signals, and reconstruct the remaining components as noise.

Step4: Use wavelet decomposition to denoise the high-frequency noise portion reconstructed in step 3) again, and use correlation coefficients to determine the effective signal. Reconstruct the low-frequency signal obtained in step 3) with the signal processed by wavelet decomposition to obtain the final denoised signal. The specific noise reduction process of COA-VMPE-WD is shown in Figure 5:

2.2.5. Evaluation Index

In order to quantitatively demonstrate the denoising effect of the above methods, we selected root mean square error (RMSE)[46], signal-to-noise ratio (SNR), mean absolute error (MAE), and cross-correlation coefficient (R) as the denoising evaluation indicators, and their calculation formulas are as follows [47]:

$$RMSE=\sqrt{{\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N{(x_i-y_i)}^2}}$$

(12)

$$SNR=10\times \lg \left({\displaystyle \frac{{\displaystyle \sum _{i=1}^N{x_i}^2}}{{\displaystyle \sum _{i=1}^N{(x_i-y_i)}^2}}}\right)$$

(13)

$$MAE={\displaystyle \frac{1}{N}}{\displaystyle \sum _{i=1}^N\left|(x_i-y_i)\right|}$$

(14)

$$R(x,y)={\displaystyle \frac{{\displaystyle \sum _{i=1}^N(x_i-\overline{x})(y_i-\overline{y})}}{\sqrt{{\displaystyle \sum _{i=1}^N{(x_i-\overline{x})}^2}}\sqrt{{\displaystyle \sum _{i=1}^N{(y_i-\overline{y})}^2}}}}$$

(15)

$i$ is the epoch time; $x$、$y$ are respectively the denoised signal and the original signal; $\overline{x}$ and $\overline{y}$ are the average values of $x$and$y$, respectively; $N$is data length [48].

3. Results

3.1. Simulation Signal Decomposition

To verify the reliability of the algorithm, we used simulated signals with known true values to simulate measured signals for experiments, and compared them with existing EMD, EEMD, CEEMDAN, and WD methods. During the experiment, EMD, EEMD, CEEMDAN, and WD methods all used the correlation coefficient method to separate noise and signals. The simulated signal consists of three periodic terms and Gaussian white noise, with a sampling interval of 1 second and 4636 sampling points. The ratio of Gaussian to white noise is 9:1. The component waveform of the simulated signal is shown in Figure 3, and its mathematical expression is [49]:

$$\left.\begin{array}{l}{y}_1=6\sin (2\pi t/600)\sin (2\pi t/250)\\ y_2=4\sin (2\pi t/400)+2\cos (2\pi t/200)\\ y_3=2\sin (2\pi t/50)+\sin (16\pi t)\cos (150\pi t)\\ \epsilon =noise\\ y=y_1+y_2+y_3+\epsilon \end{array}\right\}$$

(16)

When using COA-VMD for noise reduction, the COA optimization algorithm is first used to find the optimal parameter combination for VMD decomposition. The envelope entropy is used as the fitness function, and the fitness function value changes with the number of iterations during the COA optimization process, as shown in Figure 5. The fitness value in Figure 6 reaches its minimum at the fourth iteration, at which point[K,$\alpha$]=[8,1425] is the optimal parameter combination.

We used the optimal parameter combination obtained from COA to perform VMD decomposition on the signal, and obtained 8 IMF components as shown in Figure 7. From Figure 6, the low-frequency components are mainly concentrated in the first 2 modes.

To effectively separate low-frequency signals and high-frequency noise, it is necessary to calculate the Multi-Scale Permutation Entropy (MPE) values of each IMF component. When calculating MPE, appropriate parameters need to be set. Chen et al. (2023) set the scale factor to 12, embedding dimension to 6, and delay time to 1 [50]; For the setting of key parameters, we also take $s=12$, $m=6$, $\tau =1$ to calculate the permutation entropy mean of each IMF component under different scale factors as the final MPE value. The closer the MPE value is to 1, the greater the random volatility and higher irregularity of the time series. Through experimental analysis, we set the MPE threshold to 0.6 and consider IMF components greater than 0.6 as noise components, while those less than 0.6 are considered low-frequency signal components. The MPE values of each IMF component decomposed by VMD are shown in Figure 8 and Table 1.

From Figure 8 and Table 1, the MPE values of IMF1~IMF8 gradually increase, indicating that the random volatility of the sequence is increasing and the noise components are gradually increasing, which is consistent with the waveform description in Figure 9(a). The MPE values of IMF1~IMF2 are all less than 0.6, therefore, IMF1~IMF2 are reconstructed as low-frequency signals, IMF3~IMF8 are reconstructed as high-frequency noise, and then wavelet decomposition is used to process the high-frequency noise part. The results of COA-VMPE-WD decomposition are shown in Figure 9(b).

3.2. Evaluation of Noise Reduction Effect of Simulated Signals

To further compares the results of WD denoising after MPE judgment reconstruction, we used WD、EMD、EEMD and CEEMDAN to denoise the simulated signals and compared them with COA-VMPE-WD. Figure 10 shows the waveform comparison of the three denoising methods used. From Figure 10, the denoised signal waveforms of WD、EMD、EEMD and CEEMDAN have poor fitting effects with the original sequence, while COA-VMPE-WD, compared to WD、EMD、EEMD and CEEMDAN has a waveform that is closer to the simulated signal after denoising. The curve is smoother, avoiding the endpoint effect and mode mixing problem in the denoising process of EMD and other methods, which can extract more effective signals and achieve better denoising effect. The evaluation indicators for noise reduction using different methods are shown in Table 2.

According to Table 2, for simulated signals, COA-VMPE-WD has the best noise reduction evaluation indicators compared to WD, EMD, EEMD and CEEMDAN methods. The RMSE (0.2291 mm) of COA-VMPE-WD is reduced by 26.5% compared to the suboptimal EMD, and the MAE (0.1780 mm) is reduced by 23.6%, indicating its outstanding residual control ability; The SNR of COA-VMPE-WD (27.8397 dB) is improved by 10.7% compared to EMD, indicating effective separation of noise and signal; The R (0.9757) of COA-VMPE-WD is close to 1, verifying its stronger signal fidelity. EMD is limited by endpoint effects and modal aliasing, and its accuracy is still lower than COA-VMPE-WD. WD、EEMD and CEEMDAN are misaligned due to noise addition strategies, resulting in significant degradation of RMSE and MAE, and a sharp drop in SNR (about 13 dB), reflecting their insufficient anti-interference ability. In summary, the advantage of COA-VMPE-WD lies in the coupled COA optimized VMD, which adaptively decomposes modes through MPE entropy constraints, avoiding the mode mixing problem of traditional EMD series.

3.3. Application of COA-VMPE-WD for GNSS Time Series Denoising Analysis

We conducted experiments on the time series of 52 GNSS stations with E (East), N (North), and U (Up) components using COA-VMPE-VMD, and the processing process was consistent with the simulation signal experiment. Among them, the results of COA optimized VMD parameters in the COA-VMPE-VMD method are shown in Figure 11 and Figure 12, the K value of the U component fluctuates greatly, while the K values of the E and U components are mainly concentrated between 5 and 11, and the values of $\alpha$each component mainly fluctuate within the range of 800 to 3600. We present in the appendix the adaptive curve of FJPT station solved by COA-VMPE-VMD method in Figure A1, the modal number distribution decomposed by COA-VMPE-VMD method in Figure A2, and the IMF value distribution of FJPT station solved by COA-VMPE-VMD method in Figure A3, Figure A4 and Figure A5.

3.4. Comparative Analysis of COA-VMPE-WD

The study takes FJPT station as an example to analyze the performance of COA-VMPE-WD denoising method. Figure 13 shows the denoising effects of five methods: WD, EMD, EEMD, CEEMDAN, and COA-VMPE-WD. From Figure 13, all five methods can effectively extract nonlinear changes from the original data, and the denoised signal exhibits significant periodic changes in the U component. Compared with the other four methods, the COA-VMPE-WD denoised signal has a smoother curve, indicating that it effectively filters out high-frequency noise. However, the denoised signals of the other four methods still contain some high-frequency noise, especially in the EMD method.

Figure 14 shows the residual comparison between the denoised signal obtained using five methods and the original time series. As shown in the figure, the residual fluctuations are significant after EMD and EEMD denoising. The COA-VMPE-WD method obtains the smallest average residual between the denoised signal and the original signal, indicating that this method effectively filters out high-frequency noise and has a more significant denoising effect.

To better demonstrate the effectiveness of our method in denoising GNSS measured data, we evaluated the accuracy of 52 stations after denoising processing. Figure 15 shows the comparison of four evaluation index values of the denoised signals from 52 GNSS reference stations using the above five methods. The results of NEU component are shown in (a) and (b) of Figure 15, Figure 16 and Figure 17. The COA-VMPE-WD method has the smallest RMSE and MAE values, and the WD denoising method has the worst performance The noise reduction effect of EMD method is consistent with that of EEMD method in some stations. Overall, the difference in noise reduction effect between the two methods is not significant and relatively close. The reason for this may be due to the complex crustal tectonic movements in some stations, and further analysis is needed to determine the specific reasons; Compared with CEEMDAN, EMD, and EEMD methods, the RMSE and MAE obtained by denoising with COA-VMPE-WD method are smaller. As shown in Figure 15, Figure 16 and Figure 17, the R (close to 1) and SNR after denoising with COA-VMPE-WD method are larger. All evaluation indicators after denoising with COA-VMPE-WD method are optimal, indicating that its denoising effect is better and it can effectively extract more useful signals.

4. Discussion

4.1. Analysis of Optimal Noise Models Using Different Methods

In order to illustrate the denoising effect of the method proposed in this paper more intuitively, different noise models assume that there are differences in estimating the velocity and uncertainty of the reference station. Therefore, the study selects four types of noise models: FN+WN, FN+RW, FN+RW+WN, and PL+WN (where WN is white noise, FN is flicker noise, RW is random walk noise, and PL is power-law noise). The maximum likelihood estimation method of Hector software is used for analysis. The optimal noise models for the E, N, and U components of each reference station are determined by an improved Bayesian information criterion (BIC_tp) [4], and the noise amplitude and velocity uncertainty of each component before and after denoising in the optimal noise model are obtained. The results are shown in Table 3 to Table 4. As shown in Figure 18.

According to Table3, before denoising, the N component of the original GNSS coordinate time series is dominated by the PLWN noise model, the optimal noise models for the E component coordinate time series are dominated by PLWN and FNRWWN, and the U component is dominated by FNWN and FNRWWN. Table4 shows the optimal noise model of GNSS stations after different denoising methods. Table4 compares the noise type changes of different components (N, E, U) after processing with the original noise model (Origin) and five denoising methods (WD, EMD, EEMD, CEEMDAN, COA-VMPE-WD) On the N-component, the original noise is mainly PLWN, and all denoising methods can effectively preserve the PLWN characteristics without introducing high-order noise, indicating that the methods have good adaptability to the N-component. On the E component, the original noise includes PLWN and FNRWWN. After denoising, most methods (such as WD and EMD) simplify FNRWWN to PLWN. However, COA-VMPE-WD detects the characteristics of FNRW noise at some stations while detecting PLWN. On the U-component, the original noise is mainly FNRWWN and FNWN. After denoising, most methods (such as CEEMDAN) optimize it to PLWN, but COA-VMPE-WD still generates FNRW in FNWN, indicating that COA-VMPE-WD can detect the influence of RW. The U-component noise is more complex, and RW has the greatest impact on the U-component.

4.2. Analysis of GNSS Station Velocity And Annual Term Using Different Methods

Figure 18 shows the noise amplitude distribution of E, N, and U components before and after denoising under the main noise model. According to Table 4, it can be seen that in the E, N, and U components, after denoising by WD, EMD, EEMD, CEEMDAN, and COA-VMPE-WD methods, RW noise can be effectively detected, and the average amplitude of colored noise is also significantly reduced. Moreover, the COA-VMPE-WD method reduces the average amplitude of colored noise compared to the other four methods. Overall, the COA-VMPE-WD method achieves good denoising effect. However, in Figure 18, the U component has the highest amplitude value under the CEEMDAN denoised PLWN model. This may be due to the significant impact of RW in the noise model on the amplitude of denoised PL. When the noise model is PLWN, the amplitude of denoised FN decreases significantly but cannot be completely eliminated. The reason for this phenomenon may be that the noise characteristics of the denoised data have changed, and the noise model does not match FNWN. When using the PLWN model to solve the noise in the denoised signal, RW is transformed into PL.

Figure 19 shows the uncertainty of GNSS station velocity before and after denoising, and the average GNSS station velocity obtained from the denoised signal compared to the original data is shown in Table 5. From Figure 19, it can be seen that the E, N, and U components of the majority of reference stations are denoised using WD, EMD, WD, CEEMDAN, and COA-VMPE-WD methods, resulting in varying degrees of reduction in station velocity and velocity uncertainty for each component. The COA-VMPE-WD method reduces the average velocity of GNSS stations after denoising by 44.44%, 50.00%, 52.38%and 44.44% compared to the WD, EMD, EEMD, and CEEMDAN methods with N-component, respectively. On the E component, the COA-VMPE-WD method reduced station velocity by an average of 14.29%, 18.18%, 21.74%, and 21.74% compared to the WD, EMD, EEMD, and CEEMDAN methods. On the U component, the COA-VMPE-WD method reduced station velocity by an average of 50.00%, 59.09%, 18.18%, and 64.00% compared to the WD, EMD, EEMD, and CEEMDAN methods. To sum up, the CEEMDAN method has the most significant effect on reducing velocity uncertainty, followed by the WD method. The EMD and EEMD methods are poor and both have situations where the velocity uncertainty of some reference stations increases (such as the U component of SDCY station, the E component of XJBL station, etc.).

According to Table 5, on the N component, the original mean station velocity is 0.018mm/a, while WD and CEEMDAN maintain their original values. EMD and EEMD slightly increase (0.020mm/a-0.021 mm/a), while COA-VMPE-WD significantly decreases to 0.010 mm/a, indicating that this method has the strongest noise suppression effect on the N component. On the E component, the original mean station velocity is 0.023 mm/a, and the results of various methods are relatively close (0.018 mm/a -0.023 mm/a). Among them, COA-VMPE-WD has the largest reduction (0.018 mm/a), while other methods basically maintain the original characteristics. On the U component, there are significant differences between different methods: EMD remains at its original value, CEEMDAN increases to 0.025 mm/a, EEMD and COA-VMPE-WD decrease to 0.011 mm/a and 0.009 mm/a, respectively, reflecting the complex noise characteristics of the U component and the significant differences in processing effects between different algorithms. In summary, COA-VMPE-WD exhibits strong noise suppression ability for all components, while CEEMDAN may excessively preserve noise features for the U component. There are significant differences in the response characteristics of denoising methods for different components, and the U component has the highest dispersion in processing effect.

5. Conclusions

The study proposes a COA-VMPE-WD denoising algorithm based on the nonlinear and non-stationary characteristics of GNSS coordinate time series. In response to the problem that traditional VMD cannot adaptively determine the values of modal number and penalty factor [$K$,$\alpha$]in the denoising process, resulting in insufficient signal decomposition or over decomposition affecting the effective denoising of time series, a COA optimized VMD parameter with envelope entropy as the fitness function is proposed. Multi scale permutation entropy (MPE) is used as the screening criterion for noise and signal, and wavelet decomposition algorithm (WD) is combined to improve VMD. This article conducts noise reduction research and analysis through simulated signals and measured data from 52 GNSS stations, verifying the reliability and effectiveness of the method. It compares it with traditional WD, EMD, EEMD, and CEEMDAN method. Our main conclusions are as follows:

(1)

Simulation signal experiments show that compared to traditional EMD and WD methods, VMD method can effectively alleviate modal aliasing and endpoint effects, and has better signal feature extraction ability. In terms of noise reduction, the COA-VMPE-WD method outperforms the WD, EMD, EEMD, and CEEMDAN methods in terms of RMSE, R, and SNR evaluation metrics, effectively removing noise from the original data.

(2)

Experimental data shows that WD, EMD, EEMD, and CEEMDAN methods can remove white noise from the original coordinate time series and significantly reduce the amplitude of colored noise. Compared with WD, EMD, EEMD, and CEEMDAN methods, COA-VMPE-WD has a more significant noise reduction effect and better preserves the characteristics of the original signal.

(3)

The WD, EMD, EEMD, CEEMDAN COA-VMPE-WD methods have a significant impact on the U-component in terms of the velocity and uncertainty of the reference station, the COA-VMPE-WD method reduced station velocity by an average of 50.00%, 59.09%, 18.18%, and 64.00% compared to the WD, EMD, EEMD, and CEEMDAN methods, The noise reduction effect is higher than the other four methods. The above results verify the effectiveness and reliability of the COA-VMPE-WD denoising method. However, the method proposed in this paper reduces computational efficiency to a certain extent. Therefore, further research is needed on how to construct a simple and efficient denoising method based on VMD.