Talbot Algorithm with a Trapezoidal Wave in the 2.5D Airborne Transient Electromagnetic Method in Marine Investigations

Abstract

The airborne transient electromagnetic (TEM) method is widely used in land applications but faces challenges in marine applications due to the strong masking effects of seawater. The accuracy of the inverse Laplace transform algorithm strongly affects the applicability of the 2.5D TEM method in marine research; thus, an appropriate transmitted waveform is required. To resolve these issues, a trapezoidal wave was utilized as the waveform of the current source, owing to the substantial energy contribution of the low-frequency range. Further, the Talbot algorithm was selected for the inverse Laplace transform as it can yield an accurate response with smaller summation terms than the commonly used Gaver–Stehfest (G-S) algorithm. On this basis, a rugged seabed and a subsea tunnel were also investigated. The voltage response is obtained when the flight heights of the loop source change. The results confirmed that the Talbot algorithm with a trapezoidal waveform is more reliable and robust for complex marine geological models and is expected to provide an effective approach for marine explorations.

1. Introduction

Transient electromagnetic (TEM) methods can detect a pure secondary field after the cutoff of the primary field, and they have been widely used in land applications [1,2,3,4,5]. In recent decades, with rapid ocean development and utilization, greater attention has been given to the application of TEM methods in the marine field [6,7,8,9,10,11]. To overcome the masking effects of seawater, most studies have adopted the towed configuration [10,11,12]. However, this configuration has low work efficiency and high acquisition cost. In addition, many shallow areas are inaccessible for exploration ships because of rocks or the depth of seawater. For configurations above the sea surface, the airborne TEM method is considered a promising candidate because it is convenient and highly efficient [1,2,13,14,15,16]. In practice, experiments conducted at sea need to consider not only the complicated environment but also other factors such as manpower and costs. Thus, numerical design studies and accurate date interpretation need to be reinforced before launching a sea-going measurement. Some existing interpretation methods for land applications can be transferred directly to marine TEM methods. However, the seawater layer causes these to be different from land TEM methods in terms of sounding ability, device design, and numerical modeling. Therefore, some issues pertaining to the airborne TEM method in marine investigations need to be further studied for it to become a feasible method of analysis.

An accurate numerical method is necessary in marine TEM methods. The method of 2.5D modeling, i.e., a 2D geoelectric model with a 3D loop source, is often used because it is closer to the practical geological conditions than either 1D or 2D models, and has a lower computation cost than 3D models. Moreover, many marine activities such as surveys of seafloor topography and subsea targets, can be viewed as 2.5D models [17,18,19,20]. Analysis by 2.5D forward modeling is based on the background/anomalous-filed algorithm. In order to eliminate the time and strike components, it is usually manipulated in the Laplace and Fourier domains and calculated using the finite element method (FEM) [21,22,23,24]. The inverse Laplace and Fourier transforms are used to recover the solution in the time domain. Each step in the modeling inevitably introduces some calculation errors. As the Fourier inversion can be evaluated accurately by the sine and cosine transforms [25,26], the Laplace inversion becomes a key factor that influences the accuracy of the marine TEM response. There are many algorithms for computing the inverse Laplace transform [27,28,29], but a favorable one should balance the computational load and calculation accuracy.

In land investigations, for the airborne TEM methods with a step waveform source, the Gaver–Stehfest (G-S) algorithm is widely used to compute the inverse Laplace transform because it is a pure real number operation and can be implemented simply and easily [22,30]. However, in most airborne TEM systems, commonly used current waveforms include trapezoidal [31,32,33], half-sine [34,35], and triangular waves [15]; these waveforms have a considerable influence on the TEM response. It has been demonstrated that the G-S algorithm for a complex wave yields large errors, especially in the late observation time. As seawater is a high-conductivity medium, the masking effect causes EM waves to yield a lower penetration depth and attenuate quickly in seawater. The response of the target in the sea is low, and sampling it for analysis takes a long time. In addition, the accuracy of the G-S algorithm is closely related to the problem and the summation terms, which means it has poor robustness. Therefore, an alternative algorithm matched with the current waveform is necessary for the inverse Laplace transform in marine TEM calculations.

In this study, 2.5D forward modeling in marine TEM methods was investigated using the FEM. The trapezoidal wave was selected as the source waveform. For the inverse Laplace transform, the Talbot algorithm is selected to match the waveform in terms of accuracy and robustness. Based on such an algorithm and waveform, two examples in shallow water—a rugged seabed and a subsea tunnel—were modeled, and the effect of the loop attitude on the airborne TEM method was considered and studied.

2. Methods

2.1. Airborne TEM System and Waveform of Current Source

The airborne TEM system with a central-loop configuration is shown in Figure 1. A current source with a certain waveform is loaded in the transmitter loop and generates the primary field. When the current pulse is cut off, the eddy current in the anomaly produces a secondary field, which can cause an induced voltage response in the receiver loop [36,37,38]. In contrast to land applications, both the primary and secondary fields penetrate the sea layer. As seawater has a higher conductivity than soil and rock, it exhibits strong absorption in the EM field. According to the masking effect, the lower the frequency is, the deeper the EM wave penetrates. Therefore, the current waveform should have more power in the low-frequency range.

Common waveforms in TEM systems include square, trapezoidal, half-sine, and triangular waves, all of which have been studied deeply and applied widely in land TEM systems [39]. As stated above, the airborne TEM system for marine investigations is heavily influenced by the masking effect; thus, a lower frequency with more power is necessary for deep detection. According to the frequency spectra of these waves, it can be observed that the trapezoidal and square waves have similar distributions in lower frequencies, and they have more power than the others. However, it is difficult to yield an absolute square wave in practical systems [31,33,38]. Therefore, the trapezoidal wave was selected as the source waveform in our investigation. The trapezoidal wave can be expressed in the Laplace domain as

$$f_s=\frac{1}{t_d{s}^2}-\frac{1}{t_d{s}^2}{e}^{-s{t}_d}-\frac{1}{t_d{s}^2}{e}^{-s(t_d+t_s)}+\frac{1}{t_d{s}^2}{e}^{-s(2t_d+t_s)},$$

(1)

where $t_s$ is the duration time of the stationary period; and $t_d$ is the falling edge, which is equal to the rising edge.

2.2. Inverse Laplace Transform Algorithm

2.2.1. Vertical Component of Magnetic Field in the ($k_y$, s) Domain

For 2.5D forward modeling in the FEM, the electromagnetic field is first analyzed in the Laplace and Fourier domain, i.e., the ($k_y$, s) domain, where $k_y$ is the wavenumber of the along-strike field; and s is the Laplace transform variable. It is then recovered in the time domain. For a layered geoelectric model, the vertical component of the magnetic field under the trapezoidal waveform excitation can be written as

$${\widehat{H}}_z=f_s·H,$$

(2)

$$H=\{\begin{array}{c}G{\displaystyle {\int }_0^{\infty }\left(\frac{u_0{R}_1-u_1}{u_0{R}_1+u_1}{e}^{-2u_0{z}_0+u_{0}z}+e^{u_{0}z}\right)}\frac{k_x^2+k_y^2}{u_0{k}_x}\sin (k_{x}a)\cos (k_{x}x)d{k}_x,\mathrm{for}z<0\\ G{\displaystyle {\int }_0^{\infty }\left(\frac{u_0{R}_1-u_1}{u_0{R}_1+u_1}{e}^{-2u_0{z}_0+u_{0}z}+e^{-u_{0}z}\right)}\frac{k_x^2+k_y^2}{u_0{k}_x}\sin (k_{x}a)\cos (k_{x}x)d{k}_x,\mathrm{for}0<z<z_0\end{array},$$

(3)

$$G=\sqrt{\frac{2}{\pi }}\frac{I_0\sin (k_{y}b)}{\pi k_y}.$$

(4)

For a homogeneous half-space, R₁ = 1, and for a layered model, R₁ is defined as follows:

$$R_1=\frac{u_1+u_2\tanh (u_1{h}_1)}{u_2+u_1\tanh (u_1{h}_1)},$$

(5)

$$\begin{array}{l}{u}_0=\sqrt{k_x^2+k_y^2+k_0^2},k_0^2=\mu s\left(\epsilon s+{\sigma }_0\right),\\ u_1=\sqrt{k_x^2+k_y^2+k_1^2},k_1^2=\mu s\left(\epsilon s+{\sigma }_1\right),\\ u_2=\sqrt{k_x^2+k_y^2+k_2^2},k_2^2=\mu s\left(\epsilon s+{\sigma }_2\right),\end{array}$$

(6)

where $z_0$ is the height of the transmitter loop; $k_x$ and $k_y$ are the wavenumbers along the x-and y-directions, respectively; $h_1$ is the depth of the first layer; ${\sigma }_0$, ${\sigma }_1$, and ${\sigma }_2$ are conductivities of the air, seawater, and seabed layer, respectively; $I_0$ is the amplitude of the incident current; the current waveform $f_s$ is considered in the Laplace domain; $a$ and $b$ are the values for half of the loop side length.

In Equation (2), the field is divided into two parts; one is the current waveform $f_s$, and the other is H, which only has a relation with parameters of the transmitter instrument and the geoelectric model. In addition, the Laplace form of the trapezoidal wave is divided into four parts, as shown in Equation (1). When solving the Laplace inversion of Equation (2), the combination of H and the first part of $f_s$ can be solved; then, the other parts can be obtained by the translation property of the inverse Laplace transform.

2.2.2. Laplace Inversion Algorithm

Various numerical approaches have been applied for the Laplace inversion [27,28,29], and they have different influences on the results. Considering the computational load and accuracy, two typical numerical approaches, the G-S algorithm and the Talbot algorithm, will be applied to solve the inverse Laplace transform, respectively. The G-S algorithm is a real operation, which is easy to calculate and has been widely used in 2.5D land TEM. The Talbot algorithm is a complex operation. From the results yielded by numerical examples, it has the advantage of superior accuracy and weak dependence on problems under the same calculation conditions [29,40,41,42].

• G-S algorithm

The G-S algorithm [28,30] converts a function $F$ in the Laplace domain into the time-domain solution according to the relation:

$$f(t)=\frac{\ln 2}{t}{\displaystyle \sum _{j=1}^N{V}_{j}F\left(\frac{\ln 2}{t}j\right)},$$

(7)

and the weight coefficient $V_i$ is calculated by the following equation:

$$V_j={\left(-1\right)}^{\frac{N}{2}+j}{\displaystyle \sum _{k=INT\left(\frac{j+1}{2}\right)}^{min\left(j,N/2\right)}\frac{k^{N/2}\left(2k\right)!}{k\left(k-1\right)!\left(N/2-k\right)!\left(j-k\right)!\left(2k-j\right)!}},$$

(8)

where N is the number of terms, the value is usually a positive even number (N = 8, 10, 12, 14, 16). INT is the integral function, and min is the minimum value function.

• Talbot algorithm

The fixed version for the algorithm has been demonstrated by Abate and Whitt [40,42], and the algorithm is shown below:

$$f(t)=\frac{0.4}{t}{\displaystyle \sum _{j=0}^{N-1}\mathrm{Re}\left[c_j^{T}F(\frac{{\delta }_j}{t})\right]},$$

(9)

where ${\delta }_j$ is calculated by ${\delta }_0=0.4N$; ${\delta }_j=0.4\pi j\left[\cot (j\pi /N)+i\right]$ for $0<j\le N$; and the number of terms N is a positive integer. The corresponding coefficient is given by

$$c_0^T=\frac{1}{2}{e}^{{\delta }_0},$$

(10)

$$c_j^T=\left\{1+i\left(\frac{j\pi }{N}\right)\left[1+{\cot }^2\left(\frac{j\pi }{N}\right)\right]-i\cot \left(\frac{j\pi }{N}\right)\right\}{e}^{{\delta }_j},$$

(11)

where i is the imaginary unit, $i=\sqrt{-1}$.

3. Performance of the Laplace Inversion Algorithm

To select an appropriate Laplace inversion method for the airborne TEM system in marine investigations, the responses from the G-S and Talbot algorithms were compared when they were used for a homogeneous sea half-space model and a two-layered marine geology model, respectively. In fact, many marine activities can be viewed as 2.5D problems, and these two models are highly representative because they are often used as the background medium in 2.5D modeling, and the accuracy of their results affects the effectiveness of forward and inverse modeling. Moreover, 2.5D forward modeling can be carried out through numerical methods such as the FEM [22]. In the following, the FEM was combined with the Talbot inverse Laplace transform to calculate the TEM response in the sea model.

3.1. Model 1: Homogeneous Sea Half-Space

The two algorithms were first applied to a homogeneous sea half-space. The seawater conductivity was set to 4 S/m. To simplify the analysis, it was assumed that a central-loop source was placed on the sea surface with a transmitting square loop of 2 × 2 m, and the equivalent area of the receiving loop was 1 m² [16]. The current source carried a trapezoidal waveform, where $t_s$ = 8 ms, and $t_d$ = 1 ms, as shown in Equation (1). The pulse width $T_0$ was assumed to be 10 ms. The current amplitude $I_0$ was assumed to be 100 A. These parameters remained constant in the following study.

In Equations (7) and (9), the number of terms N in the summation was selected from 8 to 14. The voltage responses derived from the G-S and Talbot algorithms were compared with the analytic solution [43]. As shown in Figure 2a,b, the response obtained by the G-S algorithm deviated from the analytic solution for all values of N, and the relative errors were high for most observation times, especially at late times. However, as shown in Figure 2c,d, the best results were derived using the Talbot algorithm, which agreed well with the analytic solution and had the lowest relative errors (below 1% for most of the observation period). Compared with the results obtained by the G-S algorithm, the response for the Talbot algorithm was less sensitive to the value of N. Even when N had a smaller value, the Talbot algorithm had high accuracy for computing the inverse Laplace transform.

3.2. Model 2: Two-Layered Marine Geology Model

The two algorithms were applied to a two-layered marine geology model, which included a seawater layer (4 S/m), and a seabed layer (0.5 S/m). The depth of the seawater was assumed to be 20 m. The source parameters remained the same as those in Model 1. As there is no analytic solution for the layered model, the voltage response was compared with a conventional solution obtained by the Hankel transform [25], which was operated in the frequency domain.

As shown in Figure 3a,b, the voltage response obtained by the G-S algorithm deviated from the conventional solution, especially in later times. Moreover, the relative errors were larger for most of the observation times. By contrast, the Talbot algorithm worked well in all the observation times, as shown in Figure 3c,d, and the relative errors were below 1% in the majority of the observation times for all values of N. Moreover, the Talbot algorithm rarely depended on N. The G-S algorithm is often preferred in TEM land modeling owing to the minimum number of summation terms and rapid computation; however, in marine modeling, the skin effect of the seawater necessitates a long duration of observations. Thus, the response in late observation times is important for the investigation.

The following can be determined from the results above: (1) the G-S algorithm is not suitable for the trapezoidal wave source and also has poor performance on TEM applications for marine modeling, although it is simple to implement; and (2) the Talbot algorithm does not require many summation terms, and can yield an accurate response for the trapezoidal waveform regardless of where the voltage response originates from in the sea homogeneous half-space and sea-layered modeling. Therefore, the Talbot algorithm is more appropriate for computing the inverse Laplace transform in marine TEM investigations.

4. Simulation Models in Marine Investigations

After the Laplace inversion method matched with the trapezoidal wave were determined, the outstanding question was how well they can be applied in marine TEM. Thus, the voltage responses of a rugged seabed and a subsea tunnel were investigated, which are two common models in the shallow sea. Furthermore, considering that the airborne TEM method has been widely applied, the influence of the Talbot algorithm on the voltage response of the loop distance above sea level with the trapezoidal wave were also discussed.

4.1. Model 1: Subsea Tunnel

As shown in Figure 4, an immersed tunnel (similar to that of the Hong Kong-Zhuhai-Macau Bridge) was modeled. To simplify the analysis, the tunnel was regarded as a rectangle having a size of 10 m × 40 m (height and width) and a buried depth of 10 m beneath the seabed. The sea depth $h_1$ was 30 m, and the conductivities ${\sigma }_0$, ${\sigma }_1$, ${\sigma }_2$ were set to 0, 4, and 0.5 S/m, respectively. The tunnel equivalent conductivity ${\sigma }_a$ was assumed to be 10³ S/m. The loop source was moved from −100 m to 100 m on the x-axis, and the center of the rectangle was on the y-axis, which was just above the origin of the axis. The loop height $h_0$ was set to 10 m.

The number of terms N for the Talbot algorithm was first set as 8. In Figure 5, the response curves are flat with high amplitudes at early times. However, an abnormal hump becomes evident near 25 ms. As the conductivity of the tunnel is greater than that of the seabed, the field encountering the tunnel yielded stronger responses than that springing from the seabed. In addition, because the secondary field decayed slowly in seawater, the hump became gradually more apparent at later stages. Moreover, the peak of the hump appeared at the sounding point of 0 m. It can be thought that the TEM response reflects the subsea tunnel directly. Furthermore, it was shown that the Talbot algorithm matched with a trapezoidal wave source in the airborne TEM can reflect the low-resistance body distinctly.

In Figure 6a, the number of terms N is changed from 8 to 14, and the voltage response is selected at the observation time of 100 ms because the anomaly in the curve is more apparent at this time. As the value of N increased, the voltage response curve remained almost constant. When N had a smaller value, the Talbot algorithm with the trapezoidal wave in the FEM method yielded a high accuracy for the subsea tunnel model. However, with increasing N, the computation time also increased. Thus, the number of terms N was set as 8 for this model. The impact of the loop source height on this model was also analyzed; the voltage response is presented in Figure 6b. The response was also chosen at the observation time of 100 ms. As the source height $h_0$ was increased, the voltage response curve gradually became flatter. At greater height, the hump gradually vanished, and flat curves were recovered. Along with the increase of $h_0$, the field diffusing into the seawater was reduced. Due to the strong attenuation of the field in seawater, the detecting depth decreased, and the tunnel had little impact on the response. This demonstrates that the Talbot algorithm with a trapezoidal wave source has good performance for the sounding of a conductive target in the seabed within the range allowed by the detecting depth.

4.2. Model 2: Rugged Seabed

For further illustration, the seafloor topography in form of a rugged seabed is modeled in Figure 7 that had a cosine form with the period of 100 m and an amplitude of 10 m. The trench length $l_1$ was 100 m, and the trench depth $h_2$ was 20 m. The sea depth $h_1$ was 15 m. The center of the loop source was at the origin of the axis, which was also just above the peak of the hump in the rugged seabed.

Figure 8 shows the voltage response when the loop height is 10 m and the number of terms N for the Talbot algorithm is 8. Due to the rugged terrain of the seabed, different parts of the seabed produced different responses. At the location of the two trenches, two distinct humps in the curves were observed. However, the ridge part of the seabed yielded an apparent concavity in the curves. Since the conductivity of the ridge was smaller than that of the seawater in the trenches, this lead to a downward deviation between the response curves of the seawater and those of the ridge. Thus, the response generated by the ridge was smaller than that of the trenches. As time increased, the anomaly in the curves became gradually apparent. Other obvious anomalies were noted later, especially at 100 ms.

Further, the value for the number of terms N was varied from 8 to 14, and the voltage response of the observation time of 100 ms is shown in Figure 9a. When the value of N was increased, the voltage response curve remained unchanged. However, its computation time considerably increased, which was consistent with the results obtained by the subsea tunnel model. Thus, the number of terms N was set as 8, which can meet the computing needing for this model. Although the rugged seabed was reflected directly when the loop source was placed at the height of 10 m, in practical TEM investigations, the loop source has a greater height; thus, the source height $h_0$ was also analyzed for this model. In Figure 9b, when $h_0$ was less than 20 m, there were two humps and a concavity corresponding to the location of the rugged terrain in the seabed. However, if $h_0$ exceeded 20 m, the response curves became flat. In this modeling, it was noted that the rugged seabed was the anomaly. Considering the effect of the seabed terrain on the response, the maximum height of the source was different than that of the low-resistance body modeling. This demonstrates that in a TEM system a suitable detecting height of the source must be selected according to the subsea terrain and target. It was further validated that the Talbot algorithm provides favorable results for methods that require the inverse Laplace transform even when the Talbot algorithm has a smaller value of N. In addition, the marine TEM method with a trapezoidal wave source works well for complex marine modeling.

5. Conclusions

This paper focused on the algorithm for the inverse Laplace transform in 2.5D airborne TEM for marine investigation. Considering the masking effect of seawater, a trapezoidal wave was selected as the source waveform because its frequency spectrum mainly aggregates in the low-frequency range. For the inverse Laplace transform, the Talbot algorithm was found to have the capacity to match this waveform better than the common G-S algorithm owing to a higher accuracy and robustness. Based on the Talbot algorithm with the trapezoidal waveform, two typical applications—a rugged seafloor and a subsea tunnel—were modeled by using the 2.5D FEM. The voltage response for the number of terms and the height of the loop source were investigated. It was demonstrated that the Talbot algorithm with a trapezoidal waveform provides good performance for the complex marine geology model, where the suitable detecting height of the loop source is selected based on the subsea terrain and target. The airborne TEM method with the Talbot algorithm can effectively and rapidly detect an anomaly in shallow water at low cost, and therefore is expected to be a promising technology in marine investigations.

In future work, the actual geologic model, such as seafloor topography or subsea tunnels, will be investigated by the existing airborne TEM system. By comparing the measured data with the simulation results, the parameters of the simulation model will be optimized in future research. Furthermore, the forward algorithm will also be optimized. The results and data could provide a basis for the development of a compact airborne TEM system for marine investigations.