Adaptive Magnetic Anomaly Detection Method with Ensemble Empirical Mode Decomposition and Minimum Entropy Feature

Due to the fast attenuation of the magnetic field along with the distance, the magnetic anomaly generated by the remote magnetic target is usually buried in the magnetic noise. In order to improve the performance of magnetic anomaly detection (MAD) with low SNR, we propose an adaptive method of MAD with ensemble empirical mode decomposition (EEMD) and minimum entropy (ME) feature. The magnetic data is decomposed into the multiple intrinsic modal functions (IMFs) with different scales by EEMD. According to a defined criterion, the magnetic noise and magnetic signal are reconstructed based on IMFs, respectively. Entropy feature of reconstructed magnetic signal is extracted based on the probability density function (PDF) of the noise which is updated by the reconstructed magnetic noise. Compared to the traditional minimum entropy method, the entropy feature extracted by the proposed method is more obvious. The magnetic anomaly is detected whenever the entropy feature drops below the threshold. Thus, it is effective for revealing the weak magnetic anomaly by the proposed method. The measured magnetic noise is used to validate the performance of the proposed method. The results show that the detection probability of the proposed method is higher with low input SNR.


Introduction
The magnetic field generated by a magnetic target, which changes the distribution of the ambient magnetic field around it, is called the magnetic anomaly. Magnetic anomaly detection is a passive method of detecting the ferromagnetic objects. It has been used in many areas, such as unexploded ordnance detection [1][2][3], underwater pipeline or cable detection [4,5], and human medical investigation [6,7]. Due to fast attenuation of the magnetic field along with the distance, the magnetic anomaly generated by the remote target is usually buried in the magnetic noise [8,9]. Thus, we need the effective method to improve the performance of magnetic anomaly detection.
Recently, researchers have proposed several methods of magnetic anomaly detection [10][11][12][13][14][15][16][17][18][19]. These methods can be divided into two categories. One category is to improve the SNR of the anomaly signal by the noise suppression [20][21][22][23]. Zhou et al. [20] proposed a trend filtering based on empirical mode decomposition for MAD. Wavelet denoising and the Karhunen-Loeve expansion have been proposed to detect the magnetic anomaly. Liu et al. [8] proposed an adaptive coherent noise suppression method using coherence. However, those methods are not sufficient for magnetic anomaly detection when SNR < 1. For this situation, Ginzburg et al. [22] proposed a method of MAD using three orthonormal basis functions (OBF), which relies on the matched filtering. Fan et al. [24] proposed an improved method using four orthonormal basis functions, which can not only detect magnetic anomaly but also determine the orientation of the target. These methods are optimal for detecting the magnetic anomaly buried in the Gaussian white noise. In some situations, the geomagnetic noise is considered the noise with a power spectral density (PSD) of 1/f α , where 0 < α < 2. In order to handle the noise with PSD of 1/f α , Sheinker et al. [25] designed a whitening filter, which improved the performance of MAD. This method can be considered an autoregressive (AR) process and work effectively with a high order of filter. The other category is to reveal the magnetic anomaly by the transformation. The minimum entropy method transforms the measured data into entropy value to detect the anomaly [12]. Because the most common-mode geomagnetic noise can be removed by the differential signal, the method based on the information entropy of the differential signal was more appropriate for the target detection in the complex environment [17]. The above methods need the PDF of the magnetic noise. The high-order crossing method is an alternative method for spectral analysis using zero-crossing count [14]. The advantages of these methods are easy to implement and low computational complexity. However, the detection performance of these methods may be limited by low SNR of input signal. In addition, some methods based on machine learning are proposed for magnetic anomaly detection [26][27][28][29][30].
In order to improve the detection performance in the case of low SNR and the complex magnetic environment, we propose an adaptive method of MAD with ensemble empirical mode decomposition and minimum entropy (EEMD-ME) feature in this paper. Meanwhile, the similar methods with empirical mode decomposition and minimum entropy (EMD-ME) feature and complete ensemble empirical mode decomposition with adaptive noise and minimum entropy (CEEMDAN-ME) feature are also implemented. In the methods, the measured magnetic data is decomposed into multiple IMFs with different scales. According to the characteristics of the IMFs, the parameters of PDF of magnetic noise are updated in real time using the reconstructed magnetic noise. Then, the entropy feature of the reconstructed magnetic signal is calculated in a moving window with the updated PDF of the noise. Finally, the magnetic anomaly is detected whenever the entropy feature drops below the threshold. In brief, there are three advantages of the proposed method: (1) the parameters of PDF of the noise are updated in real time; (2) SNR of input signal is improved; and (3) low computational complexity is implemented in the detection. Thus, the proposed method has better robustness, higher detection probability in the weak magnetic anomaly detection.

Detection Theory
2.1. EMD, Ensemble EMD, and Complete Ensemble EMD with Adaptive Noise. Empirical mode decomposition is pioneered by Huang for adaptively processing the nonlinear and nonstationary signals [31]. It decomposes the signal into a series of IMFs with different time scales, which must satisfy two conditions: (1) the number of extreme points in IMF is equal to or not more than that of one zero-crossing point; (2) the mean of the envelope determined by the maximum and minimum values is zero. The algorithm can be described as follows: Step 1. Set k = 0 and find all the maxima (minima) points of the original signal r 0 = xðtÞ. Interpolate between maxima (minima) points to obtain the upper (lower) envelope r + k ðtÞ (r − k ðtÞ) Step 2. Calculate the mean envelope m k ðtÞ = ðr + 0 ðtÞ + r − 0 ðtÞÞ/2 Step 3. Obtain the IMF candidate c k ðtÞ = r k ðtÞ − m k ðtÞ Step 4. Determine whether the candidate c k ðtÞ satisfy two conditions: (i) Yes, save d k+1 ðtÞ as an IMF and calculate the residual r k+1 ðtÞ = xðtÞ − ∑ k i=1 c i . Do k = k + 1 and treat r k+1 ðtÞ as the input signal in Step 2 (ii) No, treat d k+1 ðtÞ as the input signal in Step 2.
In order to alleviate the "mode mixing" phenomenon in EMD, noise-assisted versions have been proposed, such as EEMD [32,33] and complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) [34][35][36]. The main idea of the noise-assisted versions is that the different finite variance white noise n a ðtÞ plus the original signal and the average of the corresponding IMFs is calculated as the final IMFs.
The algorithm of EEMD is described as follows: Step 1. Add the finite variance white noise n i ðtÞ ði = 1, ⋯,IÞ to the original signal with I times and obtain the noise-added signal: x i ðtÞ = xðtÞ + εn i a ðtÞ, where ε denotes the standard deviation of added noise amplitude and ε > 0 Step 2. Decompose the noise-added signal x i ðtÞ into IMFs c i k ðtÞ and a residual r i ðtÞ by EMD Step 3. Calculate the average of each IMF and the residual as the final IMF and residual: c k ðtÞ = 1/I∑ I i−1 c i k ðtÞ and rðtÞ = 1/I∑ I i−1 r i ðtÞ. For CEEMDAN, let E k ð•Þ be the kth mode of a given signal decomposed by EMD. The algorithm of CEEMDAN [36] is given as follows: Step 1. Add the finite variance white noise n i a ðtÞ ði = 1, ⋯,IÞ to the original signal with I times and obtain the noise-added signal: Step 5. Decompose the signal r k ðtÞ + ε k E k ðn i a ðtÞÞ by EMD to obtain the first IMF mode. Define the (k + 1)th modes as:

Journal of Sensors
Step 6. Repeat Steps 4-5 until the residual r k+1 ðtÞ satisfies the stopping criterion.
Thus, the signal xðtÞ can be described by the IMF components and a residual component as follows: where c k ðtÞ denotes the kth IMF component and rðtÞ denotes the residual component. Generally, the IMF components are sorted from high frequency to low frequency. Therefore, there will be a component indexed by h, which can be considered the division of energy distribution from noise to signal. After the hth IMF component, the dominant energy of the magnetic anomaly signal is distributed among these IMF components. The criterion based on continuous mean square error (CMSE) [34] is proposed to find the dividing point. The CMSE is defined as where N is the length of each IMF components and K is the number of IMF components. If the first significant change of CMSE is (h − 1)th IMF component, the dividing point is h. It means that the energy components of adjacent IMF components are significantly different. Thus, the dividing point can be used to split the raw signal. In the magnetic data, the components before c h are considered the main energy of the noise. The components after c h contain the main energy of the magnetic signal.

Minimum Entropy Feature of Magnetic Anomaly.
A magnetic target can be considered a magnetic dipole when the distance between the target and the magnetometer is more than three times the target largest dimension. For detecting a magnetic target, the MAD model is shown in Figure 1. The magnetic field generated by the target can be expressed as [37] where μ 0 = 4π × 10 -7 H/m is the permeability of free space M is the magnetic moment of the target, and R is the distance between the target and the magnetic sensor.
In practice, the magnetic field B measured by the magnetic sensor is consisted of the ambient magnetic field B e and magnetic field B a generated by the target when the target is present. When the target is far from the sensor, jB a j is much less than jB e j. When using the total field magnetometers to detect the target, we can obtain the scalar value of the measurement as [13,16] where u denotes the unit vector of the ambient field.
In the measurement, the magnetic anomaly generated by the target can change the magnetic noise pattern. In the information theory, entropy is used as a measurement of information. In Ref. [12], entropy is used to detect the changes in the magnetic noise pattern. It means that the entropy feature of the magnetic field will dramatically change when the ambient magnetic field contains the magnetic anomaly. Thus, entropy value can be considered a feature of the magnetic anomaly. The value of entropy can be calculated in a moving window of L samples. The expression is given as follows: where pðx n Þ denotes the probability density function of ambient magnetic field. The detection occurs whenever the entropy value drops below a threshold value.
In the ME method, a priori condition of PDF of magnetic noise is needed when using the entropy feature for the target detection. The PDF is acquired by using statistical methods based on the measured magnetic noise. However, when the environmental conditions change, it is necessary to update the PDF by the new data [17]. It is hard to acquire the new background noise in real time when the magnetometers mounted on the platform are moving. Thus, the detection performance of ME method will be affected in this situation.

Magnetic Anomaly Detection with EMD and Minimum
Entropy Feature. In order to overcome the problems in the ME method in Ref. [12], we propose an adaptive MAD method with ensemble empirical mode decomposition and minimum entropy feature in this paper. Based on EEMD, the SNR of the reconstructed signal can be improved according to the dividing point. Meanwhile, the reconstructed noise can be approximated as magnetic background noise which is considered the Gaussian distribution [12]. The reconstructed noise can be used to update the parameters of the PDF of the noise in real time. Thus, the proposed method not only overcomes the problem of updating the PDF of the ambient magnetic noise in the ME method but also improves the detection      Journal of Sensors probability of the weak magnetic anomaly. Similarly, the EMD-ME and CEEMDAN-ME methods are also implemented. The steps of the proposed method are as follows: Step 1. Decompose the original signal into a series of IMFs components and a. residual component Step 2. Reconstruct the signal and noise according to the dividing point Step 3. Update the parameters of PDF of noise based on the reconstructed noise Step 4. Calculate the entropy feature of the reconstructed signal with the PDF.
Step 5. According to the threshold, determine whether the detection occurs.

Experiment
3.1. Experiment Design. In order to reduce the influence of temporal magnetic variations on magnetic anomaly detection, we designed a one-dimensional sensor array (shown in Figure 2) to collect the magnetic background noise. The magnetic sensor was the optically pumped magnetometer with high sensitivity, of which the intrinsic noise was about The magnitude of the moment was 30 A m 2 , and the unit direction vector was (-0.65, -0.71, -0.26). The typical magnetic anomaly can be obtained by the simulation. Then, the real magnetic noise measured by sensor array was added to the simulated signal. Finally, the synthetic magnetic signal with real-word noise and the simulated anomaly signal can be obtained.

3.2.
Results. The synthetic signal of magnetic field was obtained, shown in Figure 3. It can be seen that the magnetic anomaly signal was completely buried by the noise. The SNR of the synthetic signal was about -10.47 dB. The EEMD-ME method was applied to the synthetic signal, in which the parameters of EEMD were set as follows: the standard deviation of added noise amplitude ε = 0:1, and the ensemble member I = 100. By the EEMD, the synthetic signal was decomposed into a set of IMFs, shown in Figure 4(b). Similarly, the EMD-ME and CEEMDAN-ME methods were also applied to the synthetic signal. In the CEEMDAN-ME method, the parameters were the same as that of the EEMD-ME method. The decomposition results by EMD and CEEMDAN were shown in Figures 4(a) and 4(c). From Figure 4, we can see that different IMFs reveal the different time scales of the magnetic field signal. According to formula (2), the dividing point was calculated as h = 5. Thus, IMF1-IMF5 represented the high-frequency components, in which the main energy of the noise was contained. IMF6-Residual represented the low-frequency components, in which the main energy of the signal was contained. Based on the IMFs, the reconstructed magnetic anomaly signal and the reconstructed magnetic nosie were obtained, shown in Figure 5.
After the reconstruction, an approximated PDF of magnetic noise was obtained by the reconstructed magnetic noise. The parameters of PDF of the noise were updated based on the reconstructed noise. The entropy feature of the reconstructed signal can be calculated with the PDF of the noise. The outputs of the methods are shown in It can be seen that the magnetic anomaly was obviously revealed by the proposed methods. In addition, the traditional ME method was also applied to the synthetic signal. The output of the ME method is shown in Figure 6(a). From the results, the detection performance of the proposed method was better than that of the ME method in low input SNR.
Compared to the traditional ME method, there are two advantages of the proposed method in weak magnetic anomaly detection. (1) The approximated PDF of the noise can be updated in real time by the reconstructed magnetic noise. It overcomes the problem that the PDF of magnetic noise is difficult to update in the changing movement. (2) The input SNR can be improved by the reconstruction signal based on IMFs. Compared to the SNR of the original signal with -10.47 dB, the SNR of the reconstructed signal is improved to 9.37 dB in the EEMD-ME method, 8.16 dB in the EMD-ME method, and 7.14 dB in the CEEMDAN-ME method, respectively. In addition, the execution time of updating the parameters of PDF of the noise is about 0.87 s in the EEMD-ME method, 0.07 s in the EMD-ME method and 2.51 s in the CEEMDAN-ME method, respectively. It can be considered the update of PDF in real time, especially in the EMD-ME and EEMD-ME methods.
The mentioned SNR is calculated by the following equation: where x denotes the signal, n denotes the noise, and N is the length of the signal.   In order to evaluate the effectiveness of the proposed method, the Monte Carlo simulation was carried out. In the simulation, the simulated target was moving in parallel to the x-axis with a constant velocity of 2/m, the range of X position in (-200, 200) m. The magnitude of magnetic moment of the target was 30 A m 2 . The orientation (inclination and declination of the moment) of the target was randomly varied. According to the Neyman-Pearson criterion [21], the threshold values of EEMD-ME, EMD-ME, and CEEMDAN-ME were set to 0.9991, 0.9955, and 0.9984 with the false alarm rate 3.4% based on the real magnetic noise. The closest proximity approach (CPA) from the magnetic sensors to the moving track of the target was changed to evaluate the performances of the methods. Table 1 shows the detection performance versus with different CPA.
From Table 1, the magnetic anomaly generated by the target can be detected by the EMD-ME, EEMD-ME, and CEEMDAN-ME methods, and the detection accuracy was 100% when CPA was less than 25.18 ( ffiffiffiffiffiffiffiffiffiffiffiffiffiffi P 2 y + P 2 z q ). As the CPA increased, the detection accuracy decreased. When CPA was 30.14, the detection accuracy of EEMD-ME was 99.5%, which was higher than that of the CEEMDAN-ME and EMD-ME methods. Similarly, when the CPA was increased to 35.12, the detection accuracy of EEMD-ME dropped to 90.5%. Under the same condition, the detection accuracy of CEEMDAN-ME and EMD-ME methods dropped to 79.6% and 57.4%, respectively. From the result, it is noted that the detection performance of EEMD-ME is better than of EMD-ME and CEENDAN-ME methods in weak magnetic anomaly detection. Meanwhile, the execution time of the EEMD-ME method is about 1 s. Thus, the EEMD-ME method is better in the weak magnetic anomaly detection.

Conclusions
In this paper, we propose the adaptive method with ensemble empirical mode decomposition and minimum entropy feature, which can improve the performance of magnetic anomaly detection with low SNR. In the method, the raw magnetic data is processed by EEMD, and IMFs are obtained. According to the properties of IMFs, the reconstructed magnetic noise is used to update the parameters of PDF of it. The reconstructed magnetic signal is used to extract the entropy feature of the magnetic anomaly using the updated PDF. Whenever the entropy feature drops below the threshold, the magnetic anomaly is detected. Compared to the traditional ME method, the entropy feature extracted by the proposed method is more obvious and useful to reveal the weak magnetic anomaly. The execution time of EEMD-ME method is about 1 s. Meanwhile, the EMD-ME and CEEMDAN-ME methods are implemented and used to detect the weak magnetic anomaly signal. The results show that when the CPA is increased from 30.14 to 35.12, the detection accuracy of EEMD-ME method drops from 99.5% to 90.5%. However, the detection accuracy of EMD-ME and CEEMDAN-ME methods drops from 91.0% and 97.7% to 57.4% and 79.6%, respectively. Therefore, the detection performance of EEMD-ME method is better in the weak magnetic anomaly detection.

Data Availability
The magnetic field data used to support the findings of this study were supplied by Geomagnetic Sensing and Application Lab, Harbin Engineering University, under license and so cannot be made freely available. Requests for access to these data should be made to Prof. Chong Kang, kangchongheu@163.com.

Conflicts of Interest
The authors declare that they have no conflict of interests.