Processing and Interpretation of UAV Magnetic Data: A Workflow Based on Improved Variational Mode Decomposition and Levenberg–Marquardt Algorithm

Abstract

Unmanned aerial vehicles (UAVs) have become a research hotspot in the field of magnetic exploration because of their unique advantages, e.g., low cost, high safety, and easy to operate. However, the lack of effective data processing and interpretation method limits their further deployment. In view of this situation, a complete workflow of UAV magnetic data processing and interpretation is proposed in this paper, which can be divided into two steps: (1) the improved variational mode decomposition (VMD) is applied to the original data to improve its signal-to-noise ratio as much as possible, and the decomposition modes number K is determined adaptively according to the mode characteristics; (2) the parameters of target position and magnetic moment are obtained by Euler deconvolution first, and then used as the prior information of the Levenberg–Marquardt (LM) algorithm to further improve its accuracy. Experiments are carried out to verify the effectiveness of the proposed method. Results show that the proposed method can significantly improve the quality of the original data; by combining the Euler deconvolution and LM algorithm, the horizontal positioning error can be reduced from 15.31 cm to 4.05 cm, and the depth estimation error can be reduced from 16.2 cm to 5.4 cm. Moreover, the proposed method can be used not only for the detection and location of near-surface targets, but also for the follow-up work, such as the clearance of targets (e.g., the unexploded ordnance).

1. Introduction

Unmanned aerial vehicles (UAVs) have been applied in many fields during the past decade, e.g., remote sensing, geological prospecting, and unexploded ordnance (UXO) detection [1,2,3]. Among these applications, the use of UAV for magnetic survey is a booming branch [4]. Compared with the traditional terrestrial magnetic survey, a UAV magnetic survey can cover a wider range with a higher efficiency. Besides, it is also easy to operate, low-cost, and has good safety compared with the manned aircraft magnetic survey. Moreover, a UAV magnetic survey can be conducted in areas that are difficult to access or would pose a potential hazard to operators, which means that some gaps in traditional magnetic surveys can now be studied [5,6].

Considerable progress has been made in the integration of UAV magnetic systems in the past decade, however, the low quality of UAV magnetic data and the lack of effective data interpretation methods have become the biggest obstacles to its further development. Several attempts have been made to improve the quality of UAV magnetic data, e.g., the magnetic gradiometer is towed about 2 m below the UAV in [7] to mitigate the interference generated by the platform; a similar method is also adopted in [8,9,10,11]. However, this design introduced the directional and positional error, meanwhile increasing the risk of impact to the magnetometers. Filtering methods [12,13,14] are considered to improve the quality of UAV magnetic data, however, there are three main shortcomings: (1) filter parameters need to be determined according to a priori knowledge; (2) the frequency characteristic of the target signal will change with the conditions, such as the flight speed of the platform, which means that the parameters of the filter should also be changed, and (3) once the frequencies of the target signal and noise overlap, the filtering method will be difficult to apply. Liu et al. [15] proposed an adaptive elimination method of geomagnetic background noise, but it should be noted that this method requires an additional magnetometer, which means a shorter flight time for UAV magnetic systems.

Recently, a method based on the empirical mode decomposition (EMD) and its improved form has been applied to UAV magnetic data processing [16]. Although the effectiveness of this method has been proven by experimental results, the processing workflow is somewhat complex, and meanwhile, the physical meaning of the obtained intrinsic mode function (IMF) is not obvious. Inspired by this successful case, variational mode decomposition (VMD) is introduced in this paper for the processing of UAV magnetic data. To the best of our knowledge, this is the first time that VMD has been applied to this field. VMD is a completely non-recursive decomposition method proposed by Dragomiretskiy et al. [17], and has been widely used in gear fault diagnosis, feature extraction, and seismic time-frequency analysis [18,19,20]. However, the performance of VMD is affected by the decomposition modes number K, under-decomposition or over-decomposition will occur if K is not appropriate. In [21], K is set to be consistent with that of EMD, which is obviously not convenient. Appropriate K can be obtained by using the optimization algorithm [22], but it is time-consuming and inefficient. Previous studies [23,24] have shown that the mode characteristics can be used to optimize K, and this method is also used in this paper.

Once K is determined, the original data can be processed, and the magnetic map of the survey area can be obtained. The anomalies caused by ferromagnetic targets will be reflected in the magnetic map, and the parameter estimation of these targets is the main task of data interpretation [25]. Euler deconvolution [26,27,28] is often used to estimate the position and magnetic moment of the target. Davis et al. [29] proposed a terrestrial magnetic survey workflow for the detection of UXO by using sliding window-based Euler deconvolution, however, the changing size of the sliding window makes the process complex and cumbersome. Zhou and Zhang [30] proposed a method for the location of multiple magnetic targets by combining gradient tensor inversion and K-means clustering. Simulation results show that the positioning error of this method is about 0.05 m, which is acceptable for most applications. In addition, methods based on Euler deconvolution, such as tilt-Euler deconvolution [31] and joint-gradient Euler deconvolution [32], have also been proposed to better interpret magnetic data.

The optimization method is a good choice considering the overdetermined problem to be solved here. Ding et al. [33] proposed a new method for estimating locations and moments of multiple dipole-like magnetic sources using differential evolution (DE); the results of synthetic and field tests showed that the location of magnetic source can be estimated well. A local optimization method, the Levenberg–Marquardt (LM) algorithm, is introduced in this paper to obtain a more accurate estimation of the target parameters. LM algorithm, also known as the Damped Least Squares (DLS) algorithm, is very suitable for solving the nonlinear problem with the advantage of fast convergence [34,35,36]. The result of Euler deconvolution is used as the initial value of the LM algorithm. By minimizing the misfit between the predicted and observed magnetic anomalies, the target location accuracy can be significantly improved. It should be noted that the proposed method is not superior to other methods such as DE but provides a simple and fast methodology in target parameters estimation.

All in all, the main contributions of this manuscript are as follows:

• A complete workflow based on the improved VMD, and the joint Euler deconvolution-LM algorithm is proposed for the processing and interpretation of UAV magnetic data.
• The VMD method is applied for the processing of UAV magnetic data, and the decomposition modes number K is adaptively determined according to the mode characteristics.
• The positioning accuracy of near-surface targets is significantly improved by combining Euler deconvolution and the LM algorithm without increasing too much complexity, which is very helpful for the follow-up work.

The rest of this paper is organized as follows: Section 2 presents the related theories of VMD, Euler deconvolution, and LM algorithms. The proposed workflow for the processing and interpretation of UAV magnetic data is described in detail in Section 3. Experiments are carried out, and the results are analyzed in Section 4. Section 5 is the discussion and final remarks of the proposed workflow.

2. Theoretical Background

The proposed method can be divided into two steps, namely, magnetic data processing based on the improved VMD method and magnetic data interpretation based on the combined Euler deconvolution-LM algorithm. The theory and implementation of the relevant methods are briefly introduced in this part.

2.1. Magnetic Data Processing

2.1.1. VMD Algorithm

As declared in [17], VMD decomposes the original signal $x\left(t\right)$ into K band-limited IMFs, which have the sparsity properties, the narrow-band IMF $u_k\left(t\right),k=1,2,\dots ,K$ is defined as an amplitude-modulated-frequency-modulated (AM-FM) signal, as given by:

$$u_k\left(t\right)=A_k\left(t\right)\cos \left({\varphi }_k\left(t\right)\right),$$

(1)

where $A_k\left(t\right)$ and ${\varphi }_k\left(t\right)$ are the instantaneous amplitude and phase, respectively. The instantaneous frequency ${\omega }_k\left(t\right)$ is defined as the derivative of ${\varphi }_k\left(t\right)$

$${\omega }_k\left(t\right)=d{\varphi }_k\left(t\right)/dt={\varphi }_k^{\prime }\left(t\right).$$

(2)

The Hilbert transform is applied to each mode, then the mode’s spectrum is shifted to baseband according to Fourier transform. After that, the Gaussian smoothness is used to estimate the bandwidth of each mode. The constrained variational problem of VMD is given by:

$$\begin{array}{c}\underset{\left\{u_k\right\},\left\{{\omega }_k\right\}}{\min }\left\{{\displaystyle \sum _{k=1}^K}\Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)*u_k\left(t\right)\right]\cdot e^{-j{\omega }_{k}t}{\Vert }_2^2\right\}\\ s.t.{\displaystyle \sum _{k=1}^K}{u}_k=x\left(t\right)\end{array},$$

(3)

where $\left\{u_k\right\}$ and $\left\{{\omega }_k\right\}$ are the K mode components and center frequencies, respectively. $\delta \left(t\right)$ is the Dirac distribution, $*$ refers to convolution, and ${\Vert \cdot \Vert }_2^2$ is the square of the Euclidean norm. (3) can be transformed into unconstrained problem by utilizing the quadratic penalty factor $\alpha$ and Lagrange multiplier $\lambda \left(t\right)$, the augmented Lagrangian form of (3) can be expressed as:

$$L\left(\left\{u_k\right\},\left\{{\omega }_k\right\},\lambda \right)=\alpha {\displaystyle \sum }_{k=1}^K\Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)*u_k\left(t\right)\right]\cdot e^{-j{\omega }_{k}t}{\Vert }_2^2+\Vert x\left(t\right)-{\displaystyle \sum }_{k=1}^K{u}_k{\Vert }_2^2+\langle \lambda \left(t\right),x\left(t\right)-{\displaystyle \sum }_{k=1}^K{u}_k\rangle ,$$

(4)

where $\langle a,b\rangle$ represents the inner product of $a$ and $b$. To solve this unconstrained variational problem, the alternate direction method of multipliers is applied, and the values of $u_k$, ${\omega }_k$, and $\lambda$ are updated in frequency respectively. $u_k$ can be updated using the equivalent minimization problem, defined as follows:

$$u_k^{n+1}\left(t\right)=\underset{u_k\in X}{\mathrm{argmin}}\left\{\alpha \Vert {\partial }_t\left[\left(\delta \left(t\right)+\frac{j}{\pi t}\right)*u_k\left(t\right)\right]\cdot e^{-j{\omega }_{k}t}{\Vert }_2^2+\Vert x\left(t\right)-{\displaystyle \sum }_i{u}_i\left(t\right)+{\frac{\lambda \left(t\right)}{2}\Vert }_2^2\right\},$$

(5)

where $X$ is the functional space of signal and modes, the ${\cdot }^n$ and ${\cdot }^{n+1}$ are omitted for the fixed ${\omega }_k$ and $u_{i\ne k}$.

The formula for updating $u_k$, ${\omega }_k$, and $\lambda$ are given by:

$${\widehat{u}}_k^{n+1}\left(\omega \right)=\left[\widehat{x}\left(\omega \right)-{\displaystyle \sum }_{i\ne k}{\widehat{u}}_i\left(\omega \right)+\frac{\widehat{\lambda }\left(\omega \right)}{2}\right]\cdot \frac{1}{1+2\alpha {\left(\omega -{\omega }_k\right)}^2}$$

(6)

$${\omega }_k^{n+1}=\frac{{{\displaystyle \int }}_0^{\infty }\omega {\left|{\widehat{u}}_k\left(\omega \right)\right|}^{2}d\omega }{{{\displaystyle \int }}_0^{\infty }{\left|{\widehat{u}}_k\left(\omega \right)\right|}^{2}d\omega }$$

(7)

$${\widehat{\lambda }}^{n+1}\left(\omega \right)={\widehat{\lambda }}^n\left(\omega \right)+\tau \left[\widehat{x}\left(\omega \right)-{\displaystyle \sum }_k{\widehat{u}}_k^{n+1}\left(\omega \right)\right],$$

(8)

where ${\widehat{u}}_k\left(\omega \right)$, $\widehat{x}\left(\omega \right)$, and $\widehat{\lambda }\left(\omega \right)$ are the Fourier transform of $u_k\left(t\right)$, $x\left(t\right)$, and $\lambda \left(t\right)$. $\tau$ is the parameter of noise tolerance. The above calculation is continuously implemented until the following iteration termination condition is met:

$${{\displaystyle \sum }}_k\Vert {\widehat{u}}_k^{n+1}-{\widehat{u}}_k^n{\Vert }_2^2/\Vert {\widehat{u}}_k^n{\Vert }_2^2<\epsilon ,$$

(9)

where $\epsilon$ is the tolerance parameter of the convergence criterion. The default values of $\epsilon$ and $\tau$ of the standard VMD are utilized in this study.

A key parameter affecting the performance of VMD is the selection of the decomposition modes number K. Previous studies are generally based on experience or optimization algorithms to determine K, but all these methods present drawbacks. Ref. [23] proposed an adaptive VMD method for signal processing based on mode characteristics, several indicators, e.g., permutation entropy (PE), extreme value in frequency domain, kurtosis, dominant frequency spacing between neighboring modes, and energy loss coefficient are used to determine K. Li et al. [37] put forward a frequency-aided method to solve the problem of choosing parameters for VMD and judged whether the decomposition was excessive by comparing the difference of IMF center frequencies. Inspired by these studies, an improved VMD method based on mode characteristics is utilized in the following subsection to determine K.

2.1.2. The Improved VMD Algorithm

As mentioned earlier, inappropriate K can lead to under-decomposition or over-decomposition, which both affect the performance of VMD. Here, we propose an improved VMD method which is completely data-driven and based on the mode characteristics to determine K. The method mainly consists of six steps, explained as follows:

• The starting point is to determine the following parameters: the minimum mode number $K_{min}$, the maximum mode number $K_{max}$, the threshold $\mu$, and the penalty parameter $\alpha$;
• Let $K=K_{min}$, perform VMD on the input signal $x\left(t\right)$, a set of IMFs can be obtained and recorded as $\left\{IM{F}_{i,K},i\in \left[1,K\right]\right\}$, the corresponding center frequencies are recorded as $\left\{f_{i,K},i\in \left[1,K\right]\right\}$;
• Let $K_n=K+1$, we can also obtain a set of IMFs and the corresponding center frequencies, namely, $\left\{IM{F}_{p,K_n},p\in \left[1,K_n\right]\right\}$, and $\left\{f_{p,K_n},p\in \left[1,K_n\right]\right\}$;
• Compare $\left\{f_{i,K}\right\}$ and $\left\{f_{p,K_n}\right\}$ to determine the new center frequency in $\left\{f_{p,K_n}\right\}$, and the IMF corresponding to the new center frequency is obtained as $r\left(t\right)$;
• Calculate the correlation coefficients $\left\{Coe{f}_i\right\}$ between $r\left(t\right)$ and $\left\{IM{F}_{i,K}\right\}$;
• Execute judgement, if “$\exists Coe{f}_i>\mu$” is true, it means that the decomposition is excessive, then make the optimal decomposition number $K_{best}=K$. Otherwise, continue to implement 3–5 until the criterion is met.

In this paper, the threshold $\mu$ is set to 0.4 empirically, besides, to judge whether there is under-decomposition or not, we introduce the energy loss coefficient $e_l$ which is defined as the energy ratio of the decomposition residual to the original signal [23]:

$$e_l=\frac{\Vert x\left(t\right)-{{\displaystyle \sum }}_k{u}_k{\Vert }_2^2}{{\Vert x\left(t\right)\Vert }_2^2},$$

(10)

where $u_k$ is the kth IMF, and ${\displaystyle \sum }_k{u}_k$ is the reconstructed signal. When under-decomposition occurs, there is a large error between the reconstructed signal and the original signal, which means that the $e_l$ is also relatively large. Therefore, the under-decomposition can be determined by judging whether the $e_l$ is greater than the threshold $\delta$. In this paper, $\delta$ is set to 0.1. By using the above method, a suitable K can be obtained, which can be used to process the UAV magnetic data.

2.2. Magnetic Data Interpretation

After completing the processing of the UAV magnetic data, the data interpretation, or inversion, can be carried out. The purpose of data interpretation is to infer the parameters of the target, e.g., the spatial position and magnetic moment, from the observed magnetic anomalies. Euler deconvolution has been adapted as a tool for the quick initial interpretation due to its simplicity. However, Euler deconvolution, like other geophysical inversion methods, can only provide the approximate location of the target, which cannot meet our requirements. Therefore, an optimization method based on the LM algorithm is proposed to further improve the target positioning accuracy. The Euler deconvolution and the LM algorithm are introduced in the following subsections, respectively.

2.2.1. Euler Deconvolution

The theoretical basis of Euler deconvolution is that the magnetic anomaly field generated by the target satisfies the following Euler’s homogeneity equation [26,29]:

$$\frac{\partial T}{\partial x}\left(x-x_0\right)+\frac{\partial T}{\partial y}\left(y-y_0\right)+\frac{\partial T}{\partial z}\left(z-z_0\right)=-NT,$$

(11)

where $T$ is the observed field at $\left(x,y,z\right)$ generated by the target at location $\left(x_0,y_0,z_0\right)$, $N$ denotes the structural index (SI) which varies with different source types.

Applying (11) to several neighboring data points with a prior value for N, the location of the target can be solved in the least squares sense as follows:

$$P={\left(A^{T}A\right)}^{-1}{A}^{T}b,$$

(12)

$$P=\left[\begin{array}{c}{x}_0\\ y_0\\ z_0\end{array}\right],A=\left[\begin{array}{ccc}\frac{\partial T_1}{\partial x}& \frac{\partial T_1}{\partial y}& \frac{\partial T_1}{\partial z}\\ \frac{\partial T_2}{\partial x}& \frac{\partial T_2}{\partial y}& \frac{\partial T_2}{\partial z}\\ \begin{array}{c}⋮\\ \frac{\partial T_n}{\partial x}\end{array}& \begin{array}{c}⋮\\ \frac{\partial T_n}{\partial y}\end{array}& \begin{array}{c}⋮\\ \frac{\partial T_n}{\partial z}\end{array}\end{array}\right],b=\left[\begin{array}{c}{x}_1\frac{\partial T_1}{\partial x}+y_1\frac{\partial T_1}{\partial y}+z_1\frac{\partial T_1}{\partial z}+N{T}_1\\ x_2\frac{\partial T_2}{\partial x}+y_2\frac{\partial T_2}{\partial y}+z_2\frac{\partial T_2}{\partial z}+N{T}_2\\ \begin{array}{c}⋮\\ x_n\frac{\partial T_n}{\partial x}+y_n\frac{\partial T_n}{\partial y}+z_n\frac{\partial T_n}{\partial z}+N{T}_n\end{array}\end{array}\right],$$

(13)

where $n$ is the number of data points.

Usually, a specified size of window is used to select the neighboring data points, and the observed magnetic anomalies of these data points are thought to be caused by an isolated target.

2.2.2. Implementation of the LM Algorithm

Generally, the magnetic anomaly field generated by the near-surface target can be considered as a magnetic dipole [38], which can be calculated by:

$$\mathit{B}=\frac{{\mu }_0}{4\pi }\left[\frac{3\left(\mathit{m}·\mathit{r}\right)\mathit{r}}{r^5}-\frac{\mathit{m}}{r^3}\right],$$

(14)

where ${\mu }_0$ is the permeability in vacuum, $\mathit{m}$ is the dipole moment of the target, $\mathit{r}$ is the displacement vector from the target to the measurement point, and $r=\left|\mathit{r}\right|$. The anomaly field generated by the target is usually much smaller than the geomagnetic field, hence the observed anomaly signal $\Delta B$ is the projection of the anomaly field in the direction of the geomagnetic field [39], which can be expressed as:

$$\Delta B=f\left(\mathit{r},\mathit{m}\right)=\mathit{B}\cdot \frac{{\mathit{B}}_{\mathit{e}}}{\left|{\mathit{B}}_{\mathit{e}}\right|},$$

(15)

where ${\mathit{B}}_{\mathit{e}}$ and $\left|{\mathit{B}}_{\mathit{e}}\right|$ are the geomagnetic field and its magnitude, respectively. The main task of data interpretation is to invert the position and magnetic moment of the target through the magnetic anomaly data $\left(\Delta B_i,i=1,2,\dots ,N\right)$ of multiple observation points $\left({\mathit{r}}_i,i=1,2,\dots ,N\right)$. By minimizing the misfit between the observed anomaly field $\Delta B_i$ and the predicted anomaly field $f_i$ based on (15), the objective function can be established and expressed as follows:

$${\mathit{x}}^{*}=\underset{x}{\min }g\left(\mathit{x}\right)={{\displaystyle \sum }}_{i=1}^N{\left(f_i\left(\mathit{x}\right)-\Delta B_i\right)}^2.$$

(16)

where $\mathit{x}=\left(r_x,r_y,r_z,m_x,m_y,m_z\right)$ is the parameter of the target, N is the number of observation points. In fact, (16) is a least square problem which can be solved by the LM algorithm. The following update is used in the LM method to avoid the problem of rank-defect of Jacobian matrix that may occur in the Gauss–Newton method:

$${\mathit{x}}_{k+1}={\mathit{x}}_k-{\left(\mathit{J}{\left({\mathit{x}}_k\right)}^T\mathit{J}\left({\mathit{x}}_k\right)+{\lambda }_k\mathit{I}\right)}^{-1}\mathit{J}{\left({\mathit{x}}_k\right)}^{T}r\left({\mathit{x}}_k\right),$$

(17)

where $\mathit{J}\left({\mathit{x}}_k\right)$ is the Jacobian matrix that contains the first derivatives of the residual $r\left({\mathit{x}}_k\right)$, and

$$r\left({\mathit{x}}_k\right)=\left[\begin{array}{c}{f}_1\left({\mathit{x}}_k\right)-\Delta B_1\\ \begin{array}{c}{f}_2\left({\mathit{x}}_k\right)-\Delta B_2\\ ⋮\end{array}\\ f_N\left({\mathit{x}}_k\right)-\Delta B_N\end{array}\right],$$

(18)

$\mathit{I}$ is the unit matrix, and ${\lambda }_k$ represents the search direction. If ${\lambda }_k=0$, LM algorithm becomes the Gauss–Newton method, and when ${\lambda }_k$ is large, the LM algorithm becomes gradient descent with a small step size. Therefore, the search direction of the LM method combines the Gauss–Newton method and the steepest descent method. More details about LM algorithm can be found in [34,35,36]. It is worth emphasizing that due to the use of Taylor expansion for approximation, an implicit premise is that an initial value must be provided. Here, the solution of Euler deconvolution is used as the initial value of the LM method.

3. The Proposed Workflow Using Improved VMD and Joint Euler-LM Algorithm

This paper proposes a complete workflow for the processing and interpretation of UAV magnetic data, which can be mainly divided into two parts: (1) data processing based on the improved VMD method, where the optimal decomposition modes number K is determined adaptively by mode characteristics, and (2) data interpretation conducted by Euler deconvolution and LM algorithm to obtain a more realistic position and magnetic moment parameters of the target. The proposed workflow is depicted in Figure 1. The steps of the workflow can be summarized as follows:

Step 1:	Mission planning is the first step for UAV magnetic survey, many factors, e.g., the survey task, target size, topography of the survey area and weather conditions should be carefully considered to determine the parameters of line spacing, flight altitude and speed. Then, the data collection can be carried out.
Step 2:	Data of each flight profile are processed according to the improved VMD method described in Section 2.1.2, then the magnetic map of the survey area is obtained by interpolation.
Step 3:	Euler deconvolution is used to obtain the preliminary estimation of target position and magnetic moment, and these results are then used as the initial value of the LM algorithm to obtain more accurate target parameters.
Step 4:	The result of data interpretation is finally evaluated according to the real situation, and the follow-up work (e.g., the clearance of the target) can be carried out.

4. Field Experiments and Analysis

In this section, field experiments are carried out to prove the feasibility of the proposed method. Firstly, the UAV magnetic survey system is briefly introduced. Then, the collection of UAV magnetic data is described in detail. After that, the data processing results are compared with several common schemes. Finally, the parameter estimation of the target is obtained based on the joint Euler deconvolution-LM method, and the effectiveness of the proposed method is verified according to the real situation.

4.1. UAV Magnetic Survey System

A six-rotor-based UAV magnetic survey system has been deployed for the purpose of near-surface targets detection, as shown in Figure 2. The detailed parameters of the system can be found in [16], which is provided in the Supplementary Materials. This system consists 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 of the UAV, a data acquisition module, and a power module. The two OPMs are rigidly mounted below the center of the UAV by a boom, with a vertical distance of 0.45 m. The TMI and VMI data are synchronized by the pulse per second signal with a sampling frequency of 160 Hz.

4.2. The Collection of UAV Magnetic Data

As we mentioned in Section 3, the line spacing, flight altitude and speed are determined according to a variety of factors, e.g., the task purpose, target size, topography of the survey area, and weather conditions. Considering that the survey area for our experiment is quite flat, a constant flight altitude is adopted. It should be noted that in areas with complex terrain, it is necessary to use a radar (or laser) altimeter with the function of terrain following. The data used in Section 4.3 are obtained from a maneuvering experiment with a flight altitude of 20 m above ground level (AGL), and a speed of about 4 m/s. For the near-surface targets detection application, the flight altitude and line spacing are usually determined through field experiments. For the field experiment proposed in Section 4.4, the flight altitude is set to 2 m AGL, with a line spacing of 0.75 m and a speed of about 2 m/s. The sampling interval along the profile is about 1 cm. The data acquisition module collects differential GPS and magnetic data at the same time and stores them in the SD card.

4.3. UAV Magnetic Data Processing

Data containing 6000 sampling points can be obtained after the UAV maneuvers back and forth along a preset path; the flight altitude of the UAV is about 20 m. The original data and its time-frequency spectrogram after DC component removed are shown in Figure 3.

Three features can be observed from Figure 3: (1) the signal-to-noise ratio (SNR) of the original data is low due to the influence of many kinds of noise; (2) there are two kinds of noise which are obviously different from the target signal in frequency, including the power frequency interference (energy concentrated in 50 Hz) and the interference generated by the platform, the energy of which is distributed in the frequency band of 30-40 Hz; and (3) there are some periodic signals in the frequency band of 0.5–3 Hz, which should be related to the maneuver of the platform. To make a better comparison, EMD is first applied to the original data, 10 IMFs and their corresponding spectrum are obtained, as shown in Figure 4.

As can be seen from Figure 4, the boundary effects and mode mixing are inevitable in EMD. In addition, how to reconstruct signal is a key problem because the physical meaning of the obtained IMFs is not obvious; this uncertainty is caused by the weak mathematical theory of the EMD method itself. Therefore, the VMD method is expected to achieve better results by constructing and solving variational problems.

According to Figure 3b and Figure 4, the decomposition modes number K of VMD is set from 4 to 7, and the penalty factor $\alpha$ is set to 2000. The center frequencies of IMFs obtained and the energy loss coefficient $e_l$ with different K are shown in

Table 1. The center frequencies of IMFs and the energy loss coefficient with different K.

K	Center Frequency/Hz						Energy Loss Coefficient
4	1.45	34.45	37.32				0.2534
5	1.43	18.81	34.49	37.74			0.2469
6	1.41	17.86	34.02	36.99	50.04		0.0980
7	1.41	17.63	32.37	35.09	37.27	50.04	0.0939

. The center frequency of IMF1 is close to DC, so it is no longer listed in Table 1.

It can be found that, when K < 6, the energy loss coefficient is greater than 0.1, which means that the signal has not been completely decomposed. When K increases from 6 to 7, the center frequency of the new IMF is very close to that of the previous IMFs, and the change of the energy loss coefficient is very small, which means that over-decomposition may occur at this time. To further confirm the best K, the correlation coefficients between the new IMF (decomposition number = K) and the previous IMFs (decomposition number = K − 1) are calculated according to the method proposed in Section 2.1.2, and the results are shown in

Table 2. The correlation coefficients between the new IMF and the previous IMFs.

K	New IMF	Correlation Coefficients
5	IMF 3	0.0031	0.0262	0.0474	0.0294
6	IMF 6	0.0013	0.0113	0.0783	0.1103	0.0821
7	IMF 4	0.0004	0.0030	0.0418	0.6327	0.0832	0.0039

.

When K = 7, IMF4 is a new frequency component, and the correlation coefficient between it and the corresponding IMF4 when K = 6 is 0.6327, greater than the preset threshold, which means that the over-decomposition occurs. Therefore, the optimal decomposition modes number is set to 6. The corresponding decomposition result is shown in Figure 5.

The following conclusions can be drawn from Figure 5:

• VMD can significantly reduce the number of IMFs, and each IMF has a relatively clear physical meaning.
• The center frequency of IMF1 is close to DC, which needs to be paid more attention to in the magnetic survey because the anomaly signal caused by the target is also quasi-static.
• The energy of IMF2 is distributed in the frequency band of 0.5-3 Hz, which is the interference related to the maneuver of the UAV.
• IMF3-IMF5 is the interference produced by the UAV platform, which may be related to the airborne electronic equipment.
• The center frequency of IMF6 is 50 Hz, which is the power frequency interference.

To further verify the effectiveness of the proposed VMD method, a low-pass filter [14], an improved EMD method, and a complete ensemble empirical mode decomposition with adaptive noise (CEEMDAN) method [16], are used to process the original signal shown in Figure 3a. PE and the processing time are used as evaluation indicators, i.e., the smaller the PE is, the more obvious the interference suppression effect is; the shorter the processing time is, the better the real-time performance is. Results of PE and the processing time obtained by the four methods are given in

Table 3. PE and the processing time obtained by the four different methods.

Method	Low-Pass Filter	EMD	CEEMDAN	Proposed
PE	0.4362	0.4414	0.3966	0.3091
Processing time/s	0.112	1.008	10.145	4.305

. Although the proposed method is inferior to the low-pass filter and EMD-based method in processing time, it achieves the best effect of interference suppression. In addition, compared with the CEEMDAN method, the proposed method has shorter processing time, which means that it has better real-time performance.

4.4. UAV Magnetic Data Interpretation

The field experiment was conducted in Hunan Province, China, on 12 October 2021. A cylindrical steel pipe is buried in the center of the survey area, with a depth of about 0.58 m. The pre-programmed flight profiles run along the north-south direction, with a line spacing of 0.75 m. The actual flight trajectory and the relative position of the target in the survey area are shown in Figure 6; due to the influence of crosswind, the actual line spacing is 0.2–1.2 m. The flight altitude of the UAV is 2 m, with an airspeed of about 2 m/s.

The magnetic map of the survey area obtained by the original data is shown in Figure 7a; the characteristics of the target are masked due to the low SNR of the original data. The proposed VMD-based method is applied to process the original data, then the magnetic map of the survey area is obtained by interpolating the processed data, as shown in Figure 7b. It can be clearly found that the quality of the processed data has been greatly improved compared with the original data, which lays a foundation for the further interpretation.

The magnetic anomaly area of Figure 7b is selected, then Euler deconvolution is used to obtain a preliminary estimate of the position and magnetic moment of the target. Finally, these parameters are taken as the initial values of the LM algorithm, and the further estimation of the parameters of the target is obtained.

Table 4. Results of estimated target parameters obtained by two methods.

Parameter	Real	Euler-Deconvolution	Joint Euler-LM Method
Location (m)	(21.802, 21.964, −0.580)	(21.806, 22.117, −0.742)	(21.791, 21.925, −0.526)
Magnetic moment (Am2)	-	(−0.044, 0.534, −1.113)	(−0.106, 0.630, −1.235)

shows the real position of the target and the estimated parameters obtained by the above two different methods, in which the real position of the target is obtained by the differential real-time kinematic (RTK) with centimeter-level accuracy. According to Table 4, the target position estimated by the joint Euler deconvolution-LM method is closer to the real value compared with the Euler deconvolution method, the horizontal positioning error is reduced from 15.31 cm to 4.05 cm, and the depth estimation error is reduced from 16.2 cm to 5.4 cm.

To further evaluate the convergence of the inversion model, the coefficient of determination (also known as R-squared) is introduced. By comparing the observed and predicted values of magnetic anomaly, R-squared can be defined as follows:

$$R^2=1-\frac{{{\displaystyle \sum }}_i{\left(B_{pred,i}-B_{obs,i}\right)}^2}{{{\displaystyle \sum }}_i{\left(B_{pred,i}-\widehat{B_{obs}}\right)}^2},$$

(19)

where $i\in \left[1,N\right]$, and N is the number of observed points. $B_{obs,i}$ and $B_{pred,i}$ are the observed and predicted values of magnetic anomaly at the ith point, and $\widehat{B_{obs}}$ is the mean value of $B_{obs,i}$. By combining the Euler deconvolution and the LM algorithm, the value of $R^2$ reaches 0.9736, which is significantly better than the 0.8749 based on the Euler deconvolution method. It is proved that the proposed method can obtain more realistic inversion results, which is very important for further data interpretation.

Although some global optimization algorithms such as DE can also gain accurate estimation of the target parameters, these methods usually take more time than the joint Euler deconvolution-LM method proposed in this paper. Figure 8 shows the variation of the norm of residual error with iterations when using DE and the proposed method, respectively. For the sake of fairness, the search space of DE is set to be relatively compact, and its upper and lower bounds are ${\mathit{x}}_{max}$= [23, 23, 1.2, 0.8, 0.8, 2] and ${\mathit{x}}_{min}$= [21, 21, −0.5, −0.8, −0.8, −2], respectively. The population size is set to be 60, and the mutation scale factor F and the crossover rate CR are set to be 0.6 and 0.4, respectively.

As can be seen from Figure 8, the results of the two methods are very similar; the DE method obtains a stable solution after about 80 iterations, while the proposed method completes the process in less than 6 iterations. Under the same computer configuration, the DE method spends 46.69 s to complete the task; for comparison, the proposed method takes 1.39 s, which shows excellent efficiency.

5. Final Remarks

This paper presents a complete workflow for the processing and interpretation of UAV magnetic data, which mainly includes two parts: (1) UAV magnetic data processing based on the improved VMD method, in which the decomposition modes number K is adaptively determined by mode characteristics; and (2) UAV magnetic data interpretation based on Euler deconvolution and LM algorithm, where the target parameters are first estimated by Euler deconvolution, and then used as priori information of the LM algorithm. The results of field experiments show that the proposed method can be applied to UAV magnetic data processing, and the estimated target position is closer to the real situation.

In Section 4.3, the improved VMD method is applied to the UAV magnetic data processing. It is worth noting that, when the decomposition modes number K is appropriate, the result of VMD has a clear physical meaning, which is not available in other methods. We believe that this discovery is of great significance for the UAV magnetic interference mapping and suppression and is worthy of in-depth study. Section 4.4 establishes a step-by-step data interpretation method based on Euler deconvolution and LM algorithm, which makes full use of the local optimization and the fast convergence of the LM algorithm. The horizontal positioning error of the target is reduced to less than 10 cm, which is very beneficial to the follow-up work, such as the clearance of targets.

Nevertheless, there are some drawbacks in this paper. For example, this paper only optimizes the decomposition modes number K of VMD, and better results may be obtained by considering the influence of penalty factor $\alpha$. Moreover, it should be noted that the proposed data interpretation method is only suitable for an isolated target, and if multiple targets exist or the distance between targets is very close, it is recommended to consider gradient or gradient tensor data [33,40]. In addition, due to the limited information of the collected magnetic anomaly data, more detailed parameters (e.g., shape, attitude, and material) of the target cannot be obtained. More comprehensive information about the target can be obtained by enriching sensor data, and there have been some successful attempts in this regard [41,42,43].

Our future work will focus on these issues, and mainly includes the following three parts: 1. Further optimize the VMD method to get better data processing results; 2. To improve the spatial resolution of the target through gradient measurement; 3. To obtain more abundant target parameters by combining other types of sensors (e.g., hyperspectral, electromagnetic) for data interpretation.