High-Sensitivity Real-Time Tracking System for High-Speed Pipeline Inspection Gauge

Abstract

Real-time tracking of pipeline inspection gauges (PIGs) is an important aspect of ensuring the safety of oil and gas pipeline inline inspections (ILIs). Transmitting and receiving extremely low frequency (ELF) magnetic signals is one of the preferred methods of tracking. Due to the increase in physical parameters of the pipeline including transportation speed, wall thickness and burial depth, the ELF magnetic signals received are short transient (1-second duration) and very weak (10 pT), making the existing above-ground-marker (AGM) systems difficult to operate correctly. Based on the short transient very weak characteristics of ELF signals studied with a 2-D finite-element method (FEM) simulation, a data fusion model was derived to fuse the envelope decay rates of ELF signals by a least square (LS) criterion. Then, a fast-decision-tree (FDT) method is proposed to estimate the fused envelope decay rate to output the maximized orthogonal signal power for the signal detection through a determined topology and a fast calculation process, which was demonstrated to have excellent real-time detection performance. We show that simulation and experimental results validated the effectiveness of the proposed FDT method, and describe the high-sensitivity detection and real-time implementation of a high-speed PIG tracking system, including a transmitter, a receiver, and a pair of orthogonal search coil sensors.

1. Introduction

Today, over 70 percent of oil and gas transportation relies on a total of over 3,500,000 km of pipelines in 120 countries, which is the safest, most reliable and efficient method, compared with highways, railways and marine transport, and so on [1]. To meet the rapidly growing demand for oil and gas, the construction of trunk pipelines which can operate in high-pressure, high-speed and high-capacity conditions has been increasing in recent years, giving the pipelines larger diameters, thicker pipe walls, deeper burial depths and larger transportation distances [2,3,4,5]. Due to corrosion and high-pressure loadings on the pipe wall, pipeline failures including defects, deformations, blockages, leakages and explosions occur every year, resulting in serious risks to pipeline safety [6,7,8]. In China, it was reported that there was a pipeline failure every 4 km [9], and in the USA, the economic and human losses due to the significant pipeline incidents in the past 20-year period are $7 billion and 330 fatalities along with 1292 injuries, respectively [10]. Pipeline inspection gauges (PIGs) are periodically used for pipeline in-line inspections (ILIs) to minimize these risks and maintain pipeline integrity [11,12,13,14,15,16]. However, during the pipeline ILI operation, the PIG is likely to be blocked by the deformed pipelines, while an odometer installed in the PIG can only record the PIG position offline, because the shielding effect of the steel pipe wall makes it difficult to transmit the position information from inside the pipeline to the ground through wireless communication technologies, such as Wi-Fi, Bluetooth, or 3G and 4G networks [17,18]. Thus, the pipeline ILI usually needs the use of above-ground-marker (AGM) technologies to track the position of the PIG in real time [19,20].

Existing AGM methods utilize acoustic wave, static magnetic field, extremely low frequency (ELF) magnetic fields for PIG localization and tracking, and have been widely practiced in the oil and gas industry [21,22,23,24,25,26,27]. Here, the working principle of the AGM system is shown in Figure 1. The PIG is propelled by the flowing medium and the speed is mainly determined by the pressure inside the pipeline. A receiver placed above the ground detects the tracking signals generated by the PIG moving right below, and records the PIG position simultaneously. With the help of the global position system (GPS) and general packet radio service (GPRS) technologies, all receivers along the pipeline usually with an interval of 1 km can send the current location and local time of the PIG to the monitoring center to form a tracking network, and the information will also be used to compensate for the accumulated locating error caused by the PIG odometer.

The acoustic wave method is used to detect the acoustic wave generated by the friction between the PIG and the pipe wall. In [21,22], it is used to improve the tracking precision of the underground PIG, and in [23], a micro-electro-mechanical-system (MEMS) acoustic vector sensor is found to have advantages of high sensitivity for detecting low-frequency acoustic signals. However, the acoustic signals generated by the moving PIG are easily cluttered by the environmental acoustic sources, such as moving vehicles, trains and airplanes, and it produces a high number of false alarms. The static magnetic field method is used to detect the static magnetic field generated by the permanent magnets assembled in the PIG. In [24,25], experimental results showed that this method was effective when the magnetic signals emitted by the permanent magnets were strong enough. Furthermore, it is only applicable to certain types of PIG due to the requirement to carry heavy and bulky permanent magnets, thus limiting the scope of application. To overcome the aforementioned challenges, the ELF method is introduced, which is a single-frequency magnetic signal transmission and reception method that is suitable for most types of PIG [26,27]. A small and lightweight transmitter equipped in the PIG generates the ELF magnetic field and the receiver detects the ELF magnetic field using magnetic sensors and real-time signal detection methods. Comparing to high-frequency magnetic field, the ELF magnetic field with frequencies between 22 Hz and 24 Hz has a good penetrating ability for signal transmitting through steel pipe wall and soil, and is less susceptible to environmental interferences [28,29]. In [30], the tracking and localization of the in-pipe PIG is realized by applying the equivalent magnetic dipole model and the inverse calculation of the ELF magnetic field. In [27,28,30], experimental results showed that the existing ELF tracking systems worked well for long-duration and high-magnitude ELF magnetic signals. However, due to the higher speed of the PIG during the pipeline ILI operation, the received time-domain ELF signals become a transient signal with a very short duration [29]. Due to the thicker pipe wall and deeper burial depth of pipelines, as well as the desired low power consumption requirement of the transmitter for the long working time, the magnitude of the received ELF magnetic field is very low and immersed in strong background noise, leading to a very poor signal-to-noise ratio (SNR) [31]. The increases of the speed, pipe wall thickness, and burial depth means the existing ELF tracking systems still face the challenges of detecting the short transient very weak ELF magnetic signals.

According to the signal detection theory, solving these challenges requires the use of high-sensitivity magnetic sensors with advanced signal detection methods, so that the probability of detection (POD) of the short transient very weak ELF signals can be high enough for a given low probability of false alarm (PFA). There are several high-sensitivity magnetic sensors that can be used to detect very weak magnetic signals, such as the search coil sensor, fluxgate magnetometer and SQUID sensor [32,33]. A search coil sensor can detect magnetic fields as weak as 10⁻² pT with a frequency range from 1 Hz to 1 MHz. A fluxgate magnetometer can detect magnetic fields from 10 to 10⁹ pT, while the frequency response is limited by the excitation field and response time of ferromagnetic materials. A SQUID sensor is more sensitive at low frequencies (<1 Hz) and the sensitivity is as low as 10⁻² pT. However, the need for liquid-helium coolant makes the SQUID sensor heavy and inconvenient. As for the signal detection, the generalized likelihood ratio (GLR) is considered to be one of the typical signal detection methods [34,35]. Poor et al. [36] develops the quickest change detection method based on the distribution of test statistics. Wavelet transform becomes an effective tool for joint time-frequency analysis and has a broader application in the diagnosis of mechanical faults and the evaluation of non-destructive testing (NDT) signals [37,38,39]. Marius Birsan et al. [40] presents an underwater detection of ELF electromagnetic field signals from ships using the cross wavelet transform and wavelet coherence method for high SNR situations. Qin et al. [41] present transient weak fault feature extraction based on an optimized Morlet wavelet transform and kurtosis index. Meanwhile, a maximum likelihood (ML) estimator is also developed to extract the amplitude from fault-generated transient impulses in [42].

This paper proposes a high-sensitivity real-time tracking system for high-speed PIG to detect short transient very weak ELF magnetic signals. A 2-D finite-element method (FEM) simulation was conducted to understand the relationships between the time-domain ELF signals and the physical parameters of the pipeline, including transportation speed, wall thickness, and burial depth. According to the short transient very weak characteristics of ELF signals, a data fusion model was derived to well represent the ELF signals by fusing the envelope decay rates through a least square (LS) criterion. Then, a fast-decision-tree (FDT) method is proposed to detect the ELF signals. The main differences of the FDT method are: (1) it first estimates the fused envelope decay rate to maximize the orthogonal signal power of ELF signals, and the maximized orthogonal signal power is used as the test statistic for the signal detection, which provides an excellent detection performance; and (2) the deterministic topology and fast calculation process of the FDT method avoid massive matrix inversion and multiplication, which significantly reduces the computation costs and provides excellent real-time performance. The detection performance of the proposed FDT method is discussed, with a comparison between other existing detection methods as well as theoretical bounds through Monto-Carlo simulations. In consideration of high sensitivity as well as convenience, a pair of orthogonal search coil sensors was designed to receive the very weak ELF magnetic signals. Then, an ELF tracking system including an ELF transmitter, an ELF receiver and a pair of orthogonal search coil sensors was developed and tested, and the experimental results validate the effectiveness of the proposed FDT method. The main reason for the excellent detection performance obtained is discussed through the analyses of the normalized power spectrum, estimation performance of speed and SNR. Finally, the high-sensitivity detection and real-time implementation of high-speed PIG tracking system is achieved.

2. Simulation and Signal models

2.1. 2-D FEM Simulation Studies

In order to study the characteristics of the ELF signals for further developing a targeted signal detection method as well as sensing system, simulation studies of 2-D harmonic magnetic analyses are conducted using ANSYS software, and the 2-D FEM model is shown in Figure 2. The FEM model is axisymmetric and consists of transmitter, oil, steel pipe wall, soil and air. The cylindrical transmitter supported by the PIG’s mechanical connecting structure is located at the symmetric line of the pipeline. The transmitter consists of a transmitting coil and an iron core which has a high permeability to increase the strength of the ELF magnetic field. The transmitting coil is fed a sinusoidal 4 mA (root mean square value) current with frequency of 23 Hz. The simulated speed of the transmitter ranges from 1 m/s to 15 m/s, and the ELF signals in the pipe longitudinal direction (X-axis direction) and pipe radial direction (Y-axis direction) are calculated at the point above the soil surface (location highlighted by “star” in Figure 2). The element and boundary condition of the FEM simulation is set to PLANE53 and flux-parallel condition, respectively. The material parameters including conductivity, permeability and mesh size are listed in

Table 1. Material parameters of 2-D FEM simulation model.

Property	Iron Core	Coil	Oil	Q235 Steel Pipe	Soil	Air
Conductivity (S/m)	1 × 107	5.7 × 107	0.01	2 × 106	0.005	0
Permeability (H/m)	1000 µ0 1	µ0	µ0	500 µ0	µ0	µ0
Mesh size (m)	0.01	0.01	0.1	0.02	0.2	0.2

¹µ₀is permeability of vacuum.

, and the model geometry is listed in

Table 2. Model geometry of 2-D FEM simulation model.

Property	Value
Length of pipeline (L)	45 m
Burial depth (DB)	2.5~5 m
Outer radius of pipeline (R)	254 mm
Thickness of pipe wall (T)	6~15 mm
Inner radius of pipeline	(R-T) mm
Length of coil (LC)	130 mm
Width of transmitting coil (2WC)	59 mm
Width of iron core (2WI)	25 mm
Coil turns	25,000
Transmitting current (RMS)	4 mA
Current frequency	23 Hz
Speed of transmitter (V)	1~15 m/s

.

The simulated ELF signals of the X-axis and Y-axis, when the speed is 15 m/s, pipe wall thickness is 15 mm and burial depth is 5 m, are shown in Figure 3a,b. The envelopes of the X-axis and Y-axis ELF signals with three different speeds are shown in Figure 3c,d. The sampling interval of the 2-D FEM simulation is set to 1 millisecond, and a sample equal to 0 means that the ELF transmitter passes right below the measuring point. Figure 3a,b shows that the number of samples of the ELF signals is 1000 points, corresponding to a signal duration as low as 1 second, and the amplitudes are as low as 15 pT. Figure 3c,d shows that the steepness of the envelopes is closely related to the speed of the transmitter, and the envelopes appear to be short transient as the transmitter moves faster.

The amplitudes of the ELF signals of the X-axis and Y-axis, changing with the pipe wall thickness and burial depth, are shown in Figure 4a,b. The coordinate of thickness and burial depth is set to a linear scale, while that of the magnetic flux density is set to a logarithmic scale, and the amplitudes decrease logarithmically as the thickness and burial depth increase linearly. Taking the X-axis results as an example, compared to the amplitude for thin wall thickness (6 mm) and shallow burial depth (2.5 m) which is about 1000 pT, the amplitudes for thick wall thickness (15 mm) or deep burial depth (5 m) are about 100 pT, which decreases by an order of magnitude. Then, the amplitude for thick wall thickness and deep burial depth is about 10 pT, which decreases by two orders of magnitude. The trend of Y-axis results is the same as for the X-axis. The simulation results discussed above show that under the circumstances of a high speed of 15 m/s, the conditions of a thicker pipe wall and deeper burial depth will cause the ELF signals to become short transient in the envelope and very weak in amplitude.

2.2. Data Fusion Model

The very weak ELF signals are easily interfered with by strong narrow-band background noise. The ELF signals of X-axis and Y-axis with noise can be expressed as:

$$x[n]=s_x[n]+w_x[n] n=-\frac{N_x}{2},\dots ,-1, 0, 1,\dots ,\frac{N_x}{2}$$

(1)

$$y[n]=s_y[n]+w_y[n] n=-\frac{N_y}{2},\dots ,-1, 0, 1,\dots ,\frac{N_y}{2}$$

(2)

where w_x[n] and w_y[n] are narrow-band background noise, N_x and N_y are the number of samples along X-axis and Y-axis, respectively.

Table 3. Notations used in this article.

Notations	Description
Nx	Number of samples of X-axis ELF signal
Ny	Number of samples of Y-axis ELF signal
βxy	Fused envelope decay rate
${\widehat{\beta }}_{\mathit{xy}}$	Estimated fused enveloped decay rate
βmax	Fused envelope decay rate corresponding to maximum speed of transmitter
Hx	Observation matrix of X-axis ELF signal
${\widehat{P}}_{\mathit{xy}}$	Maximized orthogonal signal power
$P_{{\mathrm{H}}_{\mathrm{X}}}$	Energy of X-axis observation vector
$P_{{\mathrm{H}}_{\mathrm{Y}}}$	Energy of Y-axis observation vector
η	Normalized signal energy
Cm	Computation cost of multiplication
Ca	Computation cost of addition
$\widehat{\gamma }$	Threshold for signal detection

summarizes the important notations and their definitions used in this article.

As shown in Figure 3a,b, the time-domain ELF signals have two main properties including envelope and oscillation. Here, a Morlet mother wavelet is used to represent the steepness of the envelope as well as the oscillation:

$$\psi (t)=e^{-\beta t^2/2}\cos \omega t$$

(3)

where β is defined as the envelope decay rate, and the exponential component is used to represent the steepness, while the cosine component is used to represent the single-frequency oscillation. There are a main lobe peak and two side lobe peaks in the X-axis signals. Since the two side lobe peaks are much smaller than the main lobe peak and are more susceptible to noise interference, we set up the X-axis signal model with the main lobe peak:

$$s_x[n]=A_x{e}^{-{\beta }_x{[n/f_s]}^2}\cos [{\omega }_{0}n/f_s+{\phi }_0]$$

(4)

taking the origin and two peak points as cut-off points, the Y-axis signal model is set by four intervals:

$$s_y[n_1,n_2]=\{\begin{array}{l}{A}_y{e}^{- {\beta }_{yl}{[n_1/f_s]}^2}\cos [{\omega }_0{n}_2/f_s+{\phi }_0] n_1=-\frac{N_y-M_y}{2},\cdots ,-1,0 n_2=-\frac{N_y}{2},\cdots ,-\frac{M_y}{2}\\ A_y{e}^{- {\beta }_{yh}{[n_1/f_s]}^2}\cos [{\omega }_0{n}_2/f_s+{\phi }_0] n_1=0,1,\cdots ,\frac{M_y}{2} n_2=-\frac{M_y}{2},\cdots ,-1,0\\ A_y{e}^{- {\beta }_{yh}{[n_1/f_s]}^2}\cos [{\omega }_0{n}_2/f_s+{\phi }_0] n_1=-\frac{M_y}{2},\cdots ,-1,0 n_2=0,1,\cdots ,\frac{M_y}{2}\\ A_y{e}^{- {\beta }_{yl}{[n_1/f_s]}^2}cos[{\omega }_0{n}_2/f_s+{\phi }_0] n_1=0,1,\cdots ,\frac{N_y-M_y}{2} n_2=\frac{M_y}{2},\dots ,\frac{N_y}{2}\end{array}$$

(5)

where n₁ is window variable of the exponential function in Y-axis signal model, n and n₂ are both window variables of the sine function in X-axis and Y-axis signal models, respectively. The normalized angular frequency is ω₀ = 2πf₀/f_s in unit of rad/sample. The ELF signal frequency f₀ is 23 Hz and the sampling frequency f_s is 1 kHz. A_x and A_y are amplitudes of the X-axis and Y-axis signals, respectively, and variable M_y is the coordinate intervals of the two peaks in the Y-axis. Variable β_x is the envelope decay rate of main lobe peak of X-axis, and β_yl and β_yh are used to represent the Y-axis envelopes deviating from the origin and facing the origin, respectively.

It is noted that the envelope decay rates (β_x, β_yl and β_yh) increase as the speed of the transmitter increases, as shown in Figure 3c,d. In order to simplify these three parameters to match the speed of the transmitter well, the fused envelope decay rate is derived here by a LS criterion. Taking the X-axis signal as an example, the LS criterion minimizes error function by estimating the parameters (A_x, φ₀ and β_x) in s_x[n]:

$$J_x={\displaystyle \sum _{n=-\frac{N_x}{2}}^{\frac{N_x}{2}}{(s_x[n]-x[n])}^2}.$$

(6)

Then, Equation (6) can be decomposed into the linear form as:

$$s_x[n]=H_x{\alpha }_x$$

(7)

where the observation matrix is:

$$H_x=[H_{x1},H_{x2}]=\left[\begin{array}{ll}{e}^{-{\beta }_x{[-\frac{N_x}{2f_s}]}^2}\cos [{\omega }_0(-\frac{N_x}{2})]& e^{-{\beta }_x{[-\frac{N_x}{2f_s}]}^2}\sin [{\omega }_0(-\frac{N_x}{2})] \\ ⋮& ⋮\\ e^{-{\beta }_x{[\frac{1}{f_s}(-1)]}^2}\cos [{\omega }_0\times (-1)]& e^{-{\beta }_x{[\frac{1}{f_s}(-1)]}^2}\sin [{\omega }_0\times (-1)]\\ e^{-{\beta }_x{[\frac{1}{f_s}\times 0]}^2}\cos [{\omega }_0\times 0]& e^{-{\beta }_x{[\frac{1}{f_s}\times 0]}^2}\sin [{\omega }_0\times 0]\\ e^{-{\beta }_x{[\frac{1}{f_s}(1)]}^2}\cos [{\omega }_0\times 1]& e^{-{\beta }_x{[\frac{1}{f_s}(1)]}^2}\sin [{\omega }_0\times 1]\\ ⋮& ⋮\\ e^{-{\beta }_x{[\frac{N_x}{2f_s}]}^2}\cos [{\omega }_0(\frac{N_x}{2})]& e^{-{\beta }_x{[\frac{N_x}{2f_s}]}^2}\sin [{\omega }_0(\frac{N_x}{2})]\end{array}\right].$$

(8)

Thus, the estimation of α_x (A_x and φ₀) is expressed as:

$${\widehat{\alpha }}_x=(H_x{}^T{H}_x{)}^{-1}{H}_x{}^{T}x.$$

(9)

Calculating the minimization of (6) is equal to calculate the maximum estimation of the X-axis signal energy, which is derived as:

$${\widehat{P}}_x=x^T{H}_x{(H_x{}^T{H}_x)}^{-1}{H}_x{}^{T}x.$$

(10)

Thus, the estimation of β_x can be estimated by maximizing (10) as:

$${\widehat{\beta }}_x=\arg \underset{{\beta }_x}{\max }[{\widehat{P}}_x=x^T{H}_x{(H_x{}^T{H}_x)}^{-1}{H}_x{}^T)x].$$

(11)

Similarly, the estimation of β_yl and β_yh are expressed as:

$$({\widehat{\beta }}_{yl},{\widehat{\beta }}_{yh})=\arg \underset{({\beta }_{yl},{\beta }_{yh})}{\max }[{\widehat{P}}_y=y^T{H}_y{(H_y{}^T{H}_y)}^{-1}{H}_y{}^T)y]$$

(12)

where H_y is transformed by s_y[n₁,n₂] and the transformation process is similar to H_x.

The purpose of the fusion is to find a relationship between these three envelope decay rates to minimize the error functions of the X-axis and Y-axis simultaneously. Assume that β_x = k_xyβ_yh = k_yβ_yl, and define the estimation of the orthogonal signal power as:

$${\widehat{P}}_{xy}=\frac{{\widehat{P}}_x}{N_x}+\frac{{\widehat{P}}_y}{N_y}.$$

(13)

Then, the relationship between k_y and normalized ${\widehat{P}}_y$ for three different speeds is shown in Figure 5a, and the relationship between k_xy and normalized ${\widehat{P}}_{\mathit{xy}}$ for three different speeds is shown in Figure 5b.

Figure 5a shows that ${\widehat{P}}_y$ reaches a maximum when k_y is around 4 in the cases of 5 m/s, 10 m/s and 15 m/s. Figure 5b shows that ${\widehat{P}}_{\mathit{xy}}$ reaches a maximum when k_xy is around 1 in the cases of 5 m/s, 10 m/s and 15 m/s. Thus, the fused envelope decay rate β_xy is defined through the linear relationship as:

$${\beta }_{xy}={\beta }_x={\beta }_{yh}=4{\beta }_{yl}.$$

(14)

The linear relationship has two main meanings: (1) the three envelope decay rates of the derived data fusion model have the similar varying trend as the speed of the transmitter changes, and the varying trend can be represented by a fixed linear relationship; and (2) the linear relationship satisfies the LS criterion by maximizing ${\widehat{P}}_x$, ${\widehat{P}}_y$ and ${\widehat{P}}_{\mathit{xy}}$, which contributes to the optimal estimation of the fused envelope decay rate.

3. Fast decision Tree Method

The data fusion model sets up the relationship between the orthogonal ELF signals and the fused envelope decay rate. According to the Equations from (10) to (14), maximizing the orthogonal signal power is equivalent to estimating the fused envelope decay rate, which also improves the ability of extracting the ELF signals from strong narrow-band background noise. Thus, the optimal estimation of the fused envelope decay rate is critical, and the maximized orthogonal signal power can be used as the test statistic for the signal detection, which is also the detection principle of the proposed FDT method. First, the fast calculation process of the proposed FDT method is derived as follows. Take the X-axis ELF signal as an example, according to the Equation (8), we have:

$$\begin{array}{l}{H}_{x1}{}^T{H}_{x2}=H_{x2}{}^T{H}_{x1}={\displaystyle \sum _{n=-\frac{N_x}{2}}^{\frac{N_x}{2}}{e}^{-{\beta }_x{[n/f_s]}^2}\cos [{\omega }_{0}n/f_s]\times e^{-{\beta }_x{[n/f_s]}^2}\sin [{\omega }_{0}n/f_s]}\\ = {\displaystyle \sum _{n=-\frac{N_x}{2}}^{\frac{N_x}{2}}{e}^{-2{\beta }_x{[n/f_s]}^2}\sin [2{\omega }_{0}n/f_s]}\approx 0\end{array}$$

(15)

Thus, simplify the Equation (10) as:

$${\widehat{P}}_x=x^T{H}_x{(H_x{}^T{H}_x)}^{-1}{H}_x{}^{T}x\approx x^T{H}_x{\left[\begin{array}{l}{H}_{x1}{}^T{H}_{x1} 0\\ 0 H_{x2}{}^T{H}_{x2}\end{array}\right]}^{-1}{H}_x{}^{T}x$$

(16)

where an approximate equation can be derived as:

$$P_{H_x}=H_{x1}{}^T{H}_{x1}\approx H_{x2}{}^T{H}_{x2}$$

(17)

where $P_{H_x}$ is defined as the energy of X-axis observation vector. Then, we have:

$${\widehat{P}}_x\approx \frac{x{H}_x{}^T{H}_x{}^{T}x}{P_{H_X}}=\frac{{(H_x{}^{T}x)}^T{H}_x{}^{T}x}{P_{H_X}}.$$

(18)

Similarly, the process of the Y-axis ELF signal is as follows:

$$P_{H_y}=H_{y1}{}^T{H}_{y1}\approx H_{y2}{}^T{H}_{y2}.$$

(19)

$${\widehat{P}}_y\approx \frac{y{H}_y{}^T{H}_y{}^{T}y}{P_{H_y}}=\frac{{(H_y{}^{T}y)}^T{H}_y{}^{T}y}{P_{H_y}}.$$

(20)

Next, the topology of the FDT method is derived as follows. The FDT method evenly divides the fused envelope decay rate into three categories at each of the tree nodes, and then calculates the normalized signal energy of the three child nodes, which is defined as:

$$\eta =\left(\frac{{\widehat{P}}_x}{P_x}+\frac{{\widehat{P}}_y}{P_y}\right)/2.$$

(21)

where P_x = x^Tx and P_y = y^Ty are received signal energy of X-axis and Y-axis signals, respectively. From Equations of (10) and (12), ${\widehat{P}}_x$ and ${\widehat{P}}_y$ are always smaller than P_x and P_y, respectively, and that η ∈ (0,1). The FDT method estimates the fused envelope decay rate through the decision rule, which is defined as determining the maximum η of the three child nodes as the next iteration direction. Here, the evolution of the FDT method and the decision rule are clearly shown in Figure 6. Variables L and n in β_L,n and η_L,n represent the tree layer and tree node. β_max corresponds to the maximum speed of the transmitter. As shown in Figure 6, the FDT method sets the maximum η as the next iteration direction, which is equivalent to selecting and estimating the fused enveloped decay rate that satisfies the LS criterion through the Equations (11) and (12).

The iteration stopping criterion is defined as follows: the FDT method calculates the η of three child nodes and then distinguishes them into three values: η_L,max, η_L,med and η_L,min. The iteration is stopped when the η meets:

$$\Delta \eta ={\eta }_{L,\max }-{\eta }_{L,\mathrm{med}}\le 0.001.$$

(22)

The varying trend of η at each FDT layer through the iteration procedure is listed in

Table 4. Iteration procedure of the FDT method.

L		V = 5 m/s, βxy = 5.5			V = 10 m/s, βxy = 21.2			V = 15 m/s, βxy = 46.9
1	${\widehat{\beta }}_{\mathit{xy}}$	8.4	25.2	42	8.4	25.2	42	8.4	25.2	43
	η	0.9839	0.8184	0.6975	0.9104	0.9679	0.9234	0.7881	0.9032	0.9303
2	${\widehat{\beta }}_{\mathit{xy}}$	2.8	8.4	14	19.6	25.2	30.8	36.4	42	47.6
	η	0.9703	0.9839	0.9288	0.9704	0.9679	0.9568	0.9261	0.9303	0.9313
3	${\widehat{\beta }}_{\mathit{xy}}$	6.53	8.4	10.27	17.73	19.6	21.47	45.73	47.6	49.47
	η	0.9946	0.9839	0.9676	0.9681	0.9704	0.9709	0.9312	0.9313	0.9310
4	${\widehat{\beta }}_{\mathit{xy}}$	5.9	6.53	7.16	20.84	21.47	22.09	46.98	47.6	48.2
	η	0.9956	0.9946	0.9920	0.9709	0.9709	0.9707	0.9513	0.9513	0.9512
Δη		0.9956 − 0.9946 = 0.001			0.9709 − 0.9709 < 0.0001			0.9513 − 0.9513 < 0.0001
Δβxy1		5.9 − 5.5 = 0.40			21.47 − 21.2 = 0.27			47.6 − 46.9 = 0.70

¹ Δβ_xy is error between β_xy (obtained from 2-D FEM simulation) and ${\widehat{\beta }}_{\mathit{xy}}$ (estimated from FDT method).

. The speed range of the FDT topology is set from 1 m/s to 15 m/s, and β_max is set to 50.

Table 4 shows in detail that the FDT method estimates β_xy by comparing and determining the maximum (which is bold) of the η in each layer. In layer 4, the iteration of the FDT meets the stopping criterion with small errors for Δη and Δβ_xy.

According to the topology and fast calculation process, once the maximum speed of the transmitter is confirmed, β_max of the FDT and β_xy of each FDT nodes are both confirmed. Meanwhile, since the parameters including N_x, f_s, ω₀ are priori knowledge, $P_{{\mathrm{H}}_{\mathrm{X}}}$ and $P_{{\mathrm{H}}_{\mathrm{Y}}}$ are also priori knowledge by the equation as:

$$(P_{H_x},P_{H_Y})=f({\beta }_{xy},N_x,f_s,{\omega }_0).$$

(23)

Therefore, the computation cost of the proposed FDT method is derived as:

$$\{\begin{array}{l}{C}_m=2N_x+3, C_a=2N_x+1 \mathrm{to} \mathrm{get} {\widehat{P}}_x \mathrm{in} (18)\\ C_m=2N_y+3, C_a=2N_y+1 \mathrm{to} \mathrm{get} {\widehat{P}}_y \mathrm{in} (20)\\ C_m=N_x+N_y, C_a=N_x+N_y+2 \mathrm{to} \mathrm{get} P_x \mathrm{and} P_y \mathrm{in} (21)\\ C_m=3N_x+3N_y+3, C_a=3N_x+3N_y+3 \mathrm{to} \mathrm{get} \eta \mathrm{in} (21)\end{array}$$

(24)

where C_m and C_a are computation costs of multiplication and addition, respectively. It is worth noting that the proposed FDT method only requires the computation costs of C_m = 3N_x + 3N_y + 3 and C_a = 3N_x + 3N_y + 3 to calculate one η at a time instead of massive matrix inversion and multiplication mentioned in (10) and (12), which is the key for ensuring a real-time performance.

Finally, the FDT method detects the ELF signals through a hypothesis test defined as:

$$\{\begin{array}{l}{H}_0:{\widehat{P}}_{xy}=\frac{{\widehat{P}}_x}{N_x}+\frac{{\widehat{P}}_y}{N_y}<\widehat{\gamma }\\ H_1:{\widehat{P}}_{xy}=\frac{{\widehat{P}}_x}{N_x}+\frac{{\widehat{P}}_y}{N_y}>\widehat{\gamma }\end{array}$$

(25)

where $\widehat{\gamma }$ is a threshold; H₀ is referred as the null hypothesis which shows nonexistence of signals and H₁ as the alternative hypothesis which shows existence of signals. The Newman-Pearson criterion is then used to build the threshold which maximizes POD for a given PFA, and the relationship between the threshold $\widehat{\gamma }$ and PFA is defined as:

$$\mathrm{PFA}=Q(\frac{\widehat{\gamma }-\widehat{\mu }}{\widehat{\sigma }})={\displaystyle {\int }_{\frac{\widehat{\gamma }-\widehat{\mu }}{\widehat{\sigma }}}^{\infty }\frac{1}{\sqrt{2\pi }}\exp (-\frac{1}{2}{t}^2)}dt$$

(26)

where $\widehat{\mu }$ and $\widehat{\sigma }$ are the mean and the variance of the test statistic due to narrow-band background noise, respectively, Q(x) is called the right-tail probability function. The threshold is calculated by (26) for a given PFA before the signal detection.

4. Performance Evaluation

4.1. Detection Methods and Theoretical Bounds

A maximum likelihood (ML)-based detection method sets the estimated amplitude of the ELF signals as the test statistic, which is based on the time-domain signals without considering the short transient characteristics of the envelopes. A discrete wavelet transform (DWT)-based detection method reconstructs ELF signal envelopes using the wavelet coefficient at the characteristic scales, using a Symlet wavelet based on the time-frequency analysis. Meanwhile, two theoretical bounds are used for comparison with the proposed FDT method’s results. The upper bound is derived from the matched filter (MF) method:

$$D_{MF}=\frac{s_{x}x}{{\sigma }_x{}^2}+\frac{s_{y}y}{{\sigma }_y{}^2}.$$

(27)

where s_x and s_y are the noiseless ELF signals and regarded as a prior known here; σ_x and σ_y are noise variances for the X-axis and Y-axis signals, respectively. The lower bound is derived from the energy detection (ED) method:

$$D_{ED}=\frac{x^{T}x}{{\sigma }_x{}^2}+\frac{y^{T}y}{{\sigma }_y{}^2}.$$

(28)

The estimation of SNR is defined as:

$$\mathrm{SNR}=10{\log }_{10}\left(\frac{s_x{}^T{s}_x+s_y{}^T{s}_y}{w_x{}^T{w}_x+w_y{}^T{w}_y}\right)\approx 10{\log }_{10}\left(\frac{{\widehat{P}}_x+{\widehat{P}}_y}{N_x{\widehat{\sigma }}_x+N_y{\widehat{\sigma }}_y}\right).$$

(29)

where w_x and w_y are the narrow-band background noise, which have the same definition as Equations (1) and (2). The mean of the narrow-band background noise can be considered to be zero, while the variances including ${\widehat{\sigma }}_x$ and ${\widehat{\sigma }}_y$ are obtained from the long-time observation. It is worth noting that ${\widehat{P}}_x$ and ${\widehat{P}}_y$ calculated from the FDT method are close to the noiseless ELF signal energy, and the estimation performance of the SNR will be discussed in Section 5.

4.2. Monte Carlo Simulation Study

The simulation parameters are chosen as follows: N_x = 1200, N_y = 2000, β_max = 50. The FDT method is compared with other detection methods using the receiver operating characteristics (ROC) in the cases of 5 m/s, 10 m/s and 15 m/s, as shown in Figure 7a–c, respectively. Meanwhile, when the PFA is set to 0.01, the result of POD versus SNR is shown in Figure 7d–f. All results are obtained using Morte Carlo simulations with number of trial equals to 10⁵ times. According to the given SNR, the simulated narrow-band background noise is generated numerically for each trial. The range of SNR of interest here is from −3 dB to 0 dB. An SNR equal to 0 dB means that the signal and the noise have the same energy, and the SNR equal to −3 dB means that the signal energy is half of the noise energy. Thus, the performance evaluation is based on comparing the detection methods’ ability to accurately extract and identify the target ELF signals that are buried in the strong narrow-band background noise, which is critical for the practical AGM system, because it also requires the same ability to detect the short transient very weak tracking signals with very poor SNR.

Figure 7a–f shows that the POD of the FDT method is higher than the PODs of the ML and DWT methods, especially when the speed is up to 15 m/s. Compared with the ML and DWT methods, the POD of the FDT method exhibits a very small decrease as the speed of the transmitter increases at a given SNR or PFA. Figure 7a–f shows that the detection performance of the FDT method is close to the MF method (upper bound), which is very attractive. Taking Figure 7c as an example, when PFA is close to 0, the PODs of the MF and FDT methods are both high, corresponding to the upper left corner of Figure 7c. Taking Figure 7f as another example, when the SNR is 0 dB, they are both more than 95% in POD and the difference between them is less than 5%. An in-depth discussion on the main reasons for the obtained results will be presented in Section 5.

4.3. Field Testing and Validation

The field testing is conducted by the newly designed and developed tracking system and 20˝ PIG by the authors’ group at Tsinghua University, as shown in Figure 8. The tracking system developed mainly includes an ELF transmitter, an ELF receiver and a pair of orthogonal search coil sensors. The transmitter generates the transmitting current which is a 23 Hz triangle wave using the H-bridge inverter controlled by MSP430 low power processor [43]. The titanium alloy shell is selected to protect the transmitter circuit for the high-pressure in-pipe situations due to the non-magnetic and lightweight properties. Figure 8c shows that the transmitter is equipped in the front of the PIG, and the PIG is pulled into the Q235 steel pipe at different speeds. Figure 8d,e shows that the ELF receiver mainly contains TMS320F28335 digital signal processor (DSP) and high-Q high-gain narrow-band filter and amplifier. The Q-factor of the designed 16th order band-pass filter is 16 (f₀ = 23 Hz), and the overall gain of the multi-stage amplifier is 100,000 (100 dB). The X-axis and Y-axis ELF signals received from the orthogonal search coil sensors are first filtered and amplified, and then converted to digital signals by 12-bit A/D converter. The DSP chips are programed by the proposed FDT methods in C programming language. The digital ELF signals of the X-axis and Y-axis, and the test statistics, are recorded in real time. Figure 8f shows the developed search coil sensors, which also uses the titanium alloy shell for the protection.

The narrow-long search coil sensor is developed to receive very weak ELF magnetic signals. The sensor is inserted by a rod of ferromagnetic material with a high magnetic permeability. The diameter of enameled wires is 0.25 mm. The structure is shown in Figure 9 and the parameters are listed in

Table 5. Parameters of search coil sensor.

D (mm)	Dc (mm)	n2	n1	lw (mm)	lc (mm)
30	20	2400	23	720	800

.

There were two groups of experiments: a high SNR situation with thin wall thickness, and a low SNR situation with thick wall thickness. Each group consisted of low speed (2 ~ 5 m/s), medium speed (5 ~ 10 m/s) and high speed (10 ~ 15 m/s), and the burial depth was 5 m. The experimental results are shown in Figure 10. The experimental parameters in Figure 10a–f, including pipe wall thickness and true speed, and the corresponding estimated parameters calculated by the receiver based on the proposed FDT method, are shown in

Table 6. Estimated parameters of experimental results.

Figure 10	Pipe Wall Thickness	True Speed	Estimated Speed	Estimated SNR
(a)	10 mm	4.7 m/s	4.9 m/s	10.35 dB
(b)	10 mm	8.0 m/s	7.6 m/s	9.45 dB
(c)	10 mm	11.1 m/s	10.1 m/s	6.95 dB
(d)	15 mm	6.1 m/s	6 m/s	4.00 dB
(e)	15 mm	8.3 m/s	7.8 m/s	0.55 dB
(f)	15 mm	12.5 m/s	13.8 m/s	−0.9 dB

. The true speed is recorded by the PIG odometer while the estimated speed is obtained from ${\widehat{\beta }}_{\mathit{xy}}$ of Table 4.

It is worth noting that when the normalized test statistic is equal to 1 and over the threshold (indicated by arrows), it represents that the transmitter is passing right below the receiver, and the receiver successfully detects the ELF signals in real time. Figure 10a–c shows that the X-axis ELF signal appears as an obvious unimodal characteristic while Y-axis ELF signal appears an obvious bimodal characteristic, which has the same characteristics as the 2-D FEM simulation. As the speed of the transmitter increases, the ELF signal envelopes appeared to be short transient and the durations of the test statistics exceeding the threshold were also shortened. Figure 10d–f shows that the ELF signals were very weak and heavily interfered by narrow-band background noise. For the case of Figure 10f shown in Table 5, the estimated SNR is even less than 0 dB. However, the normalized test statistic appears to be an obvious peak over the threshold around t = 17 s, which determines the hypothesis of the signal existence. At the same time and same measuring point, the ELF magnetic signals were also measured by Mag-03 flux-gate high-sensitivity magnetometer and found to be smaller than 10 pT. This result confirms that the developed tracking system can effectively detect the very weak ELF magnetic signals. The experiment results validate the proposed FDT method and the tracking system developed for high- sensitivity real-time PIG tracking.

5. Discussions

The proposed FDT method and tracking system developed is summarized as the flow chart shown in Figure 11. The tracking system mainly developed relies on high-sensitivity search coil sensors, a high-Q narrow-band filter and high-gain amplifier, and the proposed FDT detection method to realize real-time detection of short transient very weak ELF magnetic signals. The 2-D FEM simulation results indicate that the magnitude of the very weak ELF magnetic signals is as low as 10 pT, which is lower than one millionth of the magnitude of the Earth’s magnetic field. Thus, the orthogonal search coil sensors, which have high number of turns, high-permeability iron core, and narrow long structure, are designed and developed to measure the very weak ELF magnetic signals along the X-axis and Y-axis. Then, the narrow-band filter and amplifier effectively eliminates out-of-band noise and amplify 23 Hz target signals. Finally, as shown in Figure 10d–f, even if the ELF signals are immersed in strong narrow-band background noise, the test statistics calculated from the FDT method are over the thresholds under a low PFA when the transmitter passes right below the receiver.

The normalized power spectrum of the three different signals, including noiseless ELF signals, narrow-band noise and ELF signals with noise are shown in Figure 12. Since the frequency of narrow-band noise is concentrated at 23 Hz, it is difficult to discriminate between ELF signals and noise by frequency domain, which indicates that the frequency-domain information may not be useful for the detection. Thus, the performance of the time-frequency analysis based DWT method is lower than the time-domain characteristic based FDT method.

The probability density function (PDF) of ${\widehat{\beta }}_{\mathit{xy}}$ and estimated SNR is shown in Figure 13. Figure 13a shows a very concentrated distribution of ${\widehat{\beta }}_{\mathit{xy}}$ when the SNR is 6 dB or 3 dB, and even if the SNR is as low as 0 dB, the PDF also shows the concentrated distribution around 40~50 with a low variance (β_xy is 46.9 from Table 4). The ${\widehat{\beta }}_{\mathit{xy}}$ is also used to estimate the real-time speed of the transmitter according to the relationship established between the 2-D FEM simulation results of ELF signals and the data fusion model. Then, Table 6 shows the errors between the estimated speed and the true speed are less than 1.5 m/s. Thus, these results are the main reason that the detection performance of the FDT method is higher than the amplitude-estimation based ML method and close to the upper theoretical bound. Meanwhile, it can be seen from Table 6 that the estimated SNR can effectively distinguish high SNR situations and low SNR situations, and can reflect the decrease of SNR as the speed increases because of the short transient characteristics. The reason for the good estimation performance of the SNR can also be verified by Figure 13b.

The real-time performance of the ELF tracking system developed is based on the data fusion model, fast calculation process, and the determined topology of the FDT. The data fusion model sets up the linear relationship between the X-axis signal model and the Y-axis signal model by fusing the envelope decay rates. Thus, the fused envelope decay rate can match well with the speed of the transmitter. Then, the fast calculation process was derived to calculate the η without massive matrix inversion and multiplication, and the calculation amount was found to be low enough for the DSP chip. The most important factor is that each tree nodes of the FDT is determined when the maximum speed of the transmitter is given. It is worth noting that the matrices including $H_X$ and $H_Y$, and the values including $P_{{\mathrm{H}}_{\mathrm{X}}}$ and $P_{{\mathrm{H}}_{\mathrm{Y}}}$ used for the calculation of the η and test statistic are both prior-known information because of the fast calculation process and determined FDT topology. Therefore, the proposed FDT method is very easily programmed, stored and calculated by the embedded DSP chip in practical applications. The experiments conducted by the receiver using a TMS320F28335 DSP chip finally validate the real-time performance of the FDT method.

6. Conclusions

Real-time detection of the short transient very weak ELF magnetic signals that are buried in strong narrow-band background noise is one of the key research areas in the nondestructive testing (NDT) field to ensure the safety of pipeline ILI for the oil and gas industry. The existing AGM systems work well for high-SNR and low-speed situations, but still face the challenges brought by high-sensitivity high-speed PIG tracking. This paper proposes a novel FDT detection method, as well as an ELF tracking framework for the high-speed PIG, to detect the short transient very weak ELF magnetic signals accurately and effectively. The 2-D FEM simulation results indicate that the duration and magnitude of the ELF signals can be as low as 1 s and 10 pT, which is validated by the experimental results. According to the short transient very weak characteristic as well as the linear relationship between X-axis and Y-axis ELF signal models, the data fusion model was established, and the fused envelope decay rate was used to match the speed of the transmitter. The proposed FDT method estimates the fused envelope decay rate to calculate the test statistic in real time, which was compared with the threshold to determine the existence of the ELF signals. The simulation results show that the FDT method has a higher POD than the other existing methods, especially when the speed is faster, for example, at 15 m/s. Meanwhile, the performance of the FDT method is close to the upper theoretical bound and the difference between them is less than 5% when SNR is 0 dB. The experimental results further validate that the proposed FDT method can detect the short transient very weak ELF signals in real time. Based on these results, the proposed FDT method and ELF tracking system developed are validated to be effective in real-time tracking of high-speed PIG.