A Novel Noise Reduction Method of UAV Magnetic Survey Data Based on CEEMDAN, Permutation Entropy, Correlation Coefficient and Wavelet Threshold Denoising

Despite the increased attention that has been given to the unmanned aerial vehicle (UAV)-based magnetic survey systems in the past decade, the processing of UAV magnetic data is still a tough task. In this paper, we propose a novel noise reduction method of UAV magnetic data based on complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN), permutation entropy (PE), correlation coefficient and wavelet threshold denoising. The original signal is first decomposed into several intrinsic mode functions (IMFs) by CEEMDAN, and the PE of each IMF is calculated. Second, IMFs are divided into four categories according to the quartiles of PE, namely, noise IMFs, noise-dominant IMFs, signal-dominant IMFs, and signal IMFs. Then the noise IMFs are removed, and correlation coefficients are used to identify the real signal-dominant IMFs. Finally, the wavelet threshold denoising is applied to the real signal-dominant IMFs, the denoised signal can be obtained by combining the signal IMFs and the denoised IMFs. Both synthetic and field experiments are conducted to verify the effectiveness of the proposed method. The results show that the proposed method can eliminate the interference to a great extent, which lays a foundation for the further interpretation of UAV magnetic data.


Introduction
The past decade has seen a variety of applications conducted by unmanned aerial vehicles (UAVs) in many fields, e.g., archaeology, remote sensing, geological prospecting, and unexploded ordnance (UXO) detection [1]. Among these applications, the use of UAVs for magnetic surveys is a booming branch of research [2]. UAV magnetic surveys can cover a wider range with a higher efficiency compared with the traditional terrestrial magnetic surveys, and are also easy to operate, have a low-cost, and have a good safety profile compared with manned aircraft magnetic surveys [3]. In addition, UAV-based magnetic surveys can also be carried out in areas that are difficult to access or that would pose a potential hazard to operators (e.g., near active volcanoes), which means that some gaps in traditional magnetic surveys can now be studied [4,5].
A substantial body of research has accumulated on the integration of UAV magnetic survey systems; however, the processing of UAV magnetic data remains an open problem. Several attempts have been made to process UAV magnetic data, e.g., Malehmir et al. [6] used a median filter to process the spiky sample points, while data with severe noise were excluded. However, the quality of the collected data has not been evaluated. Walter et al. [7] investigated the periodic variations caused by the swing of magnetometers, spectral analysis and lowpass filter were applied to identify and remove the periodic signal, yet the target signal may be removed if its frequency overlaps with the periodic signal. Mu et al. [8] proposed a lowpass filter to remove the interference field; however, the cutoff frequency and the order of the filter need to be determined according to a priori knowledge. Similar methods were also used to process the multi-rotor UAV magnetic data, as described in [9,10]. Liu et al. [11] proposed an adaptive cancellation of geomagnetic background noise for magnetic anomaly detection. The system is regarded as a two-channel linear time-invariant (LTI) system, where the first sensor records the background noise as a reference, and the second sensor records the target signal and background noise at the same time. It should be noted that if the system is a single channel (i.e., a magnetometer), it will be difficult to use this method. Wang et al. [12] used higher-order statistics to suppress the interference of Gaussian colored noise in magnetotelluric data. However, if the noise is complex and no longer obeys Gaussian distribution, this method will not be able to suppress interference effectively. In addition, this method requires redundant data, which limits its further implementation. Overall, processing methods alone are not enough, since the collected data are usually non-stationary and easily affected by noise from many sources, e.g., interference generated by UAV platform, geological noise, industrial frequency interference, and instrument noise.
Empirical mode decomposition (EMD), proposed by Huang [13], is an adaptive timefrequency analysis method which is suitable for non-linear and non-stationary signals. EMD analyzes the signal according to the characteristics of the signal itself and does not need a basis function [14]. However, the boundary effects and mode mixing heavily impact the effect of EMD. To overwhelm these problems, a noise-assistant analysis method, i.e., ensemble empirical mode decomposition (EEMD), is proposed [15]. EEMD basically overcomes the mode mixing; however, two new problems, the difference in intrinsic mode function (IMF) numbers and the introduction of extra noise, have arisen. A complete ensemble empirical mode decomposition with adaptive white noise (CEEMDAN) is proposed to surmount these obstacles [16]. The CEEMDAN method can significantly reduce the reconstruction error and requires fewer iterations compared with EMD and EEMD. To date, CEEMDAN has been widely used in the field of non-linear and non-stationary signal processing, e.g., biological signal processing [17,18], wind speed forecasting [19,20], financial time series forecasting [21], gear fault diagnosis [22][23][24], underwater acoustic signal denoising [25,26], and structural damage localization and quantification [27].
To further study the characteristics of non-linear and non-stationary signals, permutation entropy (PE) and correlation coefficient (CC) are proposed to evaluate the complexity of obtained IMFs, and identify whether the IMFs require denoising, as noted in several previous studies [22,25,28,29]. In addition, wavelet threshold denoising is adopted as part of the combined method [17,25,26]. To the best of our knowledge, there have been no previous studies on UAV magnetic data denoising based on CEEMDAN. Moreover, the determination of noisy IMFs is generally based on the artificial threshold, which is not only difficult to achieve in practice, but also does not make full use of the characteristics of the signal itself.
In this paper, a novel noise reduction method for UAV magnetic data is proposed by taking advantage of CEEMDAN, PE, CC, and wavelet threshold denoising. The main contributions of the proposed method are as follows: 1.
The adaptive decomposition algorithm, i.e., CEEMDAN, is applied to multi-rotor UAV magnetic data for the first time. The original data are decomposed into a set of IMF components with different scales.

2.
The IMFs are divided into four categories, i.e., noise IMFs, noise-dominant IMFs, signal-dominant IMFs, and signal IMFs according to the quartiles of PE, which is completely determined by the characteristics of the signal itself without setting a threshold artificially.

3.
The real signal-dominant IMFs are identified using CC, while the wavelet soft threshold denoising (WSTD) is applied to further suppress noise. Simulation results show that the signal-to-noise ratio (SNR) of the signal can be improved by about 16-20 dB after denoising by means of the proposed method.
This paper is organized as follows: Section 2 presents the relevant principles of CEEMDAN, PE, CC, and wavelet threshold denoising; the proposed noise reduction algorithm for UAV magnetic data is presented in Section 3; in Sections 4 and 5, the proposed method is applied to both synthetic and real UAV magnetic data, respectively; Section 6 contains the conclusion of this paper.

Principles of EMD, EEMD, and CEEMDAN Algorithm
In this section, the mathematical principles of EMD, EEMD, and CEEMDAN are introduced, and the specific implementation steps, diagrams, and pseudocodes are given to better understand how these methods work.

EMD Algorithm
EMD can adaptively decompose the original signal s(t) into several IMFs and a residue. The IMF meets the following two conditions: (a) the number of extremum points and zero-crossing points must be equal or not exceed one, and (b) the mean value of the upper envelope formed by the local maximum points and the lower envelope formed by the local minimum points is zero. The procedure of EMD is summarized as follows [30]: Step 1: Identify all the extremum points of the original signal s(t) and define the upper and lower envelope u(t) and l(t), respectively, using a cubic spline interpolation.
Step 2: Calculate the mean envelope of the upper and lower envelope.
Step 3: Subtract the mean envelope from s(t) to obtain the first intermediate signal.
Step 4: If I 1 (t) satisfies the criteria of the IMF, then define I MF 1 = I 1 (t), otherwise treat I 1 (t) as the new signal and repeat the above procedure k times until I k+1 (t) satisfies the IMF conditions. The acquisition of IMF usually requires several iterations. To finish the iteration, the stopping criterion is defined as follows: where n is the length of the intermediate signal. The iteration will be stopped when S D < δ (in this study δ is set to 0.3).
Step 5: Let r 1 (t) = s(t) − I MF 1 , treat r 1 (t) as the new signal and repeat Step 1-4 to obtain the next IMF, until r N (t) becomes either a constant or a monotonic function. Finally, the original signal s(t) after EMD can be expressed as: where N is the number of IMFs, and r N (t) is the final residue. The pseudocode of EMD is described in detail in Algorithm 1. Figure 1 is the schematic diagram of the main steps of EMD.

Input:
The original signal x. Output: Several IMF and a residue, i.e., I MF i , (i = 1, 2, . . . , n) and r. 1: function EMD (x, ResidueThreshold, S DT ) 2: while residue> ResidueThreshold do 7: i←i+1 8: S D ← ∞ 10: while S D > S DT do 11: for j = 1→j = N do 12: [ end while 24: return IMF and residue 25: end function  The effect of EMD is easily affected by mode mixing, i.e., signals of different feature scales appear in the same IMF, or signals with the same feature scale are dispersed into different IMFs [31]. This problem not only decreases the decomposition efficiency but also degrades the subsequent denoising performance.

EEMD Algorithm
The EEMD algorithm is proposed to eliminate the mode mixing of EMD. The procedure of EEMD is summarized as follows [17,26]: Step 1: Different Gaussian white noise signals with zero mean and unit variance n i (t) are added to the original signal s(t) to obtain a set of new signals Step 2: Decompose each s i (t) by EMD to obtain I MF ik , where k = 1, 2, . . . , N denotes the number of IMFs.
Step 3: Average the I MF ik to obtain the EEMD mode The pseudocode of EEMD is described in detail in Algorithm 2.

Input:
The original signal x, the amplitude of the added Gaussian noise σ, and the number of ensemble trials m.

CEEMDAN Algorithm
Since the number of ensemble average is finite, a reconstruction error still exists in the result of EEMD. CEEMDAN can effectively overcome the mode mixing, with the reconstruction error and computational cost significantly reduced. The procedure of CEEMDAN is summarized as follows [25,26,32]: Step 1: The white noise ε 0 n i (t) is added to the original signal s(t), and the first IMF of CEEMDAN is obtained by calculating the ensemble average: where E n ( * ) is defined as the nth mode component of EMD.
Step 2: The first residual component can be obtained Step 3: Construct the new signal and decompose it by EMD. The second mode component can be obtained: Step 4: The nth residual signal and the (n+1)th IMF can be obtained according to the process of Step 3 Step 5: Repeat Step 4 until the residual signal is no longer decomposed. The original signal can be expressed as where K is the number of IMFs by CEEMDAN, and r(t) is the final residual mode.
The pseudocode of CEEMDAN is described in detail in Algorithm 3.

Input:
The original signal x, the amplitude of the added Gaussian noise σ, and the number of ensemble trials m. Output: Several IMF and a residue, i.e., I MF k , (k = 1, 2, . . . , p) and r. 1: n i ← N(σ * std(x), 1) 6: x i ← x + n i 7: end for 10: The flow chart of CEEMDAN is shown in Figure 3.

Permutation Entropy
PE was initially introduced by Bandt and Pompe [33] as a tool for measuring the complexity of time series, the advantages of PE are its simplicity, fast calculation, better robustness, and strong anti-noise ability, which make it suitable for the feature extraction of non-linear data. The specific steps of PE are summarized as follows [23,29]: Step 1: The first step in the calculation of permutation entropy requires extracting ordinal information from the time series. Given a time series X = {x 1 , x 2 , . . . , x N }, K reconstructed time series can be obtained as: where m and τ represent the embedding dimension and time delay, respectively.

Correlation Coefficient
The correlation coefficient is a dimensionless index that is widely applied in multivariate statistics to represent the relationship between two groups of variables [28]. Its value ranges from −1 to 1. The larger the absolute value of the correlation coefficient is, the stronger the correlation between the two variables is. For the two groups of variables x and y, the correlation coefficient ρ xy is defined as follows [35] ρ xy = cov(x, y) cov(x, x)·cov(y, y) , where cov(x, y) is the covariance of x and y, cov(x, x) and cov(y, y) are the variance of x and y, respectively. Therefore, the correlation coefficient can be expressed as

Wavelet Threshold Denoising
The wavelet transform is an effective time-frequency analysis tool which has the characteristic of multi-resolution and has been widely used for signal processing [36,37]. For a one-dimensional noisy signal, where s(t), f (t), and e(t) are defined as the noisy signal, real signal, and Gaussian noise signal, respectively. The specific steps of wavelet-based denoising are as follows [25]: Step 1: Proper wavelet basis function and decomposition level are selected to conduct wavelet decomposition on the noisy signal s(t).
Step 2: The thresholds are estimated according to appropriate threshold selection criteria for the high-frequency coefficients at different decomposition scales.
Step 3: The low-frequency coefficients of decomposition and the threshold high-frequency coefficients are used to reconstruct signals.
The key to wavelet threshold denoising is the selection of the wavelet basis function and the threshold function. The db4 wavelet basis function and a soft threshold method are selected in this paper.

The Proposed Method for UAV Magnetic Data Denoising
A denoising algorithm for UAV magnetic survey data based on CEEMDAN, PE, CC, and WSTD is proposed in this paper. The proposed method is based on the premise that there is a significant difference in complexity between the target signal and noise, so PE can be used to measure whether the IMF is dominated by target signal or noise. To avoid the influence of unreasonable threshold setting on subsequent processing, IMFs are divided into four categories by the quartiles of PE. The real signal-dominant IMFs are further confirmed by CC, and the noise is further suppressed by WSTD. The flow chart of the proposed method is shown in Figure 4. The specific procedures are summarized as follows: Step 1: The original signal s(t) is decomposed into several IMFs by CEEMDAN and arranged from high frequency to low frequency.
Step 2: Calculate the PE of all IMFs, the PE sequence is arranged in ascending order, and the extremum and the quartiles of the PE sequence are found, namely MI N pe , Q1, Q2, Q3, and MAX pe . Step 4: The CCs between signal (noise)-dominant IMFs and the real signal which is constituted by signal IMFs are obtained. The median of the absolute value sequence of CCs is recorded as M cc , and the IMFs corresponding to the CC greater than M cc are defined as the real signal-dominant IMFs.
Step 5: WSTD is applied to the real signal-dominant IMFs. The wavelet basis function and the decomposition level are db4 and 4, respectively.
Step 6: The denoised signal can be obtained by combining the signal IMFs and the denoised IMFs.

Acquisition of Synthetic Signal
Generally, the target signal of UAV magnetic surveys can be considered as a magnetic dipole, since the distance between sensors and the target is usually 2.5 times greater than the maximum dimension of the target [38]. The magnetic anomaly field generated by the target can be calculated by: where µ 0 is the permeability in vacuum, m is the dipole moment of the target, r is the displacement vector from the target to the measurement point, and r = |r|. As illustrated in Figure 5, the target center is located 1.  The projection of the magnetic field generated by the target in the direction of the geomagnetic field is the real signal; in addition, the geomagnetic field, UAV's interference field, equipment noise, and power frequency interference constitute the actual signal. Gaussian white noise with different SNRs is added to the real signal as simulated synthetic signals.

Evaluation of Different Denoising Methods
To clearly verify the denoising effectiveness of the proposed method, two combined noise reduction methods, i.e., EMD-PE-WSTD and EEMD-PE-WSTD, are chosen to compare with the proposed CEEMDAN-PE-CC-WSTD method. The former two methods use EMD and EEMD to decompose the original signal into several IMFs, then the noise IMFs are identified and removed using the quartile of the calculated PE. Both noise-dominant and signal-dominant IMFs are denoised using WSTD. Finally, the processed signal is obtained by combining the signal IMFs and the denoised IMFs. For EEMD and CEEMDAN-based methods, the amplitude of the added noise and the number of ensemble trials are 0.12 and 50, respectively. For the synthetic signal with a SNR of −10 dB, the decomposition results using EMD, EEMD, and CEEMDAN are shown in Figure 6. For each set of IMFs, results of PE can be obtained, as shown in Table 1.   For each decomposition method, IMFs are divided into four categories according to the corresponding quartiles of PEs, as shown in Table 2. Noise IMFs are first identified and removed. For the EMD and EEMD-based method, both noise-dominant and signaldominant IMFs (IMF3-IMF7) are denoised using WSTD. For the CEEMDAN-based method, signal IMFs constitute the real signal, and CCs of the remaining IMFs and the real signal are obtained, as shown in Table 3. The median of the absolute value of CC sequence is 0.0232, and the IMFs corresponding to the CC greater than this value are selected as the real signal-dominant IMFs, e.g., IMF7 and IMF8 are selected in this case. Then, the real signal-dominant IMFs are denoised by WSTD, and the denoised signal can be obtained by combining the denoised IMFs and the signal IMFs. Table 2. Four categories of IMFs obtained by EMD, EEMD, and CEEMDAN methods. EMD  IMF1, IMF2  IMF3, IMF4, IMF5  IMF6, IMF7  IMF8, IMF9   EEMD  IMF1, IMF2  IMF3, IMF4, IMF5  IMF6, IMF7  IMF8, IMF9   CEEMDAN  IMF1, IMF2, IMF3  IMF4, IMF5  IMF6, IMF7, IMF8  IMF9, IMF10   Table 3. Correlation coefficients between each IMF and the real signal of the CEEMDANbased method. To clearly compare different methods, SNR and root mean square error (RMSE) are used to evaluate the denoising performance. SNR shows an energy relationship between signal and noise, which is an intuitive method to evaluate the effect of the denoised signal by analyzing whether the SNR is improved. RMSE shows the difference between the denoised signal and the real target signal; the smaller the RMSE, the better the denoising effect. The formulas of SNR and RMSE are, respectively, given as follows

Mode
where s(i) is the original signal,ŝ(i) is the denoised signal, and n is the number of sampling points. Figure 7 shows the synthetic signal with −10 dB SNR and the denoised signal based on EMD-PE-WSTD, EEMD-PE-WSTD, and CEEMDAN-PE-CC-WSTD methods. The denoised results of synthetic signals with the SNR of −10 dB, −5 dB, 0 dB, and 5 dB are shown in Table 4. As can be seen from Figure 7 and Table 4, the three denoising methods all can reduce noise, and the proposed CEEMDAN-PE-CC-WSTD method has a better performance than the other two methods. There are two main reasons to explain these results: (1) all three methods can effectively suppress noise, indicating that the quartile of PE can indeed classify IMFs into four categories according to the dominance of signal and noise, and (2) the real signal-dominant IMFs are further confirmed by the median of CC, and hence the noise in the original signal is further suppressed. It should be noted that both the classification of IMFs and the confirmation of real signal-dominant IMFs are realized through the characteristics of the signal itself, and there is no need to set a threshold artificially. The use of the quartile of PE and the median of CC makes this method completely adaptive for the analysis of noisy signals.    The WSTD alone is applied to the synthetic data with different SNRs, using two different wavelet basis functions, db4 and sym4, respectively. The decomposition level is from 3 to 7. WSTD results are shown in Table 5. It can be seen that the effect of WSTD depends heavily on the selection of wavelet basis function and the number of decomposition levels. In addition, SNR of denoised signal reaches a maximum when the decomposition level increases to a certain value. The best result of WSTD reaches the level of the EMD-PE-WSTD and EEMD-PE-WSTD method. However, it is worth noting that for the method in Table 4, we do not optimize the wavelet basis function and the number of decomposition level. Therefore, as suggested by Tables 4 and 5, it is expected that our proposed method has a better performance than the WSTD method.

The Multi-Rotor UAV Magnetic Survey System
A multi-rotor-based UAV-magnetometer system was deployed for the purpose of nearsurface targets detection, where the parameters of the system can be found in [8]. However, the magnetometers were semi-rigidly mounted below the UAV, which could not meet the requirements of vertical take-off and landing (VTOL), and the potential impact risk of magnetometers also existed [8,39]. To surmount these obstacles, a magnetic survey system based on a six-rotor UAV was developed, as shown in Figure 8. This system consisted of two cesium optically pumped magnetometers (OPMs) and a fluxgate magnetometer to record the total magnetic intensity (TMI) data and the vector magnetic intensity (VMI) data, a differential GPS to provide the location information of the UAV, a data acquisition module, and a power module. The two OPMs were rigidly mounted below the center of the UAV by a boom, with a vertical distance of 0.45 m. The TMI and VMI data were synchronized by the pulse per second signal, with a sampling frequency of 160 Hz. The technical specifications of the multi-rotor magnetic survey system are given in Table 6.

Evaluation of UAV-Borne Magnetic Survey Results after Denoising
Experiments were carried out in Sichuan, China for the detection of near-surface buried targets. A segment pipeline made of steel was used as the target with a buried depth of 3 m, a length of 2 m, and a diameter of 0.15 m. A 12 m × 16 m rectangular area was selected as the survey area, and the target was buried near the center of the survey area. The pre-programmed flight profiles ran along the north-south direction, with a line spacing of 0.5 m. Once the flight profiles were programmed, the multi-rotor UAV magnetic system was able to automatically perform survey tasks, including take-offs and landings. The flight altitude was set to 3.5 m above ground level (AGL), with a flight speed of 2 m/s. Figure 9 shows the flight profiles of the UAV. The origin of the local cartesian coordinates was the starting point of the flight. Profiles above the survey area can be obtained after cutting off the undesired and curved flight data (see Figure 9). The two-dimensional magnetic map of the survey area can be obtained by interpolating the original data, as shown in Figure 10a. The characteristics of the target signal were masked due to the original signal containing a lot of random noise. The denoised data of each flight profile were obtained using EMD-PE-WSTD, EEMD-PE-WSTD, and the proposed CEEMDAN-PE-CC-WSTD method, and the corresponding results of magnetic maps are shown in Figure 10b-d, respectively.  To evaluate the quality improvement of magnetic maps after denoising, we introduce two parameters, i.e., peak signal-to-noise ratio (PSNR) and structural similarity (SSIM) [40]. The target area (East-West: −12 m to −6 m, South-North: 16 m to 25 m) of Figure 10 was selected, and the reference magnetic map of the target area can be obtained using the lowfrequency electromagnetic field simulation software, ANSYS Maxwell 19.0. Information about the geomagnetic field (e.g., inclination and declination) can be obtained according to the International Geomagnetic Reference Field (IGRF) model. The results of PSNR and SSIM of different denoising methods are shown in Table 7. As shown in Table 7, the results obtained by the proposed method have the largest PSNR and SSIM, therefore, the effectiveness of the proposed method is proved. Considering the quasi-static characteristics of the target signal, the complexity of the flight profile data can reflect the noise level on the other hand, i.e., the flight profile data with lower PE means less noise. The results of the PE of the data of each flight profile before and after noise reduction are shown in Figure 11. The average PE of the original data, denoised data using EMD-PE-WSTD, EEMD-PE-WSTD, and CEEMDAN-PE-CC-WSTD methods are 0.9904, 0.4970, 0.4679, and 0.4231, respectively. Data denoised by the proposed method have the lowest average PE among the three methods, which indicates that the complexity of the data is significantly reduced, and the noise is greatly suppressed. Figure 11. Results of PE of each flight profile data before and after noise reduction.

Conclusions
In this paper, a novel noise reduction method for multi-rotor UAV magnetic survey data based on CEEMDAN, PE, CC and WSTD is proposed. The CEEMDAN method is used to decompose the raw magnetic data into a series of IMFs with different scales. The quartile of PE is applied to divide the IMFs into four categories, i.e., the noise IMFs, noise-dominant IMFs, signal-dominant IMFs, and signal IMFs. Correlation coefficients are introduced to identify the real signal-dominant IMFs, and WSTD is applied to the real signal-dominant IMFs. Finally, the denoised signal can be obtained by combining the signal IMFs and the denoised IMFs. The proposed method is validated through experiments on both simulated synthetic signals and multi-rotor UAV magnetic survey data. The denoised data obtained by the proposed method are qualitatively and quantitatively analyzed and compared with EMD-PE-WSTD and EEMD-PE-WSTD methods. The results show that the proposed CEEMDAN-PE-CC-WSTD method can significantly suppress the noise and obtain a clearer target signal, which is very beneficial to the follow-up data interpretation. Our future work will include further verification of the proposed method via more UAV magnetic survey applications.