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.
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 and penalty term must be artificially set in the VMD of signals. If the value of 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 the value based on kurtosis [27] and energy factor [28]. In these methods, only the number of decomposed modes 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 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 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 into IMFs with different center frequencies . The variational model is correspondingly described by Equation (1). where is IMF components after decomposition, is the center frequency of each component, is the unit pulse function, is the imaginary unit, is the convolution operation, and is the partial derivative with respect to . Lagrangian multiplier and secondary penalty term are introduced to obtain the optimal solution of constrained variation as follows:
The alternation of multiplication operators is conducted to iteratively update and obtain the optimal mode component , center frequency , and Lagrangian multiplier . where is the Fourier transform, is the fidelity coefficient, and is the iteration times. Iteration is constantly updated but halted when relative error is less than convergence precision. The decomposition result is eventually outputted.
2.2. GPSO Optimized VMD
The effect of VMD depends significantly on the number of decomposed modes 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. Generate particles randomly with , where and are the geometrical location and velocity vector. So the initial generation of particle swarm is identified as(2)VMD Decomposition. Calculate the fitness of individuals in the initial population, and select the individual’s optimal value and the global optimal value 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 particle , we letwhere and are the random numbers in (0, 1), is the inertia weight, and and are the learning factors. Therefore, the th generation of particle swarm is formed as follows:
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 value to the minimum value . It varies with the iteration times of the algorithm as follows: where is the current iteration steps, is the maximum iteration steps, and and 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 , selection, crossover, and mutation operations are conducted to obtain its next-generation population , 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 ( and ) of the population based on the fitness, and calculate the optimal particles that have been found so far for each particle . Calculate the optimal particles that have been found so far for the current population .(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 to obtain a matrix : where is the embedding dimensions, is the delay time, is the number of reconstructed components, and is the components in the th column of the reconstruction matrix. The elements in each column of the reconstruction matrix are reorganized in ascending order, so that a group of symbol sequence can be obtained for each column of the matrix . The occurrence probability of each symbol sequence can be calculated. At this time, the PE of different symbol sequences for the time series can be defined in the form of Shannon entropy as .
In the above definition, the value of 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.
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 number of modal components and the penalty factor of VMD. The optimal parameters 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, and and at the minimum PE value are taken as the optimal decomposition values.
Step 2. Set the obtained and as VMD parameters and decompose the vibration signal into 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.
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 signals consisted of sinusoidal signal , frequency-modulated signals , amplitude-modulated signals , exponentially decayed sinusoidal signals , and other random noises with mean 0, standard deviation 1, and amplitude 1.4. Moreover, they had and sampling frequency 1000 Hz. These simulation signals are defined by Equation (14).
The time-frequency domain waveform of each subsignal is given in Figure 3, 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.
Prior to VMD, number of decomposed modes 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 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 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.
As shown in Figure 5, 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 the 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.
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, with 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 was decomposed into five IMF components, whose permutation entropies were calculated separately.
Based on the PE of each IMF component in Figure 7, 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.
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 Figure 8 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.
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 where indicates the signals that do not contain noise and means the denoised signals. The SNR and MSE values in the test of algorithms with simulation signals are presented in Table 2.
As revealed in Table 2, SNR was -7.2690 after adding random noise into a simulation signal . 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 signal and took the data from 1000 sampling points for analysis, as shown in Figure 9.
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 modes 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.
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 Figure 11. As revealed in the slope of PE in Figure 12, the mutation point of PE for the signal was IMF, so that the noise components of the signal 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.
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. Figure 13 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 parameters and in VMD. Based on the obtained parameters 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.
Data Availability
No data were used to support the findings of the study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was funded by the National Natural Science Foundation of China under Grants (61873275), the Foundation of Basic Strengthening Technology of the Military Science and Technology Commission (2019JCJQJJ047), and the Natural Science Foundation of Hubei Provincial of China (2017CFB377).