Aeromagnetic data processing method of multi-rotor UAV based on improved empirical mode decomposition

Abstract

The present invention relates to a multi-rotor unmanned aerial vehicle aeromagnetic data processing method based on improved empirical mode decomposition, comprising the following steps: step (1) performing empirical mode decomposition on an original signal x(t) to obtain N IMF components; step (2) calculating the correlation coefficient Coef(i) of each IMF component with the original signal x(t), i=1,2,...,N; step (3) calculating the permutation entropy value PE(i) of each IMF component, i=1,2,...,N, and recording the median value of PE(i) as Median PE ; Step (4) to judge, if PE(i)≤Median PE , then keep Coef(i) unchanged, otherwise set Coef(i) = 0, and get the updated Coef(i) sequence; Step (5) calculates the joint coefficient C c (i), i = 1, 2, ..., N, C c (i)=[(1‑PE(i))+Coef(i)] / 2, let PE(i)≤Median PE When true, the corresponding C c (i) The median of the sequence is denoted as Median cc ; Step (6) to judge, if C c (i)≥Median cc , the corresponding IMF components are retained, otherwise the corresponding IMF components are discarded; step (7) reconstructs the signal s(t) using the remaining IMF components.

Description

Technical Field

[0001] The invention belongs to the field of geophysical aeromagnetic detection, and in particular relates to a multi-rotor unmanned aerial vehicle aeromagnetic data processing method. Background Art

[0002] In recent years, drone-based aeromagnetic detection systems have developed rapidly and attracted widespread attention. Compared with traditional aerial magnetic survey or ground magnetic survey technology, drone aeromagnetic systems have the advantages of high efficiency, good safety, low cost, and high data resolution. Most of the existing drone aeromagnetic detection systems are based on helicopters or fixed-wing aircraft. For example, the Chinese patent "A drone aeromagnetic detection system and method" (patent number CN104122597A) has problems such as large interference or the need for a runway in actual application scenarios. Multi-rotor drones have the advantages of low interference, no need for take-off and landing runways, and simple operation. In recent years, they have been used in fields such as geophysical exploration. For example, in 2020, C. Walter et al. used the DJI S900 multi-rotor drone to conduct mineral exploration in Ontario, Canada. See C. Walter, A. Braun, and G. Fotopoulos, “High-resolution unmanned aerial vehicle aeromagnetic surveys for mineral exploration targets,” Geophys. Prospect., vol. 68, no. 1, pp. 334–349, 2020.

[0003] Due to technical and other reasons, most literature does not describe the processing methods of multi-rotor UAV aeromagnetic data. The lack of effective processing and interpretation of data has become an obstacle to the promotion of the application of multi-rotor UAV aeromagnetic systems.

[0004] At present, there is no literature that uses improved empirical mode decomposition for multi-rotor UAV aeromagnetic data processing. Pan Guangzhen et al. determined the signal and noise components based on the Mahalanobis distance and used an optimization algorithm to determine the threshold of the noise component. The algorithm is relatively complex and not conducive to real-time processing (Pan Guangzhen, Wang Feng, Sun Yanqing. Application of improved EEMD algorithm in ECG signal denoising [J]. Data Acquisition and Processing, 2018, 33(4): 646-653.); Wang Yujing et al. used kurtosis value and correlation coefficient to determine the IMF component used to reconstruct the signal (Wang Yujing, Kang Shouqiang, Zhang Yun, Liu Xue, Jiang Yicheng, Mikulovich VI. Rolling bearing state identification method based on ensemble empirical mode decomposition sensitive intrinsic mode function selection algorithm [J]. Journal of Electronics and Information Technology, 2014, 36(3): 595-600.), which performs poorly at low signal-to-noise ratio. Summary of the invention

[0005] In order to solve the above technical problems, the present invention proposes a multi-rotor UAV aeromagnetic data processing method based on improved empirical mode decomposition. The multi-rotor UAV aeromagnetic data processing method based on improved empirical mode decomposition proposed in the present invention comprehensively considers the correlation coefficient and the permutation entropy value. The algorithm is completely driven by data and performs well even at a low signal-to-noise ratio.

[0006] The object of the present invention is to provide a multi-rotor UAV aeromagnetic data processing method based on improved empirical mode decomposition to solve the problems existing in the above-mentioned background technology.

[0007] To achieve the above object, the technical solution provided by the present invention is a multi-rotor UAV aeromagnetic data processing method based on improved empirical mode decomposition, comprising the following steps:

[0008] (1) Perform empirical mode decomposition on the original signal x(t) to obtain N IMF components;

[0009] (2) Calculate the correlation coefficient Coef(i) of each IMF component and the original signal x(t), i = 1, 2, ..., N;

[0010] (3) Calculate the permutation entropy value PE(i) of each IMF component, i = 1, 2, ..., N, and record the median of PE(i) as Median PE ;

[0011] (4) Make a judgment, if PE(i)≤Median PE , then keep Coef(i) unchanged, otherwise set Coef(i)=0 to obtain the updated Coef(i) sequence;

[0012] (5) Calculate the combined coefficient C c (i),i=1,2,…,N,C c (i)=[(1-PE(i))+Coef(i)] / 2, put “PE(i)≤Median PE " is true when the corresponding C c (i) The median of the sequence is denoted as Median cc ;

[0013] (6) Make a judgment, if C c (i)≥Median cc , then the corresponding IMF component is retained, otherwise the corresponding IMF component is discarded;

[0014] (7) Reconstruct the signal s(t) using the remaining IMF components.

[0015] The original signal described in step (1) specifically refers to the magnetic sensor data signal, including the optical pump magnetometer data signal and the fluxgate magnetometer data signal. The step of performing empirical mode decomposition on the original signal described in step (1) is to find all the extreme points of the original signal, and use cubic spline interpolation to fit the upper and lower envelopes; calculate the mean of the upper and lower envelopes, and subtract the mean envelope from the original signal to obtain the intermediate signal; if the intermediate signal meets the IMF condition, then the intermediate signal is an IMF component, if not, repeat the above process; obtain the IMF 1 (t), the IMF is subtracted from the original signal 1 (t), as the new original signal, repeat the above process until the empirical mode decomposition of the original signal is completed. The intermediate signal satisfies the IMF condition if: (a) in the entire data segment, the number of extreme points and the number of zero-crossing points must be equal or the difference does not exceed 1, and (b) at any time, the upper and lower envelopes are locally symmetrical relative to the time axis.

[0016] The correlation coefficient in step (2) refers to the Pearson correlation coefficient, which is defined as Where cov(IMF i (t),x(t)) refers to the covariance between the ith IMF component and the original signal, σ(IMF i (t)) and σ(x(t)) refer to the variance of the ith IMF component and the original signal, respectively.

[0017] The correlation coefficient Coef(i) and permutation entropy value PE(i) in steps (2) and (3) are normalized respectively.

[0018] Beneficial effects:

[0019] In view of the deficiencies of the prior art, the present invention combines the correlation coefficient and the permutation entropy value to screen the IMF components dominated by the effective components in the original signal. The algorithm is completely data-driven, does not require prior knowledge, and does not require manual setting of thresholds. It has the characteristics of simple implementation. In the case of low signal-to-noise ratio, the processing result can also significantly reduce the influence of noise and has good robustness. BRIEF DESCRIPTION OF THE DRAWINGS

[0020] Figure 1 This is a flow chart of a multi-rotor UAV aeromagnetic data processing method based on improved empirical mode decomposition according to the present invention;

[0021] Figure 2 It is the original data signal with high signal-to-noise ratio;

[0022] Figure 3 The eight IMF components are obtained by empirical mode decomposition of the original data signal with high signal-to-noise ratio;

[0023] Figure 4 The signal is reconstructed by using the last two IMF components for the original data signal with high signal-to-noise ratio;

[0024] Figure 5 It is the original data signal with low signal-to-noise ratio;

[0025] Figure 6 The eight IMF components are obtained by empirical mode decomposition of the low signal-to-noise ratio original data signal;

[0026] Figure 7 The signal reconstructed by directly using the IMF1 and IMF8 components screened by the joint coefficients;

[0027] Figure 8 This is the signal reconstructed from the last two IMF components screened out by the method proposed in the present invention. DETAILED DESCRIPTION

[0028] The technical scheme of the present invention is described in detail below in conjunction with the embodiments and drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work belong to the scope of protection of the present invention.

[0029] The basic idea of ​​the present invention is to use the permutation entropy value to adaptively adjust the correlation coefficient between the IMF component and the original signal, and to select the IMF component dominated by the useful signal by constructing the joint coefficient to reconstruct the signal. The complete process is shown in Figure 1 .

[0030] Example 1

[0031] In this example, the processing of aeromagnetic data of a multi-rotor UAV under high signal-to-noise ratio is considered. For the original signal x(t), empirical mode decomposition is performed on it to obtain:

[0032]

[0033] Among them, c i (t) is the IMF (intrinsic mode function), which contains the local characteristics of the original signal at different time scales, r(t) is the remainder, and N is the number of IMF components obtained by performing empirical mode decomposition on the original signal.

[0034] The correlation coefficient Coef(i) between each IMF and the original signal x(t) is obtained by formula (2):

[0035]

[0036] Among them, cov(IMF i(t),x(t)) is the covariance between the ith IMF component and the original signal, σ(IMF i (t)) and σ(x(t)) are the variances of the ith IMF component and the original signal, respectively.

[0037] Permutation entropy can be used to measure the complexity of time series data. Its mathematical expression is: for a time series {X i , i=1,2,…,P} to reconstruct the phase space and obtain:

[0038]

[0039] Wherein, m is the embedding dimension, τ is the delay time, K=P-(m-1)τ, each row in the matrix Y is a reconstruction component, and there are K reconstruction components in total.

[0040] Rearrange each reconstruction component in ascending order of value, and for the j-th reconstruction component:

[0041] X[i+(j 1 -1)τ]≤X[i+(j 2 -1)τ]≤…≤X[i+(j m -1)τ] (4)

[0042] Among them, j 1 ,j 2 ,…,j m Indicates the index of the column where each element in the reconstructed component is located.

[0043] Reconstructing each row of the matrix yields a sequence of symbols:

[0044] S(l)=(j 1 ,j 2 ,…,j m ),l=1,2,…,k,k≤m! (5)

[0045] There are m! different symbol sequences in the m-dimensional phase space mapping. If the probability of each symbol sequence appearing is recorded as P 1 ,P 2 ,…,P K According to the definition of information entropy, the time series X i The permutation entropy of is defined as:

[0046]

[0047] Considering that the maximum value of permutation entropy is ln(m!), the permutation entropy value is normalized:

[0048]

[0049] The size of the permutation entropy value indicates the degree of randomness of the time series: the smaller the entropy value, the simpler and more regular the time series; conversely, the more complex and random the time series is.

[0050] When conducting aeromagnetic detection with a multi-rotor drone, the data with target survey lines can be considered as a high signal-to-noise ratio situation. The raw data with high signal-to-noise ratio is as follows: Figure 2 As shown, the original data signal is subjected to empirical mode decomposition to obtain 8 IMF components, such as Figure 3 As shown, in order to compare the effects of the present invention, the correlation coefficients and permutation entropies of the IMF components and the original signals are first calculated according to equations (2) and (6). The results are shown in Table 1.

[0051] Table 1 Correlation coefficients and permutation entropy of each IMF component of the original signal with the original signal at high signal-to-noise ratio

[0052] IMF 1 2 3 4 5 6 7 8 Coef(i) 0.3297 0.0103 0.0163 0.0422 -0.0022 0.0066 0.5921 0.5149 PE(i) 0.9974 0.9863 0.9765 0.9815 0.8248 0.6549 0.4375 0.4223

[0053] Two different methods are discussed. Method A uses formula (8) to directly calculate the joint coefficient C of each IMF component. c (i) The IMF components corresponding to the first 25% of the joint coefficients arranged from large to small are taken to reconstruct the signal; Method B according to the present invention Figure 1 The process shown calculates the joint coefficient C of each IMF component c (i) Take the IMF components that meet the conditions to reconstruct the signal:

[0054]

[0055] The joint coefficients of each IMF component obtained by two different methods are shown in Table 2.

[0056] Table 2 Joint coefficients of IMF components obtained by two different methods at high signal-to-noise ratio

[0057] IMF 1 2 3 4 5 6 7 8 <![CDATA[C c (i)-A]]> 0.1662 0.0120 0.0199 0.0304 0.0865 0.1759 0.5773 0.5463 <![CDATA[C c (i)-B]]> 0.0013 0.0069 0.0117 0.0092 0.0865 0.1759 0.5773 0.5463

[0058] It can be seen from Table 2 that when the signal-to-noise ratio is high, whether the joint coefficient is calculated directly or the joint coefficient is calculated using the method proposed in the present invention, the results are consistent. The last two IMF components are retained and used to reconstruct the signal. The results are as follows: Figure 4 As shown, it can be seen that the noise in the original signal is significantly suppressed.

[0059] Example 2

[0060] This example considers the processing of aeromagnetic data of a multi-rotor UAV under low signal-to-noise ratio conditions. The original data signal is as follows: Figure 5 As shown. Performing empirical mode decomposition on the original signal, 8 IMF components are obtained. The results are as follows Figure 6According to equation (2) and equation (6), the correlation coefficients of each IMF component and the original signal and their permutation entropy are obtained. The results are shown in Table 3.

[0061] Table 3 Correlation coefficients and permutation entropy of each IMF component of the original signal with the original signal at low signal-to-noise ratio

[0062] IMF 1 2 3 4 5 6 7 8 Coef(i) 0.6150 0.0286 0.0055 0.0143 0.0405 0.0152 0.0430 0.7675 PE(i) 0.9973 0.9992 0.9511 0.7978 0.7428 0.5972 0.4907 0.0010

[0063] According to the steps in Example 1, the results of directly calculating the combination coefficient and calculating the combination coefficient according to the method proposed by the present invention are respectively obtained, as shown in Table 4.

[0064] Table 4 Joint coefficients of IMF components obtained by two different methods at low signal-to-noise ratio

[0065] IMF 1 2 3 4 5 6 7 8 <![CDATA[C c (i)-A]]> 0.3089 0.0147 0.0272 0.1082 0.1489 0.2090 0.2762 0.8823 <![CDATA[C c (i)-B]]> 0.0014 0.0001 0.0245 0.1011 0.1489 0.2090 0.2762 0.8823

[0066] When the signal-to-noise ratio is low, if the IMF components corresponding to the first 25% of the joint coefficients are directly taken from large to small, IMF1 and IMF8 are selected, but the main component of IMF1 is 50Hz power frequency noise, which is not an effective signal component. The reconstructed signal result is shown in Figure 7 According to the method proposed by the present invention, IMF7 and IMF8 are selected for reconstructing the signal, and the results are shown in Figure 8 When the signal-to-noise ratio is low, the original signal is seriously contaminated by noise, resulting in a certain similarity between IMF1 and the original signal. At this time, if the joint coefficient is directly calculated, it is easy to select the IMF component of the noise-dominated mode. The algorithm of the present invention uses the permutation entropy value to reweight the correlation coefficient. In the entire algorithm process, there is no need to manually set parameters such as thresholds. The entire process is completely driven by data. Figure 8 and Figure 7 The comparison results show that the method proposed in the present invention also has a good processing effect on low signal-to-noise ratio raw data.

[0067] Matters not covered by the present invention are known technologies.

[0068] The above embodiments are only for illustrating the technical concept and features of the present invention, and their purpose is to enable people familiar with the technology to understand the content of the present invention and implement it accordingly, and they cannot be used to limit the protection scope of the present invention. Any equivalent changes or modifications made according to the spirit of the present invention should be included in the protection scope of the present invention.

Claims

1. A multi-rotor UAV aeromagnetic data processing method based on improved empirical mode decomposition, characterized in that: The following steps are involved: Step (1) The original signal Perform empirical mode decomposition and obtain IMF components; Step (2) Calculate the IMF components and the original signal Correlation coefficient ; Step (3) Calculate the permutation entropy value of each IMF component ,Will The median value of ; Step (4) is to make a judgment. If , then keep unchanged, otherwise , get the updated sequence; Step (5) Calculate the joint coefficient , ,Will When true, the corresponding The median of the sequence is denoted by ; Step (6) is to judge if , then the corresponding IMF component is retained, otherwise the corresponding IMF component is discarded; Step (7) Reconstruct the signal using the remaining IMF components ; The original signal in step (1) refers to a magnetic sensor data signal, including an optical pump magnetometer data signal and a fluxgate magnetometer data signal.

2. The method for processing aeromagnetic data of a multi-rotor UAV based on improved empirical mode decomposition according to claim 1, characterized in that: The step of performing empirical mode decomposition on the original signal in step (1) is as follows: Step (1.1), find all extreme points of the original signal, and use cubic spline interpolation to fit the upper and lower envelopes; Step (1.2), find the mean of the upper and lower envelopes, and subtract the mean envelope from the original signal to get the intermediate signal; Step (1.3), if the intermediate signal satisfies the IMF condition, then the intermediate signal is an IMF component. If not, repeat the above process for it. The intermediate signal satisfies the IMF condition if: (a) in the entire data segment, the number of extreme points and the number of zero-crossing points must be equal or the difference does not exceed 1, (b) at any time, the upper and lower envelopes are locally symmetrical with respect to the time axis; Step (1.4), get the first IMF component, denoted as , using the original signal minus As the new original signal, repeat the above process until the empirical mode decomposition of the original signal is completed.

3. The method for processing aeromagnetic data of a multi-rotor UAV based on improved empirical mode decomposition according to claim 1, characterized in that: The correlation coefficient in step (2) refers to the Pearson correlation coefficient, which is defined as: ， in, Refers to the The covariance of the IMF components and the original signal, and Respectively refer to The variance of the IMF components and the original signal; the correlation coefficient of steps (2) and (3) and permutation entropy Normalize them separately.

Images