# A Vibration Signal Denoising Method of Marine Atomic Gravimeter Based on Improved Variational Mode Decomposition

Authors: Wenbin Gong; An Li; An Liao; Hao Che; Chunfu Huang; Fangjun Qin
Journal: Wireless Communications and Mobile Computing
 (2022)

---

## Abstract

As a high-precision gravity measuring device, a marine atomic gravimeter is highly sensitive to vibration signals. Accurate measurement and analysis of vibration signal is the primary condition to realize vibration compensation and vibration suppression. Denoising plays a crucial role in the processing of these vibration signals. The vibration signals of a marine gravimeter contain numerous nonlinear and nonstationary components. In this paper, a vibration signal denoising method of marine atomic gravimeter based on improved variational mode decomposition (VMD) was put forward to effectively suppress the noise. An improved genetic particle swarm optimization (GPSO) was first adopted for the parametric optimization of VMD by taking minimum permutation entropy (PE) as fitness function and adaptively determining the optimal parameters of VMD. PE was then utilized to calculate the proportion of noise-containing components in the intrinsic mode function (IMF) components obtained by VMD. The components were classified into noise and signal components by searching for the mutation points of two adjacent IMF permutation entropies. On this basis, noise components were denoised by Savitzky-Golay (SG) filter. In the end, the denoised components were reconstructed with the signal components to generate denoised vibration signals. To verify the effectiveness, the proposed method was applied in denoising, simulated and measured vibration signals of a marine atomic gravimeter, and compared with Daubechies (db) wavelet, Symlets (sym) wavelet, and empirical mode decomposition (EMD). The results showed that the proposed method could effectively remove the noise from nonlinear vibration signals and retain the authentic and useful information, so that it was able to provide the supporting data for gravity compensation of marine atomic gravimeter.


---

## Body

## 1. Introduction

As a high-precision absolute gravimeter, an atomic gravimeter may be applied in inertial navigation and earth gravity field measurement [1]. In terms of marine measurement, Bidel et al. [2] employed an atomic gravimeter in marine gravity measurement for the first time. An atomic gravimeter was attached to a gyroscope stabilized platform, so that it could always measure the actual virtual component of gravity field regardless of waves and sways. The platform overcame effectively the poor verticality of a ship-borne atomic gravimeter in measurement. Based on an inertial stabilized platform, Cheng et al. [3] conducted an absolute measurement experiment with a ship-borne atomic gravimeter when the ship was moored, and Li et al. [4] performed a further lake test of an atomic gravimeter in absolute gravity measurement.During measurement, an atomic gravimeter is severely affected by the vibration of Raman retroreflection mirror since its measurement precision and reliability are highly restricted by vibration noise. For this reason, vibration noise isolation and attenuation are crucial to obtaining accurate atomic interference phase and realizing accurate gravity field information detection. In order to guarantee the measurement precision of an atomic gravity, vibration isolation techniques [5, 6] and vibration compensation methods [7–9] are often employed for vibration noise suppression. The noise in the marine environment is very complex and vulnerable to impulse noise, so it is difficult to process the signal [10]. In the structure of a ship, noise is mainly caused by the vibration of main engine, diesel engine, main propulsion system, main propeller, and other devices. The noise caused by crew activities intermittently is also coupled into the signals of main vibration noise. These vibration signals are nonlinear and nonstationary [11–14]. For the purpose of vibration noise suppression in measurement with an atomic gravimeter, accurate measurement and analysis of vibration signals lays a basis for vibration compensation and suppression. Normally, an accelerometer is utilized to gather vibration signals, but these vibration signals are inevitably mingled with ambient noise and circuit noise in the process. Denoising analysis must be therefore conducted with the gathered vibration data to extract the actual vibration information of vibration sources, so as to provide the reliable data for vibration isolation design and vibration compensation [15].Wavelet transform denoising is a denoising method suitable for nonstationary signals. It has been widely applied in engineering, but its final effect is affected by the selection of threshold and wavelet basis function [16]. In 1998, Huang et al. put forth an analysis method for nonlinear and nonstationary signals based on Hilbert-Huang transform, that is, empirical mode decomposition (EMD), which could adaptively decompose a signal into a number of intrinsic mode functions (IMFs) with actual physical implications. The method had been widely applied in denoising nonlinear and nonstationary signals. Kopsinis and McLaughlin combined wavelet decomposition with EMD to denoise signals and used different thresholds for IMFs in filtering and reconstruction to realize signal denoising [17]. Rezaee and Osguei [18] made an improvement to EMD by introducing a new parameter to obtain a new local mean. In this way, they enhanced the precision and efficiency of EMD and effectively applied it in the analysis of vibration signals. Nevertheless, the application of EMD was also troubled by end effect and mode mixing especially when the signal to noise ratio (SNR) was low. In order to resolve this problem, some improvements of EMD have been explored including ensemble empirical mode decomposition (EEMD), improved complete ensemble empirical mode decomposition (ICEEMD), and partial ensemble empirical mode decomposition (PEEMD) [16, 19, 20]. Nevertheless, these methods can inhibit mode mixing to some extent but intrinsically extract local extremum and interpolate envelopes. They are still empirical and lack a solid mathematic basis.In order to effectively inhibit the mode mixing of EMD, Dragomiretski and Zosso [21] put forward variational mode decomposition (VMD) in 2014. Based on the three-dimensional variational constraint theory, this algorithm estimated multiple modes simultaneously by virtue of nonrecursion and improved the computational efficiency while guaranteeing the integrity of features. Hence, it could satisfactorily resist noise and reduce mode mixing. With its solid mathematic basis for the decomposition of nonstationary signals, VMD has been applied in earthquake time-frequency analysis, signal filter denoising, and ground vibration attenuation, which is a sufficient proof of its effectiveness and superiority in signal decomposition [22–24]. Nevertheless, number of decomposed modes K and penalty term α must be artificially set in the VMD of signals. If the value of K is set too high, overdecomposition may be caused to generate false components. If it is set too low, under decomposition occurs and results in the mixing of modes close to frequency. Moreover, the penalty factor α also affects the extraction of single-component modes. If it is set too large, the bandwidth will be narrower for single modes, causing to intercept the effective frequency components outside bandwidth. If it is set too small, the bandwidth will be wider, and the two adjacent modes will share the center frequency and result in information redundancy [25, 26].In order to adaptively decompose signals, some methods have been developed to determine theK value based on kurtosis [27] and energy factor [28]. In these methods, only the number of decomposed modes K is optimized while the penalty factor α is ignored. Therefore, optimal decomposition cannot be achieved with these methods. Along with the emergence of intelligent optimization algorithms, attempts have been made to apply some optimization algorithms in the optimization of VMD parameters, and satisfying results have been obtained [24, 29, 30]. Zhou et al. [31] put forward the particle swarm optimization (PSO) to optimize the VMD parameters. In this method, they used mean permutation entropy (MPE) as its fitness function and determined the optimal combination of K and α by searching for the minimum of the fitness function. It was an efficient search algorithm because of its fast convergence while requiring the setting of fewer parameters. Nevertheless, the PSO is troubled by premature convergence and faces slower convergence in the late stage since population diversity disappears in the searched space. Meanwhile, it cannot be further optimized after reaching a certain precision of convergence, so that its final precision is not good. Based on genetic algorithm (GA), Kumar et al. took kernel-based mutual information (KEMI) as fitness function to find out the optimal parameters K and α of VMD [29]. The GA algorithm has strong global search capability and avoids local optimum, but it may be easily affected by such problems as premature convergence, numerous computations, slow convergence, and poor stability.To solve the above problems, a vibration signal denoising method combining improved VMD parameter optimization algorithm and Savitzky-Golay (SG) filter is proposed in this paper. Firstly, the improved genetic particle swarm optimization (GPSO) is applied to VMD, and the permutation entropy (PE) is used as the fitness function to optimize the VMD parameters and decompose the noisy vibration signal. Secondly, the noise proportion in IMF is calculated by PE, and the IMF component is divided into noise component and the signal component. SG filter was applied for denoising noisy components. Finally, the denoised component and signal component are reconstructed to obtain the denoised vibration signal. The method is applied to the analysis of vibration simulation signals and measured signals.


## 2. Theoretical Background

### 2.1. Variational Mode Decomposition

As for the structural variation of signals in VMD, Wiener filter is introduced to solve variation by virtue of iterative computation. Each input signal is decomposed intoK IMFs with different center frequencies ωk. The variational model is correspondingly described by Equation (1).
(1)minuk,ωk∑k=1K∂tδt+jπt∗ukte−jωkt22s.t.∑k=1Kuk=f,where uk=u1,⋯,uk is K IMF components after decomposition, ωk=ω1,⋯,ωk is the center frequency of each component, δt is the unit pulse function, j is the imaginary unit, ∗ is the convolution operation, and ∂t is the partial derivative with respect to t. Lagrangian multiplier λ and secondary penalty term α are introduced to obtain the optimal solution of constrained variation as follows:
(2)Luk,ωk,λ=α∑k∂tδt+jπt∗ukte−jωkt22+ft−∑iuit+λt222+λt,ft−∑kukt.The alternation of multiplication operators is conducted to iteratively updateuk,ωk,λ and obtain the optimal mode component uk, center frequency ωk, and Lagrangian multiplier λ.
(3)u^kn+1ω=f^ω−∑i≠ku^iω+λ^ω/21+2αω−ωk2,(4)ωkn+1=∫0∞ωu∧ω2dω∫0∞u^ωdω,(5)λ∧n+1ω=λ∧nω+τf^ω−∑ku∧n+1ω,where ∧ is the Fourier transform, τ is the fidelity coefficient, and n is the iteration times. Iteration is constantly updated but halted when relative error is less than convergence precision. The decomposition result uk is eventually outputted.
(6)∑ku^kn+1ω−u^knω22u^kn22<ε.


### 2.2. GPSO Optimized VMD

The effect of VMD depends significantly on the number of decomposed modesK and penalty term α. Considering their limitations in VMD, GA and PSO are combined in this paper. In this way, genetic computation is adopted to obtain the optimal solution of PSO. Under specific conditions, genetic operations including reproduction, crossover, and mutation are carried out for particles and their displacement rate, so that the particles can be alienated from local optimal to obtain the global optimal solution. The GPSO has quicker convergence and better grouping quality than traditional PSO. The combined algorithm makes full use of the advantages of both GA and PSO, so as to guarantee the quick convergence to global optimal solution. The process flow is presented in Figure 1. The specific steps are given as follows:
(1)
Randomly initialize the particles in a population and set their corresponding parameters. Generaten particles randomly with Xi=pi,vii=1,2,⋯,n, where pi and viare the geometrical location and velocity vector. So the initial generation of particle swarm t=0 is identified as(7)X0=X10,X20,⋯,Xn0=p10,v10,p20,v20,⋯,pn0,vn0.(2)
VMD Decomposition. Calculate the fitness of individuals in the initial population, and select the individual’s optimal value Ppbt and the global optimal value Pgbt in the particle swarm. In this paper, PE is taken as the fitness function, and minimum PE is employed to determine optimal solution. By calculating the fitness function, the complexity of the signal is obtained from the PE value. The more complex the signal is, the greater the calculated PE value is, vice versa. After the vibration signal is decomposed by VMD, if there are many noise components included in the IMF component, the higher the complexity of the signal is, the greater the PE value is. If a few noise components are included in the IMF component, the stronger the regularity of the signal, the simpler the signal, and the lower the PE value.(3)
Update the position and velocity of each particle. For each particleXit=pit,vit, we let(8)pit+1=pit+vit+1,(9)vit+1=C1vit+C2r10,1Ppbt−Pit+C3r20,1Pgbt−Pit,where r10,1 and r20,1 are the random numbers in (0, 1), C1 is the inertia weight, and C2 and C3 are the learning factors. Therefore, the t+1th generation of particle swarm is formed as follows:
(10)Xt+1=X1t+1,X2t+1,⋯,Xnt+1=p1t+1,v1t+1,p2t+1,v2t+1,⋯,pnt+1,vnt+1.Figure 1
Process flow of VMD parametric optimization by GPSO.Larger inertia weight helps get out of local minimum point and facilitates global search, while smaller inertia weight is conductive to precise local search in the current region of search and helpful to the convergence of algorithm. For this reason, a linearly changing weight can be employed. An inertia weight decreases linearly from the maximum valueCmax to the minimum value Cmin. It varies with the iteration times of the algorithm as follows:
(11)C1=Cmax−t∗Cmax−Cmintmax,where t is the current iteration steps, tmax is the maximum iteration steps, and Cmax and Cmin are the maximum and minimum inertia weights, respectively, which are set to 1.2 and 0.6 in this paper.
(4)
Crossover Operation. Apply a crossover operator to a population, and switch some chromosomes in each selected pair of individuals at a probability to generate new individuals. The crossover probability is set to 0.8 in this paper.(5)
Mutation Operation. Apply a mutation operator to the population. Change any or some genes of the selected individuals to other allele(s) at a probability. For a population Pt, selection, crossover, and mutation operations are conducted to obtain its next-generation population Pt+1, whose fitness is calculated. The fitness is then used for sequencing. These genetic operations will be repeated. The crossover probability is set to 0.3 in this paper.(6)
Recalculate the fitness of particles in the new population, update the optimal solution (Ppbt and Ggbt) of the population based on the fitness, and calculate the optimal particles Xpbt=Ppbt,vpbt that have been found so far for each particle i. Calculate the optimal particles Xgbt=Pgbt,vgbt that have been found so far for the current population Xt.(7)
Determine whether the conditions for the end of iteration are satisfied. If not, return to Step 3.


### 2.3. Fitness Function

When the GPSO is employed to optimize the VMD parameters, a fitness function must be determined to evaluate its optimization results. PE is a mean entropy parameter to measure the complexity of one-dimensional time series, which can be used to detect dynamic mutation and time series randomness [32]. Phase space reconstruction is carried out for a set of time series Xi,i=1,2,⋯,N to obtain a matrix Y:
(12)Y=x1x1+τ⋯x1+d−1τx2x2+τ⋯x2+d−1τxjxj+τ⋯xj+d−1τ⋮⋮⋮xKxK+τ⋯xK+d−1τ,j=1,2,⋯,K,where d is the embedding dimensions, τ is the delay time, K is the number of reconstructed components, and xj is the components in the Jth column of the reconstruction matrix. The elements in each column of the reconstruction matrix Y are reorganized in ascending order, so that a group of symbol sequence Sl=j1,j2,⋯,jd can be obtained for each column of the matrix Y. The occurrence probability of each symbol sequence P1,P2,⋯,Pk can be calculated. At this time, the PE Hp of K different symbol sequences for the time series Xi can be defined in the form of Shannon entropy as Hpd=−∑j=1kPjlnPj.In the above definition, the value ofHp represents how time series is stochastic. A smaller value implies simpler time series, while a larger value leads to more complex and stochastic time series. Hence, minimum PE is taken as a fitness function in this paper and then solved to determine optimal decomposition parameters.


### 2.4. Savitzky-Golay Filter

SG filter, a method put forth by Savitzky and Golay, is widely applied in denoising the nonstationary signals containing non-Gaussian noise [33]. In the method, univariate P-order polynomials are adopted to fit the fixed length neighborhood of each data point in the selected data. Least squares criterion is followed to determine the polynomial coefficients by minimizing fitting error, so as to obtain the optimal fitted value of the data point, which is the value obtained after denoising. In this way, signals are denoised. At the time of denoising, the SG filter method can effectively retain the variation information of signals.
(13)min∑j=−mmYj−yj2,s.t.Yi=c0+c1i+c2i2+⋯+cpip.


## 2.1. Variational Mode Decomposition

As for the structural variation of signals in VMD, Wiener filter is introduced to solve variation by virtue of iterative computation. Each input signal is decomposed intoK IMFs with different center frequencies ωk. The variational model is correspondingly described by Equation (1).
(1)minuk,ωk∑k=1K∂tδt+jπt∗ukte−jωkt22s.t.∑k=1Kuk=f,where uk=u1,⋯,uk is K IMF components after decomposition, ωk=ω1,⋯,ωk is the center frequency of each component, δt is the unit pulse function, j is the imaginary unit, ∗ is the convolution operation, and ∂t is the partial derivative with respect to t. Lagrangian multiplier λ and secondary penalty term α are introduced to obtain the optimal solution of constrained variation as follows:
(2)Luk,ωk,λ=α∑k∂tδt+jπt∗ukte−jωkt22+ft−∑iuit+λt222+λt,ft−∑kukt.The alternation of multiplication operators is conducted to iteratively updateuk,ωk,λ and obtain the optimal mode component uk, center frequency ωk, and Lagrangian multiplier λ.
(3)u^kn+1ω=f^ω−∑i≠ku^iω+λ^ω/21+2αω−ωk2,(4)ωkn+1=∫0∞ωu∧ω2dω∫0∞u^ωdω,(5)λ∧n+1ω=λ∧nω+τf^ω−∑ku∧n+1ω,where ∧ is the Fourier transform, τ is the fidelity coefficient, and n is the iteration times. Iteration is constantly updated but halted when relative error is less than convergence precision. The decomposition result uk is eventually outputted.
(6)∑ku^kn+1ω−u^knω22u^kn22<ε.


## 2.2. GPSO Optimized VMD

The effect of VMD depends significantly on the number of decomposed modesK and penalty term α. Considering their limitations in VMD, GA and PSO are combined in this paper. In this way, genetic computation is adopted to obtain the optimal solution of PSO. Under specific conditions, genetic operations including reproduction, crossover, and mutation are carried out for particles and their displacement rate, so that the particles can be alienated from local optimal to obtain the global optimal solution. The GPSO has quicker convergence and better grouping quality than traditional PSO. The combined algorithm makes full use of the advantages of both GA and PSO, so as to guarantee the quick convergence to global optimal solution. The process flow is presented in Figure 1. The specific steps are given as follows:
(1)
Randomly initialize the particles in a population and set their corresponding parameters. Generaten particles randomly with Xi=pi,vii=1,2,⋯,n, where pi and viare the geometrical location and velocity vector. So the initial generation of particle swarm t=0 is identified as(7)X0=X10,X20,⋯,Xn0=p10,v10,p20,v20,⋯,pn0,vn0.(2)
VMD Decomposition. Calculate the fitness of individuals in the initial population, and select the individual’s optimal value Ppbt and the global optimal value Pgbt in the particle swarm. In this paper, PE is taken as the fitness function, and minimum PE is employed to determine optimal solution. By calculating the fitness function, the complexity of the signal is obtained from the PE value. The more complex the signal is, the greater the calculated PE value is, vice versa. After the vibration signal is decomposed by VMD, if there are many noise components included in the IMF component, the higher the complexity of the signal is, the greater the PE value is. If a few noise components are included in the IMF component, the stronger the regularity of the signal, the simpler the signal, and the lower the PE value.(3)
Update the position and velocity of each particle. For each particleXit=pit,vit, we let(8)pit+1=pit+vit+1,(9)vit+1=C1vit+C2r10,1Ppbt−Pit+C3r20,1Pgbt−Pit,where r10,1 and r20,1 are the random numbers in (0, 1), C1 is the inertia weight, and C2 and C3 are the learning factors. Therefore, the t+1th generation of particle swarm is formed as follows:
(10)Xt+1=X1t+1,X2t+1,⋯,Xnt+1=p1t+1,v1t+1,p2t+1,v2t+1,⋯,pnt+1,vnt+1.Figure 1
Process flow of VMD parametric optimization by GPSO.Larger inertia weight helps get out of local minimum point and facilitates global search, while smaller inertia weight is conductive to precise local search in the current region of search and helpful to the convergence of algorithm. For this reason, a linearly changing weight can be employed. An inertia weight decreases linearly from the maximum valueCmax to the minimum value Cmin. It varies with the iteration times of the algorithm as follows:
(11)C1=Cmax−t∗Cmax−Cmintmax,where t is the current iteration steps, tmax is the maximum iteration steps, and Cmax and Cmin are the maximum and minimum inertia weights, respectively, which are set to 1.2 and 0.6 in this paper.
(4)
Crossover Operation. Apply a crossover operator to a population, and switch some chromosomes in each selected pair of individuals at a probability to generate new individuals. The crossover probability is set to 0.8 in this paper.(5)
Mutation Operation. Apply a mutation operator to the population. Change any or some genes of the selected individuals to other allele(s) at a probability. For a population Pt, selection, crossover, and mutation operations are conducted to obtain its next-generation population Pt+1, whose fitness is calculated. The fitness is then used for sequencing. These genetic operations will be repeated. The crossover probability is set to 0.3 in this paper.(6)
Recalculate the fitness of particles in the new population, update the optimal solution (Ppbt and Ggbt) of the population based on the fitness, and calculate the optimal particles Xpbt=Ppbt,vpbt that have been found so far for each particle i. Calculate the optimal particles Xgbt=Pgbt,vgbt that have been found so far for the current population Xt.(7)
Determine whether the conditions for the end of iteration are satisfied. If not, return to Step 3.


## 2.3. Fitness Function

When the GPSO is employed to optimize the VMD parameters, a fitness function must be determined to evaluate its optimization results. PE is a mean entropy parameter to measure the complexity of one-dimensional time series, which can be used to detect dynamic mutation and time series randomness [32]. Phase space reconstruction is carried out for a set of time series Xi,i=1,2,⋯,N to obtain a matrix Y:
(12)Y=x1x1+τ⋯x1+d−1τx2x2+τ⋯x2+d−1τxjxj+τ⋯xj+d−1τ⋮⋮⋮xKxK+τ⋯xK+d−1τ,j=1,2,⋯,K,where d is the embedding dimensions, τ is the delay time, K is the number of reconstructed components, and xj is the components in the Jth column of the reconstruction matrix. The elements in each column of the reconstruction matrix Y are reorganized in ascending order, so that a group of symbol sequence Sl=j1,j2,⋯,jd can be obtained for each column of the matrix Y. The occurrence probability of each symbol sequence P1,P2,⋯,Pk can be calculated. At this time, the PE Hp of K different symbol sequences for the time series Xi can be defined in the form of Shannon entropy as Hpd=−∑j=1kPjlnPj.In the above definition, the value ofHp represents how time series is stochastic. A smaller value implies simpler time series, while a larger value leads to more complex and stochastic time series. Hence, minimum PE is taken as a fitness function in this paper and then solved to determine optimal decomposition parameters.


## 2.4. Savitzky-Golay Filter

SG filter, a method put forth by Savitzky and Golay, is widely applied in denoising the nonstationary signals containing non-Gaussian noise [33]. In the method, univariate P-order polynomials are adopted to fit the fixed length neighborhood of each data point in the selected data. Least squares criterion is followed to determine the polynomial coefficients by minimizing fitting error, so as to obtain the optimal fitted value of the data point, which is the value obtained after denoising. In this way, signals are denoised. At the time of denoising, the SG filter method can effectively retain the variation information of signals.
(13)min∑j=−mmYj−yj2,s.t.Yi=c0+c1i+c2i2+⋯+cpip.


## 3. The Proposed Vibration Signal Denoising Method

In this paper, a vibration signal denoising method combining improved VMD parameter optimization algorithm and SG filter is proposed. This method adopts the improved heuristic method GPSO and takes PE as the fitness function to automatically determine the numberK of modal components and the penalty factor α of VMD. The optimal parameters K and α are used to perform the VMD of noise-containing vibration signals and obtain a number of IMFs. Subsequently, PE is employed to calculate the portion of noise in decomposed signal components. The components are classified into noise and signal components by searching for the mutation points of two adjacent IMF permutation entropies. Noise components are denoised by virtue of SG filter. In the end, the denoised components are reconstructed with signal components to eventually obtain denoised vibration signals. Through decomposition and reconstruction, the main information of the signal is retained and a large amount of noise contained in the signal is eliminated. This method can be adaptively select optimization parameters and noise components. We do not directly abandon noise-containing components but optimize them to prevent overfiltering from causing signal distortion. The signals reconstructed after denoising contain more signal information. The denoising process of the proposed algorithm is given in Figure 2. The exact procedures of the proposed algorithm can be expressed as follows:Step 1.
Taking PE as the fitness function, GPSO-VMD is applied to decompose the vibration signal, andK and α at the minimum PE value are taken as the optimal decomposition values.Step 2.
Set the obtainedK and α as VMD parameters and decompose the vibration signal into K IMFs.Step 3.
Calculate the PE of IMF obtained by decomposition. The signal-dominated IMF and noise-dominated IMF are distinguished by the mutation point of PE of two adjacent IMF.Step 4.
Apply SG filter to denoise the noise component.Step 5.
Reconstruct the IMFs dominated by the signal and the IMF component after denoising to obtain the final denoised signal.Figure 2
Process flow of the proposed denoising method.


## 4. Simulation

### 4.1. Construction of Simulation Signals

The vibration signals measured by an atomic gravimeter in a ship environment are very complicated, nonlinear, and nonstationary. These signals actually contain lots of unpredictable disturbance noises. In order to verify the effectiveness of the proposed method, vibration simulation signals were designed with nonlinear and nonstationary features. The simulation signalsxt consisted of sinusoidal signal x1t, frequency-modulated signals x2t, amplitude-modulated signals x3t, exponentially decayed sinusoidal signals x4t, and other random noises with mean 0, standard deviation 1, and amplitude 1.4. Moreover, they had t=0,0.001 and sampling frequency 1000 Hz. These simulation signals are defined by Equation (14).
(14)x1t=0.6sin20πt+π/6,x2t=0.4cos60πt+sin10πt,x3t=1+0.3cos20πt⋅sin400πt,x4t=2e−5t⋅sin600πt,x=x1+x2+x3+x4+1.4⋅randn.The time-frequency domain waveform of each subsignal is given in Figure3, and the mixed signal is presented in Figure 4. As revealed in the frequency domain waveform of mixed signals, noises had higher power than signals, so that signals were submerged in a highly noisy environment. Moreover, noises were evenly distributed in the entire frequency domain of signals, which makes it very difficult to accurately extract feature signals.Figure 3
Time-frequency waveform of each subsignal in simulation signals.Figure 4
Time-frequency waveform of simulation signals.Prior to VMD, number of decomposed modesK and penalty term α should be properly selected. The improved GPSO was adopted to optimize the VMD algorithm. Thus, we set population size 50, crossover 0.8, and mutation probability 0.3. The number of decomposed modes K was set in the range [2, 10], while the penalty term α was set in the range [200, 3000]. Minimum PE was adopted as the fitness function. The fitness function convergence curve of GPSO is shown in Figure 5. Convergence was achieved at the time of the second iteration, when the optimal VMD parameters K and α were 5.1 and 2219.4, respectively, and rounded to 5 and 2219 since they must be integral. In order to prove its superiority, the GPSO was compared with GA and PSO algorithms, respectively. All optimization methods employed minimum PE as fitness function and had the same population size and maximum iteration times. The experiment used Windows 10 operating system, Intel Corei7-8750H and matlab2019a for simulation. The convergence results of the fitness function for these three algorithms are presented in Figure 5 and summed up in Table 1.Figure 5
Comparison of fitness function convergence curves.Table 1
Comparison of optimization results.
MethodBest parameter (α/K)Time (s)GA-VMD2077.6/5.52262PSO-VMD2359.8/4.41805GPSO-VMD2219.4/5.11627As shown in Figure5, the GPSO, GA, and PSO converged after the second, fifth, and seventh iteration, respectively. The GPSO had the lowest fitness after convergence and achieved the fastest convergence among them, which proves the strong global search capability and fast convergence rate.


### 4.2. Denoising Analysis of Simulation Signals

Noise simulation signals were decomposed in VMD based on theK and α values obtained by GPSO, so as to gain a number of IMF components as shown in Figure 6. Obviously, signals could be effectively decomposed in VMD. The decomposition results were greatly consistent with simulation signals, proving the efficacy of the parameters obtained by GPSO.Figure 6
Decomposition results of simulation signals by VMD.The IMF components obtained by VMD were arranged from low frequency to high frequency. Noises were mainly concentrated in high-frequency components, but there were still some effective IMF components. If they were simply abandoned, some effective signals would be lost. For this reason, effective components must be separated from noise-containing components, so as to denoise the latter. For this purpose, PE of each IMF component was calculated separately and used to determine the proportion of noise-containing components in the IMF components. This was achieved by calculating the mutation point in the PE of two adjacent components, that is,β=maxHPIMFl+1−HPIMFl with l for mutated IMF components. The noise-containing components were denoised by SG filter and then reconstructed with the effective IMF components to finally obtain denoised signals. A simulation signal xt was decomposed into five IMF components, whose permutation entropies were calculated separately.Based on the PE of each IMF component in Figure7, the difference between the permutation entropies of adjacent IMF components was 0.00154, 0.00271, 0.00252, and 0.00018, respectively. The difference between IMF2 and IMF3 was the largest. Hence, IMF3 was a mutation point, which helped identify IMF3-IMF5 as high-frequency noise-containing components. These high-frequency noise-containing components were treated by SG filter to obtain a signal. In this paper, the SG filter parameters are set as polynomial order 3 and data frame length 41. This signal was reconstructed with other signal components including IMF1 and IMF2 to obtain denoised vibration signal.Figure 7
PE of each IMF component.We also compared the proposed method with the classic wavelet denoising and empirical mode decomposition (EMD) to verify its effectiveness. Daubechies (db) wavelet and Symlets (sym) wavelet with good orthogonality in the wavelet transform denoising were selected to denoise simulation signals. The green and cyan signals in Figure8 show the results of wavelet transform denoising with five layers of soft threshold by db5 and sym6 wavelets, respectively. After analyzing these results, it was found that the denoising by db5 and sym6 wavelets might achieve the desired effect of denoising but filtered lots of useful high-frequency information, resulting in information distortion. The yellow signal in Figure 8 shows the result of denoising by EMD. In order to highlight the advancement of the proposed VMD algorithm, the denoising by EMD was performed in the same way as the proposed algorithm. A noise was first decomposed by EMD to obtain a number of IMF components. The PE of each component was calculated to find out the mutation point. The SG filtering was carried out for the IMF components in front of the mutation point. At last, the denoised signal was obtained through reconstruction. Based on the denoising results presented in Figure 8, EMD was more effective than db5 and symb6 wavelets in terms of denoising but still troubled by severe information distortion, causing the excessive loss of useful information. The denoising results of the proposed VMD algorithm are shown in the red signal in Figure 8. Obviously, the proposed algorithm could effectively retain the information of signals, achieve good denoising effect, and achieve high consistency regardless the denoised signal or the original signal.Figure 8
Simulation signals denoised by different algorithms.In order to quantitatively analyze the denoising effect of the proposed method, we took signal-to-noise ratio (SNR) and Mean Square Error (MSE) as the indicators to evaluate the denoising effect with denoised signals and noise-containing signals [34]. The greater SNR, the better effect of denoising. The lower MSE, the poorer effect of denoising. The SNR and MSE values are expressed by
(15)SNR=10lg∑i=1nxi2∑inyi−xi2,(16)MSE=∑i=1n1nxi−yi2,where xi indicates the signals that do not contain noise and yi means the denoised signals. The SNR and MSE values in the test of algorithms with simulation signals are presented in Table 2.Table 2
Denoising results with different algorithms.
MethodsSNR (dB)MSEOriginal signal-7.26902.0352db5-4.48652.1325sym6-4.48102.1436EMD-PE-SG-0.42771.1413The proposed method1.25170.7204As revealed in Table2, SNR was -7.2690 after adding random noise into a simulation signal xt. It was evident that the signal had lower power than noise, so that it was entirely submerged by the latter, making it very difficult to perform denoising. After being denoised by db5, sym6, EMD, and the proposed method, the proposed method achieved SNR 1.2517, which was the largest among them. Moreover, the proposed method realized the lowest MSE, revealing the most remarkable improvement by the proposed method. The proposed algorithm proves that denoising can be significantly achieved even in a highly noisy environment.


## 4.1. Construction of Simulation Signals

The vibration signals measured by an atomic gravimeter in a ship environment are very complicated, nonlinear, and nonstationary. These signals actually contain lots of unpredictable disturbance noises. In order to verify the effectiveness of the proposed method, vibration simulation signals were designed with nonlinear and nonstationary features. The simulation signalsxt consisted of sinusoidal signal x1t, frequency-modulated signals x2t, amplitude-modulated signals x3t, exponentially decayed sinusoidal signals x4t, and other random noises with mean 0, standard deviation 1, and amplitude 1.4. Moreover, they had t=0,0.001 and sampling frequency 1000 Hz. These simulation signals are defined by Equation (14).
(14)x1t=0.6sin20πt+π/6,x2t=0.4cos60πt+sin10πt,x3t=1+0.3cos20πt⋅sin400πt,x4t=2e−5t⋅sin600πt,x=x1+x2+x3+x4+1.4⋅randn.The time-frequency domain waveform of each subsignal is given in Figure3, and the mixed signal is presented in Figure 4. As revealed in the frequency domain waveform of mixed signals, noises had higher power than signals, so that signals were submerged in a highly noisy environment. Moreover, noises were evenly distributed in the entire frequency domain of signals, which makes it very difficult to accurately extract feature signals.Figure 3
Time-frequency waveform of each subsignal in simulation signals.Figure 4
Time-frequency waveform of simulation signals.Prior to VMD, number of decomposed modesK and penalty term α should be properly selected. The improved GPSO was adopted to optimize the VMD algorithm. Thus, we set population size 50, crossover 0.8, and mutation probability 0.3. The number of decomposed modes K was set in the range [2, 10], while the penalty term α was set in the range [200, 3000]. Minimum PE was adopted as the fitness function. The fitness function convergence curve of GPSO is shown in Figure 5. Convergence was achieved at the time of the second iteration, when the optimal VMD parameters K and α were 5.1 and 2219.4, respectively, and rounded to 5 and 2219 since they must be integral. In order to prove its superiority, the GPSO was compared with GA and PSO algorithms, respectively. All optimization methods employed minimum PE as fitness function and had the same population size and maximum iteration times. The experiment used Windows 10 operating system, Intel Corei7-8750H and matlab2019a for simulation. The convergence results of the fitness function for these three algorithms are presented in Figure 5 and summed up in Table 1.Figure 5
Comparison of fitness function convergence curves.Table 1
Comparison of optimization results.
MethodBest parameter (α/K)Time (s)GA-VMD2077.6/5.52262PSO-VMD2359.8/4.41805GPSO-VMD2219.4/5.11627As shown in Figure5, the GPSO, GA, and PSO converged after the second, fifth, and seventh iteration, respectively. The GPSO had the lowest fitness after convergence and achieved the fastest convergence among them, which proves the strong global search capability and fast convergence rate.


## 4.2. Denoising Analysis of Simulation Signals

Noise simulation signals were decomposed in VMD based on theK and α values obtained by GPSO, so as to gain a number of IMF components as shown in Figure 6. Obviously, signals could be effectively decomposed in VMD. The decomposition results were greatly consistent with simulation signals, proving the efficacy of the parameters obtained by GPSO.Figure 6
Decomposition results of simulation signals by VMD.The IMF components obtained by VMD were arranged from low frequency to high frequency. Noises were mainly concentrated in high-frequency components, but there were still some effective IMF components. If they were simply abandoned, some effective signals would be lost. For this reason, effective components must be separated from noise-containing components, so as to denoise the latter. For this purpose, PE of each IMF component was calculated separately and used to determine the proportion of noise-containing components in the IMF components. This was achieved by calculating the mutation point in the PE of two adjacent components, that is,β=maxHPIMFl+1−HPIMFl with l for mutated IMF components. The noise-containing components were denoised by SG filter and then reconstructed with the effective IMF components to finally obtain denoised signals. A simulation signal xt was decomposed into five IMF components, whose permutation entropies were calculated separately.Based on the PE of each IMF component in Figure7, the difference between the permutation entropies of adjacent IMF components was 0.00154, 0.00271, 0.00252, and 0.00018, respectively. The difference between IMF2 and IMF3 was the largest. Hence, IMF3 was a mutation point, which helped identify IMF3-IMF5 as high-frequency noise-containing components. These high-frequency noise-containing components were treated by SG filter to obtain a signal. In this paper, the SG filter parameters are set as polynomial order 3 and data frame length 41. This signal was reconstructed with other signal components including IMF1 and IMF2 to obtain denoised vibration signal.Figure 7
PE of each IMF component.We also compared the proposed method with the classic wavelet denoising and empirical mode decomposition (EMD) to verify its effectiveness. Daubechies (db) wavelet and Symlets (sym) wavelet with good orthogonality in the wavelet transform denoising were selected to denoise simulation signals. The green and cyan signals in Figure8 show the results of wavelet transform denoising with five layers of soft threshold by db5 and sym6 wavelets, respectively. After analyzing these results, it was found that the denoising by db5 and sym6 wavelets might achieve the desired effect of denoising but filtered lots of useful high-frequency information, resulting in information distortion. The yellow signal in Figure 8 shows the result of denoising by EMD. In order to highlight the advancement of the proposed VMD algorithm, the denoising by EMD was performed in the same way as the proposed algorithm. A noise was first decomposed by EMD to obtain a number of IMF components. The PE of each component was calculated to find out the mutation point. The SG filtering was carried out for the IMF components in front of the mutation point. At last, the denoised signal was obtained through reconstruction. Based on the denoising results presented in Figure 8, EMD was more effective than db5 and symb6 wavelets in terms of denoising but still troubled by severe information distortion, causing the excessive loss of useful information. The denoising results of the proposed VMD algorithm are shown in the red signal in Figure 8. Obviously, the proposed algorithm could effectively retain the information of signals, achieve good denoising effect, and achieve high consistency regardless the denoised signal or the original signal.Figure 8
Simulation signals denoised by different algorithms.In order to quantitatively analyze the denoising effect of the proposed method, we took signal-to-noise ratio (SNR) and Mean Square Error (MSE) as the indicators to evaluate the denoising effect with denoised signals and noise-containing signals [34]. The greater SNR, the better effect of denoising. The lower MSE, the poorer effect of denoising. The SNR and MSE values are expressed by
(15)SNR=10lg∑i=1nxi2∑inyi−xi2,(16)MSE=∑i=1n1nxi−yi2,where xi indicates the signals that do not contain noise and yi means the denoised signals. The SNR and MSE values in the test of algorithms with simulation signals are presented in Table 2.Table 2
Denoising results with different algorithms.
MethodsSNR (dB)MSEOriginal signal-7.26902.0352db5-4.48652.1325sym6-4.48102.1436EMD-PE-SG-0.42771.1413The proposed method1.25170.7204As revealed in Table2, SNR was -7.2690 after adding random noise into a simulation signal xt. It was evident that the signal had lower power than noise, so that it was entirely submerged by the latter, making it very difficult to perform denoising. After being denoised by db5, sym6, EMD, and the proposed method, the proposed method achieved SNR 1.2517, which was the largest among them. Moreover, the proposed method realized the lowest MSE, revealing the most remarkable improvement by the proposed method. The proposed algorithm proves that denoising can be significantly achieved even in a highly noisy environment.


## 5. Vibration Signal Denoising of Marine Atomic Gravimeter

### 5.1. Data Collection

The vibration data of an atomic gravimeter was collected by navigation test. The test platform was composed of an atomic gravimeter, an inertial stabilized platform, and a vibration measuring device. The collection of vibration data was performed by a collecting unit formed by a data collector, an accelerometer, and a computer. The accelerometer was attached to the atomic gravimeter placed on the inertial stabilized platform. The accelerometer converted the vibration information into analog voltage output. The data collector performed the analog-to-digital conversion of analog signals collected by sensors and transferred the digital information to the computer. Data collection software was installed in the computer to store and process the received digital information. We selected the vibration signalxsignal and took the data from 1000 sampling points for analysis, as shown in Figure 9.Figure 9
Time domain waveform of vibration signals of a marine atomic gravimeter.


### 5.2. Implementation of Denoising

The GPSO was employed to optimize the VMD parameters in the proposed method, so as to obtain the optimal number of decomposed modesK and penalty term α. We set population size 50, crossover probability 0.8, mutation probability 0.3, number of decomposed modes in [2, 10], and penalty term in [200, 3000]. PE was used as the fitness function. At the third iteration, the fitness function has completed convergence, and the minimum PE value is 0.69137. The corresponding optimization results were [6, 2301]. The convergence curve of the fitness function is shown in Figure 10.Figure 10
Convergence curve of vibration signal fitness function.Subsequently, the parameters optimized by GPSO were used in the decomposition by VMD to obtain six IMF components. Based on these components, it was found that noise signals were concentrated in IMF5 and IMF6.The PE of each IMF component was calculated. Then the maximum difference between the permutation entropies of two adjacent IMF components was calculated to find out a mutation point and classify IMF components into effective signal IMF components and noise-containing IMF components as shown in Figure11. As revealed in the slope of PE in Figure 12, the mutation point of PE for the signal xsignal was IMF, so that the noise components of the signal xsignal were IMF5-IMF6. The SG filtering was conducted for noise components. The filtered IMF components were reconstructed with the effective IMF components to obtain the denoised vibration signal as shown in Figure 13.Figure 11
Decomposition by VMD.Figure 12
PE of each IMF for the vibration signal.Figure 13
Comparison of denoised vibration signals.The vibration signal of a marine atomic gravimeter was measured during navigation, so that it was impossible to obtain a noise-free original vibration signal. For this reason, SNR and MSE could not be used to quantitatively analyze and compare the proposed method with other algorithms. Figure13 shows the results of wavelet transform denoising by five layers of soft threshold with db5 and sym6 wavelets and the effects of denoising by EMD-PE-SG, respectively. As revealed in Figure 13, EMD-PE-SG achieved better denoising than db5 and sym6 wavelets and could retain more actual information. However, it still lost too much useful information compared with the proposed method, so it resulted in severe loss and distortion of signal information.


## 5.1. Data Collection

The vibration data of an atomic gravimeter was collected by navigation test. The test platform was composed of an atomic gravimeter, an inertial stabilized platform, and a vibration measuring device. The collection of vibration data was performed by a collecting unit formed by a data collector, an accelerometer, and a computer. The accelerometer was attached to the atomic gravimeter placed on the inertial stabilized platform. The accelerometer converted the vibration information into analog voltage output. The data collector performed the analog-to-digital conversion of analog signals collected by sensors and transferred the digital information to the computer. Data collection software was installed in the computer to store and process the received digital information. We selected the vibration signalxsignal and took the data from 1000 sampling points for analysis, as shown in Figure 9.Figure 9
Time domain waveform of vibration signals of a marine atomic gravimeter.


## 5.2. Implementation of Denoising

The GPSO was employed to optimize the VMD parameters in the proposed method, so as to obtain the optimal number of decomposed modesK and penalty term α. We set population size 50, crossover probability 0.8, mutation probability 0.3, number of decomposed modes in [2, 10], and penalty term in [200, 3000]. PE was used as the fitness function. At the third iteration, the fitness function has completed convergence, and the minimum PE value is 0.69137. The corresponding optimization results were [6, 2301]. The convergence curve of the fitness function is shown in Figure 10.Figure 10
Convergence curve of vibration signal fitness function.Subsequently, the parameters optimized by GPSO were used in the decomposition by VMD to obtain six IMF components. Based on these components, it was found that noise signals were concentrated in IMF5 and IMF6.The PE of each IMF component was calculated. Then the maximum difference between the permutation entropies of two adjacent IMF components was calculated to find out a mutation point and classify IMF components into effective signal IMF components and noise-containing IMF components as shown in Figure11. As revealed in the slope of PE in Figure 12, the mutation point of PE for the signal xsignal was IMF, so that the noise components of the signal xsignal were IMF5-IMF6. The SG filtering was conducted for noise components. The filtered IMF components were reconstructed with the effective IMF components to obtain the denoised vibration signal as shown in Figure 13.Figure 11
Decomposition by VMD.Figure 12
PE of each IMF for the vibration signal.Figure 13
Comparison of denoised vibration signals.The vibration signal of a marine atomic gravimeter was measured during navigation, so that it was impossible to obtain a noise-free original vibration signal. For this reason, SNR and MSE could not be used to quantitatively analyze and compare the proposed method with other algorithms. Figure13 shows the results of wavelet transform denoising by five layers of soft threshold with db5 and sym6 wavelets and the effects of denoising by EMD-PE-SG, respectively. As revealed in Figure 13, EMD-PE-SG achieved better denoising than db5 and sym6 wavelets and could retain more actual information. However, it still lost too much useful information compared with the proposed method, so it resulted in severe loss and distortion of signal information.


## 6. Conclusion

A vibration signal denoising method based on improved VMD is put forward in this paper. In this method, an improved GPSO based on PSO and GA is first adopted for the parametric optimization of VMD. Minimum PE is taken as the fitness function to adaptively search for the optimal parametersK and α in VMD. Based on the obtained parameters K and α, a noise-contained signal is decomposed into a number of IMF components. PE is utilized to calculate the proportion of noise-containing components in the signal components obtained by decomposition. A mutation point is found with PE of adjacent IMF components to classify these components into noise and signal components. The SG filter is carried out to denoise these noise components. At last, the denoised components are reconstructed with signal components to generate the denoised vibration signal. The reconstructed signal contains more physical information. The proposed method makes use of signal mode decomposition to adaptively extract noise but does not use any fixed priori threshold. In order to demonstrate its effectiveness, the proposed method is applied in denoising with the vibration data collected by vibration and measured by a marine atomic gravimeter. As proved in the test, the proposed method can effectively separate noise from vibration signals and achieve great denoising. Its potential has been demonstrated in filtering noise and improving the quality of vibration data to provide the supporting data for the vibration compensation of a marine atomic gravimeter.



---

*Source: 1744689-2022-03-27.xml*