# Three-Dimensional Localization of Buried Polyethylene Pipes Using Acoustic Method

William Xerri, Gineth Saracco, Alessandra Ribodetti, Laurent Zomero, Philippe Picon

► To cite this version:

William Xerri, Gineth Saracco, Alessandra Ribodetti, Laurent Zomero, Philippe Picon. Three-Dimensional Localization of Buried Polyethylene Pipes Using Acoustic Method. Sensors, 2022, Sensors (Advances in Non-Destructive Testing Methods), 22 (23), pp.9433.1-9433.23. 10.3390/s22239433. hal-03883870

HAL Id: hal-03883870

https://hal.science/hal-03883870

Submitted on 4 Dec 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
1 Article

2 3D Localization of Buried Polyethylene Pipes using Acoustic Method

4 William Xerri ¹,³,*, Gineth Saracco ¹, Alessandra Ribodetti ², Laurent Zomero ³, Philippe Picon ³

5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
1
CNRS-UMR7330 CEREGE, AMU, CdF, IRD, Europole de l'Arbois, BP 80, 13545 Aix-en-Provence Cedex 4,
France; xerri@cerege.fr; ginet@cerege.fr

2 IRD-UR082 CNRS-UMR7329 GEOAZUR, Campus Azur Université et Observatoire de la côte d'azur, 250 rue Albert Einstein CS 10269, 06905 Sophia Antipolis, France; Alessandra.ribodetti@geoazur.unice.fr

³ MADE-SA, 167 Impasse de la Garrigue, 83210 La Farlède, France; l.zomero@made-sa.com; p.picon@made-sa.com

* Correspondence: xerri@cerege.fr

**Abstract:** Localization of buried polyethylene pipes is an important issue for network managers. This study focuses on an acoustic method, which consists of vibrating the pipe and observing the signal with a receiver placed on the ground surface. This method provides an estimate of the path of the pipe but gives no information on the depth. We developed a multi-sensor method based on the principle of vibrating the pipe, which allows estimating the depth while being non-invasive and non-destructive and without, a priori, information on the propagation medium. These sensors are positioned perpendicular to the pipe. We developed a new estimator to estimate the depth and the propagation velocity in the medium which is an important variable in our problem. This estimator is based on the MUSIC algorithm and is adapted to our choice of modeling. In this paper, two models of travel times in typical situations are presented. The first one represents the case where all sensors can be placed inside the trench (on the ground surface) in which the pipe is buried. The second one represents the case where sensors are placed inside and outside the trench. These travel time models aim to provide a fast result to allow the method to be used by field agents. They are compared with a full wavefield modeling by finite differences.

**Keywords:** acoustic method; buried polyethylene pipe; MUSIC algorithm; propagation time modelling; signal processing; full wavefield

29 Academic Editor: Firstname
30 Lastname

Received: date 31

Accepted: date 32

Published: date 33

**Publisher's Note:** MDPI stay¹⁴ neutral with regard to jurisdiction¹⁵ claims in published maps¹⁶ institutional affiliations. 37

<figure><img src="image_2.png" /></figure>

38
39
Copyright: © 2022 by the author⁴⁰
Submitted for possible open access⁴¹
publication under the terms of the Creative Commons⁴²
conditions of the Creative Commons⁴³
Attribution (CC BY) license⁴⁴
(https://creativecommons.org/licenses/by/4.0/).

1. Introduction

Different methods exist for locating buried pipes [1] and, more generally, for investigating near-surface structures [2-7]. The choice of method to locate buried pipes depends on the context. The use of tracer wires is becoming increasingly widespread. This method consists of burying an electric wire with the pipe, which allows, by passing a current through the wire, the creation of an electromagnetic field and, therefore, the use of electromagnetic methods. Electromagnetic methods provide very good results for locating pipes when they are applicable. If there is no tracer wire, then two possible cases can be discerned. Either the pipe is metallic, or it is not. In the case of a metal pipe, an electric current can be injected into the pipe; electromagnetic methods, which work very well, can, therefore, be used. In the case of non-metallic pipes, two types of processes can be distinguished in a non-exhaustive way: i) Ground-Penetrating Radar (GPR) [8, 9] is a versatile tool for locating buried infrastructures; however, this type of tool often needs to be calibrated according to the type of soil to which it is applied. Moreover, it is difficult to differentiate a water pipe from a gas pipe or a buried electrical cable. ii) Acoustic methods are applicable in the case of non-conductive pipes [1, 10-12]. According to [1], we can refer to several types of acoustic methods: seismic wave methods [13], point
vibration measurements [14] and pipe excitation methods [15]. The first two categories
can use a specific source placed on the ground surface for probing. These methods do not
differentiate between the pipe of interest and other pipes present. What distinguishes the
last category is its ability vibrate the pipe, and this vibration will then be diffused into the
ground.

This study focuses on the pipe excitation methods that have the advantage of
discriminating the pipe of interest in a dense urban environment with noise.

The principle of this is to inject an acoustic signal (called acoustic signature) into the
pipe. A method that has been used for many years in the field uses a receiver (geophone)
placed on the ground surface, which measures the vertical vibrations of the ground.
Through successive measurements of the energy of the received signal, this method can
estimate the passage of the pipe and follow its path. However, this method does not
provide any information on the pipe depth. The GasTracker tool, based on this method,
was developed by the company MADE-SA. In this study, a GasTracker is used to
estimate, experimentally on the ground surface, the X and Y position (at ± 10 cm) of the
pipe in order to position the sensor network in the area where the pipe is located.

The objective of this study is to implement a multi-sensor method, based on the
differences in travel times between sensors, to estimate the depth in addition to the
passage of the pipe [10, 11].

The aim of this work is to estimate the depth with an accuracy below 10 cm with a
non-destructive and non-invasive method and without any a priori information on the
propagation medium (e.g., characteristics of the soil, which can change completely from
one application area to another).

Nevertheless, some information on how the pipe is buried are known. The pipe is
buried between 0.4 m and 1.5 m. It is buried in a trench whose width can vary from
approximately 30 cm to over 1 m. Therefore, two different media can be distinguished,
the inside of the trench and the outside of the trench (vertical stratifications). Moreover,
in this problem the depth is at the meter scale.

An important aim of this study is to define a model of the travel time of the signal
that provides the fastest result, in order to directly estimate the depth in the field. In the
travel time modelings presented in this paper, we consider the first arrival time of the
signals. Only P-waves and converted S-waves are considered. Geophones used are single
component sensors. They measure the vertical component of the received vibrations.

The travel time modelings that we propose are confronted using a simulation with
much more complete propagation modeling, that of the complete wavefield with finite
differences. This comparison aims to verify the coherence of the estimated quantities, in
particular, the velocity.

A cylindrical scattering of the vibration induced by the pipe is considered [16-19]. To
ensure that all sensors observe the same section of the vibrating pipe, they are placed
perpendicular to the pipe. This is why the problem is represented in a plane orthogonal
to the passage of the pipe. This section can be approximated by a point source.

The variables to be estimated in this problem are the position of this source, its
lateral position Sx and especially its depth Sz, but also the average velocity of
propagation in the ground. An average velocity is assigned to each medium of
propagation.

In a first part, the sensors are assumed to be all positioned inside the trench.
Therefore, a modelling of the travel times of the signal between the pipe and the receivers
by considering a weakly heterogeneous medium and by assigning a mean propagation
velocity is presented. In order to ensure an acceptable estimation accuracy, an estimator
based on the MUSIC (Multiple Signal Classification) algorithm adapted to this modelling
of travel times is developed. This estimator allows us to evaluate variables of interest
(position of the source and the velocity). Measurements are carried out on a test area in
order to work on real data. In the second part, the sensors are positioned inside and
outside the trench. As an extension of the method developed in the first part, the
modelling is adapted to this case to fit a larger number of situations. The different compaction of the soil on either side of the trench is taken into account. This requires considering a velocity specific to each of these environments. This second modelling of travel times is validated within the framework of the problem, initially through numerical simulations, then through comparison with finite differences, and finally by comparison with real data.

## 2. Modelling with a single propagation medium (M1)

### 2.1. Propagation time modelling M1

In this section, all sensors are considered located inside the trench. Since no information on the characteristics of the propagation medium is known, the first approach is to consider the propagation medium as weakly heterogeneous and to assign an average propagation velocity. The problem is represented in a plane orthogonal to the passage of the pipe (Figure 1). The pipe section has been represented using a point source S, the receivers $R_i$ were placed on the ground surface perpendicular to the pipe. An average velocity $V_0$ has been assigned to the propagation medium. The position of the source S ($S_x$; $S_z$) and the velocity of the wave $V_0$ were unknown. The only known parameters were the coordinates of the receivers $R_i$ ($R_{ix}$; $R_{iz} = 0$).

<figure><img src="image_3.png" /><figcaption>Figure 1. Scheme of the Model M1: case of a single propagation medium.</figcaption></figure>

We defined $\theta_{M1}$ as the vector of variables to be estimated in the case of Model M1.

$$ \theta_{M1} = [S_x \ S_z \ V_0]^T \qquad (1) $$

Here, the symbol $[.]^T$ is the transposed operator, $S_x$ the plumb of the pipe, $S_z$ the pipe depth and $V_0$ is the average propagation velocity in the medium.

We note $\tau_i(\theta_{M1})$ the wave travel time between the source S and the receiver $R_i$ for Model M1.

$$ \tau_i(\theta_{M1}) = \frac{|SR_i|}{V_0} \qquad (2) $$

We note $\tau_{1i}(\theta_{M1})$ the relative delay time between sensors $R_i$ and $R_i$ for Model M1.

$$ \tau_{1i}(\theta_{M1}) = \frac{|SR_i| - |SR_1|}{V_0} \qquad (3) $$

Relative delay times are considered because, in this problem, the emission time of the source S is unknown, but the receivers are triggered synchronously.

### 2.2. Validation of travel time estimate of M1 through comparison with a finite difference modelling of the full wavefield

In this section, the P-wave first arrival travel time estimation proposed in this paper is validated through comparison with full wavefield modelling using finite differences.
The acoustic and isotropic version of the code developed by Operto and al. is used [19]
from the original formulation proposed by Jean Virieux [20, 21].

A comparison of the travel time estimate of M1 through the finite difference
modelling of the full P-wavefield is presented here, which is computationally fast and
accurate in the presence of flat topography, as in the test site analyzed in this study.

The following case study has been performed: 30 sensors were aligned and spaced
0.01 m apart; sensor 1 (R₁) was placed just above the source (Sₓ = 0); the rest of the sensors
were placed increasingly far away from the source; the depth of the source (Sz) was 1.5 m,
and the average velocity of propagation (V₀) was 500 m/s.

To ensure simplicity in the simulation, the origin of the reference frame was
translated, and the source was placed at the top and the sensors at the bottom, but they
were symmetrical. In the simulation presented here, the source was placed at S (0.35; 0),
and the sensors at Ri (0.35 + (i-1) 0.01; 1.5), but this corresponds well to the situation
described above. In Figure 2, we observe the physical model used for the simulation. The
source used was a Ricker function, also called the "mexican hat"; this function is a second
derivative of a Gaussian function.

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Physical velocity model in the case of the propagation of the full wavefield by finite differences with one propagation medium. S indicates the source at the surface. R₁ is the first receiver.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Comparison of travel time estimate of Model M1 (dashed blue line) with the full P-wavefield propagation model using finite differences.</figcaption></figure>
The results of this simulation are presented in Figure 3. The vertical axis represents the spatial position of the sensors, and the horizontal axis shows the recorded times. The simulation of the received signals using finite differences is shown in black and, the travel time curve calculated using the Model M1 is shown in blue. At the scale of the problem, the Model M1 is in agreement with the more complex full wavefield model using finite differences.

## 2.3. Cramer-Rao Bound from Model M1

The Cramer-Rao Bound (CRB) was calculated to obtain information on the estimation accuracy [27, 28]. The CRB represents the smallest possible standard deviations of all unbiased estimates of the model variables. First, the CRB was computed, and second, simulations using the CRB were performed.

### 2.3.1. Calculation of the Cramer-Rao Bound from Model M1

The CRB was calculated according to the Fisher information matrix, which is denoted by F.

$$ \mathrm{CRB}(\theta_{\mathrm{M}1}) = F^{-1}(\theta_{\mathrm{M}1}/\tau_{1i}(\theta_{\mathrm{M}1})) \quad (4) $$

where /.. represented the known operator. The Fisher information matrix is expressed

$$ F(\theta_{\mathrm{M}1}/\tau_{1i}(\theta_{\mathrm{M}1})) = \sum_{i=2}^{N} \frac{1}{\mathrm{Var}(\tau_{1i}(\theta_{\mathrm{M}1}))} \nabla_{\theta_{\mathrm{M}1}}(\tau_{1i}(\theta_{\mathrm{M}1})) \nabla_{\theta_{\mathrm{M}1}}^{\mathrm{T}}(\tau_{1i}(\theta_{\mathrm{M}1})) \quad (5) $$

where N is the number of sensors, Var(τ₁ᵢ(θₘ₁)) is the variance of the relative delay time between sensor 1 and i, and ∇θₘ₁[...] is the gradient operator as a function of θₘ₁. Let us focus on the calculation of the gradient of τ₁ᵢ(θₘ₁). The expression of this gradient can be decomposed from Equation (3),

$$ \nabla_{\theta_{\mathrm{M}1}}(\tau_{1i}(\theta_{\mathrm{M}1})) = \nabla_{\theta_{\mathrm{M}1}}\left(\frac{|\mathrm{SR}_i|}{V_0}\right) - \nabla_{\theta_{\mathrm{M}1}}\left(\frac{|\mathrm{SR}_1|}{V_0}\right) \quad (6) $$

By developing the calculation, the expression becomes

$$ \nabla_{\theta_{\mathrm{M}1}}(\tau_{1i}(\theta_{\mathrm{M}1})) = \begin{bmatrix} \frac{S_X - R_{iX}}{V_0|\mathrm{SR}_i|} & -\frac{S_X}{V_0|\mathrm{SR}_1|} \\ \frac{S_Z - R_{iZ}}{V_0|\mathrm{SR}_i|} & -\frac{S_Z}{V_0|\mathrm{SR}_1|} \\ \frac{|\mathrm{SR}_1| - |\mathrm{SR}_i|}{V_0^2} & 0 \end{bmatrix} \quad (7) $$

Using this result in Equation (5), the Fisher information matrix can be calculated, and by calculating its inverse, the CRB is obtained, as shown in Equation (4).

### 2.3.2. Numerical simulation using the Cramer-Rao Bound

The figures presented in this section are the results of the numerical simulations. These results are not exhaustive and are used to give an idea of the impact of the variables of the problem on the accuracy of the depth estimation. To obtain these curves, one variable was varied and the others were fixed. The impact on the possible accuracy of depth estimation was observed.

Figure 4.a shows the evolution of the CRB with changes in depth as a function of the error on the relative delay times. The situation considered for this simulation is as follows: five sensors were used, the distance between sensors was 0.2 m, the plumb of the
pipe ($S_x$) was placed at 0 m, and the average propagation velocity ($V_0$) was 500 m/s. Curves for two different depths ($S_z$), 0.4 m and 1 m, were observed. An accuracy of the relative delay times in a microsecond range was needed to obtain an accuracy of 0.1 m in the depth. In the worst case, an accuracy of the order of 0.1 µs was needed.

Figure 4.b represents the evolution of the CRB with changes in depth as a function of the error in the propagation velocity. The situation considered for this simulation is as follows: the distance between the sensors was 0.2 m, the plumb of the pipe ($S_x$) was 0 m, and the pipe depth ($S_z$) was 0.4 m. The curves for different average propagation velocities ($V_0$) were observed. A propagation velocity accuracy of the 10% range was required to obtain an error of less than 0.1 m in the depth.

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Depth error boundary as a function of the error of the other variables: (a) depth error as a function of the error in the relative delay times; (b) depth error as a function of the error in the propagation velocity.</figcaption></figure>

## 2.4. Adaptation of the MUSIC algorithm for Model M1

In the previous section, a relative delay time accuracy of the order of $10^{-6}$ s has been defined as a necessary condition. The choice to focus on the MUSIC algorithm (Multiple Signal Classification) [29-35] for antenna processing, the so-called 'high resolution algorithm' became interesting.

### 2.4.1. Presentation of the MUSIC algorithm adapted to our problem

The MUSIC algorithm is usually used to discern different sources and their direction of arrival through an antenna array [29]. It can be adapted in the near field to estimate the distance of sources in addition to their direction of arrival [30-35]. In this paper, the aim was to estimate the location of a single source $S$ in the near field and also to estimate the propagation velocity $V_0$ according to Model M1.

The algorithm can be divided into several steps as follows:

1. Estimate the variance-covariance matrix of the system from the signals received by sensors.

2. Decompose the variance-covariance matrix of the system into eigenvalues and eigenvectors.

3. Definition of the noise subspace, denoted by $U_b$, with the eigenvectors corresponding to the smallest eigenvalues.

4. Construction of a family of vectors, denoted by 'a', parametrized by the variables we want to estimate, $S$ and $V_0$ (Equation (1)). This family of vectors is constructed from the modelling of relative delay times (Equation (3)).
$$a(S, V_0) = \left[ 1 \ e^{-j2\pi f_0 \frac{|SR_2|-|SR_1|}{V_0}} \ \dots \ e^{-j2\pi f_0 \frac{|SR_N|-|SR_1|}{V_0}} \right]^T \quad (8)$$

where $f_0$ is the signal frequency, and $N$ is the number of sensors.

5. Knowing that the signal subspace and the noise subspace are orthogonal, the project of a $(S,V_0)$ on the noise subspace must be at the minimum for the values of $S$ and $V_0$ corresponding best to the received signal. It is traditional to take the inverse of this projection and to look for the values of $S$ and $V_0$ that maximize this criterion, which is denoted by $C_{\text{music}}$.

$$C_{\text{MUSIC}}(S, V_0) = \frac{1}{a(S, V_0)^H U_b U_b^H a(S, V_0)} \quad (9)$$

where $[.]^H$ is the conjugate transposed operator.

### 2.4.2. Test of estimator using numerical simulation

In this part, the algorithm worked in the simulation of an ideal case is verified. The statistics of the estimator (mean and variance) using the Monte Carlo method are presented. All estimates made with 1,000 runs of the noise and noise is applied to the propagation times. For these simulations, a situation has been fixed in which all the variables are known. The signals received are simulated as a function of these variables, and then these ones were used to run the algorithm.

The following situation is considered: a depth ($Sz$) of 0.7 m, a propagation velocity ($V_0$) of 500 m/s, and the plumb of the pipe ($Sx$) at 0 m. Five sensors were placed at 0.2 m, the first one vertically above the pipe.

In Table 1, the statistics of the results of the estimator, using the Monte Carlo method, are presented. 1,000 runs of noise were performed. The same situation as before was considered, and the theoretical propagation times between the source $S$ and sensors were noised with a white Gaussian noise of standard deviation $\sigma_{\text{noise}}$. The propagation times are in the millisecond range and the propagation time differences with sensor 1 vary between $10^{-5}$ and $10^{-4}$ seconds.

<table><caption>**Table 1.** Statistics of the MUSIC estimator adapted to Model M1 using the Monte Carlo method (for 1,000 runs, 5 sensors spaced 0.2 apart). Noise is applied to the propagation times which are in the order of a millisecond.</caption><thead><tr><th></th><th></th><th>Mean</th><th>Standard deviation</th><th>True Value</th></tr></thead><tbody><tr><td rowspan="3">Numerical Simulation 1<br>with 5 sensors and σ<sub>noise</sub> = 1.10<sup>-7</sup></td><td>Depth Sz (m)</td><td>0.7056</td><td>0.0486</td><td>0.7</td></tr><tr><td>Average velocity V<sub>0</sub> (m)</td><td>494</td><td>17</td><td>500</td></tr><tr><td>Plumb of the pipe S<sub>x</sub> (m)</td><td>0.0047</td><td>0.0065</td><td>0</td></tr><tr><td rowspan="3">Numerical Simulation 2<br>with 6 sensors and σ<sub>noise</sub> = 5.10<sup>-7</sup></td><td>Depth Sz (m)</td><td>0.6658</td><td>0.0317</td><td>0.7</td></tr><tr><td>Average velocity V<sub>0</sub> (m)</td><td>499</td><td>7</td><td>500</td></tr><tr><td>Plumb of the pipe S<sub>x</sub> (m)</td><td>0.0172</td><td>0.0109</td><td>0</td></tr><tr><td rowspan="3">Numerical Simulation 3<br>with 7 sensors and σ<sub>noise</sub> = 1.10<sup>-6</sup></td><td>Depth Sz (m)</td><td>0.7386</td><td>0.0466</td><td>0.7</td></tr><tr><td>Average velocity V<sub>0</sub> (m)</td><td>501</td><td>9</td><td>500</td></tr><tr><td>Plumb of the pipe S<sub>x</sub> (m)</td><td>0.0325</td><td>0.0315</td><td>0</td></tr></tbody></table>

For this set of variables, which is representative, the observed accuracy is lower than the desired accuracy of 0.1 m. It is interesting to note that increasing the number of sensors reduces the accuracy required on the relative delay times.

### 2.5. *Experimental measurements*

Experimental measurements were performed on a semi-controlled test area in order to work on real data.
2.5.1. Experimental set-up

The measurement chain presented in this section is not innovative; this
measurement chain was designed to be adapted to experimental needs. In particular, it
was necessary to be able to control the emission and to change the type of signal emitted
at any time. Moreover, in the reception chain, the filter should not be so fine that it could
not adapt to different types of signals that could be emitted.

The emission chain (Figure 5) was composed of a computer to control the emitted signal; an amplifier, which received the signal from the computer sound card; and finally, a loudspeaker. The loudspeaker was fixed at one end of the pipe and emitted the signal inside the pipe.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Experimental set-up of the acoustic emission chain: (a) diagram of the main components; (b) photograph of the material with the battery to power the amplifier.</figcaption></figure>

Experiments were carried out in an anechoic room to characterize the loudspeaker (CNRS-LMA-Marseille, C. Pinhède). A white noise is sent in command to the loudspeaker and the emitted signal is measured with a microphone placed just in front of the loudspeaker output. This experimentation allows to obtain the transfer function of the loudspeaker. The modulus and phase of its transfer function are shown in Figure 6.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Transfer function of the loudspeaker estimated by emitting a white noise: (a) modulus; (b) phase.</figcaption></figure>

The responses of the loudspeaker to the two types of signals used in this study are
presented.

The first signal of interest is a Ricker function. The time response of the loudspeaker to a Ricker function control signal is presented in Figure 7 and the frequency response in Figure 8. The properties of this type of signal are used in the calculation of the propagator in the simulation of the propagation of the complete wave field by finite differences. The real measurements with a Ricker function emission are to compare the trend of the delay times with the simulations.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Example of the loudspeaker response (Ricker function): (a) theoretical source signal (red line); (b) experimental source signal sent by the loudspeaker through the pipe (blue line) recorded during the controlled experiment.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Fourier transform of the Ricker function: (a) theoretical Ricker spectrum (red line) showing the central frequency of the theoretical source around 500 Hz; (b) experimental Ricker spectrum (loudspeaker output; blue line) showing a central frequency around 500 Hz but distorted in relation to the command.</figcaption></figure>

The second signal of interest is a monochromatic signal at 500Hz, a burst signal. It is used to estimate the pipe depth. The time response of the loudspeaker to a burst control signal is presented in Figure 9 and the frequency response in Figure 10.

The loudspeaker distorts the burst less than the Ricker function, because the Ricker function is richer in frequency.

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Example of the loudspeaker response (burst of 100 ms at 500 Hz): (a) theoretical source signal; (b) experimental source signal sent by the loudspeaker through the pipe.</figcaption></figure>

318
319
320

321
322
323
324
325
326
327
328
329
330

331
332
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Fourier transform of the burst of 100 ms at 500 Hz: (a) theoretical burst spectrum; (b) experimental burst spectrum (loudspeaker output).</figcaption></figure>

The acoustic reception chain is presented in Figure 11. The sensors are geophones that measure vibrations in their vertical axis. They are cylindrical with a diameter of 2.54 cm. During the measurements, a fatty substance was added between the sensors and the ground in order to obtain a better coupling. Geophones receive the signals that are filtered and amplified by electronic cards. At the input of electronic cards, signals are filtered using a high-pass filter with a cut-off frequency of 100 Hz, and the output, using a low-pass filter with a cut-off frequency of 100 kHz. This wide analog filtering allows us to digitally refine the filtering. The amplification is adjustable between 0 dB and 112 dB, which allow the amplification to be adjusted in the field. At the output of the electronic cards, the signals are digitized by the acquisition card; it is an analog-to-digital converter. The digitized signals are recovered on a computer.

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Experimental set-up of the reception chain: (a) scheme of the main components; (b) photograph of the material.</figcaption></figure>

Figure 12 shows an example of measurements taken on the semi-controlled test area. The trench and the sensors positioned perpendicular to the pipe passage can be observed. Sensors are spaced 5 cm apart.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Example of positioning of the experimental device for signal acquisition (reception chain). The black arrow indicates the acquisition design composed of a line of five sensors and the support (graduated ruler) where the sensors are moved during the experiment. The green arrow indicates the trench where a pipe is buried. The box contains the electronic cards presented in Figure 11 and, on the left of the box, the acquisition card connected to the computer to digitalize the signals received from the buried pipe is shown.</figcaption></figure>

### 2.5.2. Experimental results

Experimental measurements were performed on a semi-controlled test area. In this test area, the position of the pipe was known, and it was a polyethylene pipe, but the characteristics of the propagation medium were unknown. The propagation velocity is unknown.

To take the measurements, a line of five sensors was positioned perpendicular to the pipe route. To focus on the depth, it was assumed that the plumb of the pipe ($S_x$) was known, and sensor 1 was placed just above the pipe. Sensors were spaced 0.2 m apart. The signal emitted in the pipe was a monochromatic signal at 500 Hz.

On Figure 13, an example of signals measured by sensors is presented. The signal emitted by the loudspeaker is a burst of 100 ms at 500 Hz repeated every second. Figure 13.a shows the signals received by the sensors before digital filtering and Figure 13.b shows the same signals digitally filtered between 480 Hz and 520 Hz. Sensor 1 is the closest to the pipe and sensor 5 the farthest. The sensors farthest from the pipe receive the most attenuated signals.

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Example of signals measured by sensors (burst emission of 100 ms at 500 Hz): (a) signals received by sensors before digital filtering; (b) signals digitally filtered between 482 Hz and 520 Hz.</figcaption></figure>

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377
Results on real signals are presented, for which depth estimation with the desired accuracy of 0.1 m (Table 2 and Table 3) is reached. For a depth of 0.42 m, an accuracy of 0.03 m (Table 2) is obtained. For a depth of 0.7 m, an accuracy of 0.05 (Table 3) is obtained.

<table><caption>**Table 2.** First example of estimation obtained with the MUSIC algorithm adapted to the problem and to Model M1 (5 sensors spaced 0.2 m apart).</caption><thead><tr><td></td><td>Estimate value</td><td>Reference value</td></tr></thead><tbody><tr><td>Depth (Sz) (m)</td><td>0.39</td><td>0.42</td></tr><tr><td>Average velocity (V₀) (m/s)</td><td>360</td><td>unknown</td></tr></tbody></table>

<table><caption>**Table 3.** Second example of estimation obtained with the MUSIC algorithm adapted to the problem and to Model M1 (5 sensors spaced 0.2 m apart).</caption><thead><tr><td></td><td>Estimate value</td><td>Reference value</td></tr></thead><tbody><tr><td>Depth (Sz) (m)</td><td>0.75</td><td>0.70</td></tr><tr><td>Average velocity (V₀) (m/s)</td><td>540</td><td>unknown</td></tr></tbody></table>

Figure 14.a shows measured signals (filtered between 480 Hz and 520 Hz) on the test area that match with Model M1 (a zoom of Figure 13.b). With Model M1, the more receivers are far from the source S, the more they should receive in a delayed way the information of the sound wave propagating in the pipe. It is on real signals of this type that the results, presented in Tables 2 and 3, are obtained.

Figure 14.b shows measured signals on the test area that not match with Model M1. A phenomenon can be discerned several times, sensors farthest from the source S receive the signal before the closest sensors. This phenomenon shows the role played by the vertical discontinuities (the trench), not included in Model M1.

It is needed to evolve the model to cover a greater number of situations. This is why an evolution of this model is presented in part 3. Model M1 is sufficient when the transition is smooth or the medium weakly heterogeneous.

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Example of real received signals digitally filtered between 480 Hz and 520 Hz: (a) corresponding to the Model M1; (b) not corresponding to the Model M1.</figcaption></figure>
### 3. Modelling with two propagation media (M2)

In this part, an evolution of the model is presented. Now sensors are placed inside and outside the trench. A vertically stratified change in medium is considered, which represents the trench in which the pipe is buried. First, the modelling of the propagation times is presented, and then the calculation of the intermediate variables is presented more in detail.

#### 3.1. Propagation time modelling M2

##### 3.1.1. Presentation of the propagation time Model M2

In the problem modelling, a change of medium was added to take into account the trench (Figure 15). The pipe section was still represented using a point source S and the sensors $R_i$ were arranged on the ground surface in the same way as before. An average velocity $V_0$ was assigned to the propagation medium inside the trench and an average velocity $V_1$ outside the trench. For the sensors outside the trench, new intermediate variables $P_i$ appeared, which represented the interface points between the two media. For each $R_i$ outside the trench, there was a corresponding $P_i$.

<figure><img src="image_17.png" /><figcaption>Figure 15. Scheme of the problem modelling in the case of two propagation media: Model M2.</figcaption></figure>

We defined $\theta_{M2}$ as the vector of variables to be estimated in the case of Model M2.

$$ \theta_{M2} = [S_x \ S_z \ V_0 \ V_1]^T \qquad (10) $$

We noted $\tau_i(\theta_{M2})$ the wave travel time between the source S and the receiver $R_i$ for Model M2. Two cases can be discerned:

* If $R_i$ is inside the trench, it is returned to the single propagation medium case of Model M1 (Equation (2))

* If $R_i$ is outside the trench, then

$$ \tau_i(\theta_{M2}) = \frac{|S P_i|}{V_0} + \frac{|P_i R_i|}{V_1} \qquad (11) $$

In the following, only the case of $R_i$ outside the trench is considered because the case $R_i$ inside the trench was already covered in part 2.

Relative delay times are also considered because the emission time of the source S is unknown. We note $\tau_{1i}(\theta_{M2})$ the relative delay time between the sensors $R_i$ and $R_i$ for Model M2.

$$ \tau_{1i}(\theta_{M2}) = \frac{|S P_i| - |S R_i|}{V_0} + \frac{|P_i R_i|}{V_1} \qquad (12) $$
The modelling proposed in this paper depended on the interface point Pᵢ, which was unknown. We had to express Pᵢ as a function of θₘ₂.

### 3.1.2. Analytical expression of the interface point Pᵢ.

The variable τ₁ᵢ(θₘ₂) (Equation (12)) depended on the interface point Pᵢ (Pᵢₓ; Pᵢᵧ). It was assumed that the X coordinate Pᵢₓ was known a priori (i.e., we knew the position of the medium change). However, we had no a priori information on the Z coordinate Pᵢᵧ. We focused on the analytical expression of Pᵢₓ.

The Snell–Descartes law was considered at the interface between the two media. We can then write

$$ \frac{S_z - P_{iz}}{|S P_i|} = \frac{P_{iz}}{|P_i R_i|} \frac{V_0}{V_1} \quad (13) $$

After squaring and developing, the following fourth order equation was obtained:

$$ \begin{aligned} & P_{iz}^4 [V_1^2 - V_0^2] \\ & + P_{iz}^3 [-2R_{iz}V_1^2 - 2S_zV_1^2 + 2S_zV_0^2] \\ & + P_{iz}^2 [S_z^2V_1^2 + R_{iz}^2V_1^2 + V_1^2(R_{ix} - P_{ix})^2 + 4R_{iz}S_zV_1^2 - S_z^2V_0^2 - V_0^2(P_{ix} - S_x)^2] \\ & + P_{iz} [-2R_{iz}S_z^2V_1^2 - 2S_zR_{iz}^2V_1^2 - 2S_zV_1^2(R_{ix} - P_{ix})^2] \\ & + [R_{iz}^2S_z^2V_1^2 + S_z^2V_1^2(R_{ix} - P_{ix})^2] = 0 \end{aligned} \quad (14) $$

We propose

$$ P_{iz}^4 m_1 + P_{iz}^3 m_2 + P_{iz}^2 m_3 + P_{iz} m_4 + m_5 = 0 \quad (15) $$

After solving Equation (15), we obtained

$$ P_{iz} = \frac{\sqrt{A} - \sqrt{A - 2(g + A + \frac{u}{\sqrt{A}})}}{2} - \frac{m_2}{4m_1} \quad (16) $$

with

$$ A = \left( \frac{-K + \sqrt{K^2 + \frac{4J^3}{27}}}{2} \right)^{1/3} - \frac{J}{3 \left( \frac{-K + \sqrt{K^2 + \frac{4J^3}{27}}}{2} \right)^{1/3}} - \frac{2g}{3}, \quad (17) $$

$$ J = -\frac{g^2}{3} - 4w, \quad (18) $$

$$ K = \frac{8wg}{3} - \frac{2g^3}{27} - u^2, \quad (19) $$

$$ g = \frac{m_3}{m_1} - \frac{3m_2^2}{8m_1^2}, \quad (20) $$

$$ u = \frac{m_4}{m_1} - \frac{m_2 m_3}{2m_1^2} + \frac{m_2^3}{8m_1^3}, \quad (21) $$

$$ w = \frac{m_5}{m_1} - \frac{m_2 m_4}{4m_1^2} + \frac{m_2^2 m_3}{16m_1^3} - \frac{3m_2^4}{256m_1^4} \quad (22) $$
$m_1$ to $m_5$ depended on $\theta_{M2}$ (Equation (10)), the position of sensors $R_i$ (known) and the position of the medium change $P_{ix}$ (known). The expression of the interface point depth $P_{iz}$ depended on $m_i$, so we expressed $P_{iz}$ as a function of variables of interest $\theta_{M2}$. The theoretical delay times with the Model M2 from $\theta_{M2}$ can be estimated following this approach.

## 3.2. Validation of travel time estimate of M2 through comparison with a finite difference modelling of the full wavefield

Here, the travel times of the Model M2 are compared with the full P-wavefield propagation modelling using finite differences, in the same way as in Section 2.2. The simulations were performed for the same case study in order to compare the results. The following set-up was considered: sensor 1 is the origin of the reference frame; the point source $S$ is at a depth of $S_z = 0.7$ m; sensor 1 is in line with the source $S_x = 0$ m; sensors are spatially distributed every 0.05 m and move away from the source; the position of the change in medium is at $P_x = 0.1$ m; the average propagation velocity inside the trench is $V_0 = 300$ m/s and outside the trench is $V_1 = 600$ m/s. These values were proposed for the propagation velocities because they were close to those that would be estimated for the real data in the test site.

To ensure simplicity in the simulation, the origin of the reference frame was translated and the source was placed at the top and the sensors at the bottom, but they were symmetrical. In the simulation presented here, the source was placed at $S$ (0.75; 0) and the sensors at $R_i$ (0.75 + (i-1)0.05; 0.7), but this corresponds well to the situation described above. In Figure 16, the physical model used for simulation, and the wavefront propagation in the case of the full wavefield propagation modelling using finite differences are observed. Each panel is taken at a different time and the evolution of the wavefront over time is observed. In this simulation, the source used was a Ricker function (Figure 7.a, Figure 8.a).

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Physical velocity model and wavefront in the case of the full P-wavefield propagation modelling using the finite differences method.</figcaption></figure>

The results of this simulation are presented in Figure 17. The abscissa represents the spatial position of the sensors, and the ordinate indicates the recorded time. The seismograms represent the simulation of the received signals using finite differences, and the travel time curve calculated from Model M2 is shown in blue. At the scale of the problem, Model M2 was in agreement with the more complex full wavefield model using finite differences. This allowed us to obtain a first validation of the travel time estimate of Model M2 through numerical simulation.
<figure><img src="image_19.png" /><figcaption>**Figure 17.** Seismograms obtained by full wavefield propagation using the finite differences method in the physical velocity model presented in Figure 16. The horizontal axis indicates the position of the sensors. The red arrow indicates the sensors located at the plumb of the source. The theoretical travel times calculated in the same model (M2) by the method proposed in section 3.1 are plotted with blue dashed line to assess the good agreement. The transmitted signal is considered as the same emitted Ricker function (Figure 7.a).</figcaption></figure>

The same simulation is also carried out but using as source signal the signal actually transmitted by the loudspeaker (Figure 7.b and Figure 8.b). On Figure 18, a small shift between the travel times calculated with Model M2 and the simulations of the signals received using finite differences can be observed, but the trend of the curve remains correct.

<figure><img src="image_20.png" /><figcaption>**Figure 18.** Seismograms obtained by full wavefield propagation using the finite differences method in the physical velocity model presented in Figure 16. The horizontal axis indicates the position of the sensors. The red arrow indicates the sensors located at the plumb of the source. The theoretical travel times calculated in the same model (M2) by the method proposed in section 3.1 are plotted with blue dashed line to assess the good agreement. We consider that the emitted signal was transformed by the loudspeaker (Figure 7.b).</figcaption></figure>

### 3.3. Cramer-Rao Bound from Model M2

In this section, as in Section 2.3, the Cramer-Rao Bound (CRB) was calculated to obtain information on the estimation accuracy.

#### 3.3.1. Calculation of the Cramer-Rao Bound from Model M2

The CRB was calculated according to the Fisher information matrix, which is denoted by F.

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512
$$ \mathrm{CRB}(\theta_{\mathrm{M2}}) = F^{-1}(\theta_{\mathrm{M2}}/\tau_{1i}(\theta_{\mathrm{M2}})) \qquad (23) $$

where /./ is the known operator.

The Fisher information matrix is expressed

$$ F(\theta_{\mathrm{M2}}/\tau_{1i}(\theta_{\mathrm{M2}})) = \sum_{i=2}^{N} \frac{1}{\mathrm{Var}(\tau_{1i}(\theta_{\mathrm{M2}}))} \nabla_{\theta_{\mathrm{M2}}} (\tau_{1i}(\theta_{\mathrm{M2}})) \nabla_{\theta_{\mathrm{M2}}}^{\mathrm{T}} (\tau_{1i}(\theta_{\mathrm{M2}})) \qquad (24) $$

where N is the number of sensors, Var(τ₁ᵢ(θₘ₂)) is the variance of the relative delay time between sensor 1 and i, and ∇θₘ₂[·] is the gradient operator as a function of θₘ₂.

The gradient of τ₁ᵢ(θₘ₂) is calculated from Equation (12),

$$ \nabla_{\theta_{\mathrm{M2}}} (\tau_{1i}(\theta_{\mathrm{M2}})) = \nabla_{\theta_{\mathrm{M2}}} \left( \frac{|S P_i|}{V_0} \right) - \nabla_{\theta_{\mathrm{M2}}} \left( \frac{|S R_i|}{V_0} \right) + \nabla_{\theta_{\mathrm{M2}}} \left( \frac{|P_i R_i|}{V_1} \right) \qquad (25) $$

The calculation of the gradient can be decomposed as follows

$$ \nabla_{\theta_{\mathrm{M2}}} \left( \frac{|SP_i|}{V_0} \right) = \left[ \frac{(S_X - P_{1X}) - (S_Z - P_{1Z}) \frac{\partial}{\partial S_X} P_{1Z}}{V_0 |SP_i|} \quad \frac{(S_Z - P_{1Z}) (1 - \frac{\partial}{\partial S_Z} P_{1Z})}{V_0 |SP_i|} \quad \frac{-V_0 (S_Z - P_{1Z}) \frac{\partial}{\partial V_0} P_{1Z} - |SP_i|^2}{V_0^2 |SP_i|} \quad \frac{-(S_Z - P_{1Z}) \frac{\partial}{\partial V_1} P_{1Z}}{V_0 |SP_i|} \right]^T \qquad (26) $$

$$ \nabla_{\theta_{\mathrm{M2}}} \left( \frac{|SR_1|}{V_0} \right) = \left[ \frac{S_X - R_{1X}}{V_0 |S R_1|} \quad \frac{S_Z - R_{1Z}}{V_0 |S R_1|} \quad \frac{-|S R_1|}{V_0^2} \quad 0 \right]^T \qquad (27) $$

$$ \nabla_{\theta_{\mathrm{M2}}} \left( \frac{|P_i R_i|}{V_1} \right) = \left[ \frac{(P_{1Z} - R_{1Z}) \frac{\partial}{\partial S_X} P_{1Z}}{V_1 |P_i R_i|} \quad \frac{(P_{1Z} - R_{1Z}) \frac{\partial}{\partial S_Z} P_{1Z}}{V_1 |P_i R_i|} \quad \frac{(P_{1Z} - R_{1Z}) \frac{\partial}{\partial V_0} P_{1Z}}{V_1 |P_i R_i|} \quad \frac{V_1 (P_{1Z} - R_{1Z}) \frac{\partial}{\partial V_1} P_{1Z} - |P_i R_i|^2}{V_1^2 |P_i R_i|} \right]^T \qquad (28) $$

All these expressions depended on the gradient of Pᵢz. The gradient of Pᵢz, ∇θₘ₂Pᵢz, needed to be calculated to be able to calculate the CRB. This gradient calculation is presented in Appendix A to avoid overloading the text.

### 3.3.2. Numerical simulation using the Cramer-Rao Bound from Model M2

The figures presented in this section are the result of numerical simulations. These results are not exhaustive and are used to give an idea of the impact of the variables of the problem on the accuracy of the depth estimation.

<figure><img src="image_21.png" /><figcaption>**Figure 19.** Depth error as a function of the error in the relative delay times: (a) for several numbers of sensors; (b) for several depths.</figcaption></figure>

Figure 19 shows the evolution of the CRB on the depth as a function of the error on the relative delay times. For these simulations, the following situation is considered: the distance between sensors is 0.2 m; the plumb of the pipe (Sₓ) is at 0 m; the average propagation velocity inside the trench (V₀) is 300 m/s; the average propagation velocity outside the trench (V₁) is 600 m/s; and the position of the change in medium (Pₓ) is at 0.15 m.
For Figure 19.a, a depth of 0.7 m is considered, and the curves for different numbers of sensors are observed. The accuracy constraint on the relative delay times is relaxed by increasing the number of sensors.

For Figure 19.b, seven sensors are considered, and the curves for different depths are observed. The estimation of relative delay times should be more accurate as the depth increases. For example, here, with seven sensors, to obtain a precision of 0.1 m for the depth, the precision of the relative delay times is needed to be in a 0.1 µs range.

### 3.4. Validation of the travel time Model M2 on real data

To compare travel time modelling with real data, an experiment was performed on the semi-controlled test area in a case similar to the simulation in Section 3.2. The same signal as in the simulation, a Ricker, was emitted. Sensors were placed every 0.05 m perpendicular to the pipe passage (as in Figure 15). Since we did not have a large enough number of sensors, several successive measurements were performed by leaving one sensor fixed and moving the others to mesh the space and act as if there was a large number of sensors. The time bases were recalculated with respect to the fixed sensor of each measurement to act as if only one measurement had been realized with many sensors. This added an error in the travel times, but here, we were only interested in the trend of the evolution of the relative delay times in order to validate Model M2.

Sensor 1 was taken as a reference, which was the closest to the source S, to calculate the relative delay times. Figure 20 shows the evolution of the relative delay times between sensors 1 and i in black. These relative delay times were estimated using cross-correlation between the signal received by sensor 1 and that of sensor i.

<figure><img src="image_22.png" /><figcaption>**Figure 20.** Comparison of the relative delay times between the real data and Model M2.</figcaption></figure>

The propagation velocities were unknown, a priori, for the semi-controlled test area. The velocities $V_0$ and $V_1$ of Model M2 can be estimated, a posteriori, by knowing the position of the source S in the test area. After determining all the variables in the test zone a posteriori, the relative delay times of Model M2 (blue curve of Figure 20) were calculated. At the beginning, the relative delay times increased until the break that marked the change in environment. Then, they started to increase again. In the simulation (blue curve), the position of the change in medium was fixed at 0.1 m. On the real data (black curve), the change in medium, highlighted by the break of the curve, was instead between 0.15 and 0.2 m. The Model M2 is representative of reality. It is logical and is in agreement with the observation of real data.

### 3.5. Depth estimation from Model M2

#### 3.5.1. Numerical simulation with MUSIC algorithm adapted to the Model M2

534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558

559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
In this section, numerical simulations are performed to qualify the MUSIC estimator adapted to the Model M2 (as in Section 2.4.2). The statistics of the estimator using the Monte Carlo method (1,000 runs of noise) are presented. Consider the following situation 1: a depth ($Sz$) of 0.7 m, a propagation velocity ($V_0$) of 300 m/s, a propagation velocity ($V_1$) of 600 m/s, the plumb of the pipe ($S_x$) at 0 m, and the position of the change in medium ($P_x$) at 0.15 m. Seven sensors were spaced 0.2 m apart, the first one directly at the plumb of the pipe.

In Table 4, the statistics of the results of the estimator using the Monte Carlo method are presented. 1,000 runs of noise were performed. The theoretical propagation times between the source S and sensors were noised with a white Gaussian noise of standard deviation $\sigma_{\text{noise}}$. The propagation times are in the millisecond range and the propagation time differences with sensor 1 vary between $10^{-5}$ and $10^{-3}$ seconds. In this situation, an accuracy in a microsecond range is required on the delay times to get an accuracy of less than 0.1 m on the depth. These results confirm those presented in Section 3.3.2.

<table><caption>**Table 4.** Statistics of the MUSIC estimator adapted to Model M2 using the Monte Carlo method in situation 1 (7 sensors spaced 0.2 m apart). Noise is applied to the propagation times which are in the order of a millisecond</caption><thead><tr><th></th><th></th><th>Mean</th><th>Standard deviation</th><th>True Value</th></tr></thead><tbody><tr><td rowspan="3">Numerical Simulation 1<br>with σ<sub>noise</sub> = 5.10<sup>-6</sup></td><td>Depth Sz (m)</td><td>0.7294</td><td>0.0849</td><td>0.7</td></tr><tr><td>Average velocity V<sub>0</sub> (m)</td><td>293</td><td>19</td><td>300</td></tr><tr><td>Average velocity V<sub>1</sub> (m)</td><td>592</td><td>41</td><td>600</td></tr><tr><td rowspan="3">Numerical Simulation 2<br>with σ<sub>noise</sub> = 5.10<sup>-5</sup></td><td>Depth Sz (m)</td><td>0.6666</td><td>0.2246</td><td>0.7</td></tr><tr><td>Average velocity V<sub>0</sub> (m)</td><td>261</td><td>58</td><td>300</td></tr><tr><td>Average velocity V<sub>1</sub> (m)</td><td>638</td><td>128</td><td>600</td></tr></tbody></table>

Now, the case with seventeen sensors spaced 0.05 m apart is presented. Consider the following situation 2: a depth ($Sz$) of 0.7 m, a propagation velocity ($V_0$) of 300 m/s, a propagation velocity ($V_1$) of 600 m/s, the plumb of the pipe ($S_x$) at 0 m, and the position of the change in medium ($P_x$) at 0.15 m. Seventeen sensors were spaced 0.05 m apart, the first one directly at the plumb of the pipe.

<table><caption>**Table 5.** Statistics of the MUSIC estimator adapted to Model M2 using the Monte Carlo method in situation 2 (17 sensors spaced 0.05 m apart). Noise is applied to the propagation times which are in the order of a millisecond</caption><thead><tr><th></th><th></th><th>Mean</th><th>Standard deviation</th><th>True Value</th></tr></thead><tbody><tr><td rowspan="3">Numerical Simulation 1<br>with σ<sub>noise</sub> = 5.10<sup>-6</sup></td><td>Depth Sz (m)</td><td>0.6604</td><td>0.0704</td><td>0.7</td></tr><tr><td>Average velocity V<sub>0</sub> (m)</td><td>293</td><td>18</td><td>300</td></tr><tr><td>Average velocity V<sub>1</sub> (m)</td><td>611</td><td>79</td><td>600</td></tr></tbody></table>

The results in Table 5 show that an error of the order of microseconds on delay times allows an accuracy of 0.1 m on the depth.

### 3.5.2. Estimation on real data

At least seven sensors are required (section 3.3.2.) to use the MUSIC algorithm adapted to the Model M2. As mentioned in section 3.4, we only had five sensors. The sensors were moved to perform successive measurements. The different measurements are then not synchronized. Therefore, our estimator, based on MUSIC, cannot be used on these measurements.

However, results evaluated from the delay times between the sensors are presented in Figure 20 (black curve). These delay times obtained after retiming, introducing additional errors.
The criterion presented here is the inverse of the squared error between the estimated delay times and the delay times calculated by the Model M2.

The values of the criterion are obtained by varying the variables to be estimated: the depth ($Sz$), the velocity inside the trench ($V_0$) and the velocity outside the trench ($V_1$). The lateral position of the pipe ($S_x$) and the position of the medium change ($P_x$) are considered known. Indeed, $S_x$ estimated with the GasTracker tool and $P_x$ from the temporal break observed in Figure 20.

The Figure 21 is a display of the obtained criterion. The estimation results are presented in Table 6. The error on the depth is less than 0.1 m. In addition, the velocity estimates are really small, which could be due to the large error on the delay times. The parameters $S_x$ and $P_x$ are fixed by an a priori estimate. It would be interesting to also vary these parameters around their estimated position.

<table><caption>**Table 6.** Estimation obtained with least squares from the Model M2.</caption><thead><tr><td></td><td>Estimate value</td><td>Reference value</td></tr></thead><tbody><tr><td>Depth ($Sz$) (m)</td><td>0.64</td><td>0.70</td></tr><tr><td>Average velocity inside the trench ($V_0$) (m/s)</td><td>140</td><td>unknown</td></tr><tr><td>Average velocity outside the trench ($V_1$) (m/s)</td><td>250</td><td>unknown</td></tr></tbody></table>

<figure><img src="image_23.png" /><figcaption>**Figure 21.** The inverse of the squared error between the estimated delay times and the delay times calculated by the Model M2. For this display the velocity outside the trench $V_1$ is fixed at the estimated value.</figcaption></figure>

These results obtained with least squares are encouraging but we would like to have a more discriminating criterion. The next step of our work will be to test the MUSIC algorithm on real synchronous data in order to have a criterion with a better resolution power.

## 4. Conclusion and future perspectives

We developed a method using several sensors positioned perpendicular to the passage of the pipe. This allows us to model the problem in a plane orthogonal to the passage of the pipe. The vibrating section of the pipe is represented by a point source.
In a first part, an initial model of the problem was defined by considering a single propagation medium (M1), i.e all sensors are located at the ground surface inside the trench. This Model M1 was validated, on the scale of the problem, by comparing it to a more complete modelling (full wavefield modelling by finite differences). The Cramer-Rao bound was calculated to obtain theoretical information on the accuracy reached by the variables, in particular on the relative delay times. The error on the positioning of the sensors (e.g. line of sensors not perpendicular to the pipe) induces an error on the travel times. This is why, in these simulations and theoretical studies, the error on the depth estimation was quantified as a function of the error on the delay times. The MUSIC algorithm was adapted to Model M1, and then the algorithm was tested on real data. This allowed us to question the model and thus to advance it. The results obtained from the M1 model show that the depth estimates reach the desired accuracy of 0.1m.

In the second part, a second model of the problem was defined by considering two vertically stratified propagation media (M2), i.e sensors are located at the ground surface inside and outside the trench, according to the information obtained from experiments onto a test area. This model was validated in two steps, first through numerical simulation by comparing it to a more complete model, then by comparing it to real data. The Cramer-Rao bound was also calculated to obtain theoretical information on the accuracy reached by the variables. The MUSIC algorithm was adapted to this Model M2 and was tested in simulation.

The Model M2 aims to cover the case of a marked change of medium (inside/outside the trench). The results obtained with the Cramér-Rao bound show that 7 or 8 sensors would be needed to obtain a satisfactory depth estimate. Not having real synchronous measurements realized with 8 sensors, a least squares criterion between the delay times estimated from the real signals and those given by the Model M2 was established. The results obtained on real data reach the desired accuracy.

In our future work, with more sensors, we will test the estimator based on MUSIC algorithm adapted to Model M2 on real data. Depending on the feedback from these experiments, we will be able to continue to evolve the model by potentially considering other propagation media, such as the layer of sand surrounding the pipe or the layer of tar covering the ground. This aims to produce a model closer to real situations that provides the fastest possible result to be able to apply the method in the field in real time.

**Author Contributions:** Conceptualization, W.X. and G.S.; methodology, W.X.; software, W.X. and A.R.; validation, W.X., G.S., A.R., P.P. and L.Z.; formal analysis, W.X.; investigation, W.X.; resources, P.P. and L.Z.; data curation, W.X.; writing—original draft preparation, W.X.; writing—review and editing, W.X., G.S. and A.R.; visualization, W.X. and A.R.; supervision, G.S.; project administration, not applicable; funding acquisition, not applicable. All authors have read and agreed to the published version of the manuscript.

**Funding:** This project is supported by the CIFRE contract between MADE-SA and CNRS-UMR7330, CEREGE, from the ANRT French program of CNRS (Agence National de la Recherche Technique). This project is partially supported by RICE (Research & Innovation Center of Energy), GRTgaz.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** We would like to thank the Laboratory of Mechanics and Acoustics (LMA - CNRS), Marseille, France. We also express our sincere gratitude to Cédric Pinhède, research engineer at LMA, who allowed us to access the anechoic chamber of LMA and to qualify our loudspeaker.

We gratefully acknowledge S. Operto for providing the wave propagation finite difference modelling code used in this study.

**Conflicts of Interest:** The authors declare no conflicts of interest.
# Appendix A

## Calculation of the gradient of Pᵢz

In this section, the gradient of Pᵢz is expressed with respect to θₘ₂. The notations from Section 3.1.2 are used in the calculation of Pᵢz.

From Equation (16),

$$ \nabla_{\theta_{M2}} P_{iz} = \frac{1}{2} \left( \frac{1}{2\sqrt{A}} \nabla_{\theta_{M2}}(A) - \frac{\nabla_{\theta_{M2}}(A) - 2\left[\nabla_{\theta_{M2}}(g) + \nabla_{\theta_{M2}}(A) + \nabla_{\theta_{M2}}\left(\frac{u}{\sqrt{A}}\right)\right]}{2\sqrt{A - 2\left(g + A + \frac{u}{\sqrt{A}}\right)}} \right) - \frac{m_1 \nabla_{\theta_{M2}}(m_2) - m_2 \nabla_{\theta_{M2}}(m_1)}{4 m_1^2} \quad (A. 1) $$

The gradients of all intermediate variables are calculated (A, J, K, g, u, w, m₁, m₂, m₃, m₄, m₅) to express ∇θₘ₂Pᵢz. The variable A (Equation (17)) is decomposed into three terms to ensure that the calculation was still readable.

$$ A = A_1 + A_2 + A_3 \quad (A. 2) $$

with

$$ A_1 = \left( \frac{-K + \sqrt{K^2 + \frac{4J^3}{27}}}{2} \right)^{1/3}, \quad A_2 = - \frac{J}{3 \left( \frac{-K + \sqrt{K^2 + \frac{4J^3}{27}}}{2} \right)^{1/3}}, \quad A_3 = - \frac{2g}{3} \quad (A. 3) $$

The gradient of each term is calculated.

$$ \nabla_{\theta_{M2}}(A_1) = \frac{-\nabla_{\theta_{M2}}(K) + \frac{2K\nabla_{\theta_{M2}}(K) + \frac{12J^2}{27}}{2\sqrt{K^2 + \frac{4J^3}{27}}} \nabla_{\theta_{M2}}(J)}{6 \left( \frac{-K + \sqrt{K^2 + \frac{4J^3}{27}}}{2} \right)^{2/3}} \quad (A. 4) $$

$$ \nabla_{\theta_{M2}}(A_2) = -\frac{\nabla_{\theta_{M2}}(J) \left( \frac{-K + \sqrt{K^2 + \frac{4J^3}{27}}}{2} \right)^{1/3} - J\nabla_{\theta_{M2}}(A_1)}{\left( \frac{-K + \sqrt{K^2 + \frac{4J^3}{27}}}{2} \right)^{2/3}} \quad (A. 5) $$

$$ \nabla_{\theta_{M2}}(A_3) = -\frac{2}{3}\nabla_{\theta_{M2}}(g) \quad (A. 6) $$

Once the gradient of A was calculated, the gradients of J and K can be calculated (Equation (18) and (19)).

$$ \nabla_{\theta_{M2}}(J) = -\frac{2g}{3} \nabla_{\theta_{M2}}(g) - 4\nabla_{\theta_{M2}}(w) \quad (A. 7) $$

$$ \nabla_{\theta_{M2}}(K) = \frac{8}{3} \left[ g \nabla_{\theta_{M2}}(w) + w \nabla_{\theta_{M2}}(g) \right] - \frac{2g^2}{9} \nabla_{\theta_{M2}}(g) - 2u \nabla_{\theta_{M2}}(u) \quad (A. 8) $$
The calculation of the gradients of g, u and w can be realized (Equation (20), (21) and (22)).

$$ \nabla_{\theta_{M2}}(g) = \frac{m_1 \nabla_{\theta_{M2}}(m_3) - m_3 \nabla_{\theta_{M2}}(m_1)}{m_1^2} - \frac{3m_2 m_1 \nabla_{\theta_{M2}}(m_2) - m_2^2 \nabla_{\theta_{M2}}(m_1)}{m_1^3} \quad (A. 9) $$

$$ \nabla_{\theta_{M2}}(u) = \frac{m_1 \nabla_{\theta_{M2}}(m_4) - m_4 \nabla_{\theta_{M2}}(m_1)}{m_1^2} - \frac{1}{2} \frac{[m_2 \nabla_{\theta_{M2}}(m_3) + m_3 \nabla_{\theta_{M2}}(m_2)]m_1 - 2m_3 m_2 \nabla_{\theta_{M2}}(m_1)}{m_1^3} + \frac{3m_1 m_2^2 \nabla_{\theta_{M2}}(m_2) - m_2^3 \nabla_{\theta_{M2}}(m_1)}{m_1^4} \quad (A. 10) $$

$$ \nabla_{\theta_{M2}}(w) = \begin{cases} \frac{m_1 \nabla_{\theta_{M2}}(m_5) - m_5 \nabla_{\theta_{M2}}(m_1)}{m_1^2} - \frac{1}{4} \frac{[m_2 \nabla_{\theta_{M2}}(m_4) + m_4 \nabla_{\theta_{M2}}(m_2)]m_1 - 2m_2 m_4 \nabla_{\theta_{M2}}(m_1)}{m_1^3} & (A. 11) \\ + \frac{1}{16} \frac{[m_2^2 \nabla_{\theta_{M2}}(m_3) + 2m_3 m_2 \nabla_{\theta_{M2}}(m_2)]m_1 - 3m_3 m_2^2 \nabla_{\theta_{M2}}(m_1)}{m_1^4} - \frac{3}{64} \frac{m_1 m_2^3 \nabla_{\theta_{M2}}(m_2) - m_2^4 \nabla_{\theta_{M2}}(m_1)}{m_1^5} & \end{cases} $$

Then, the gradients from $m_1$ to $m_5$, on which all other gradients depended, were calculated.

$$ \nabla_{\theta_{M2}}(m_1) = [0 \quad 0 \quad -2V_0 \quad 2V_1]^T \quad (A. 12) $$

$$ \nabla_{\theta_{M2}}(m_2) = [0 \quad 2[V_0^2 - V_1^2] \quad 4S_Z V_0 \quad -4V_1[R_{iZ} + S_Z]]^T \quad (A. 13) $$

$$ \nabla_{\theta_{M2}}(m_3) = [2V_0^2(P_{iX} - S_X)] \quad 2S_Z V_1^2 + 4R_{iZ} V_1^2 - 2S_Z V_0^2 \quad -2V_0[S_Z^2 + (P_{iX} - S_X)^2] \quad 2V_1[S_Z^2 + R_{iZ}^2 + (R_{iX} - P_{iX})^2 + 4R_{iZ} S_Z]^T \quad (A. 14) $$

$$ \nabla_{\theta_{M2}}(m_4) = [0 \quad -2V_1^2[2R_{iZ}S_Z + R_{iZ}^2 + (R_{iX} - P_{iX})^2] \quad 0 \quad -4V_1S_Z[R_{iZ}S_Z + R_{iZ}^2 + (R_{iX} - P_{iX})^2]]^T \quad (A. 15) $$

$$ \nabla_{\theta_{M2}}(m_5) = [0 \quad 2S_Z V_1^2 [R_{iZ}^2 + (R_{iX} - P_{iX})^2] \quad 0 \quad 2S_Z^2 V_1 [R_{iZ}^2 + (R_{iX} - P_{iX})^2]]^T \quad (A. 16) $$

The gradients of $m_1$ to $m_5$ were calculated and expressed according to the gradients of $P_{iZ}$.

## References

1. Liu, Y.; Habibi, D.; Chai, D.; Wang, X.; Chen, H.; Gao, Y.; Li, S. A Comprehensive Review of Acoustic Methods for Locating Underground Pipelines. Applied sciences **2020**, 10, 1031, DOI. 10.3390/app10031031

2. Jeong, H.S.; Arboleda, C.A.; Abraham, D.M.; Halpin, D.W.; Bernold, L.E. Imaging and Locating Buried Utilities. Joint Transportation Research Program **2003**, DOI. 10.5703/12882843131237

3. Harmankaya, U.; Kaslilar, A.; Thorbecke, J.; Wapenaar, K.; Draganov, D. Locating near-surface scatterers using non-physical scattered waves resulting from seismic interferometry. Journal of Applied Geophysics **2013**, 91, 66-81, DOI. 10.1016/j.jappgeo.2013.02.004

4. Filippi, C.; Leparoux, D.; Grandjean, G.; Bitri, A.; Côte, P. New robust observables on Rayleigh waves affected by an underground cavity: from numerical to experimental modelling. Geophysical Journal International **2019**, 218, 1903-1918, DOI. 10.1093/gji/ggz256

5. Muggleton, M.J.; Rustighi, E. Mapping the Underworld: recent developments in vibro-acoustic techniques to locate buried infrastructure. Géotechnique Letters **2013**, 3, 137-141, DOI. 10.1680/geolett.13.00032

6. Hao, T.; Rogers, C.D.F.; Metje, N.; Chapman, D.N.; Muggleton, J.M.; Foo, K.Y.; Wang, P.; Pennock, S.R.; Atkins, P.R.; Swingler, S.G.; and al. Condition assessment of the buried utility service infrastructure. Tunnelling and Underground Space Technology **2012**, 28, 331-344, DOI. 10.1016/j.tust.2011.10.011

7. Dou, Q.; Wei, L.; Magee, D.R.; Atkins, P.R.; Chapman, D.N.; Curioni, G.; Goddard, K.F.; Hayati, F.; Jenks, H.; Metje, N.; Muggleton, J.; Pennock, S.R.; Rustighi, E.; Swingler, S.G.; Rogers, C.D.F.; Cohn, A.G. 3D Buried Utility Location Using A Marching-Cross-Section Algorithm for Multi-Sensor Data Fusion. Sensors **2016**, 16, 1827, DOI. 10.3390/s16111827

8. Li, J.; Guo, T.; Leung, H.; Xu, H.; Liu, L.; Wang, B.; Liu, Y. Locating Underground Pipe Using Wideband Chaotic Ground Penetrating Radar. Sensors **2019**, 19, 2913, DOI. 10.3390/s19132913

9. Pettinelli, E.; di Matteo, A.; Mattei, E.; Crocco, L.; Soldovieri, F.; Redman, J.D.; Annan, A.P. GPR Response from Buried Pipes: Measurement on Field Site and Tomographic Reconstructions. IEEE Transactions Geoscience and Remote Sensing **2009**, 47, 2639-2645, DOI. 10.1109/TGRS.2009.2018301

10. Xerri, W.; Saracco, G.; Gassier, G.; Zomero, L.; Picon, P. Preliminary Acoustic Study of 3D Localization of Buried Polyethylene Pipe. July 28-30 2021, 48th Annual Review of Progress in Quantitative Nondestructive Evaluation **2022**, DOI.10.1115/QNDE2021-74945

11. Xerri, W.; Saracco, G.; Ribodetti, A.; Zomero, L.; Picon, P. Estimating the depth of buried polyethylene pipe by comparing propagation models. In Proceedings of French Congress of Acoustics (CFA), April 11-15 2022, Marseille, France.
736 12. Liu, Y.; Habibi, D.; Chai, D.; Wang, X.; Chen, H. A Numerical Study of Axisymmetric Wave Propagation in Buried Fluid-Filled Pipes for Optimizing the Vibro-Acoustic Technique When Locating Gas Pipelines. *Energies* **2019**, *12*, 3707, DOI. 10.3390/en12193707

739 13. Papandreou, B.; Rustighi, E.; Brennan, M.J. Study into the Feasibility of Using Acoustic Techniques to Locate Buried Objects; ISVR Technical Memorandum, Institute of Sound and Vibration Research: Southampton, UK, 2008.

740 14. Muggleton, J.M.; Brennan, M.J.; Rogers, C.D.F. Point vibration measurements for the detection of shallow-buried objects. Tunn. Undergr. Space Technol. 2014, 39, 27–33, DOI. 10.1016/j.tust.2012.02.006

741 15. Katz, A.; Karaa, F.; Niver, E. Innovative and Effective Techniques for Locating Underground Conduits; New Jersey Institute of Technology: Newark, NJ, USA, 2011.

745 16. McFadden, J.A. Radial Vibrations of Thick-Walled Hollow Cylinders. *Journal of the Acoustical Society of America* **1954**, *26*, 714, DOI. 10.1121/1.1907405

747 17. Gazis, D.C. Exact Analysis of the Plane-Strain Vibrations of Thick-Walled Hollow Cylinders. *Journal of the Acoustical Society of America* **1958**, *30*, 786–794, DOI. 10.1121/1.1909761

749 18. Sinha, B.K.; Plona, T.J.; Kostek, S.; Chang, S.K. Axisymmetric wave propagation in fluid-loaded cylindrical shells. I: Theory. *Journal of the Acoustical Society of America* **1992**, *92*, 1132–1143, DOI. 10.1121/1.404040

751 19. Operto, S.; Virieux, J.; Ribodetti, A.; Anderson, J.E. Finite-difference frequency-domain modeling of viscoacoustic wave propagation in 2D tilted transversely isotropic media. *Geophysics* **2009**, *74*, T75-T95, DOI. 10.1190/1.3157243

753 20. Virieux, J. P-SV wave propagation in heterogeneous media: Velocity-stress finite-difference method. *Geophysics* **1986**, *51*, 889-901, DOI. 10.1190/1.1442147

755 21. Virieux, J. SH wave propagation in heterogeneous media: Velocity-stress finite-difference method. *Geophysics* **1984**, *49*, 1933-1957, DOI. 10.1190/1.1441605

757 22. Gardner, G. H. F.; Gardner, L. W.; Gregory, A. R. Formation velocity and density - The diagnostic basics for stratigraphic traps. *Geophysics* **1974**, *39*, 770-780, DOI. 10.1190/1.1440465

759 23. Forgues, E.; Lambaré, G. Parameterization study for acoustic and elastic ray+Born inversion. *Journal of Seismic Exploration* **1997**, *6*, 253-278

761 24. Virieux, J.; Operto, S. An overview of full-waveform inversion in exploration geophysics. *Geophysics* **2009**, *74*, 127-152, DOI. 10.1190/1.3238367

763 25. Pirli, M. Heiner Igel: Computational seismology: a practical introduction. *Journal of Seismology* **2017**, *21*, 567-570, DOI. 10.1007/s10950-017-9662-4

765 26. Berenger, J. P. A perfectly matched layer for the absorption of electromagnetic waves. *Journal of Computational Physics* **1994**, *114*, 185-200, DOI. 10.1006/jcph.1994.1159

767 27. Stoica, P.; Nehorai, A. MUSIC, maximum likelihood, and Cramer-Rao bound. *IEEE Transactions on Acoustics, Speech, and Signal Processing* **1989**, *37*, 720-741, DOI. 10.1109/29.17564.

769 28. Scharf, L. L.; McWhorter, L. T. Geometry of the Cramer-Rao bound. *Signal Processing* **1993**, *31*, 301-311, DOI. 10.1016/0165-1684(93)90088-R

771 29. Schmidt, R.O. Multiple Emitter Location and Signal Parameter Estimation. *IEEE Transactions on Antennas and Propagation* **1986**, *34*, 276-280, DOI. 10.1109/TAP.1986.1143830.

773 30. Wang, B.; Zhao, Y.; Liu, J. Mixed-Order MUSIC Algorithm for Localization of Far-Field and Near-Field Sources. *IEEE Signal processing letters* **2013**, *20*, 311-314, DOI. 10.1109/LSP.2013.2245503.

775 31. Fossati, C.; Bourennane, S.; Marot, J. Chapter 3: Localization of Buried Objects in Sediment Using High Resolution Array Processing Methods. In *Underwater Acoustics*, Bourennane, S. (Ed.); InTech, Croatia, 2012; 41-58, DOI. 10.5772/33366

777 32. Han, D.; Fossati, C.; Bourennane, S.; Saidi, Z. Bearing and Range Estimation Algorithm for Buried Object in Underwater Acoustics. *Research Letters in Signal Processing* **2009**, DOI. 10.1155/2009/257564

779 33. Saidi, Z.; Bourennane, S. Cumulant-Based Coherent Signal Subspace Method for Bearing and Range Estimation. *EURASIP Journal on Advances in Signal Processing* **2007**, 084576 (2006), DOI. 10.1155/2007/84576

781 34. Sahin, A.; Miller, E. L. Object-based localization of buried objects using high-resolution array processing techniques. In proceedings of SPIE, 31 May 1996, DOI. 10.1117/12.241244

783 35. Choi, J. W.; Kim, Y. H. Spherical beam-forming and MUSIC methods for the estimation of location and strength of spherical sound sources. *Mechanical Systems and Signal Processing* **1995**, *9*, 569-588, DOI. 10.1006/mssp.1995.0043