<table><caption>Table 4.2: Average estimation errors and standard deviations on the lines parameters obtained during the last measurement step for the two scenarios in consideration. The errors are evaluated using 100 repetitions.</caption><thead><tr><td>Scenario</td><td>Range error [mm]</td><td>Angle error [degree]</td></tr></thead><tbody><tr><td>Scenario 1</td><td>3.007 ± 0.098</td><td>0.234 ± 0.0004</td></tr><tr><td>Scenario 2</td><td>10.766 ± 22.921</td><td>0.206 ± 0.134</td></tr></tbody></table>

over 100 runs are provided in Fig.4.9. On this plate, the echo detection employed by our previous method is not efficient, as it does not consider the wave dispersion effect, whereas the propagation distances are larger. This induces large misdetection rates and poor results. Hence, we display only the results of our new approach. Despite the slower convergence caused by the larger surface, and the slightly higher localization error, our method still provides precise estimates of the trajectory and plate geometry. This result indicates that our approach still works on surfaces sufficiently large to be used for realistic applications. The underlying prerequisites are a wave propagation model and filter parameters that conveniently fit the acoustic measurements and on-the-field noisy conditions. Naturally, one may also expect longer convergence times when the plate surface increases, as the echo detection is expected to be efficient mostly for short ranges as shown in Fig.4.9.a).

As a final evaluation, we determine the average mapping errors and standard deviations over 100 runs obtained during the final step for a lawn-mower path (Scenario 1) and a random walk (Scenario 2) on plate 1. Fig. 4.2 presents the results. It can be noticed that the overall results are relatively poorer for the random walk. This illustrates that the estimation accuracy also strongly depends on the robot path which shall be optimized for optimal reconstruction.

### 4.4.3 Experiment with a magnetic crawler

We evaluate our approach using a real robotic system acquiring ultrasonic data on a large plate to demonstrate its feasibility for a practical robotic inspection task. In our experimental setup, a magnetic crawler is moving on the surface of the 1700x1000x6mm steel plate, in an indoor laboratory environment. It is equipped with an onboard accelerometer to provide accurate heading measurements, and with a single transducer that is mounted in its head for simultaneous acoustic emission and reception.

This true pulse-echo setup (contrary to the pseudo-pulse echo setup with a pair of emitter/receiver used in the previous experiments) requires using an electric circuit for protecting the acquisition device during the high-voltage excitation. The prototype circuit designed to achieve that purpose is provided in Fig. 4.10. The principle is to protect the acquisition device during the excitation at high amplitude, otherwise, the acquisition device could be damaged. This is achieved by the use of Zener diodes which are placed in parallel with the acquisition device, as they can saturate the voltage level at their ends. The remaining of the excitation power is dissipated by the resistor. In this circuit, the commutation diodes placed in parallel, and in series with the generator, act as an interrupter: when the generator is no longer providing an excitation (i.e when the voltage V and intensity go to zero), it is "disconnected" by the diodes from the rest of the circuit to avoid short circuit of the echo measurements.

Despite the fact that echo signals should be preserved by the circuit, this setup leads to importance disturbance of the signals. Thus, care is taken in designing a sigmoid filter window and a band pass filter to obtain signals that can be used by our FastSLAM algorithm.
<figure><img src="image_36.png" /><figcaption>Figure 4.10: Schematic of designed electrical circuit for a true pulse-echo setup. The diodes in parallel are 1N4148 commutation diodes, while the diodes in series are Zener diodes.</figcaption></figure>

To acquire the measurements, the robot is driven along a trajectory on the plate positioned nearly vertically, and the acquisition of one measurement is triggered each time the sensor has traveled 5 cm. Also, a pump is activated to continuously bring water through the tether at the interface between the sensor and the plate for sufficient coupling. Due to variations of the contact condition between the transducer and the plate surface, the varying coupling level maintained by the water supply and the electric circuit, the signal-to-noise ratio in the measurements is relatively poorer than when the acquisition is carefully performed by hand. We will see that, in these conditions, localization and mapping can still be achieved.

In Fig. 4.11, we show qualitative results obtained with our FastSLAM method during different steps of the considered scenario. An accurate localization system such as a laser total station could be used to precisely assess the localization error by tracking the position of the sensor. However, such a system was not available at the time of the experiment. Also, snapshots of the camera recording of the experiment are provided for each corresponding step. It can be seen on the figure that, during Step 1 which is the initialization step, a square estimate of the plate geometry is returned, as only a single measurement has been integrated. During Step 2, the robot has moved to the bottom part of the plate. Yet, the mapping error appears seemingly important, and the uncertainty on the robot position an orientation is large. The reason is that the rotational noise is probably overestimated in the transition model, and the map estimate, which is still erroneous, cannot be used for correction. This issue starts to resolve itself during Step 3, where the position of the particles in the reconstructed plate appears coherent with the reality. This remains the case until the final step. Overall, this illustrates the feasibility of using ultrasonic echoes in practice to achieve localization and mapping on a metal plate with a real robotic platform, and based on the proposed FastSLAM approach.
<figure><img src="image_37.png" /><figcaption>Step 1: Camera snapshot</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>Step 1: SLAM results</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>Step 2: Camera snapshot</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>Step 2: SLAM results</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>Step 3: Camera snapshot</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>Step 3: SLAM results</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>Step 4: Camera snapshot</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>Step 4: SLAM results</figcaption></figure>

Figure 4.11: Qualitative SLAM results obtained with a magnetic crawler, and using the proposed FastSLAM approach, during different steps of the experiment.
4.5 Chapter summary

We have designed a new FastSLAM approach to achieve Simultaneous Localization and Mapping on metal plates by relying on the reflections of an ultrasonic excitation on the plate boundaries. Contrary to the proof-of-concept approach, this new method relies on a wave propagation model from the Rayleigh-Lamb theory and beamforming maps for estimating the plate edges location. Also, this method has the particularity of not requiring explicit echo detection and association, which are challenging problems, especially when considering ultrasonic data.

Experiments carried out on an undamaged and a damaged metal plate in a laboratory envi-
ronment demonstrate that this new approach achieves better results in terms of localization and
mapping accuracies and robustness, with less algorithmic complexity. More specifically, it was
demonstrated that (sub-) centimeter precision can be achieved in localization and mapping in
moderate odometry noise conditions. Yet, the proposed approach is restricted to recovering rect-
angular plate geometries. Finally, the approach has been validated using a real robotic platform,
in a laboratory environment. Despite the low quality of the signals due to noise, the qualitative
results illustrate the feasibility of ultrasonic-based on-plate SLAM with the proposed approach.
# Chapter 5

## Guided wave-based mapping with combined representations

The ability to reliably map acoustic scatterers using ultrasonic echoes is key to successful long-range robotic inspection. By relying on beamforming to locate the boundaries of a metal plate, the feasibility of on-plate ultrasonic SLAM has been demonstrated in the previous chapter. However, such an approach has some limitations: although it is relatively simple, this approach does not tackle a well-known issue of DAS beamforming, which is signal interference. Furthermore, the resulting point-estimate map lacks of uncertainty assessment. As a result, the confidence one can put in the map estimate may be unclear, especially in a realistic and noisy environment.

In this chapter, we introduce a more elaborate mapping framework based on combined map representations. More specifically, we seek to recover two different (grid and feature-based) maps of the same metal plate environment. We rely on the dense and probabilistic grid map for uncertainty assessment, and subsequently leverage the property of consistency between the two representations to mitigate interference. As the proposed framework relies on the detection of the closest edge reflector, probabilistic models are introduced to ensure robustness with respect to echo detection error. The approach is extensively evaluated in simulation against a more standard mapping method, and it is experimentally validated in a laboratory environment.

In the second part of this chapter, the designed mapping framework is integrated within a frontier exploration strategy with the aim to automatically drive an autonomous robot towards informative measurement locations. We determine that such a strategy may result in faster plate reconstruction and coverage compared to standard exploration strategies.

Last, as the reconstructed plate geometries were, so far, limited to rectangular shapes, we see how a partially reconstructed grid map can be used to extend feature-based plate boundary mapping to recover more arbitrary polygonal shapes from the ultrasonic echoes. The approach is extensively evaluated in simulation and in a real-world setting.

### 5.1 Motivation and background

The experimental results from the previous chapter have shown that delay-and-sum beamforming can be efficient for mapping a rectangular metal plate using UGWs. This method is relatively simple to implement, and it can seamlessly run online, within a FastSLAM framework. However, these benefits do not come without disadvantages and limitations. First, DAS beamforming is a non-adaptive method (i.e the way how different measurements are fused together for mapping does not depend on the data itself), and it does not tackle the detrimental effect of signal
interference (that will be further explained in the next section), whereas it is a well-known issue of DAS beamforming. Second, beamforming only provides point map estimates. As a result, due to the lack of uncertainty assessment and the time-inconsistent map estimates provided in the early steps, the approach is not suitable for mapping several metal plates on a large structure, which may require robust detection of plate boundary-crossing. Indeed, detection failures may result in mapping inaccuracies, as the present DAS approach requires that the measurement positions lie on the same metal panel. Besides, the lack of uncertainty assessment makes the approach unsuitable for active sensing, whereas it has been shown that the robot trajectory affects the estimation results. Last, this mapping method can only recover rectangular geometries. Although it may not be an issue for mapping metal panels on ship hulls and storage tanks, the extension to arbitrary shapes can be beneficial for application to more general structures, or even broader, to a wider range of acoustic SLAM problems.

We will see in this chapter that, surprisingly, the use of combined map representations may help tackle the three aforementioned limitations. We will aim at recovering grid-based and feature-based spatial representations jointly, while still limiting our case-study to mapping metal plate boundaries. In this context, the feature-based map will be represented by a set of two-dimensional lines, each line accounting for a potential boundary. The grid will consist of a dense representation of the plate environment that can account, all along the mapping process, both for areas identified as inside or outside the current metal plate, and areas whose state is still unknown, making it suitable e.g. for detecting a change of plate, or for use in a later active-sensing strategy. The construction of the grid is based on the Occupancy Grid (OG) framework, widely used in robotics, that will be briefly recalled.

Maintaining two different representations of map estimates brings its own challenges. As they should represent the same environment, they *must* always remain consistent in some sense. Yet, this consistency property shall be leveraged: on the one hand, the grid map is more spatially informative than the feature-based map, which can help to filter out candidate plate edges that are no longer relevant, and mitigate the detrimental effect of interference. On the other hand, line features estimates, when associated with individual echoes, can be useful to estimate areas outside the plate for more robust boundary detection. In this chapter, we will attempt to establish a mathematical framework for building such joint maps. It is tested in simulation, with measurements acquired experimentally by hand, and with a magnetic crawler.

Besides, the past experiments have highlighted the strong dependency of the reconstruction results with the robot's trajectory. Hence, a legitimate question would be: how the robot should navigate to maximize the information gain on the map? Can it be achieved algorithmically? Based on the availability of the probabilistic grid map, we will see that a frontier exploration strategy appears to be a reasonable option, as the robot is steered towards unexplored areas while aiming for a maximum mapping coverage. Also, one may easily comprehend that such a strategy is relevant for robotic inspection as, after operation, all the structure surface must be inspected (or "explored").

Finally, to alleviate the need of the rectangular geometry assumption, we will see that leveraging a consistency constraint between the line feature-based map with a (partially reconstructed) grid map can help to recover arbitrary polygonal shapes, without prior knowledge on the number of boundaries to estimate. Although its practical benefit for our target application may be limited, this new approach may provide an elegant solution to the challenging geometry reconstruction problem using omnidirectional acoustic waves. Indeed, through minimization of a single cost function, an arbitrary shape could be recovered without any explicit echo detection and echo association, resulting in little algorithmic complexity and sensitivity to noise.
## 5.1.1 Related works

So far, the mapping methods presented earlier are exclusively feature-based, as they can recover metal plate boundaries as line features using DAS beamforming. Yet, grid-based map representations are widely used for robotic applications, especially for exploration [94, 33, 102], as, contrary to feature-based maps, they can explicitly account for explored and unexplored areas.

The estimation of grid maps using UGWs and a moving sensor has been investigated in related works such as [90, 25, 91]. However, the proposed methods rely on directional GWs, which necessitates constant alignment of the sensor directivity with the vector perpendicular to the boundary. In other words, these works rely on substantial prior knowledge on the map for robot navigation. In [62], omnidirectional GWs are used to autonomously map, in a pitch-catch setup via a frontier exploration strategy, the inside of a metal panel by relying on a grid representation. However, the multiple echoes within the signals are not exploited in the mapping procedure, and the presence of echo measurement errors is not considered.

In the robotics literature, several works have demonstrated the potential of combining grid and feature-based maps using conventional directional sensors (e.g. camera, lidar, sonar) [5, 68, 101]. One essential motivation behind this approach (along with the cumulation of both map sparsity and denseness) is that sparse map features can be used to integrate dependency between the measurements to build the grid map, whereas they are almost systematically assumed independent in the conventional OG framework. In the proposed framework, the consistency between the two maps is leveraged to mitigate the detrimental effect of interference. In addition to the use of DAS beamforming, our method is built on the detection of the single echo due to reflection on the closest edge. While the corresponding detection errors can be a source of mapping inaccuracies, our probabilistic formulation enables our approach to be robust to such. In the literature, however, such detection issues are often overlooked by recent works addressing the similar problem of room shape reconstruction from acoustic echoes in the air [51, 69].

## 5.1.2 Binary Bayesian filtering

The binary Bayesian filter is conventional for estimating occupancy grid maps. In this representation, each grid element, which represents a subdivision of the space, is a binary random variable. Either the corresponding space is occupied by some object or obstacle (hence it is attributed a value of 1), or it is empty (and has a value 0). The elements of the discretized ground truth grid (referred to as $O(\mathbf{x})$) are either 0 or 1. A grid estimate $\hat{O}(\mathbf{x})$, however, contains the probability that each grid element is occupied. Hence, each element can take any value between 0 and 1, while a value of 0.5 means that the state is truly uncertain. Fig. 5.1 shows an example of an occupancy grid map with coarse resolution for illustration.

The core of the binary Bayesian filter for recovering a grid map estimate from external observations is the conditional independence assumption between all grid elements $o$ of $O$, and between all observations $\mathbf{z}_{0:t}$. This enables a recursive implementation of the filter for updating the state belief on each grid element: $b(o) = p(o = 1|\mathbf{x}_{0:t}, \mathbf{z}_{0:t})$, that will be written $b(o) = p(o|\mathbf{x}_{0:t}, \mathbf{z}_{0:t})$ for notation simplicity. Furthermore, the conventional filter updates log-odds ratios:

$$l_t = \log \frac{p(o|\mathbf{x}_{0:t}, \mathbf{z}_{0:t})}{1 - p(o|\mathbf{x}_{0:t}, \mathbf{z}_{0:t})},$$

instead of the belief probability, as it leads to numerically stable updates using simple additions, as shown below. One may note that $l_t$ can take any value between $-\infty$ and $+\infty$. Besides, the
<figure><img src="image_45.png" /><figcaption>Figure 5.1: An example of an occupancy grid with coarse resolution. It illustrates the typical effect of integrating a single range measurement to an obstacle into an initial uncertain map.</figcaption></figure>

occupancy probability can easily be determined from the latter with the formula:

$$b(o) = 1 - \frac{1}{1 + \exp l_t}.$$

Here, we recall how the recursive filtering is derived using inverse measurement models $p(o|\mathbf{z})$, as presented in [94]. To do so, we start from the Bayes’ formula:

$$\begin{aligned} p(o|\mathbf{x}_{0:t}, \mathbf{z}_{0:t}) &= \frac{p(\mathbf{z}_t|o, \mathbf{x}_{0:t}, \mathbf{z}_{0:t-1})p(o|\mathbf{x}_{0:t}, \mathbf{z}_{0:t-1})}{p(\mathbf{z}_t|\mathbf{x}_{0:t}, \mathbf{z}_{0:t-1})} \\ &= \frac{p(\mathbf{z}_t|o)p(o|\mathbf{x}_{0:t-1}, \mathbf{z}_{0:t-1})}{p(\mathbf{z}_t|\mathbf{x}_{0:t}, \mathbf{z}_{0:t-1})}, \end{aligned}$$

where the assumption of conditional independence between the observations is used. Now, after applying the Bayes rule to the measurement model:

$$p(\mathbf{z}_t|o) = \frac{p(o|\mathbf{z}_t)p(\mathbf{z}_t)}{p(o)},$$

we obtain:

$$p(o|\mathbf{x}_{0:t}, \mathbf{z}_{0:t}) = \frac{p(o|\mathbf{z}_t)p(\mathbf{z}_t)p(o|\mathbf{x}_{0:t-1}, \mathbf{z}_{0:t-1})}{p(o)p(\mathbf{z}_t|\mathbf{x}_{0:t}, \mathbf{z}_{0:t-1})}.$$

By applying the same reasoning for the probability of non-occupancy, one can obtain:

$$p(\neg o|\mathbf{x}_{0:t}, \mathbf{z}_{0:t}) = \frac{p(\neg o|\mathbf{z}_t)p(\mathbf{z}_t)p(\neg o|\mathbf{x}_{0:t-1}, \mathbf{z}_{0:t-1})}{p(\neg o)p(\mathbf{z}_t|\mathbf{x}_{0:t}, \mathbf{z}_{0:t-1})}.$$

Dividing the two previous equations leads to cancellation of probabilities that are difficult to calculate:

$$\begin{aligned} \frac{p(o|\mathbf{x}_{0:t}, \mathbf{z}_{0:t})}{p(\neg o|\mathbf{x}_{0:t}, \mathbf{z}_{0:t})} &= \frac{p(o|\mathbf{z}_t)p(o|\mathbf{x}_{0:t-1}, \mathbf{z}_{0:t-1})p(\neg o)}{p(\neg o|\mathbf{z}_t)p(\neg o|\mathbf{x}_{0:t-1}, \mathbf{z}_{0:t-1})p(o)}, \\ &= \frac{p(o|\mathbf{z}_t)}{1 - p(o|\mathbf{z}_t)} \frac{p(o|\mathbf{x}_{0:t-1}, \mathbf{z}_{0:t-1})}{1 - p(o|\mathbf{x}_{0:t-1}, \mathbf{z}_{0:t-1})} \frac{1 - p(o)}{p(o)}. \end{aligned}$$
After taking the logarithm of this expression, we obtain the recursive update formula of the Bayesian binary filter:

$$l_t = \log \frac{p(o|\mathbf{z}_t)}{1 - p(o|\mathbf{z}_t)} - \log \frac{p(o)}{1 - p(o)} + l_{t-1}.$$

## 5.2 Grid and feature-based mapping using ultrasonic guided waves

Our approach for mapping with combined representations builds on feature-based maps estimated via beamforming to construct inverse sensor models in an online fashion, so that a grid representation of a metal plate can be recovered within the OG framework. Simultaneously, the beamforming maps are filtered with binary masks to maintain consistency with the grid map, and at the same time, to mitigate the detrimental effect of interference.

### 5.2.1 Feature-based mapping via beamforming

As a reminder, let's briefly recall the feature-based mapping approach before highlighting some of its limitations. When a set of measurements have been acquired along a trajectory {$x_i, y_i$}$_{i=0..t}$, DAS beamforming is applied to recover the geometry of a metal plate from envelope signals $e_i$ obtained through correlation of the measurement $z_i$ with a wave propagation model $\hat{g}$. This yields 2D beamforming maps that are computed with:

$$\mathcal{L}_t(r, \alpha) = \sum_{i=1}^{t} e_i(d_i(r, \alpha)) = \sum_{i=0}^{t} e_i(|x_i \cdot \cos \alpha + y_i \cdot \sin \alpha - r|), \quad (5.1)$$

where $(r, \alpha)$ represent the polar coordinates of a candidate line boundary with equation: $x \cdot \cos \alpha + y \cdot \sin \alpha - r = 0$ expressed in a fixed frame with origin $O$, while $\mathcal{L}_t(r, \alpha)$ indicates the likelihood of the line existence. Next, a feature-based estimate of the map $\hat{\mathbf{M}}_t = \{r_l, \alpha_l\}_{l=1..4}$ is recovered by extracting the local maxima of $\mathcal{L}_t$ through the resolution of the optimization problem:

$$\hat{\mathbf{M}}_t = \arg \max_{\mathbf{M}} \mathcal{L}_t(\mathbf{M}) = \arg \max_{\mathbf{M}} \sum_{l=1}^{4} \mathcal{L}_t(r_l, \alpha_l),$$

where $\mathbf{M}$ is restricted to be a rectangle. Note that an efficient resolution approach is proposed in the previous chapter by taking this constraint into account. Despite its efficiency in high signal-to-noise (SNR) conditions, the geometry retrieval can often be ambiguous due to interference, as an infinity of lines can equally account for one range measurement. Put in other words, assuming a high correlation $e_i(r_0)$, the application of DAS beamforming results in all lines $(r, \alpha)$ satisfying $|x_i \cdot \cos \alpha + y_i \cdot \sin \alpha - r| = r_0$ to have increased likelihood $\mathcal{L}_t(r, \alpha)$, whereas it is likely that the echo that explains the high value of $e_i(r_0)$ originates from reflection on a single line. In addition, the effect of interference is further exacerbated by the presence of multiple echoes and high-order reflections that are not considered in Eq. (5.1). Fig.5.2a illustrates the interference effect in a simulated scenario, where a mobile unit is driven to acquire acoustic measurements at predetermined positions marked as red points. As most of the measurements are acquired on a line parallel to an edge, there is an ambiguity on which side the real edge is. This results in the boundary on the right not being correctly identified from the beamforming map shown in Fig.5.2b.

The intuition behind the developed approach, which uses combined map representations, is that, when relying on detecting the earliest echo within the measurements, which indicates the distance to the closest edge, some information on the spatial environment can be directly
<figure><img src="image_46.png" /><figcaption>Figure 5.2: Illustration of the interference effect. (a) shows mapping results obtained in a simulated scenario, where one edge is not correctly identified from the beamforming map shown in (b). The origin O of the reference frame is indicated along with the (r, α) coordinate of the left line. The ground truth plate outline is represented by the blue rectangle.</figcaption></figure>

integrated into the mapping results. First, it can be inferred that the disk centered at the sensor position whose radius equals the detected range contains only points on the current plate. Second, the lines that cross the total area identified as “on the plate” can be filtered out on the beamforming map to mitigate interference. This principle is illustrated in Fig. 5.3, where all the edges are now correctly identified after applying the mask filter to the beamforming map. The main issue of this approach is that the actual lines may be filtered as well when the measured ranges are overestimated.

<figure><img src="image_47.png" /><figcaption>Figure 5.3: Illustration of beamforming map filtering with the first acoustic echo for interference rejection. a) shows the mapping results obtained after applying, to the beamforming map in Fig. 5.2b, a binary mask that filters lines which can no longer exist, as they cross the area identified as "on the plate". The on-plate area identified with the ground-truth echo measurements is represented in (a) as the lighter area. The origin O of the reference frame is indicated along with the (r, α) coordinate of the left line.</figcaption></figure>

### 5.2.2 Estimation problem formulation

To obtain a more informative and time-consistent representation of the environment compared to a feature-based mapping only, our approach is designed to recover a grid representation of the 2D space,

$$O(\mathbf{x}) = \begin{cases} 0 & \text{if } \mathbf{x} = [x, y] \text{ is a point of the plate} \\ 1 & \text{otherwise} \end{cases},$$
by constructing an estimate $\hat{O}_t$ with posteriors $\hat{O}_t(\mathbf{x}) = p(O(\mathbf{x}) = 1 | \mathbf{x}_{0:t}, z_{0:t})$ based on the detection of the closest edge. However, as the single-use of the detected range cannot be sufficient to recover areas outside the plate, the probability model shall rely on the feature-based map determined via beamforming to estimate the direction of arrival (DoA) of the detected echo. Besides, the introduced representation can be leveraged to filter out lines on $\mathcal{L}_t(r, \alpha)$ that are no longer relevant, because they cross the area identified as on the plate with $\hat{O}_t$, with the aim to mitigate interference. Hence, we formally introduce the binary mask filter associated with the ground-truth grid $O(\mathbf{x})$,

$$Q(r, \alpha) = \begin{cases} 0 & \text{if } \exists \mathbf{x} \text{ of the line } (r, \alpha) \text{ s.t. } O(\mathbf{x}) = 0 \\ 1 & \text{otherwise} \end{cases}, \qquad (5.2)$$

so that the feature-based map $\hat{\mathbf{M}}_t$ can be retrieved from the masked beamforming map, where the mask is to be estimated as well: $\hat{\mathcal{L}}_t(r, \alpha) = \hat{Q}_t(r, \alpha) \mathcal{L}_t(r, \alpha)$. It is critical that such filtering is based on a probabilistic approach to avoid erroneously filtering lines. In practice, echo misdetection may occur due to noise and the propagation model's inaccuracy.

### 5.2.3 Binary filtering with adaptive inverse sensor models

The standard OG framework presented earlier is convenient for estimating grids of binary values. It relies on the essential assumptions that the measurements are conditionally independent given the grid state and that the states of the individual grid cells are conditionally independent given the measurements [94]. These assumptions, which are not always realistic, yield a computationally efficient algorithm that can run in real-time. Here, the filter involves an inverse measurement model that will be estimated from all past acoustic data based on the feature-based map retrieved via beamforming: $p(O(\mathbf{x})|\mathbf{x}_t, e_t, \hat{\mathbf{M}}_t)$. More specifically, as $\hat{\mathbf{M}}_t$, which represents the parameters of the inverse model, is updated online, the interpretation of a measurement $e_t$ can vary as well based on DoA estimates. In contrast, the standard OG framework would only involve $p(O(\mathbf{x})|\mathbf{x}_t, e_t)$. Eventually, the grid can be recursively estimated with the standard Bayesian binary filter with log-odds ratios: $l_t = \log \frac{\hat{O}_t(\mathbf{x})}{1-\hat{O}_t(\mathbf{x})}$ (and the assumption of a uniform prior):

$$l_t = l_{t-1} + \log \frac{p(O(\mathbf{x})|\mathbf{x}_t, e_t, \hat{\mathbf{M}}_t)}{1 - p(O(\mathbf{x})|\mathbf{x}_t, e_t, \hat{\mathbf{M}}_t)}, \qquad (5.3)$$

with initial cell values $\forall \mathbf{x}, \hat{O}_0(\mathbf{x}) = 1/2$. For computational efficiency, the same approach is used to build $\hat{Q}_t$ with a second inverse model $p(\hat{Q}_t(r, \alpha)|\mathbf{x}_t, e_t, \hat{\mathbf{M}}_t)$. Overall, the major challenge is specifying the two inverse models to determine the mapping strategy fully. In addition, this strategy shall be designed to maintain sufficient consistency between the grid $\hat{O}_t$ and the mask $\hat{Q}_t$, as $\hat{Q}_t(r, \alpha)$ should formally depend on the full grid $\hat{O}_t$ as in Eq. (5.2). However, in our approach, it is computed independently of the later given measurements and beamforming results.

One may observe that this approach is not perfectly rigorous mathematically speaking. Indeed, the use of Eq. (5.3) is based on the assumption of conditional independence between measurements (see how the binary filter formula is derived in the previous section 5.1), whereas $\hat{\mathbf{M}}_t$ is used as a conditional variable to reintroduce dependence with past observations. Yet, this mathematically inexact trick enables to leverage the simplicity of the binary filter with log-odd ratios for our problem. Another remark is that, here, we only consider the mapping problem, as the trajectory $\mathbf{x}_{0:t}$ is assumed known. Yet, our new mapping framework may seamlessly be integrated within a SLAM framework such as FastSLAM, as in Chapter 4.
### 5.2.4 Inverse model estimation in the Cartesian space

We first start the construction of a mixture model $p(O(\mathbf{x})|\mathbf{x}_t, e_t, \hat{\mathbf{M}}_t)$ through the modeling of imperfect echo measurements. Let $\rho_t$ be the range to the closest edge retrieved from $e_t$. We assume that the echo detection can completely fail with a probability $\gamma$ to construct a model robust to measurement error. In that event, the cell distribution shall be updated with a uniform distribution. With a probability $1 - \gamma$, the range is detected with a precision consistent with the wavelength. This is accounted for by introducing the model $p(\rho|e_t)$ that we model with a truncated Gaussian distribution $p(\rho|e_t) = N_{[0,d_{\max}]}(\rho|\rho_t, \sigma^2)$ of support $[0, d_{\max}]$,

$$p(\rho|e_t) = \begin{cases} \frac{1}{\eta\sqrt{2\pi}\sigma^2} \exp\left\{-\frac{(\rho-\rho_t)^2}{2\sigma^2}\right\}, & \text{if } \rho \in [0, d_{\max}] \\ 0 & \text{otherwise} \end{cases},$$

where $\eta$ is the normalizing factor and $d_{\max}$ is the maximum detectable range that is introduced due to the finite observation window of the signals. The variable $\rho$ denotes possible values for the range. Hence, the inverse model can first be expressed as a mixture and using a convolution:

$$p(O(\mathbf{x})|\mathbf{x}_t, e_t, \hat{\mathbf{M}}_t) = \frac{\gamma}{2} + (1 - \gamma) \int p(O(\mathbf{x})|\mathbf{x}_t, \hat{\mathbf{M}}_t, \rho)p(\rho|e_t) d\rho. \quad (5.4)$$

To build the model $p(O(\mathbf{x})|\mathbf{x}_t, \hat{\mathbf{M}}_t, \rho)$ which would state how the probability is updated for a given detected range $\rho$, we consider two scenarios. In the first one, the DoAs provided by $\hat{\mathbf{M}}_t$ would not be taken into account because the detected range would not be consistent with the expected range: i.e the map $\hat{\mathbf{M}}_t$ is assumed to be erroneous as it does not contain any line that is consistent with the range measurement, given that echo misdetection is already modeled with the probability $\gamma$. To account for this scenario, we introduce a model: $p_1(O(\mathbf{x})|\mathbf{x}_t, \rho) = h(|\mathbf{x}_t - \mathbf{x}| - \rho)/2$ where $h$ denotes the Heaviside's function, $h = 1_{\rho>0}$. Hence, only the points closer than the presumed range $\rho$ can be marked as "on the plate", while those further away would have unchanged state belief distribution. Let us first assume that $p(O(\mathbf{x})|\mathbf{x}_t, \hat{\mathbf{M}}_t, \rho)$ is represented only by this model. As it is convolved in Eq. (5.4) with the truncated Gaussian distribution, the distribution updates are weighted with the range likelihoods. The convolution can seamlessly be expressed with the error function:

$$\int h(|\mathbf{x}_t - \mathbf{x}| - \rho) N_{[0,d_{\max}]}(\rho|\rho_t, \sigma^2) d\rho = \frac{\Phi_t(|\mathbf{x}_t - \mathbf{x}|) - \Phi_t(0)}{\Phi_t(d_{\max}) - \Phi_t(0)}; \Phi_t(\rho) = \erf\left(\frac{\rho - \rho_t}{\sigma\sqrt{2}}\right).$$

The mathematical derivation can be found in the appendix at the end of this thesis. With the current model, however, only points that are believed to be on the current plate are being recovered, and there is no mechanism for estimating areas outside the plate. To remedy that, we introduce a second model that will be used to integrate the information from $\hat{\mathbf{M}}_t$. Let's consider that the retrieval of a range $\rho$ might originate from a reflection of the incident wave on an edge with polar coordinates $(\hat{r}, \hat{\alpha})$ from $\hat{\mathbf{M}}_t$. We introduce a second model to update the cells in the specific direction $\hat{\alpha}$ and beyond the range $\rho$,

$$p_2(O(\mathbf{x})|\mathbf{x}_t, (\hat{r}, \hat{\alpha}), \rho) = \frac{1}{2}\beta(\hat{r}, \hat{\alpha}, \rho, \mathbf{x}_t) \times h((x - x_t) \cos \hat{\alpha} + (y - y_t) \sin \hat{\alpha} - \rho),$$

so that the cell "out-of-plate" probability can also increase. As the range-line association is not necessarily correct, we have introduced an association coefficient $\beta$ between 0 and 1 that will make the probability update effective only when the considered range $\rho$ matches with that expected
from a reflection on the edge ($\hat{r}, \hat{\alpha}$), and with Gaussian distribution to model uncertainty on the line parameters:

$$ \beta(\hat{r}, \hat{\alpha}, \rho, \mathbf{x}_t) = \exp \left\{ -\frac{1}{2b^2} \left( |x_t \cdot \cos \hat{\alpha} + y_t \cdot \sin \hat{\alpha} - \hat{r}| - \rho \right)^2 \right\}. $$

For simplicity, the two aforementioned sub-models are then assembled together through averaging, yielding for a single line of the feature-based map:

$$ p(O(\mathbf{x})|\mathbf{x}_t, (\hat{r}, \hat{\alpha}), \rho) = \frac{p_1(O(\mathbf{x})|\mathbf{x}_t, \rho) + p_2(O(\mathbf{x})|\mathbf{x}_t, (\hat{r}, \hat{\alpha}), \rho)}{2}. $$

This will result in Eq. (5.4) in the convolution of the model $p_2$ with the truncated Gaussian which can also be expressed with the error function:

$$ \int \beta(\hat{r}, \hat{\alpha}, \rho, \mathbf{x}_t)h((x-x_t)\cos \hat{\alpha} + (y-y_t)\sin \hat{\alpha} - \rho) \times \mathcal{N}_{[0,d_{\max}]}(\rho|\rho_t, \sigma^2)d\rho \\ = \frac{\Gamma_{t,\hat{r},\hat{\alpha}}((x-x_t)\cos \hat{\alpha} + (y-y_t)\sin \hat{\alpha}) - \Gamma_{t,\hat{r},\hat{\alpha}}(0)}{\Phi(d_{\max}) - \Phi(0)} $$

where we have:

$$ \Gamma_{t,\hat{r},\hat{\alpha}}(\rho) = \frac{\lambda_t(\hat{r},\hat{\alpha}) \cdot \sigma_{\text{new}}}{\sigma} \times \text{erf}\left(\frac{\rho - \mu_{\text{new}}}{\sigma_{\text{new}}\sqrt{2}}\right), \quad \mu_{\text{new}} = \frac{b^2\rho_t + \sigma^2 d_t(\hat{r},\hat{\alpha})}{b^2 + \sigma^2} \\ \sigma_{\text{new}} = \frac{\sigma b}{\sqrt{b^2 + \sigma^2}} \text{ and } \lambda_t(\hat{r},\hat{\alpha}) = \exp\left(-\frac{(\mu_{\text{new}} - \rho_t)^2}{2\sigma^2} - \frac{(\mu_{\text{new}} - d_t(\hat{r},\hat{\alpha}))^2}{2b^2}\right). $$

How this expression is derived can be found in the appendix at the end of the thesis.

Yet, we have so far considered that only one line is available, whereas the full feature-based map $\hat{\mathbf{M}}_t$ possesses four lines. Hence, for updating the state of the cell $\hat{O}_t(\mathbf{x})$, we take the maximum value of the convolution terms over the lines, so that the probability that $\mathbf{x}$ is outside the plate can be increased as soon as there is at least one line estimate that is consistent with the range measurement. Furthermore, we introduce a multiplicative coefficient:

$$ C(\mathbf{x}, \mathbf{x}_t, \hat{r}_l, \hat{\alpha}_l) = \exp\{-\zeta|(y-y_t)\cos(\hat{\alpha}_l) - (x-x_t)\sin(\hat{\alpha}_l) - \hat{r}_l|\}, $$

so that the update can be made effective locally around the axis passing through $\mathbf{x}$ and its projection on the line ($\hat{r}_l, \hat{\alpha}_l$). All brought together, the full inverse model in the Cartesian space can be expressed with:

$$ p(O(\mathbf{x})|\mathbf{x}_t, e_t, \hat{\mathbf{M}}_t) = \frac{\gamma}{2} + \frac{1-\gamma}{4} \times \frac{\Phi_t(|\mathbf{x}_t - \mathbf{x}|) - \Phi_t(0)}{\Phi_t(d_{\max}) - \Phi_t(0)} \\ + \frac{1-\alpha}{4} \max_{(\hat{r}_l, \hat{\alpha}_l) \in \hat{\mathbf{M}}_t} \left\{ C(\mathbf{x}, \mathbf{x}_t, \hat{r}_l, \hat{\alpha}_l) \times \frac{\Gamma_{t, \hat{r}_l, \hat{\alpha}_l}((x-x_t)\cos \hat{\alpha}_l + (y-y_t)\sin \hat{\alpha}_l) - \Gamma_{t, \hat{r}_l, \hat{\alpha}_l}(0)}{\Phi_t(d_{\max}) - \Phi_t(0)} \right\}. $$

The first term of the sum accounts for echo misdetection. The second one updates on-plate points with the detected echo up to a certain precision. More specifically, given a detected range, the area situated at a lower distance will have a decreased probability of being outside the plate. The last term recovers areas outside the plate based on the estimated feature-based map, which provides DoA estimates: for every line ($\hat{r}, \hat{\alpha}$) of $\hat{\mathbf{M}}$, the points further away from the sensors than the detected range, and following a direction towards the line, have increased probability of being
<figure><img src="image_48.png" /><figcaption>Figure 5.4: Illustration of the inverse model in the Cartesian space for different values of detected range (which is equal to the radius of the light disk surface). The sensor position is represented by the red point, and the considered line estimate is in green.</figcaption></figure>

outside-the-plate when the range measurement matches the expected one. Besides, the update is made effective based the multiplicative coefficient $C(\mathbf{x}, \mathbf{x}_t, \hat{r}, \hat{\alpha})$. It restricts the updated grid cells to limit mapping errors when the orientation of the line estimate is not stabilized.

To illustrate the effect of the inverse model, we show, in Fig. 5.4, the grid $p(O|\mathbf{x}_t, e_t, (\hat{r}, \hat{\alpha}))$ for different values of detected range $\rho_t$, and for a fixed and single line estimate $(\hat{r}, \hat{\alpha})$. In Fig. 5.4.a, the detected range value matches the distance to the line estimate. Thus, the probability to lie outside the plate, assigned to the points in the direction towards the line, is substantially increased based on the coefficient $C(\mathbf{x}, \mathbf{x}_t, \hat{r}, \hat{\alpha})$. In Fig. 5.4.b, The probability increase is still noticeable, but it is weaker compared to the previous situation, as the match between the two distances is slightly worst. Last, when the difference between the measured and expected distances is too large as in Fig. 5.4.c, there are no cells that have increased probability of being outside the plate due to the lack of consistent DoA estimate.

### 5.2.5 Inverse model estimation in the line boundary space

As the mask $\hat{Q}_t$ is determined independently of the grid $\hat{O}_t$ with another binary filter, it is necessary to ensure sufficient consistency between the two grids, so that the filtered lines are consistent with the estimated on-plate area. We propose an inverse model in the line space which maintains exact agreement in range and approximate consistency in angle from the former inverse model:

$$p(Q(r, \alpha)|\mathbf{x}_t, e_t, \hat{\mathbf{M}}_t) = \frac{\gamma}{2} + \frac{1-\gamma}{4} \times \frac{\Phi_t(d_t(r, \alpha)) - \Phi_t(0)}{\Phi_t(d_{\max}) - \Phi_t(0)} \\ + \frac{1-\gamma}{4} \max_{(\hat{r}_l, \hat{\alpha}_l) \in \hat{\mathbf{M}}_t} \left\{ C'(\alpha, \hat{\alpha}_l) \times \frac{\Gamma_{t, \hat{r}_l, \hat{\alpha}_l}(|x_t \cdot \cos \hat{\alpha}_l + y_t \cdot \sin \hat{\alpha}_l - \hat{r}_l|) - \Gamma_{t, \hat{r}_l, \hat{\alpha}_l}(0)}{\Phi_t(d_{\max}) - \Phi_t(0)} \right\}$$

where $C'(\alpha, \hat{\alpha}_l) = \exp\{-\mu|\alpha - \hat{\alpha}_l|\}$. The exact consistency in range is indeed maintained, as any line in the $(r, \alpha)$ space is uniquely mapped to the orthogonal projection of the robot position on the line at each measurement step. Such an approach is convenient, as the overall algorithm is computationally efficient, and can run in real-time. Finally, we decide to filter out, on the beamforming maps, only the lines whose nonexistence is certain enough. Thus, the features $\hat{\mathbf{M}}_t$ are retrieved from $\hat{L}_t(r, \alpha) = \mathbf{1}_{\hat{Q}_t(r, \alpha)>\tau} \cdot L_t(r, \alpha)$ where $\tau$ is chosen below 0.5.
### 5.2.6 Algorithm implementation

To summarize, we provide in Algorithm 3, the overall pipeline for implementation of our mapping approach with combined map representations.

**Algorithm 3:** Mapping_with\_combined\_representations(x_t, e_t, O_{t-1}, Q_{t-1}, L_{t-1})

**Data:** Position at current time $x_t$, current envelope measurement $e_t$, grid map $O_{t-1}$, grid mask filter $Q_{t-1}$, beamforming map $L_{t-1}$.

**Result:** Grid map $\hat{O}_t$, mask filter $\hat{Q}_t$ and geometry estimate $\hat{M}_t$ at current time $t$.

// Updates the beamforming map and retrieves a feature-based map estimate ;

$\mathcal{L}_t(r, \alpha) = \mathcal{L}_{t-1}(r, \alpha) + e_t(|x_t \cos \alpha + y_t \sin \alpha - r|), \forall (r, \alpha);$

$\hat{\mathcal{L}}_t(r, \alpha) = \mathbf{1}_{\hat{Q}_{t-1}(r, \alpha) > \tau} \cdot \mathcal{L}_t(r, \alpha) \forall (r, \alpha);$

$\hat{\mathbf{M}}_t = \arg \max_{\mathbf{M}} \sum_{(r_l, \alpha_l) \in \mathbf{M}} \hat{\mathcal{L}}_t(r_l, \alpha_l);$

// Updates the grid map and grid filter with the new inverse models ;

Determines $\rho_t$ as the position of the first local maxima in $e_t$.

$$l_t(\mathbf{x}) = \log \frac{\hat{O}_{t-1}(\mathbf{x})}{1 - \hat{O}_{t-1}(\mathbf{x})} + \log \frac{p(O(\mathbf{x})|\mathbf{x}_t, e_t, \hat{\mathbf{M}}_t)}{1 - p(O(\mathbf{x})|\mathbf{x}_t, e_t, \hat{\mathbf{M}}_t)} \forall \mathbf{x};$$

$$\tilde{l}_t(r, \alpha) = \log \frac{\hat{Q}_{t-1}(r, \alpha)}{1 - \hat{Q}_{t-1}(r, \alpha)} + \log \frac{p(Q(r, \alpha)|\mathbf{x}_t, e_t, \hat{\mathbf{M}}_t)}{1 - p(Q(r, \alpha)|\mathbf{x}_t, e_t, \hat{\mathbf{M}}_t)} \forall (r, \alpha);$$

$$\hat{O}_t(\mathbf{x}) = 1 - \frac{1}{1 + \exp l_t(\mathbf{x})} \forall \mathbf{x}; \quad \hat{Q}_t(r, \alpha) = 1 - \frac{1}{1 + \exp l_t(r, \alpha)} \forall (r, \alpha);$$

return $\hat{O}_t, \hat{Q}_t, \hat{\mathbf{M}}_t$.

## 5.3 Numerical and experimental results

We assess the performance of our approach in different scenarios. In the first one, we use the ultrasonic measurements acquired on our aluminum plate of dimensions 600x450x6mm following the acquisition procedure presented in Chapter 3. As we only tackle the mapping problem here, only ground truth sensor positions will be used. In the second scenario, the data are obtained in simulation for plates of varying sizes and sensor trajectories. In the last scenario, we use the ultrasonic measurements acquired on our steel plate of dimensions 1700x1000x6mm, following the acquisition procedure presented in the previous chapter. For all the trials, the parameters of the filter are: $\gamma = 0.7$, $\sigma = 1$cm, $b = 2$cm, $\zeta = 3$, $\mu = 2$ and $\tau = 0.2$. The values of $\sigma$ and $b$ are consistent with the acoustic wavelength, approximately 2cm in our setups. The parameter $\gamma$ is set with a high value to limit the effect of one measurement update and increase the algorithm's robustness w.r.t. echo misdetection. The detection of the range to the closest edge is achieved with a standard peak detector applied to the correlation signals. The resulting performance of our approach is compared with that of a more standard mapping strategy [62] as a baseline for comparison. The latter is slightly adjusted to better match our case-study: as in [62], all the points within a given range measurement $\rho_t$ are directly marked as "on the plate". The lines, on the beamforming map, that cross the identified area are then filtered out, and a geometry estimate is then extracted from the processed beamforming result.

### 5.3.1 Results with hand-picked measurements

Fig. 5.5 shows the qualitative mapping results obtained with our approach at different steps of the scenario. The range to the closest edge is retrieved during the first step, but its orientation
<figure><img src="image_49.png" /><figcaption>Figure 5.5: Mapping results achieved on the aluminum plate during Steps 1, 42, 108 and 216. Each sub-figure depicts, on the left, the sensor trajectory along with the detected range to the closest edge represented by the circles, and on the right, the filtered beamforming maps during the same step, where the red rectangles indicate the retrieved edge estimates.</figcaption></figure>

is ambiguous, resulting in a squared plate estimate. During Step 42, an essential portion of the plate has been correctly identified as “on the plate”. Besides, the top and bottom sides of the plate have been marked as “outside the plate” based on the edge estimates retrieved from the filtered beamforming map. Nevertheless, the left edge has not been recovered because of interference that caused a line further away to be selected, as it cannot be filtered with our approach. During Step 108, all the edges are correctly retrieved from the filtered beamforming map. However, the area outside the plate on the left side has not yet been recovered because the edge estimation was incorrect when the mobile unit was close to it. This is achieved after simulating backwards the same trajectory, where the mapping is entirely achieved, as shown in Step 216. Ultimately, despite an average error of 4.3cm on range detection, the mean estimation errors on the line parameters are 0.5cm for the ranges $r_l$, and are inferior to one degree for the orientations $\alpha_l, l = 1...4$. Furthermore, the two map representations can be considered coherent with each other, as only 7% of the plate points have been incorrectly identified as outside it.

Hence, the results demonstrate the efficiency of our approach for mapping a plate structure with combined representations, where the grid can account for plate edges that are recovered at any time of the mapping process. However, the detected edges often lie at the limit of the filtered area. To assess the robustness of our approach, we evaluate the final average mapping errors over 10 repetitions of the lawnmower trajectory for varying noise levels on the echo detection. This is achieved by adding Gaussian noise to the ground truth range, determined based on the known acquisition positions. We do the same experience by using the baseline approach. Fig. 5.6 shows the results. It can be noticed that, with our approach, a range estimation error close to 1 cm and an orientation error of a few degrees at most are maintained for noise levels comparable to the wavelength. This is not the case with the baseline approach, where the estimation errors increase rapidly with the noise level, which illustrates its lack of robustness to map the inside of the plate with error-prone measurements.
<figure><img src="image_50.png" /><figcaption>Figure 5.6: Final mapping errors w.r.t. the noise level on the echo detection assessed over 10 repetitions of the lawnmower trajectory. The errors are determined both for our approach and for the baseline one. The left plot shows the range errors, while the orientation errors are shown on the right plot. The solid lines represent the average values, while the upper and lower bounds of the colored areas are situated respectively above and below one standard deviation to the mean values.</figcaption></figure>

## 5.3.2 Simulation results

To demonstrate that our method can generalize to various setups, we run it in simulation using data from three plates of different dimensions, and for a noise level of 5cm on echo detection. For comparison, we also repeat the process with the baseline approach, where all the points within the noisy range measurement are directly marked as "on the plate", and the lines, on the beamforming map, that cross the identified area are directly filtered out. We provide the resulting average errors and standard deviations assessed over 50 repetitions of random paths (at constant velocity) in Table 5.1. It can be observed that, for the three scenarios, our approach achieves a higher precision than the baseline approach.

Overall, the results demonstrate the effectiveness of our method for mapping a plate while filtering beamforming maps to reduce the effect of interference (without completely alleviating it), even in the presence of echo detection errors.

<table><caption>Table 5.1: Average estimation errors and standard deviations on the lines parameters obtained in simulation, using our approach (in **bold**) and the baseline method.</caption><thead><tr><th>Plate size [mm]</th><th>Model</th><th>Range error [cm]</th><th>Angle error [degree]</th></tr></thead><tbody><tr><td rowspan="2">600×450×6</td><td>Baseline</td><td>19.38 ± 4.83</td><td>33.46 ± 10.10</td></tr><tr><td>Ours</td><td><b>0.17 ± 0.09</b></td><td><b>0.44 ± 0.37</b></td></tr><tr><td rowspan="2">1000×600×6</td><td>Baseline</td><td>31.73 ± 14.72</td><td>13.62 ± 15.10</td></tr><tr><td>Ours</td><td><b>0.12 ± 0.11</b></td><td><b>0.37 ± 0.25</b></td></tr><tr><td rowspan="2">2000×1000×6</td><td>Baseline</td><td>69.62 ± 17.90</td><td>0.58 ± 0.89</td></tr><tr><td>Ours</td><td><b>1.59 ± 3.59</b></td><td><b>0.32 ± 0.2</b></td></tr></tbody></table>

## 5.3.3 Results with a magnetic crawler

We evaluate our approach with a real robotic system and a large plate to demonstrate its viability for a practical robotic inspection task. In our experimental setup shown in Fig. 5.7a, the robot is equipped with an onboard accelerometer to provide accurate heading measurements, and with a single transducer for both acoustic emission and reception. Also, a pump is activated to continuously bring water through the tether at the interface between the sensor and the plate
<figure><img src="image_51.png" /><figcaption>Figure 5.7: Combined feature-based and grid-based mapping with a robotic platform relying on UGWs. (a) shows the experimental setup. (b), (c), (d) and (e) present the mapping results at different steps along the robot transect. The robot pose (red arrow) and trajectory (red line) are also represented.</figcaption></figure>

for sufficient coupling. During the acquisition, the robot is driven along a trajectory on the plate positioned nearly vertically, and the acquisition of one measurement is triggered each time the sensor has traveled 5 cm. For simplicity, we test our approach by assuming flawless odometry data, even though it is imperfect, and there may be slight errors in the kinematic model. However, the mapping procedure is expected to integrate seamlessly within a SLAM framework. Figures 5.7b, 5.7c, 5.7d and 5.7e show the qualitative mapping results obtained with our approach at different steps of the robot transect. It can be observed that, at the early stage, the size of the plate is almost recovered as the estimated dimensions are 1798 × 1012mm during step (c), and the estimated orientation w.r.t. the robot is consistent with the reality. Besides, the inside and outside of the plate are correctly recovered, except for the regions near the corners, incorrectly marked as outside the plate. This slight inconsistency might be caused by the slight odometry errors and the increased complexity of the acoustic signals. Overall, these results demonstrate the applicability of our approach with a real robotic system. It may be useful for mapping a large structure made of several plates, as plate boundary detection is made robust, or by an active-sensing strategy for an optimal and autonomous mapping.

One of the limitations of our approach is that the parameters may be tuned carefully. Oth-
erwise, proper lines may eventually be filtered on the beamforming map. In that situation, it
is preferable to increase the noise parameters so that more measurements would be necessary
before a line is filtered out. Also, our approach is evaluated on an isolated plate in a laboratory
environment, while real metal structures present complex structural features such as welds or
stiffeners. The constraint of rectangular shapes is not a limitation for mapping storage tanks or
ship hulls, as they are almost entirely made of rectangular panels. However, extensions of our
approach are to be investigated to work with more general structures. It is expected it would
only require an extension of the optimization method presented in Chapter 4 to retrieve the
edges from the beamforming maps. At the same time, the design of the Bayesian framework
could be left unchanged.

5.4 Towards ultrasound-based autonomous exploration

The algorithms presented so far were dedicated to localization and mapping purposes: following a *predetermined* sensor trajectory for collecting the ultrasonic measurements, the localization and mapping performance of the different methods was evaluated. Yet, if one seeks to maximize the reconstruction performance, the question of *how* the robot should move is also central. For example, some robot trajectories may result in better map reconstruction than others if the measurements positions are such that the retrieved quantity of information is maximal. In our case study, automatic path planning and path following can be beneficial to the development of a fully automated inspection agent that can meet reconstruction optimality criteria. One could aim at maximization of the reconstruction quality of a large structure for reliable identification of potential defects, or minimize the reconstruction time to reduce inspection costs.

As a first step towards such an application, we introduce, in this section, an active-sensing
strategy based on frontier exploration [102] to enable the autonomous reconstruction of the ge-
ometry of a metal surface by a mobile robot relying on ultrasonic echoes. The method builds
on our framework for mapping with combined representations: during the mapping process, the
grid and feature-based map estimates of the unknown metal surface are leveraged to estimate the
information gain provided by measurements made at potential next positions. Following a util-
ity function that balances information gain and travel cost, the sensor is automatically directed
towards the best frontier point that lies in the boundary between "explored" and "unexplored"
areas. The developed method is evaluated in simulation, and is compared with multiple explo-
ration strategies, essentially closest and random frontier point selection. Finally, an experiment
using a mobile robot equipped with a co-localized pair of transducers is used to validate the
viability of the approach. It is acknowledged that what follows is a mere summary of the method
presented in [77].

5.4.1 Related works

Autonomous robotic exploration to map unknown environments is a major research field in robotics, aiming at algorithmic path planing and path following solutions to maximize information gain and minimize costs. Frontier-based approaches yield common solutions to the exploration problem, and are extensively used to map indoor environments [102, 27]. By moving towards new frontiers between mapped and unknown areas, a mobile robot can extend its map into new territory until the entire environment has been explored. The key to effective frontier exploration is the selection of target frontier points. In the original method [102], the closest frontier point is selected as the next target. In most cases, exploration strategies select the next best frontier by evaluating candidate locations according to different criteria, resulting in various extensions of the basic frontier-based exploration strategy [45, 63, 27]. In [98], for example, a utility function balancing the travel cost with the information gain is used.

There are currently different methods to evaluate the information gain: one is to directly
assess the undetected space size in the nearby region of the target frontier point [4]. Another
one is to calculate the information entropy to assess the uncertainty on map prediction [10, 94].
Yet, to the best of our knowledge, the use of such approaches for mapping using omnidirectional
UGWs in a pulse-echo setup has not been thoroughly studied in the literature.
A Lamb wave-based exploration strategy for mapping a metal plate is proposed in [62] using a pitch catch setup. Yet, the transducers are manually positioned in the desired (and sometimes distant) positions, which does not leverage constant signal acquisition over the continuous robot displacement. Furthermore, the mapping method only relies on detection of the first reflection on the plate boundaries, and it is subject to gridlock situations which are not predictable.

### 5.4.2 Overview of the frontier exploration approach

Our approach relies on the mapping framework with combined representations presented earlier. During Step t of the mapping process, the grid map estimate $\hat{O}_t$ is used to determine the $K_t$ frontier points $\mathbf{p}_{t,k}, k = 1..K_t$, that lie in the frontier between explored and unexplored areas on the plate. Along with the feature-based map estimate $\hat{\mathbf{M}}_t$ obtained via DAS beamforming, the grid $\hat{O}_t$ is used to evaluate the cumulative entropy:

$$I(\mathbf{p}_{t,k}, \hat{\mathbf{M}}_t, \hat{O}_t) = \sum_{\mathbf{p} \in V(\mathbf{p}_{t,k}, \hat{\mathbf{M}}_t)} E(\mathbf{p}),$$

of all the points in the vicinity $V(\mathbf{p}_{t,k}, \hat{\mathbf{M}}_t)$ of $\mathbf{p}_{t,k}$, that are expected to lie on the plate. $e(\mathbf{p})$ is the well known entropy over the state probability of a point $\mathbf{p}$ of the grid:

$$E(\mathbf{p}) = -\hat{O}_t(\mathbf{p}) \log \hat{O}_t(\mathbf{p}) - (1 - \hat{O}_t(\mathbf{p})) \log(1 - \hat{O}_t(\mathbf{p})).$$

The overall utility of a frontier point is defined as a weighted sum between the resulting entropy value and the distance between the sensor and the frontier point:

$$U(\mathbf{p}_{t,k}, \hat{\mathbf{M}}_t, \hat{O}_t) = I(\mathbf{p}_{t,k}, \hat{\mathbf{M}}_t, \hat{O}_t) - \lambda \cdot \|\mathbf{p}_{t,k} - \mathbf{x}_t\|,$$

where $\lambda$ is a positive parameter, and $\mathbf{x}_t$ the robot position at time t. Hence, by selecting the frontier point with maximum utility as next best point to visit, a moving robot will be driven towards nearby locations where the reconstruction uncertainty is still important.

During the whole mapping process, the feature-based and grid-based estimation procedures are being executed. Thus, as the maps, frontier points and their respective utility vary accordingly, the selection of the best frontier point to move to is repeated at every measurement step. Finally, a frontier point is always attributed a utility of 0 when it is considered closer to a boundary than a fixed distance, so that other frontiers can also be explored. The frontier exploration is considered complete when there are no more frontier points to evaluate regarding their distance to the edge. An illustration of a mapping output provided by the proposed approach is presented in Fig. 5.8.

To summarize, the different steps of one exploration step are the following:

* Frontier points are extracted from the grid map.

* The utility of each frontier point is evaluated.

* The candidate pose with the highest utility is selected as the next goal.

* The robot navigates to the target position.
<figure><img src="image_52.png" /><figcaption>Figure 5.8: Illustration of a mapping output provided by our approach. The blue dot represents the sensor position and the blue line is its overall trajectory. The red rectangle represents the estimate $\hat{M}_t$ obtained via DAS beamforming. The grid $\hat{O}_t$ is accounted for with the color of the cells. The valid (resp. invalid) frontier points corresponds to the red (resp. black) dots. The green point indicates the selected next best point.</figcaption></figure>

### 5.4.3 Overview of the exploration results

To assess the performance of our approach, we perform several tests conducted in simulation and with a real robot. During each experiment, as the mapping time increases, we keep track of the distance traveled by the robot, along with the ratio of covered surface of the plate as a performance criterion. The latter is defined as the ratio between the plate surface correctly labeled as "on the plate" and the entire plate surface. Contrary to the estimation errors on the plate geometry parameters, it is a good indicator of how spread the measurements positions are. Furthermore, it is relevant to aim at complete measurement coverage during a robotic inspection operation, as the entire surface of the structure is to be inspected. For both the simulated and real-world experiments, we will consider a flat rectangular metallic plate with dimensions 1700x1000x6mm.

First, we run our approach in a simulated scenario. In Fig. 5.9, we show the mapping outputs at different steps, along with the sensor trajectory resulting from the application of our exploration strategy. One can observe that, initially, the map is very uncertain. Thus, the selection of frontier point is essentially random. During Step 41, after a seemingly erratic trajectory due to high uncertainty on the map, three edges are correctly recovered. Yet, the mapping error on the grid is relatively important due to the right edge that is incorrectly recovered. Following the selected frontier points, the robot is driven to the right part of the plate, until all the points have been marked as invalid. This results in the correction of the errors in the grid map. During the last step, the sensor has followed the remaining frontier points until there is no valid point left, leading to termination of the exploration process with a seemingly low mapping error.

For quantitative comparison with baselines, we now also evaluate the coverage performance of three other exploration strategies: the classic frontier approach where the closest frontier point is systematically selected as next goal location, a random approach where the next frontier point is sampled from a uniform distribution, and our method using the true plate's geometry, which is expected to lead to the best coverage performance.

Each exploration algorithm is run for the same simulation setup: the metal plate has dimensions 1700x1000x6mm, and the sensor is initially placed at 20 cm from two boundaries, in a corner. To assess the repeatability of each exploration strategy, we evaluate the means along with the 10% and 90% quantiles of the mapping coverage. For better visualization, we only present coverage results up to 83% of mean coverage value. The results of the experiments can be seen in Fig. 5.10. It can be observed that the proposed approach outperforms the random
<figure><img src="image_53.png" /><figcaption>Figure 5.9: Mapping outputs obtained in simulation during different steps of the mapping process. The blue dot represents the sensor position, and the blue line is the sensor trajectory. The red rectangle represents the feature-based map estimate obtained via DAS beamforming. The grid $\hat{O}_t$ is accounted for with the color of the cells. The valid (resp. invalid) frontier points correspond to the red (resp. black) dots. The green point indicates the selected next best point.</figcaption></figure>

approach (which is the least efficient) and the standard closest point strategy, as higher mapping coverage values are obtained for lower traveled distances. The reason is that it takes both into account information gain and travel cost. Naturally, knowing the ground truth plate geometry leads to better results. It represents an upper bound on our approach performance.

Finally, the performance of our approach is experimentally validated using a real robotic platform. We use a turtlebot equipped with a structure for holding two transducers and maintaining them in contact with the surface, as depicted in Fig. 5.11. A layer of water is placed on the plate surface as a coupling medium. For signal acquisition, we use two-tone bursts of a sinusoidal wave at 100kHz as the excitation and a sampling frequency of 1,25 MHz for acquisition.

During all the experiments, we use a camera which provides a top view recording of each experiment. The tags depicted in the experimental setup shown in Fig. 5.11 are used to track the robot position on the plate using the AR tags tracking package available on ROS<fnref n="2" />. This system is used here for ground truth positioning. Simultaneous localization, mapping and exploration is left for future work. Also, to evaluate our exploration system independently to the quality of the controller, the movement between waypoints is implemented with the joystick. Yet, the destination points are automatically generated by our algorithm.

In the next, we compare the performance of the closest frontier point selection and random exploration strategies with the proposed approach. We perform five repetitions of each algorithm, each time placing the sensors at the same starting position corresponding to $[x_0, y_0] = [25cm, 45cm]$ expressed in the plate frame. The results are provided in Fig. 5.12, where we represent the evolution of plate coverage up to a value of 91% for the mean coverage. One may observe that the coverage percentage increases relatively fast for the proposed method, as an average coverage value of 83% is reached with a displacement of 2.38m. We remark that the random approach outperforms our method at the beginning, but plateaus around 80% coverage due to its inability to select informative waypoints to further improve the mapping results.

<footnote n="2">²http://wiki.ros.org/ar_track_alvar</footnote>
<figure><img src="image_54.png" /><figcaption>Figure 5.10: Evolution of the plate surface coverage for the four exploration strategies with respect to the traveled distance in the simulations. The discontinuous lines represent the coverage mean value over 50 repetitions of each method, while the upper and lower bounds of the colored areas represent respectively the 90% and 10% quantiles.</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>Figure 5.11: Top view of the experimental setup. A turtlebot platform is used along with AR tag markers for ground truth positioning.</figcaption></figure>
<figure><img src="image_56.png" /><figcaption>Figure 5.12: Evolution of the plate surface coverage for the three exploration strategies with respect to the traveled distance during the experiments. The solid lines represent the coverage mean value over 5 repetitions of each method, while the upper and lower bounds of the colored areas represent respectively the 90% and 10% quantiles.</figcaption></figure>

The closest point selection method is outperformed by the two others due to inefficiency of the greedy strategy. Overall, the results illustrate the potential of an active-sensing mapping system for automated inspection. Future works shall be directed towards implementation and testing in more realistic scenarios to further make possible the emergence of such an application.

## 5.5 Polygonal shapes reconstruction from acoustic echoes

In the previous sections, the mapping problem was limited to recovering the geometry of rect-
angular panels using UGWs. By leveraging DAS beamforming and the rectangular constraint,
one can efficiently recover the position of the plate boundaries, as the explicit detection of all
the echoes within the signals and their association with potential line boundaries is not required.
As seen in Chapter 4, this allows for integration of the approach within a FastSLAM framework
that can seamlessly run online. Yet, it is restricted to rectangular geometries.

This assumption may not be restrictive if one is only interested in recovering the geometry of
surfaces made of metal panels, and to achieve SLAM on structures like ship hulls or storage tanks.
However, reducing the number of assumptions on the map may be beneficial to a potential real-
world inspection application. More specifically, defects, which may also act as acoustic scatterers,
could be integrated in the map. Thus, they could be detected using the same DAS beamforming
approach. In this situation, however, the rectangular constraint can no longer hold for recovering
them. In addition, DAS beamforming could also be used for mapping more general structures,
with more arbitrary geometries.

The objective of this last section is to propose a more general approach based on beamforming
along with combined map representations to recover arbitrary 2-dimensional polygonal geome-
tries from acoustic echoes. It applies DAS beamforming in the geometrical boundary space, and
relies on the *single* detection of the closest boundary to map the inner surface of the geometric
shape, using a grid representation. The feature-based geometry estimate is then constructed
via minimization of a regularized cost function using a simulated annealing optimizer to select
boundaries, on the beamforming results, consistent with the mapped inner space. We also de-
sign a boundary rejection criterion to recover their exact number based only on a specified upper bound. We test our approach in simulation for various geometries and actual conditions using experimental ultrasonic measurements acquired on a metal plate sample. The results show that it can successfully recover the geometric shape of the reverberant environment, the exact number of boundaries, and has the same performance as the baseline method, presented earlier, whose use is limited to recovering rectangular shapes.

### 5.5.1 Related works

The use of ultrasonic echoes for recovering polygonal shapes such as a room geometry, and using a single mobile unit, has been studied in various works that have been mentioned earlier [51, 53, 69]. However, these are only tested in simulation. They rely on essential assumptions, such as first-order echoes are directly available, the echo-wall association is being flawlessly known, or the number of walls to recover is being known a priori. Due to noisy measurements and inaccurate signal models, these assumptions may not always hold for a real-world scenario, resulting in low accuracy reconstruction, as illustrated by the experiments in Chapter 4.

When deployed in real-world scenarios, other acoustic methods for mapping indoor environments rely on directivity of the signal acquisition [79]. Another setup, which is more computationally involving, is to use an array of sensors for DoA estimation [80, 61]. Yet, the mentioned works only rely on the earliest signal in the measurements, whereas late reflections also convey information on the environment. Overall, the efficient mapping of arbitrary polygonal shapes from omnidirectional acoustic echoes is still to be established.

### 5.5.2 Problem statement

Let's first re-introduce the considered problem. We consider a setup where a single (collocated) emitter/receiver acoustic device is driven to acquire acoustic measurements $z_{0:T}$ at positions $x_{0:T}$ on a bounded metal surface, where $x_i = [x_i, y_i, \theta_i]^T$ contains the two-dimensional position and orientation of the robot. At the $i$-th scanning position, a pulse $s(t)$ is used to emit an omnidirectional ultrasonic wave. The received signal $z_i(t)$ contains the reflections of the incident wave on the plate boundaries. We also assume that, at each position $x_i$, an accurate measurement (or estimate) of the range $\rho_i$ to the closest boundary is available by processing the acoustic data. This measurement may only be subject to little noise in practice, as it may be one of the most energetic echoes. Based on all range measurements $\rho_{0:T}$, the surface $A_T$ is obtained, as a grid map, using the unions of all disks centered at sensor positions $x_i$, and whose radius equal to the detected range $\rho_i$. By construction, this area cannot be crossed by a boundary. For simplicity of the presentation, we assume that noise-free odometry information (i.e linear and rotational displacement between the positions of the successive measurements) is available. Yet, the proposed method may be integrated in a FasSLAM algorithm as in Chapter 4 for localization.

Our objective is to recover the shape $M$ of the plate that we now assume polygonal, instead of the much more limiting assumption of a rectangular geometry that has been made so far, using ultrasonic measurements. One may note that not only we only address the mapping problem here, but we assume little noise on the measurement $\rho_i$. To alleviate the need of the latter assumption, one may integrate the proposed method within the mapping framework with combined representations presented in the first part of this chapter. Indeed, it has been observed that the mapping results were robust w.r.t echo detection errors.
<figure><img src="image_57.png" /><figcaption>Figure 5.13: Schematic of the problem setup and notations. An example of shape to estimate is represented by the red polygon. A candidate geometry estimate is represented by the set of green lines, whose $(r, \alpha)$ coordinates are only provided for one of them. The area A formed with a single range measurement $\rho$ is represented by the lighter area. The inner surfaces of A and $\hat{M}$ are delimited with the discontinuous dashed lines.</figcaption></figure>

### 5.5.3 A regularized cost for polygonal shape reconstruction

We propose to estimate $\mathbf{M}$ using the (now very classic) DAS beamforming approach, and only adapt the cost function whose minimization provides a feature-based geometry estimate. This is achieved by adding a regularization term to the initial cost function for maintaining consistency between the feature-based map estimate and a grid representation of the plate inner surface. The latter is obtained using the range measurements $\rho_i$.

The geometry $\mathbf{M}$ to recover is assumed to be a 2D convex polygon that can be represented by its linear boundaries: $\mathbf{M} = \{r_l, \alpha_l\}_{l=1..L}$ where the parameters $(r_l, \alpha_l)$ define the line equation in the 2D plane, expressed in a fixed frame whose origin will be chosen as the initial robot position, and the $x$ axis will be aligned with the initial robot heading. The number $L$ of boundaries is assumed a priori unknown. The adopted notations are precised in Fig. 5.13.

First, the envelope measurements $e_i$ are classically obtained by correlating the measurements $z_i$ with a wave propagation model. Next, the beamforming map to assess the likelihood of existence for any line boundary $(r, \alpha)$ is determined with:

$$ \mathcal{L}_T(r, \alpha) = \sum_{i=0}^{T} e_i (|x_i \cdot \cos \alpha + y_i \cdot \sin \alpha - r|), $$

Based on the former mapping approach, only rectangular shapes could be recovered by retrieving exactly four local maxima of $\mathcal{L}_T$ through the maximization of their cumulative likelihood. We propose to extend this approach to recover, in principle, arbitrary 2D polygonal shapes. The estimated geometry has the form: $\hat{\mathbf{M}} = \{\hat{r}_k, \hat{\alpha}_k\}_{k=1..K}$, where the maximum number of boundaries $K$ is such that $K \ge L$. This is a very minimal hypothesis, as $K$ can be set much higher than the unknown $L$.

Next, for a given robot trajectory $\mathbf{x}_{0:T}$, we determine the inner surface $A_T$ using the range measurements $\rho_{0:T}$. Based on the beamforming map $\mathcal{L}_T$ that we normalize by its maximum value, we introduce the cost function:

$$ J(\hat{\mathbf{M}}) = - \left( \sum_{(r_k, \alpha_k) \in \hat{\mathbf{M}}} \mathcal{L}_T(r_k, \alpha_k) + \lambda \frac{\hat{\mathbf{M}} \cap A_T}{\hat{\mathbf{M}} \cup A_T} \right), \quad (5.5) $$
where the first term aggregates the likelihood of presence of the estimated boundaries, and the second term is the well known Intersection over Union (IoU), widely used for computer vision purposes [75]. In Eq. (5.5), the surfaces $\hat{\mathbf{M}}$ and $A_T$ are as described in Fig.5.13. The IoU term (that we will simply write as IoU($\hat{\mathbf{M}}$, $A_T$)) is always between 0 and 1. It is zero when $\hat{\mathbf{M}}$ and $A_T$ have an empty intersection, and it is 1 when the two surfaces are equal. $\lambda$ is a positive parameter that balances the importance between the boundaries likelihood and IoU terms. Our geometry estimate is then defined as:

$$ \hat{\mathbf{M}}_T = \arg \min_{\hat{\mathbf{M}}} J(\hat{\mathbf{M}}). $$

Hence, the minimization of the cost criterion leads to maximization of the cumulative likelihood of the estimated boundaries and maintains consistency between the estimated inner surface of $\hat{\mathbf{M}}$ with that of $A_T$. In other words, the IoU in Eq. (5.5) acts as a regularization term for selecting likely boundaries on the beamforming results.

### 5.5.4 Shape reconstruction with simulated annealing

An optimization strategy is required to recover a geometry estimate by minimizing the cost. Simulated annealing (SA) [49] is a relatively simple and efficient approach for estimating a global minimum of a non-linear cost function. This metaheuristic is inspired by metallurgy, where, to form a perfect crystal (which corresponds to the state of minimal energy), a pure liquid substance is slowly cooled. A random perturbation is applied to the current state at each iteration of the optimization process. The perturbation is systematically accepted if the energy is decreased. It is accepted according to a Boltzmann probability distribution when the energy is increased to escape local minima. The temperature, which is a hyperparameter that controls the acceptance rate, is decreased slightly after each iteration, so that the probability of increasing the energy goes to zero. Even though simulated annealing offers no guarantees that the optimum will be found, this strategy can lead to satisfying solutions in a relatively restricted number of iterations. In the following, we detail the different steps of our custom implementation for mapping arbitrary polygonal geometries from acoustic echoes.

#### Initialization

First, the geometry estimate is randomly initialized with $\hat{\mathbf{M}}_T = \{\hat{r}_k, \hat{\alpha}_k\}_{k=1..K}$, where the $\hat{r}_k$ and $\hat{\alpha}_k$ are sampled from independent uniform distributions.

#### State disturbance

At iteration $p$ of the optimization process, each boundary of the estimate is iteratively disturbed with a random process. For the $k$-th boundary, we sample:

$$ \begin{cases} \tilde{r}_k = \max(\hat{r}_k + \Delta_r \cdot \chi_r^3, 0) \\ \tilde{\alpha}_k = \hat{\alpha}_k + \Delta\alpha \cdot \chi_\alpha^3 \mod 2\pi \end{cases}, $$

where $\Delta_r$ and $\Delta_\alpha$ are positive parameters, and $\chi_r$ and $\chi_\alpha$ are independently sampled from uniform distributions: $U([-1, 1])$. They are raised to the power 3 so that small variations are more likely. Still, large variations (i.e in the order of $\Delta_r$ and $\Delta_\alpha$) can be possible to explore regions of the search space that are far occasionally from the current state. Furthermore, to speed up the optimizer convergence, we determine:

$$ (\tilde{r}_k, \tilde{\alpha}_k) \leftarrow \arg \max_{(r, \alpha) \in V(\tilde{r}_k, \tilde{\alpha}_k)} \mathcal{L}_T(r, \alpha), $$
where $V(\tilde{r}_k, \tilde{\alpha}_k)$ denotes the set of boundaries that lie in the neighborhood of $(\tilde{r}_k, \tilde{\alpha}_k)$ that is to be appropriately defined. The purpose of this step is to only select $(\tilde{r}_k, \tilde{\alpha}_k)$ that are local maxima, which is a property that is expected for the true boundaries. The former edge $(\hat{r}_k, \hat{\alpha}_k)$ in $\hat{\mathbf{M}}_T$ is then substituted to the sampled one, yielding a new geometry $\tilde{\mathbf{M}}_T$.

### Simulated annealing acceptance process

The new candidate map is kept as the current state following an SA acceptance process. First, we determine the cost variation $\Delta J = J(\tilde{\mathbf{M}}_T) - J(\hat{\mathbf{M}}_T)$. Let's $\chi$ be a sample from $U([0, 1])$. The candidate geometry is kept only if we have:

$$ \chi < \exp \left( -\frac{\Delta J}{T(p)} \right), $$

where $T(p)$ is the temperature at iteration $p$ of the optimizer. Hence, if the cost is decreased, the new map is systematically kept as the new current state $\hat{\mathbf{M}}_T$. Otherwise, it is kept according to the Boltzmann distribution to escape local minima.

### Edge filtering

As it is expected that the edges number parameter $K$ is overestimated, some of them may not contribute significantly to the definition of the recovered geometry once the optimizer has converged. To derive a proper geometry with a correct number of boundaries, we apply a rejection criterion to remove non-informative edges. We iterate over each boundary $(\hat{r}_k, \hat{\alpha}_k)$ in $\hat{\mathbf{M}}_T$, and we remove it to yield a new candidate geometry $\tilde{\mathbf{M}}_T$. We then assess whether:

$$ \left| \frac{\text{IoU}(\tilde{\mathbf{M}}_T, A_T) - \text{IoU}(\hat{\mathbf{M}}_T, A_T)}{\text{IoU}(\hat{\mathbf{M}}_T, A_T)} \right| < \epsilon, $$

where we will take $\epsilon = 10^{-2}$. If the IoU variation is lower than the threshold, then the $k$-th boundary can be removed, as its contribution to the definition of the geometric surface is considered negligible. Then, the geometry $\tilde{\mathbf{M}}_T$ can be taken as the new candidate $\hat{\mathbf{M}}_T$, and the process is iterated on the next boundary. Also, the rejection criterion is not applied if the candidate map $\hat{\mathbf{M}}_T$ has 3 edges (which is a lower bound for a polygonal shape).

### 5.5.5 Simulation results

We perform different experiments to assess the efficiency of our method for recovering arbitrary 2-dimensional polygonal geometries from acoustic echoes acquired over a mobile unit trajectory. We first consider three simulated scenarios, where the geometry to be recovered has either 3, 4 or 5 boundaries. The robot paths are randomly generated with successive linear and rotational displacements. The ultrasonic measurements acquired over the path are simulated using the Image source model, and the ultrasonic model $g$ introduced in Chapter 2. It should be noted that such a data model leads to over-energetic high-order reflections, as acoustic absorption from the boundaries is not taken into account. The excitation $s(t)$ is a two tone bursts sinusoidal wave at 100 kHz. We add Gaussian noise on the measurements to maintain a fixed signal-to-noise ratio SNR = 0. The considered geometries and paths are provided in Fig. 5.14.

We run 10 repetitions of the optimizer for each scenario to assess the repeatability of our approach. The optimizer is set with parameters $\lambda = 3, T(p) = 5/p^{0.8}$, and $V(r, \alpha)$ is chosen as
<figure><img src="image_58.png" /><figcaption>Figure 5.14: The different simulated setups along with an instance of mapping results. On the left of each sub-figure, the outline of the ground-truth geometry is provided in discontinuous red, the green lines account for the estimated edges, and the robot path and acquisition positions are represented in black. The corresponding beamforming maps are provided on the right, where the red rectangles indicate the retrieved edges. The edge coordinates are expressed in the fixed frame, whose origin is the initial sensor position *O*, and the *x*-axis is aligned with the initial robot heading.</figcaption></figure>
<figure><img src="image_59.png" /><figcaption>Figure 5.15: Evolution of the performance criteria for the 3 scenarios, and for 10 repetitions of the optimizer. The top plots show, in blue, the evolution of the costs. The middle plots depict the estimated number of edges after application of the edge rejection criterion. The bottom plots show the IoU between estimated and ground-truth surfaces. The median values are in red.</figcaption></figure>

the set of edges with line coordinates within the rectangular area centered in (r, α), and with size 5.4cm × 48 °. Also, we use $K = 10$ as the maximum number of edges. Fig. 5.15 shows the evolution of the cost, as well as the number of estimated edges after use of the edge rejection criterion during 300 optimizer iterations, and for the 10 repetitions. We also display the evolution of IoU($\hat{M}$, M), which is the IoU between the estimated inner surface (based on the estimated edges) and the ground truth surface M. Hence, an IoU of 1 translates into a perfect geometry estimation. The results demonstrate that, for all the considered scenarios, our approach can seamlessly recover, after convergence, both the actual number of boundaries and the geometric shape of the simulated plates.

### 5.5.6 Experimental results

To validate our approach, we run it using the experimental ultrasonic data acquired on our 600x450x6mm rectangular aluminium metal plate by moving, by hand, two nearly collocated contact piezoelectric transducers for signal emission and reception. The optimizer is set with the *exact* same parameters than those used for the simulation tests to highlight the robustness of our approach. Also, we do not use ground-truth ranges $ρ_{0:T}$, but apply a standard peak detection method on the correlation signals to estimate them. This only leads to a slight estimation error due to high SNR. The considered simulated robot trajectory is depicted in Fig. 5.16a.

Fig. 5.16 shows the results and the evolution of the performance criteria for 10 repetitions of the optimization process. For comparison, we also display the IoU achieved with the baseline method introduced in Chapter 4, which does not require iterative optimization, but is restricted to rectangular shapes. We can see that the estimated number of edges converges seamlessly to 4, while the IoU goes to 1 for all the repetitions. Thus, our method achieves the performance of the baseline approach, despite not being provided with the ground-truth number of edges nor the rectangular constraint. Overall, the results demonstrate the potential of our approach for recovering polygonal shapes from acoustic echoes using a single mobile omnidirectional emitter/receiver device and beamforming.

Overall, we have provided an extension of the mapping approach first introduced in Chapter 4 to recover the geometry of unknown 2D polygonal structures from acoustic echoes, and using a single omnidirectional mobile sensor. It combines delay-and-sum beamforming in the geometric
<figure><img src="image_60.png" /><figcaption>(a) Mapping setup and results for the experimental scenario.</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>(b) Evolution of performance criteria during the optimization.</figcaption></figure>

Figure 5.16: Mapping results for the experimental setup. (a, left) shows the simulated sensor path (in black), the plate outline (discontinuous red) and estimated edges (green). The corresponding beamforming map is provided on the right, where the red rectangles indicate the retrieved edges. (b) shows the evolution of the loss, the number of estimated edges and the value of IoU($\hat{M}$, M) for 10 repetitions of the optimizer. The black horizontal line on the IoU plot is the value achieved with the baseline method. The median values are in red.
boundary space and a simulated annealing optimizer to retrieve boundary estimates from the beamforming results. Contrary to previous works, our approach does not rely on solid hypotheses on echo detection, identification and association, and can recover, in principle, arbitrary polygonal shapes. In future works, the approach shall be evaluated in more realistic environments. Furthermore, the optimization efficiency shall be improved, so that it can be run online. It shall be integrated within a SLAM framework, as the sensor positions are to be estimated in practice, and its extension to recover 3D volumes may also be investigated.

**5.6 Chapter summary**

In this chapter, we essentially demonstrate the potential benefits of using combined representations along with UGWs measurements acquired by a mobile robot for mapping a metal structure. The core of this approach is the simultaneous recovery of a grid and a feature-based representation of a same metal plate environment. This allows maintaining both denseness and sparsity on the map belief, while leveraging the consistency properties between the maps.

First, we propose a Bayesian framework for building such joint maps. By leveraging the simplicity of DAS beamforming and the Bayesian binary filter, the proposed framework can fully and consistently map a rectangular metal plate. Compared to the feature-based only beamforming mapping method, the proposed approach is less sensitive to interference, it provides better uncertainty representation, while being robust to detection errors of the closest boundary.

In the second part of this chapter, we leverage such mapping framework for autonomous exploration, with the aim to automatically drive an autonomous robot for efficient mapping and coverage of a rectangular metal plate surface. The proposed algorithm, which relies on frontier exploration and on uncertainty assessment on the grid map belief, is evaluated both in simulation and in the real world. The results show that it is reasonable for the target application.

Finally, we propose to use a grid map of the inner surface of a plate to recover arbitrary polygonal shapes from ultrasonic echoes. Hence, we alleviate the need of the rectangular assumption on the feature-based map. The proposed approach, which introduces a regularization term on the cost criterion to retrieve line features on the DAS beamforming results, has proven to be efficient both in simulation and in a real-world setup to recover more arbitrary polygonal geometries, without prior knowledge of the number of line features to estimate. It is envisioned that the proposed frameworks can support the emergence of a true long-range robotic (and potentially automated) inspection operation.
# Chapter 6

## Graph-based Simultaneous Localization and Mapping

In this chapter, we propose a non-probabilistic GraphSLAM framework to solve the full SLAM problem using UGWs measurements. Contrary to online SLAM, which has been considered so far in this thesis, GraphSLAM can use all the measurements all at once to infer the trajectory and map estimates, usually leading to much more accurate results. We first elaborate an observation model that can leverage multi-path multi-order acoustic echoes for localization, while DAS beamforming is used for mapping the plate boundaries. The two methods are subsequently merged within a nonlinear least squares optimizer to solve the full SLAM problem, without breaking the natural sparsity of the SLAM problem, in first approximation. We test our approach in real-world settings, on metal plates of varying size and material. The results demonstrate the improved localization accuracy of the multi-order echo observation model compared to a second model, from the previous chapters, that relies solely on first-order echoes. We also show that the entire robot trajectory and plate boundaries can be simultaneously and reliably recovered in two different scenarios, even when odometry noise induces severe drift.

### 6.1 Motivation and background

In general, a SLAM problem can come at two different tastes. When considering an online SLAM approach, the belief on the map and robot trajectory is often incrementally updated when a new measurement is made available. Solving the full SLAM problem, however, involves high-dimensional batch optimization to simultaneously recover the entire robot trajectory and map features based on all the observations. A conventional solution is provided by GraphSLAM [35]. The latter usually involves optimization over a sparse graph, where the unknown variables (robot positions and map features) are the nodes, and the edges between them represent information constraints that are induced either by odometry or measurements inputs. The objective is to find the unknown variables that minimize a cost criterion that is defined as a weighted sum of the individual constraint terms.

Unlike online filtering techniques, such as Kalman filtering or particle filtering [94], GraphSLAM can use the information contained in *all* the measurements to infer the target variables. This approach is typically known to lead to more accurate results, but at the cost of a higher computational effort. A promising prospect is that acoustic GraphSLAM, through a complete optimization over all the measurements, could provide by itself accurate robot localization, structure geometry mapping, and inspection results. Nevertheless, applying GraphSLAM techniques
to UGWs measurements is challenging due to the inherent complexity of ultrasonic measurements
that present wave dispersion effects and multi-echoes characteristics.

To make one step towards long-range robotic inspection of large metal structures, we propose, in this chapter, a non-probabilistic GraphSLAM approach to simultaneously recover the robot trajectory and individual plate geometry from all odometry and acoustic measurements. Our objective is to enable the resolution of the full SLAM problem for accurate results, even in challenging acquisition conditions, as expected in a real-world inspection scenario. We introduce an observation model that can leverage the information contained in high-order reflections from the metal plate boundaries, and we rely on DAS beamforming for plate boundaries localization. The two methods are subsequently integrated within a Non-linear Least Squares (NLS) optimizer, which iteratively updates the trajectory and map estimates after linearization of a measurement likelihood cost criterion [48]. We show experimentally that the designed multi-order echo observation model leads to better localization accuracy when compared to the simpler first-order echo only model, introduced in the previous chapters. Furthermore, we demonstrate that the designed GraphSLAM approach can successfully solve the full SLAM problem in two different experimental setups, even when the drift induced by odometry noise is severe. Last, a study on the sensitivity of the SLAM results to the cost parameters is carried out.

6.1.1 Multi-order echoes signal processing

While traditional acoustic localization and mapping methods rely solely on the reception of the direct-path incident signal from the source, the additional integration of multi-path scattering signals is a good option for more accurate and more robust results. In UGW-based NDE, for example, the integration of multi-path signals can offer better defect detection capabilities [38, 105, 106]. Indeed, all multi-path echoes contain further information on the environment. Hence, using it (instead of treating it like noise) may lead to better results. However, due to the complexity of the signals, appropriate (and often more complex) models and algorithms are needed.

By doing a parallel with the former example, the localization methods presented earlier rely
solely on detection of the first-order reflections on metal plate boundaries - or equivalently, they
rely on direct incident signals originating from fictive first-order image sources. We now propose
an observation model that can also leverage the high-order reflections in UGW measurements
for more accurate and robust localization and mapping results on metal panels.

6.1.2 Related works on GraphSLAM

GraphSLAM is an approach that is now well studied [35]. In the literature, numerous Graph-SLAM methods and solvers have been developed, with increasing computational efficiency [55, 21, 47, 46]. These techniques have gained such maturity that they can seamlessly adapt to various sensor modalities (Laser, computer vision, IMU...) and SLAM problems. The chore of state-of-the-art methods is to successively linearize the models to find a minimum of the cost criterion through non-linear least-squares. Furthermore, by leveraging the sparsity of the resulting Jacobian matrix, GraphSLAM can be used to recover increasingly large trajectories and maps, that can contain up to millions of features.

In the literature, acoustic SLAM using a single omnidirectional emitter/receiver device is the subject of various studies to enable, for example, simultaneous sensor localization and room geometry reconstruction from acoustic echoes due to reflections on the walls. Batch optimization methods are proposed in [52, 53] to solve the full SLAM problem. However, this is typically a
challenging task, as the acoustic data contain range-only information, and due to the multiple unlabeled echoes present in the measurements. For these reasons, these works are only tested in simulation, they don't explicitly address the problem of echo detection, and rely on essential assumptions, such as only first-order echoes are being observed, or the echo-map feature associations are flawlessly known. In [16], a graph-based SLAM approach that also integrates beamforming for DOA estimation is proposed for underwater localization with sound beacons. However, it requires a sensor array, and only direct incident signals are used.

## 6.2 A graph-based framework for ultrasonic SLAM

In the following, we present a novel GraphSLAM approach to solve the full SLAM problem by relying on ultrasonic measurements. We first introduce an observation model that can leverage information from high-order acoustic echoes for accurate localization. The latter approach is integrated, along with DAS beamforming, into a NLS optimizer.

### 6.2.1 Localization with a multi-order echo observation model

In a first step, we design an observation model that can leverage high-order reflections from UGWs measurements for localization only on a metal panel. Assuming a known plate geometry $\mathbf{M}$, we model, for a candidate sensor position $\mathbf{x}$, the expected envelope signal as a mixture of Gaussian functions:

$$ \hat{e}(\mathbf{x}, r) = \sum_{\mathbf{p} \in \mathcal{I}(\mathbf{x}, \mathbf{M})} \mu_{\mathbf{p}} \cdot \exp \left\{ - \left( \frac{r - \|\mathbf{p} - \mathbf{x}\|}{\sigma} \right)^2 \right\}, \quad (6.1) $$

where $\mathcal{I}(\mathbf{x}, \mathbf{M})$ is the set of image source positions when the emitter is situated in $\mathbf{x}$, and the parameter value of $\sigma$ can be chosen as the signal wavelength. As the amplitudes $\mu_{\mathbf{p}}$ may not be known a priori, we propose to estimate them by minimization of the mean-squared error between the modeled signal $\hat{e}(\mathbf{x}, r)$ and the envelope measurement $e_i(r)$: we first discretize the different variables: $\hat{e}(\mathbf{x}) = [\hat{e}(\mathbf{x}, r_1), \hat{e}(\mathbf{x}, r_2), ..., \hat{e}(\mathbf{x}, r_M)]^T$ and $\mathbf{e}_i = [e_i(r_1), e_i(r_2), ..., e_i(r_M)]^T$, $M \in \mathbb{N}$, and define a "steering" vector function:

$$ \mathbf{a}(\mathbf{x}, \mathbf{p}) = \left[ \exp \left\{ - \left( \frac{r_1 - \|\mathbf{p} - \mathbf{x}\|}{\sigma} \right)^2 \right\}, \dots, \exp \left\{ - \left( \frac{r_M - \|\mathbf{p} - \mathbf{x}\|}{\sigma} \right)^2 \right\} \right]^T $$

Let's call $N(\mathbf{x})$ the number of image sources $\mathbf{p}_n$, $n = 1..N(\mathbf{x})$, in $\mathcal{I}(\mathbf{x}, \mathbf{M})$ whose signal can be observed in $\mathbf{e}_i$, as they are such that $\|\mathbf{p}_n - \mathbf{x}\| < r_M$. After definition of the $M \times N(\mathbf{x})$ matrix:

$$ \mathbf{A}(\mathbf{x}, \mathbf{M}) = [\mathbf{a}(\mathbf{x}, \mathbf{p}_1), \mathbf{a}(\mathbf{x}, \mathbf{p}_2), \dots, \mathbf{a}(\mathbf{x}, \mathbf{p}_{N(\mathbf{x})})], $$

the expected envelope measurement obtained at a candidate position $\mathbf{x}$ is:

$$ \hat{e}(\mathbf{x}, \mathbf{M}, \mathbf{e}_i) = \mathbf{A}(\mathbf{x}, \mathbf{M}) \cdot \hat{\mu}(\mathbf{x}, \mathbf{M}, \mathbf{e}_i), $$

where $\hat{\mu}(\mathbf{x}, \mathbf{M}, \mathbf{e}_i)$ is the linear least square estimate:

$$ \hat{\mu}(\mathbf{x}, \mathbf{M}, \mathbf{e}_i) = [\mathbf{A}^T(\mathbf{x}, \mathbf{M})\mathbf{A}(\mathbf{x}, \mathbf{M})]^{-1} \mathbf{A}^T(\mathbf{x}, \mathbf{M}) \cdot \mathbf{e}_i. \quad (6.2) $$
<figure><img src="image_62.png" /><figcaption>Figure 6.1: Positions of the image sources (red dots) for a rectangular plate with dimensions $w \times h$ (blue), and a fixed emitter position $[x, y]$ (blue dot).</figcaption></figure>

The positions of the image sources required to determine $\hat{e}$ can be analytically derived for simple geometries based on the image source algorithm. For a rectangular panel with dimensions $w \times h$ and a source at $\mathbf{x} = [x, y]$ expressed in the plate frame, we have [2]:

$$ \mathcal{I}(\mathbf{x}, \mathbf{M}) = \{[2nw \pm x, 2mh \pm y]^T/(n, m) \in \mathbb{Z}^2/(0, 0)\}. $$

Fig. 6.1 illustrates the regular pattern formed by the image source positions, given a fixed real source. Fig. 6.2 illustrates how the envelope signal model can match the measurements. Fig. 6.2a shows an experimental ultrasonic measurement obtained on our aluminum plate, Fig. 6.2b is the corresponding envelope signal obtained after correlation with a wave propagation model expressed in the spatial domain. In Fig. 6.2c, the latter is superimposed with the expected envelope signal, which is calculated with Eq. (6.2) based on the ground-truth sensor position and plate size. It can be seen that there is a good match between the measured and expected envelopes, as there is an appropriate correspondence for all the positions of the local maxima.

Finally, to localize the robot on a metal plate given an envelope measurement $\mathbf{e}$, one may simply use the mean squared error between expected and measured envelope signals as a likelihood function of the sensor position:

$$ l(\mathbf{x}|\mathbf{e}) \propto -\|\hat{\mathbf{e}}(\mathbf{x}, \mathbf{M}, \mathbf{e}) - \mathbf{e}\|^2. $$

The position of the maximum of this (highly non-linear) function yields a Maximum Likelihood (ML) robot position estimate. This will be further illustrated in the results section.

### 6.2.2 Mapping only with beamforming

Within the GraphSLAM framework, we propose to address the mapping problem using the classic DAS beamforming approach. As GraphSLAM necessitates the resolution of a highly non-linear problem, the integration of DAS beamforming may efficiently reduce the optimization effort. Indeed, beamforming may be used to integrate prior knowledge on how the map can be inferred given a candidate trajectory. Assuming a known $\mathbf{x}_{0:T}$, an estimate $\hat{\mathbf{M}}(\mathbf{x}_{0:T})$ of $\mathbf{M} = \{\mathbf{m}_l\}_{l=1..4} = \{r_l, \alpha_l\}_{l=1..4}$ is obtained by solving:

$$ \hat{\mathbf{M}}(\mathbf{x}_{0:T}) = \arg \max_{\mathbf{M}} \sum_{(r, \alpha) \in \mathbf{M}} \mathcal{L}_{\mathbf{x}_{0:T}}(|x_i \cdot \cos \alpha + y_i \cdot \sin \alpha - r|), $$
<figure><img src="image_63.png" /><figcaption>(c) The envelope signal shown in (b) superimposed with the expected envelope (in red) after matching the model in Eq. (6.1) with the envelope signal. A value σ = 2cm is used.</figcaption></figure>

Figure 6.2: Illustration of the envelope signal obtained from an experimental ultrasonic mea-
surement acquired on a rectangular metal panel, and comparison with the expected envelope
calculated with the matched Gaussian mixture model.
<figure><img src="image_64.png" /><figcaption>Figure 6.3: Graph representation of the (sparse) SLAM problem, where the variables to estimate are the positions $\mathbf{x}_i$ and the map $\mathbf{M}$, the inputs are the odometry data $\mathbf{u}_i$ and measurements $\mathbf{z}_i$. The constraints (induced by the external measurements) between the variables are represented by the graph edges, which are related to the $l$ functions.</figcaption></figure>

where $\mathcal{L}_{\mathbf{x}_{0:T}}$ is the usual DAS beamforming map:

$$ \mathcal{L}_{\mathbf{x}_{0:T}}(r, \alpha) = \sum_{i=1}^{T} e_i(|x_i \cdot \cos \alpha + y_i \cdot \sin \alpha - r|), $$

and $\mathbf{M}$ is restricted to be a rectangle. It should be noted that we restrict ourselves to the rectangular case for simplicity, even though we have proposed a methodology, in the previous chapter, to extend the geometry estimation to more arbitrary polygonal shapes.

### 6.2.3 The GraphSLAM cost criterion

In this section, we formulate the full SLAM problem using ultrasonic measurements as the minimization of a GraphSLAM criterion. GraphSLAM aims at finding a maximum a posteriori estimate of the full posterior distribution over the map and robot trajectory, given odometry and ultrasonic measurements. This amounts to minimizing the negative logarithm of the distribution that we model as a sum of non-linear mean-squared error terms:

$$ J(\mathbf{x}_{0:T}, \mathbf{M}) = \sum_{i=0}^{T-1} \|\mathbf{x}_{i+1} - f_{\text{dyn}}(\mathbf{x}_i, \mathbf{u}_i)\|^2 + \lambda \cdot \sum_{i=0}^{T} \|\hat{\mathbf{e}}(\mathbf{x}_i, \mathbf{M}, \mathbf{e}_i) - \mathbf{e}_i\|^2, \quad (6.3) $$

where $\lambda$ is a parameter that is used to balance the importance between odometry and ultrasonic measurements in the cost, $f_{\text{dyn}}$ is the transition model, and $\mathbf{u}_i$ are the odometry inputs. Overall, the introduced cost consists of two parts: the first sum accounts for constraints due to odometry data, and the second one accounts for the constraints due to measurements. Hence, the cost minimization simultaneously leads to a trajectory coherent with odometry inputs, and maximizes the likelihood of sensor presence based on accurate reconstruction of the envelope measurements. The corresponding graph of the SLAM problem is represented in Fig. 6.3. One may observe that it is a sparse graph. Due to the Markov property, the total number of edges grows linearly with the number of pose nodes, which leads to a sparsely connected graph. If the latter was fully connected, the increase of the number of edges would be quadratic with the number of nodes.
<figure><img src="image_65.png" /><figcaption>Figure 6.4: Graph representation of the SLAM problem after integration of DAS beamforming. The variables to estimate are the positions $\mathbf{x}_i$, and the map $\mathbf{M}$ is hidden as it can directly be inferred from all $\mathbf{x}_i$ using DAS beamforming. For simplicity of the graph representation, an edge is drawn to account for the existence of a constraint between two variables.</figcaption></figure>

In our case-study, the minimization of $J$ through search in the joint localization and map space can be quite inefficient. To reduce the problem complexity, the cost is only evaluated for candidate trajectories $\mathbf{x}_{0:T}$ and the corresponding map obtained by DAS beamforming $\hat{\mathbf{M}}(\mathbf{x}_{0:T})$. This enables to leverage the efficiency of DAS beamforming to integrate prior knowledge on how the map can be inferred given a trajectory. Thus, in the next, we will omit the map in the cost that will be reduced to:

$$J(\mathbf{x}_{0:T}) \triangleq J(\mathbf{x}_{0:T}, \hat{\mathbf{M}}(\mathbf{x}_{0:T})). \qquad (6.4)$$

In this situation, however, the map $\mathbf{M}$ becomes a "hidden" variable of the problem. As a result, the graph is no longer sparse, as the map introduces direct dependency between all the poses. As the graph of the SLAM problem has become fully-connected, as represented in Fig. 6.4, the resolution of the full SLAM problem can no longer be computationally efficient.

### 6.2.4 Sparse non-linear least squares minimization

Finding a minimum of Eq. (6.3) can be achieved through NLS minimization, which is conventional to solve GraphSLAM problems [35]. The principle is to iteratively update an initial guess of the $\{\text{trajectory} + \text{map}\}$ variables based on a quadratic approximation of the non-linear cost.

Let's gather the unknowns in one variable $\mathbf{X} = (\mathbf{x}_0, \mathbf{x}_1, \mathbf{x}_2, ..., \mathbf{x}_T)$. At each iteration $p+1$ of the optimizer, the previous estimate $\mathbf{X}^{(p)}$ is updated with:

$$\mathbf{X}^{(p+1)} = \mathbf{X}^{(p)} + \Delta\mathbf{X}^{(p+1)}. \qquad (6.5)$$

The increment $\Delta\mathbf{X}^{(p+1)}$ is found by minimization of the quadratic approximate of the cost criterion around $\mathbf{X}^{(p)}$. Based on the formulation of the cost in Eq. (6.4), the graph of the SLAM problem is no longer sparse due to the dependency between all the poses introduced by DAS beamforming. Thus, the determination of the increment may be computationally involving.

To reintroduce sparsity in the SLAM problem while leveraging DAS beamforming at the same time, we assume that the variation of the map induced by a "small" increment $\Delta\mathbf{X}^{(p+1)}$ in trajectory is locally negligible. Hence, only during calculation of the increment, the map can be assumed constant. This results in a sparse graph with the same structure as in Fig. 6.3, except that the map $\mathbf{M}$ is not to be inferred anymore. As $J$ is a sum of squared errors, under the assumption of constant map, the cost can be rewritten as:

$$J(\mathbf{X}) = \sum_{i,j} \mathbf{E}_{i,j}^{T}(\mathbf{x}_i, \mathbf{x}_j) \mathbf{\Sigma}_{i,j} \mathbf{E}_{i,j}(\mathbf{x}_i, \mathbf{x}_j) \triangleq \mathbf{E}^{T}(\mathbf{X}) \mathbf{\Sigma} \mathbf{E}(\mathbf{X}), \qquad (6.6)$$
where the $E_{i,j}$ are vector functions that assess how much the constraint relating $x_i$ and $x_j$ is satisfied. $E$ is a stacked vector containing all the $E_{i,j}$, and $\Sigma$ is a block-diagonal matrix with block elements $\Sigma_{i,j}$. The $\Sigma_{i,j}$ are matrix parameters of the cost function, which can be deduced from Eq. (6.4). Due to sparsity of the SLAM problem, there is only a few indices $(i, j)$ for which it is non-zero. In the case of an odometry constraint, $\Sigma_{i,j} \neq 0$ when $j = i + 1$, leading to:

$$E_{i,i+1}(\mathbf{x}_i, \mathbf{x}_{i+1}) = \mathbf{x}_{i+1} - f_{\text{dyn}}(\mathbf{x}_i, \mathbf{u}_i).$$

For a constraint related to an ultrasonic measurement, we have $\Sigma_{i,j} \neq 0$ when $i = j$, and:

$$E_{i,i}(\mathbf{x}_i, \mathbf{x}_i) = \hat{e}(\mathbf{x}, \hat{\mathbf{M}}, \mathbf{e}_i) - \mathbf{e}_i.$$

Note that there is no dependency in the map, as $\hat{\mathbf{M}}$ is assumed constant. To determine $\Delta\mathbf{X}^{(p+1)}$ at iteration $p+1$ of the NLS optimizer, the full vector function $E$ is linearized around the current trajectory candidate, and using a constant map $\hat{\mathbf{M}}^{(p)}$:

$$E(\mathbf{X}^{(p)} + \Delta\mathbf{X}) \approx E^{(p)} + \mathbf{J}^{(p)}\Delta\mathbf{X}, \qquad (6.7)$$

where $\mathbf{J}^{(p)} = \mathbf{J}(\mathbf{X}^{(p)})$ is the Jacobian matrix of $E$ in $\mathbf{X} = \mathbf{X}^{(p)}$, and $E^{(p)} = E(\mathbf{X}^{(p)})$ for notation simplicity. One should note that, due to the complexity of the cost, the Jacobian matrices $\mathbf{J}^{(p)}$ should be numerically approximated. Substituting Eq. (6.7) in Eq. (6.6) yields:

$$\begin{aligned} J(\mathbf{X}^{(p)} + \Delta\mathbf{X}) &\approx [\mathbf{E}^{(p)} + \mathbf{J}^{(p)}\Delta\mathbf{X}]^T \mathbf{\Sigma} [\mathbf{E}^{(p)} + \mathbf{J}^{(p)}\Delta\mathbf{X}], \\ &\approx (\mathbf{E}^{(p)})^T \mathbf{\Sigma}\mathbf{E}^{(p)} + 2(\mathbf{E}^{(p)})^T \mathbf{\Sigma}\mathbf{J}^{(p)}\Delta\mathbf{X} + \Delta\mathbf{X}^T (\mathbf{J}^{(p)})^T \mathbf{\Sigma}\mathbf{J}^{(p)}\Delta\mathbf{X}, \\ &\approx C + 2\mathbf{B}^T\Delta\mathbf{X} + \Delta\mathbf{X}^T\mathbf{H}\Delta\mathbf{X}, \end{aligned}$$

which is the quadratic approximation of the cost. The first-order optimality condition is:

$$\frac{\partial}{\partial \Delta \mathbf{X}} J(\mathbf{X}^{(p)} + \Delta \mathbf{X}) \approx 2\mathbf{B} + 2\mathbf{H}\Delta \mathbf{X} = 0.$$

Hence, the increment $\Delta\mathbf{X}^{(p+1)}$ is found by solving the linear system:

$$\mathbf{H}\Delta\mathbf{X} = -\mathbf{B}. \qquad (6.8)$$

Due to multiple zero components of $\Omega$, $\mathbf{H}$ is a sparse matrix. This property can be leveraged by conventional linear numerical solvers through Cholesky factorization of $\mathbf{H}$, for efficient resolution of Eq. (6.8) [35]. After application of Eq. (6.5) to update the trajectory, the map is updated with DAS beamforming: $\hat{\mathbf{M}}^{(p+1)} = \hat{\mathbf{M}}(\mathbf{X}^{(p+1)})$. Then, the iteration process is repeated until a termination criterion is fulfilled.

## 6.3 Experimental results

We assess the efficiency of our GraphSLAM approach on two sets of acoustic measurements acquired on our two rectangular metal panels. We use the dataset acquired on the 600x450x6mm aluminum plate for Scenario 1, and the UGWs measurements acquired on the 1700x1000x6mm steel plate for Scenario 2.
<figure><img src="image_66.png" /><figcaption>Figure 6.5: Likelihood maps of sensor presence obtained using measurements acquired on the aluminum plate (a), (b), and on the steel plate (c), (d). The true sensor position is indicated by the red mark. The outline of each plot accounts for the plate boundaries, and the color maps only spans the range from the minimum to the maximum likelihood value.</figcaption></figure>

### 6.3.1 Comparison of the observation models performance

To assess the benefit of using a multi-order echo observation model for localization, we first compare the localization accuracy achieved both with the proposed multi-order echo model, and the simpler model used in the previous chapters, that only relies on first-order reflections. We assess the likelihood of sensor presence at position x for each model with:

$$l_{1^{\text{st}}\text{order}}(\mathbf{x}|\mathbf{e}_i) = \sum_{(r,\alpha) \in \mathbf{M}} e_i(|x \cdot \cos \alpha + y \cdot \sin \alpha - r|)$$

for the first-order model, and:

$$l_{\text{multi-order}}(\mathbf{x}|\mathbf{e}_i) = -||\hat{\mathbf{e}}(\mathbf{x}, \mathbf{M}, \mathbf{e}_i) - \mathbf{e}_i||^2$$

for the multi-order echo model. Likelihood maps are determined for each model, using different envelope measurements $e_i$, and based on a known plate geometry M. A likelihood maximum is expected at the true sensor position. We show, in Fig. 6.5, several instances of likelihood maps obtained for each model, and that illustrate the pattern tendencies based on measurements acquired on the two metal panels. It can be noted that the likelihood maps present the same two axis symmetries than the plates due to the ambiguity of the echo measurements, which cannot be disambiguated without an external localization system.

We first qualitatively assess the performance of the observation models. Fig. 6.5a depicts a scenario where the two models point to the correct sensor position. In Fig.6.5b, only the first-order model does not point to the correct sensor position. This is due to interference, that can wrongly attribute one energetic reflection in the envelope signal to reflections from a different
<figure><img src="image_67.png" /><figcaption>Figure 6.6: Average localization errors calculated using the two considered models. The position error is averaged over all the measurements acquired on the aluminium plate (top) and steel plate (bottom) with respect to $\gamma$, yielding the solid lines. The upper and lower bounds of the colored areas are situated respectively above and below one standard deviation from the average.</figcaption></figure>

boundary, leading to incorrect localization. Furthermore, highly energetic high-order reflections also contribute to localization errors because the model treats them as first-order ones. With the multi-order echo model, however, the likelihood maximum can be identified at the ground truth sensor position. The same observations hold for Fig. 6.5c, which illustrates the likelihood maps for a measurement acquired on the large steel plate. This illustrates that the multi-order approach is better suited for localization. Fig. 6.5d illustrates the typical detrimental effect of overfitting with the multi-echo model, which occurs when one echo level is relatively high compared to the others in the envelope measurement. Due to the estimation of the coefficients $\mu_p$ in Eq. (6.1) by matching expected and measured envelope signals, the former will present substantial similarity with the measurement for any position for which the echo is expected (either as a first or a high-order echo). This results in higher localization ambiguity. Yet, this phenomenon is observable in the dataset for only a very limited number of measurements. We will see later if it detrimentally affects the GraphSLAM results.

To compare the model's performances quantitatively, we assess the average localization error calculated over the whole measurement sets, for each of the two metal plates. The initial likelihood maps are first re-scaled to values between 0 and 1, and are passed through an exponential mapping function using a multiplicative scaling parameter $\gamma$. Next, the average localization errors are assessed, after considering the resulting maps as (non-normalized) distributions over the sensor position on the plate, with:

$$ \Delta_i \mathbf{x}_{1^{\text{st}}\text{order}}(\gamma) = \frac{\sum_{\mathbf{x}} \|\mathbf{x} - \mathbf{x}_i^{\text{true}}\| \cdot \exp \gamma l_{1^{\text{st}}\text{order}}^{\text{scaled}}(\mathbf{x}|\mathbf{e}_i)}{\sum_{\mathbf{x}} \exp \gamma l_{1^{\text{st}}\text{order}}^{\text{scaled}}(\mathbf{x}|\mathbf{e}_i)}, $$

$$ \Delta_i \mathbf{x}_{\text{multi-order}}(\gamma) = \frac{\sum_{\mathbf{x}} \|\mathbf{x} - \mathbf{x}_i^{\text{true}}\| \cdot \exp \gamma l_{\text{multi-order}}^{\text{scaled}}(\mathbf{x}|\mathbf{e}_i)}{\sum_{\mathbf{x}} \exp \gamma l_{\text{multi-order}}^{\text{scaled}}(\mathbf{x}|\mathbf{e}_i)}, $$
<figure><img src="image_68.png" /><figcaption>Figure 6.7: Ground truth trajectories and dead-reckoning trajectories used for Scenario 1 (a) and Scenario 2 (b). The black rectangle represents the actual outline of the metal plate.</figcaption></figure>

where $\mathbf{x}_i^{\text{true}}$ is the ground truth sensor position, and $\gamma$ is a positive scaling parameter that balances the ambiguity level due to multiple likelihood maxima. Indeed, for a low $\gamma$, different likelihood maxima can contribute to the localization error, whereas the error reduces to that achieved by the maximum-a-posteriori position estimate when $\gamma$ goes to infinity.

Fig. 6.6 shows the evolution of the average localization errors and standard deviations calculated for the two models, and w.r.t. the scaling parameter $\gamma$. The errors are evaluated for both the data acquired on the aluminium plate and the data acquired on the steel plate. It can be observed that, in the two scenarios, and no matter the value of $\gamma$, the multi-order model outperforms the baseline model in terms of localization accuracy.

Overall, these results support the benefit of using our approach that can leverage multi-order echoes for localization, compared to first-order echo only models. The following subsection will illustrate that, when integrated within the proposed GraphSLAM approach, the multi-order model can naturally lead to more accurate results than those presented in Fig. 6.6.

### 6.3.2 Experimental localization and mapping results

We now assess the performance of the designed GraphSLAM approach on the experimental data for solving the full SLAM problem. We consider lawn-mower transects as simulated robot trajectories, and use the theoretic displacement between two measurements as noise-free odometry. The resulting linear and rotational displacements are corrupted with additive Gaussian noise: $\mathbf{u}_i = [\tilde{\Delta}_{r,i}, \tilde{\Delta}_{\theta,i}]^T$, with $\tilde{\Delta}_{r,i} \sim \mathcal{N}(\Delta_{r,i}, 0.006^2)$ and $\tilde{\Delta}_{\theta,i} \sim \mathcal{N}(\Delta_{\theta,i}, 0.04^2)$. Fig. 6.7 shows the ground truth trajectories used for the two considered scenarios, along with the dead-reckoning trajectory obtained based on the generated noisy odometry data. One may observe, especially for Scenario 1, that the drift induced by odometry noise is relatively more important than the drift in the experiments of Chapter 4.

To consider a baseline for comparison, we run the FastSLAM algorithm presented in Chapter 4, in the situation corresponding to Scenario 1. The localization and mapping results obtained at different steps of the simulated robot trajectory are presented in Fig. 6.8. One can observe the overall poor localization and mapping results. Due to its inability to leverage future measurements to adjust incorrect past position estimates (and particle sampling following an incorrect map), FastSLAM is unable to recover the trajectory and metal plate geometry in that scenario
<figure><img src="image_69.png" /><figcaption>Figure 6.8: Localization and mapping results obtained with FastSLAM for the Scenario 1. The figure represents the trajectories estimated by all the particles (red lines), along with the map retrieved by the particle with the highest weight (green lines) during Steps 1, 6, 18, 49, 88 and 108 of the simulated lawn-mowing trajectory. The true outline of the plate and true sensor positions correspond to the black rectangle and blue dot, respectively.</figcaption></figure>
<figure><img src="image_70.png" /><figcaption>Figure 6.9: Representation of the SLAM results during iteration Steps 1, 15, 29, and 41 of the NLS optimization process for Scenario 1. The true outline of the plate is represented by the black rectangle. For better comparison of the reconstructed trajectory with the ground truth, they are expressed in the estimated (resp. ground truth) plate frame, while the origin of the estimated plate geometry is fixed at [0, 0].</figcaption></figure>

where the odometry drift is severe.

We now run our GraphSLAM approach on the data of Scenario 1. We use the implementation of the NLS optimizer<fnref n="3" /> provided by the Scipy library, on Python. To leverage the sparsity of the GraphSLAM problem, we also provide the structure of the Jacobian matrix (i.e the positions of the non-zero components). This greatly improves the computation of the Jacobian matrix following the procedure in [20]. The matrix is numerically approximated using the default finite difference "2-points" method. The optimizer is initialized with the dead-reckoning trajectory, and the map is initialized with the plate geometry obtained with DAS beamforming. We then run the optimizer, considering a value of $\lambda = 0.03^2$ for the cost function parameter, and $\sigma = 2$ cm for the envelope signal model parameter. The iterative optimization terminates when either the relative change of the cost function, the relative change of the independent variables, or the change in the gradient vector become inferior to their corresponding (default) value.

With our implementation, the average execution time for one NLS iteration is 2.30 seconds. As several tens of iterations might be necessary to reach a satisfying solution, the approach

<footnote n="3">³Documentation available at: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.least_squares.html.</footnote>
<figure><img src="image_71.png" /><figcaption>Figure 6.10: Evolution of the loss (left), along with the localization and mapping errors (right) during the NLS optimization process for Scenario 1. The loss evaluated at ground truth is represented as a horizontal line.</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>Figure 6.11: Evolution of the loss (left), along with the localization and mapping errors (right) during the NLS optimization process for Scenario 2. The loss evaluated at ground truth is represented as a horizontal line.</figcaption></figure>

cannot be run online. Yet, one may expect a total execution time of 1 or 2 minutes, which appears, for now, to be reasonable for an offline SLAM approach. Fig. 6.9 illustrates the SLAM results at different steps of the optimization process. It can be observed in Fig. 6.9a that initially, the plate size estimated via DAS beamforming is accurate, but the localization error is seemingly important. Step by step, the estimated trajectory converges to the ground truth one. During Step 41, both the map and trajectory have been fully recovered.

In Fig. 6.10, we show the evolution, during 50 iterations of the NLS optimization process, of the mapping error, as well as the evolution of the localization error. The former is defined as the average error on plate dimensions, while the latter is evaluated as the mean distance between ground truth and estimated trajectories expressed in the ground truth and estimated plate frames, respectively. One may first observe that, after convergence of the optimizer, a loss value lower than the loss evaluated at ground truth is achieved, which illustrates the efficiency of our implementation of NLS. One may also note that the mapping error does not significantly vary, it eventually reaches an error of 1 cm. This is due to the efficiency of DAS beamforming for recovering the plate size despite important localization error. Besides, the localization error, which is 6 cm initially, is reduced to 1 cm after convergence.

To determine if our approach still performs well on measurements acquired on a larger plate surface, and for a different plate material, we perform the same experiment using the data acquired on the steel plate (Scenario 2). Fig. 6.11 shows the evolution of the loss along with the localization and mapping errors during the NLS optimization. It can still be observed that, after convergence, the loss is lower than that achieved at ground truth, which illustrates the
<figure><img src="image_73.png" /><figcaption>Figure 6.12: Representation of the SLAM results during iteration step of the NLS optimization process for Scenario 2. The true outline of the plate is represented by the black rectangle. For better comparison of the reconstructed trajectory with the ground truth, they are expressed in the estimated (resp. ground truth) plate frame, while the origin of the estimated plate geometry is fixed at [0, 0].</figcaption></figure>

efficiency of the optimization process. Furthermore, the average localization and mapping errors drop from 3 and 5 cm to around 2 cm. One may note that the decrease of the mapping error is more noticeable in this second scenario. Fig. 6.12 shows the SLAM results in this scenario, where one can observe the satisfying matching between estimated map and trajectory with the ground truths. Overall, these results demonstrate the potential of our GraphSLAM approach that leverages a multi-order echo model and a NLS optimizer to achieve SLAM with ultrasonic guided waves on a metal plate.

### 6.3.3 A study on the sensitivity of GraphSLAM to the loss parameter

In this section, we evaluate the sensitivity of the SLAM results to the cost function parameter. To this end, we run the NLS optimizer on the data from Scenario 1 for different values $\lambda \in [0.001^2, 0.005^2, 0.02^2, 0.05^2, 0.1^2]$. In Fig 6.13, we display the evolution of the difference between the loss achieved by the optimizer, and the loss at ground truth, along with the evolution of the localization and mapping errors, for each value of $\lambda$.

One may first notice that for $\lambda = 0.05^2$ and $\lambda = 0.1^2$, the cost difference is not negative after convergence, which results in an inefficient optimization. The optimizer, which iteratively approximates the nonlinear cost criterion, becomes stuck in a local minimum, resulting in poor localization results. A value of $\lambda = 0.02^2$ yields the best localization performance. For both $\lambda = 0.001^2$ and $\lambda = 0.005^2$, the cost difference is already negative at the start of the optimization process, which indicates that the confidence in the ultrasonic measurements is underestimated. In the first case, however, there is no sensible decrease of the localization error, whereas in the second one, the localization error is less than 1 cm, after convergence. There is little variation in the mapping error, which remains lower than 1 cm in all the scenarios. This can be explained by the robustness of DAS beamforming to localization error, especially when the trajectory is initialized with dead-reckoning path.

These observations illustrate the sensitivity of the NLS approach to the cost parameter value, which is critical to accurate SLAM results. As it may not be clear which value to choose
<figure><img src="image_74.png" /><figcaption>Figure 6.13: Evolution of loss along with the localization and mapping errors during the NLS optimization in Scenario 1, and for different cost parameters $\lambda$. The top left graph shows the evolution of the difference between the achieved loss and the loss at ground truth. A close view of the same graph is provided on the bottom left plot. The middle plot shows the evolution of the localization errors, while the mapping errors are provided on the right plot.</figcaption></figure>

beforehand, one may experiment with different $\lambda$, and determine a posteriori an interval that leads to a consistent and stable trajectory (likely to be close to $\lambda \in [0.01^2, 0.05^2]$ in our scenario).

### 6.3.4 Discussions

The approach we proposed to solve the full SLAM problem using ultrasonic measurements can accommodate multiple high-order reflections without performing explicit echo extraction and association, which are known to be challenging tasks.

Nevertheless, the proposed approach has several drawbacks: the proposed GraphSLAM approach can lead to an increase in performance when compared to FastSLAM. Yet, the experiments have illustrated its sensitivity to the parameter value of the cost criterion, which needs to be carefully tuned for accurate results. Besides, the uncertainty estimation on the SLAM result has not been investigated in the present study, whereas indicators on the confidence one can have on the solution reached by the optimizer can be crucial, especially in noisy environments. Also, the proposed approach is computationally involved. In the experiments, the sensor positions and map features are in relatively low numbers, but the method cannot be run online. In contrast, state-of-the-art GraphSLAM methods, which are computationally efficient, can cope with tens of thousands of poses and up to millions of map features. Finally, our approach, as presented here, is limited to plates with rectangular geometries. This is not an issue for an inspection operation on ship hulls or storage tanks, as they are almost entirely made of rectangular panels. However, the approach could be extended to work with more general geometries by relying on the more general mapping framework presented in the previous chapter.

To further assess the viability of our approach for the target application of long-range robotic inspection, the former shall be evaluated in more realistic conditions. Indeed, a lower signal-to-noise ratio is expected on an actual structure due to larger plates, background environmental noise, the low reflectivity of weld joints, or the presence of coating. Apart from application to robotic inspection, we envision that our approach can be used for other acoustic SLAM problems, especially in environments that present 3-dimensional reflective surfaces, through adaptation of the distance function and the map feature space.
## 6.4 Chapter summary

In this chapter, we elaborate a non-probabilistic graph-based framework to solve the full SLAM problem for an inspection robot relying on ultrasonic guided waves on a metal structure. We introduce an observation model that can leverage high-order echoes within the measurements for localization, and we rely on DAS beamforming to localize the plate boundaries. We integrate these elements within a custom implementation of a NLS framework to optimize a GraphSLAM cost criterion. Experimental results performed in a laboratory environment show that the designed observation model yields better localization accuracy than our initial first-order echo model, and it can yield more accurate results, close to sub-centimeter precision, when compared to FastSLAM, even in the presence of severe odometry drift.

In future works, the computational efficiency of the approach shall be improved to work with more significant trajectories and map features, and the uncertainty assessment on the optimizer result must be investigated. Also, it shall be tested in more realistic conditions to assess its viability for inspection on actual structures, in outdoor environments. Finally, other map features such as defects should be included in the map for a complete inspection task.
# Chapter 7

## Learning the propagation model for
Lamb wave-based mapping

The inspection of sizeable plate-based metal structures such as storage tanks or marine vessel hulls is a significant stake in the industry, which necessitates reliable and time-efficient solutions. Although Lamb waves have been identified as a promising solution for long-range non-destructive testing, and despite the substantial progress made in autonomous navigation and environment sensing, a Lamb-wave-based robotic system for extensive structure monitoring is still lacking. Following our previous work on ultrasonic Simultaneous Localization and Mapping, we introduce a method to achieve plate geometry inference without prior knowledge of the material propagation properties, which may be lacking during a practical inspection task in challenging and outdoor environments. Our approach combines focalization to adjust the propagation model parameters and DAS beamforming to infer the plate boundaries location by relying directly on acoustic measurements acquired along the mobile unit trajectory. For each candidate model, the focusing ability of the corresponding beamformer is assessed over high-pass filtered beamforming maps to further improve the robustness of the plate geometry estimates. We then recover the optimal space-domain beamformer through a simulated annealing optimization process. We evaluate our method on three sets of experimental data acquired in different conditions, and show that accurate plate geometry inference can be achieved without any prior propagation model. Finally, the results show that the optimal beamformer outperforms the beamformer resulting from the predetermined propagation model in non-nominal acquisition conditions.

### 7.1 Introduction and background

The development of an accurate ultrasound-based localization and mapping system is key to the emergence of a true long-range robotic inspection operation. Indeed, on a large metal structure, a moving robot, equipped with an omnidirectional emitter/receiver ultrasonic sensor, may be able to detect structure-bound landmarks, such as weld joints or stiffeners, and subsequently use them for accurate localization, within a SLAM framework. After determination of the structural elements along with the measurement positions, acoustic tomography could be achieved to identify small acoustic scatters that could be potential defects.

To make one step towards such a challenging application, we have first proposed a SLAM framework for localization and metal plate geometry reconstruction. The core of our methodology is DAS beamforming integrated in a FastSLAM algorithm (for the online estimation problem), or within a GraphSLAM framework (for the offline batch estimation problem). The experimental
results demonstrated the efficiency of the proposed methods on two different metal panels, in a laboratory environment. Yet, a key element is necessary for the successful achievement of ultrasonic SLAM: the prior knowledge of an accurate wave propagation model.

Prior knowledge of the physical properties of isotropic metal structures may be sufficient for acoustic localization and mapping in well controlled environments. However, the hypothesis that the propagation model is known a priori may not be realistic for a practical inspection task on a large metal structure due to a wide variety of external perturbations that can significantly affect the acoustic signals, and due to a lack of knowledge of the structure physical condition (which is the reason why it is being inspected). Potential perturbations may include, but may not be limited to, temperature variations, which can substantially affect the acoustic measurements [12], the pressure due to the neighboring metal plates which are welded altogether, moisture, the varying thickness of the coating, or the effect of the coupling (namely water in the case of acoustic inspection with a robotic system). Consequently, adaptive methods that automatically calibrate the propagation model by relying directly on data acquired on the field may be necessary to achieve accurate localization and mapping results, without any human intervention.

In the next, we propose an approach that combines acoustic focalization and DAS beam-
forming to map a plate-like structure by a mobile unit, so that accurate mapping results can be
achieved during an inspection task without prior knowledge of the propagation properties of the
material. The proposed approach combines focalization to adjust the parameters of propagation
models that are derived from the Helmholtz equation, and DAS beamforming to localize the plate
edges. The focalization capability of a candidate beamformer is assessed over high-pass filtered
beamforming maps, so that the detrimental effects of interference and high-order reflections are
lowered. A Simulated Annealing (SA) optimization process is then implemented so that the
optimal beamformer can be recovered based directly on data. Results acquired on experimental
data in three different scenarios show that the proposed approach is effective. Furthermore,
mapping the plate geometry with the optimal beamformer is found to be more efficient than
mapping with a predetermined propagation model in non-nominal acquisition conditions. For
the sake of simplicity, the localization problem is not addressed here, but the previous chapters
provide indications on how localization could be considered as well.

7.1.1 Related works

Recent works in the literature either address the problem of extracting the dispersion charac-
teristics of materials from ultrasonic data [104, 13, 18, 14], or investigate acoustic localization
techniques that can also recover the propagation properties [83, 56, 103]. However, simultaneous
acoustic localization and propagation characterization from pulse-echo ultrasonic measurements
acquired by a mobile unit has not been thoroughly studied. In the literature of beamforming,
adaptive methods have been elaborated to compensate for inaccurate knowledge of the prop-
erties of the propagation media for underwater acoustic source localization purposes [81], or
more generally, for Matched Field Processing [7]. Beamforming has been combined with focal-
ization [19], an approach that considers the propagation environment as an acoustic lens, and
which seeks, through an optimization process, to adjust the propagation parameters so that the
focalization capability of beamforming is maximized. It has also been studied for joint acoustic
signal separation and source localization [57]. The potential of focalization has been successfully
demonstrated in simulation for underwater acoustic localization, and may provide an interesting
approach for Lamb waves. Yet, this has not yet been demonstrated.
## 7.2 Proposed method

In what follows, we propose an approach for mapping the geometry of a metal panel without prior knowledge of the propagation model. We first briefly recall the propagation models based on approximate solutions of the Helmholtz equation. Next, we summarize the mapping strategy to localize the plate boundaries based on a known propagation model and DAS beamforming. These elements will be the basis of an optimal beamforming process, that can adjust the propagation model parameters through maximization of the focusing ability of beamforming.

### 7.2.1 Lamb wave-based mapping

Following our general setup, we are considering a mobile platform equipped with an emitter/receiver pair of piezoelectric transducers nearly co-located, and moving on a metal surface. All along a trajectory $x_{0:T}$ (assumed to be flawlessly known, as we are only interested in the mapping problem), an excitation signal $s$ is used to generate a Lamb wave in the material. Ultrasonic measurements $z_i(t)$, which contain the reflections on the plate boundaries, are subsequently collected. As for all the chapters, it will be assumed that the excitation signal is chosen adequately so that the A0 mode is predominant (while the propagation of S0 is negligible). The objective is to infer, by relying on an isotropic propagation model of the Lamb wave, the plate geometry and the sensors' acquisition positions in the plate frame.

Following the approach used so far, a wave propagation model is first determined based on prior knowledge of the plate material. The latter is recalled here for convenience. The wave dispersion relation $k(\omega)$ for the A0 mode is calculated by numerically solving the corresponding Rayleigh Lamb equation:

$$ \frac{\tanh qh}{\tanh ph} = -\frac{(k^2 - q^2)^2}{4k^2qp}, $$

where $h$ denotes the half-thickness of the plate, $k = 2\pi/\lambda$ is the wavenumber, $\lambda$ is the wavelength, and $p$ and $q$ are defined with:

$$ p^2 = \frac{\omega^2}{c_L^2} - k^2; \quad q^2 = \frac{\omega^2}{c_T^2} - k^2. $$

In the above equation, $c_L$ and $c_T$ are respectively the longitudinal and transverse velocities related to the plate material. Their numerical values can be found in lookup tables [96]. Next, based on the wave dispersion relation, a model of the transfer function is derived with:

$$ \hat{g}(r, \omega) \approx e^{-jk(\omega)r} / \sqrt{k(\omega)r}, \qquad (7.1) $$

where $r$ is the propagation distance. To determine a plate geometry estimate from all the measurements, we first determine the envelope signals $e_i$:

$$ e_i(r) = |z'_i(r) + j\mathcal{H}(z'_i(r))|, $$

where $z'_i$ is obtained through correlation of the measurement with the wave propagation model:

$$ z'_{i}(r) = \frac{\langle z_{i}(t), \hat{z}(r, t) \rangle}{\sqrt{\langle z_{i}(t), z_{i}(t) \rangle \langle \hat{z}(r, t), \hat{z}(r, t) \rangle}}, $$

where $\hat{z}(r, t) = \hat{g}(2r, t) * s(t)$ is the expected signal for the incident wave reflecting at a distance $r$ from the transducers. Next, an estimate of the map $\mathbf{M} = \{r_l, \alpha_l\}_{l=1..4}$ is obtained by applying
DAS beamforming to the envelope measurements:

$$ \mathcal{L}_T(r, \alpha) = \sum_{i=0}^{T} e_i(|x_i \cdot \cos \alpha + y_i \cdot \sin \alpha - r|), \qquad (7.2) $$

and subsequently solving the following optimization problem:

$$ \hat{\mathbf{M}} = \arg \max_{\mathbf{M}} \mathcal{L}_T(\mathbf{M}) = \arg \max_{\mathbf{M}} \sum_{l=1}^{4} \mathcal{L}_T(r_l, \alpha_l) \qquad (7.3) $$

where $\mathbf{M}$ is restricted to be a rectangle. We have proposed to solve the problem efficiently by taking that constraint into account. First, we determine the most likely line boundary with:

$$ (\hat{r}_1, \hat{\alpha}_1) = \arg \max_{r, \alpha} \mathcal{L}_T(r, \alpha). $$

Next, assuming that $\hat{\alpha}_1$ provides the most reliable estimation of the plate orientation w.r.t. the robot, the determination of the other lines reduces to solving independent and straightforward one-dimensional optimization problems:

$$ \hat{\alpha}_l = \hat{\alpha}_1 + \frac{\pi(l-1)}{2} ; \quad \hat{r}_l = \arg \max_{r} \mathcal{L}_T(r, \hat{\alpha}_l) $$

for $l = 2, 3, 4$. The experiments from the previous chapters have shown that this approach can give accurate results in laboratory conditions. Yet, one of the major limitations of this approach is that it relies on prior knowledge of an accurate propagation model $\hat{g}$ to obtain precise results. However, such a hypothesis may not be realistic for a practical inspection task in challenging outdoor environments, where the structure state is truly unknown.

### 7.2.2 Optimal beamforming for model learning

In this section, we present an adaptive method to recover a metal plate geometry without assuming prior knowledge of the propagation model. It is based on focalization in the parameter space (i.e. the propagation model is adjusted) and beamforming for localization of the plate boundaries. First, the parametrization of the propagation models using solutions of the Helmholtz equation is presented. Next, we introduce and apply a simple high-pass filter to the beamforming maps to limit the detrimental effect of interference and high-order reflections. A loss function is then designed to assess the focusing capability of a candidate beamformer which should maximize spatial coherency (i.e. the energy that is focused at the geometry estimate on the beamforming map) in the case when the propagation model is appropriate. Finally, an optimizer based on simulated annealing is presented to recover optimal propagation parameters by minimizing the loss in a limited number of iterations. The efficiency of such an approach has been successfully demonstrated for underwater source localization purposes [19, 57].

Based on the hypothesis of linear and isotropic propagation, we aim at recovering both a propagation model $\hat{g}(r, t)$ and the plate geometry expressed as a set of lines $\hat{\mathbf{M}} = \{\hat{r}_l, \hat{\alpha}_l\}_{l=1..4}$. For candidate propagation models, we keep relying on the (approximate) solutions of the Helmholtz equation. Hence, estimating the propagation model reduces to the estimation of only three parameters which are the longitudinal velocity $c_L$, the transverse velocity $c_T$ and the plate thickness $d$, and that we will gather in the variable $\Omega = \{c_L, c_T, d\}$. This choice is convenient as we constrain the propagation models to be physically plausible while limiting the search space
for the propagation model parameters to a low dimension, which will facilitate the optimization process. As the complete state of the structure may be unknown during an inspection operation, the calibration of these wave parameters may be relevant. Indeed, the plate thickness may not be known accurately all over the surface, as the structural integrity is unknown. Besides, the effects of variation of temperature, pressure, moisture may be interpreted as variations of the effective velocities. For more complex sources of disturbances such as the effect of coupling or structure irregularities, we are not seeking to strictly compensate for their effect with such a simple model, yet sufficiently enough to predict the plate boundaries location accurately. Next, given candidate values Ω for the model parameters, the dispersion equation for the A0 mode is numerically solved, and is used to infer the propagation model $\hat{g}_{\Omega}(r, t)$ with Eq. (7.1).

Compared to plate geometry reconstruction with a known propagation model, inferring the plate geometry and propagation parameters simultaneously requires the determination of a larger number of unknowns. Consequently, an appropriate loss function that ideally prevents irregular cost surfaces with many local minima is needed to facilitate the optimization process.

Here, we introduce a loss function to assess the ability of a candidate beamformer (related to candidate parameter values Ω) to focus the energy of the wave packets contained in the measurements at the plate geometry estimate, so that minimizing the loss function w.r.t. the model parameters improves the spatial coherency achieved with the beamformer. Given candidate parameters Ω, the beamforming map $\mathcal{L}_{\Omega}(r, \alpha)$ is constructed as in Eq. (7.2), allowing a plate geometry estimate $\hat{M}_{\Omega}$ to be retrieved from it with the same optimization process as in Eq. (7.3). We then evaluate the total energy focused at the estimate $\hat{M}_{\Omega}$ over the beamforming map that has been high-pass filtered to limit the effect of interference and high-order reflections. The filtered map value at each line $(r, \alpha)$ is simply defined by the difference between its initial energy value and the minimum of energy in its vicinity that is to be appropriately defined:

$$ \tilde{\mathcal{L}}_{\Omega}(r, \alpha) = \mathcal{L}_{\Omega}(r, \alpha) - \min_{(r', \alpha') \in V(r, \alpha)} \mathcal{L}_{\Omega}(r', \alpha'). \quad (7.4) $$

where $V(r, \alpha)$ refers to the set of lines in the vicinity of $(r, \alpha)$. Compared to the standard $\mathcal{L}_{\Omega}(r, \alpha)$, the filtered map $\tilde{\mathcal{L}}_{\Omega}$, while being inexpensive to compute, has the advantage to present fewer areas where the intensity is high but homogeneous (such a situation occurs due to the combination of high-order reflections and interference). Hence, it is more compatible with Eq. (7.3) for the determination of the location of the edges. Thus, using Eq. (7.4) may appropriately filter "fuzzy" areas where the energy is spread homogeneously (i.e. not focused at a single point) due to interference on the beamforming maps, and may isolate correct intensity peaks. This will be illustrated next with experimental data.

To simultaneously recover propagation parameters and the plate geometry, performing a joint search in both the propagation model space and the geometry space would be computationally expansive. Instead, we rely on an optimal beamforming formulation: the loss is only evaluated over candidate model parameters Ω, and for the corresponding geometry estimate $\hat{M}_{\Omega}$ retrieved from the high-pass filtered beamforming map $\tilde{\mathcal{L}}_{\Omega}$ with Eq. (7.3). The loss value $l(\Omega)$ that we will seek to minimize is then defined as the opposite of the sum of the intensity levels evaluated at the retrieved edges:

$$ l(\Omega) = - \sum_{(r, \alpha) \in \hat{M}_{\Omega}} \tilde{\mathcal{L}}_{\Omega}(r, \alpha), \quad (7.5) $$

so that the energy focused at the plate geometry estimate can be maximized, while the detrimental effects of interference and high-order reflections can be minimized. The geometry estimate provided by the optimal beamformer is expected to match closely the ground truth geometry
due to the maximum of spatial coherency, as long as measurements have been acquired on a sufficient portion of the plate surface. Besides, to highlight the benefit of using high-pass filtered beamforming maps, propagation model selection based on the loss evaluated using the regular beamforming map $\mathcal{L}_{\Omega}$ will also be carried out in the following sections, and the mapping results will be compared.

### 7.2.3 Optimization with simulated annealing

An optimization process is needed to adjust the wave propagation parameters (longitudinal and transverse velocities and the plate thickness) through the minimization of the loss function. Due to the implicit definition of the propagation models, first-order methods such as gradient descent are impractical to use for our problem.

Simulated annealing is efficient for estimating a global minimum of a cost function, even when the number of unknowns is large [57]. This simple metaheuristic has been successfully used in Chapter 5 for recovering polyogonal shapes, and it is a good candidate optimizer for the new problem at hand. For Lamb-wave based mapping without a prior propagation model, the parameter search is reduced to the optimization over the value of $\Omega$, as the geometry is directly deduced using Eq. (7.3) when $\Omega$ is fixed. Our implementation of the optimizer is similar to that used in articles related to the underwater localization problem [19]. The major difference lies in the cost function design. First, we initialize the propagation parameters with uniform distributions within predefined intervals:

$$ 
\begin{aligned}
c_{L}^{(0)} &\sim \mathcal{U}([c_{\min}, c_{\max}]); \quad c_{T}^{(0)} \sim \mathcal{U}([c_{\min}, c_{L}^{(0)} - m]) \\
d^{(0)} &\sim \mathcal{U}([d_{\min}, d_{\max}]).
\end{aligned}
$$

The value of $c_T$ is drawn below the value of $c_L$, with a small margin $m > 0$, as it cannot physically be higher (i.e there are no positive real-valued solutions to the Rayleigh-Lamb equation). During each iteration of the algorithm, the parameters are randomly disturbed using the following perturbations:

$$ 
\begin{aligned}
c_{L}^{(p+1)} &= \min \left\{ \max\{c_{L}^{(p)} + \Delta c \cdot \chi_{0}^{3}, c_{\min}\}, c_{\max} \right\} \\
c_{T}^{(p+1)} &= \min \left\{ \max\{c_{T}^{(p)} + \Delta c \cdot \chi_{1}^{3}, c_{\min}\}, c_{L}^{(p+1)} - m \right\} \\
d^{(p+1)} &= \min \left\{ \max\{d^{(p)} + \Delta d \cdot \chi_{2}^{3}, d_{\min}\}, d_{\max} \right\}.
\end{aligned}
$$

Using the min-max formulation enables to leverage prior information by constraining the parameter values within plausible intervals. Similarly to the initialization step, $c_T^{(p+1)}$ is restricted to be lower than $c_L^{(p+1)}$ with the same margin $m$. $\chi_0, \chi_1, \chi_2$ are independent random values between -1 and 1, and are drawn from uniform distributions: $\chi_{1,2,3} \sim \mathcal{U}([-1, 1])$. They are raised to power 3 so that small variations are more likely. Still, large variations (i.e in the order of $\Delta c$ for the velocities, and $\Delta d$ for the thickness) can be possible to explore regions of the search space that are far occasionally from the current $\Omega^{(p)}$.

Next, for the new parameter values $\Omega^{(p+1)}$, the dispersion relation $k_{\Omega^{(p+1)}}(\omega)$ is determined by solving the Raleigh-Lamb equation. As the relation is computed numerically at each iteration of the optimization process, the computational load can be quite demanding. A trick to save computational time is to determine the wavenumber only for a restricted set of pulsation values $\omega_1, ..., \omega_K$ and to use linear interpolation to determine the wavenumber for other frequencies. The acoustic transfer function $\hat{g}_{\Omega^{(p+1)}}(r, \omega)$ is subsequently determined with the Helmholtz model given in Eq. (7.1). The beamforming map $\mathcal{L}_{\Omega^{(p+1)}}$ is then computed, and a plate geometry
estimate $\hat{M}_{\Omega(p+1)}$ is retrieved from it, allowing the determination of the loss value $l(\Omega^{(p+1)})$ with Eq. (7.5).

The variation of energy between the previous and disturbed parameters is $\Delta E = l(\Omega^{(p+1)}) - l(\Omega^{(p)})$. In the context of simulated annealing, the disturbed parameter value is not systematically kept. It is the case only when:

$$ \chi < \exp \left\{ -\frac{\Delta E}{\gamma T(p)} \right\} $$

where $\gamma$ is a strictly positive scaling parameter, $T(p)$ is the temperature at the iteration $p$, and $\chi$ is drawn randomly and uniformly between 0 and 1. Consequently, when the energy is decreased, the disturbed parameters are systematically kept. Otherwise, the acceptance rate is given by the Boltzmann distribution, which yields lower acceptance rates for more significant increases of energy. The temperature parameter is often chosen to decrease inverse logarithmically. To enable fast convergence of simulated annealing, we will decrease the temperature inverse linearly. Trials and errors are used to determine appropriate parameters for the optimizer: the values that empirically demonstrate a lower likelihood for the optimizer to be stuck in local maxima while maintaining a sufficient convergence speed are retained.

As it has already been highlighted, the cost function may be heavy to assess, in particular when the number of considered measurements is high. Indeed, each iteration of the optimizer requires solving the Rayleigh-Lamb equation for a set of frequency values, then computing the beamforming map using all the measurements, and recovering a plate geometry estimate with the optimization. For a robotic application, a few seconds might be needed to achieve convergence with around 100 measurements and using a linear decrease of the temperature. Thus, our method cannot be used in real time. Yet, this is not an issue, as it could be considered, during a practical robotic inspection task, to stop the robot for a few seconds to run the optimizer occasionally, and restrict the maximum number of acoustic measurements used to perform the optimization.

## 7.3 Experimental results

This section illustrates the efficiency of our approach for mapping a rectangular metal panel using ultrasonic guided waves in the three aforementioned scenarios. The benefit of filtering the beamforming maps is first discussed and highlighted based on the experimental data. Next, the correlation between the designed loss function and the reconstruction error is numerically assessed. The optimizer designed in the previous section is run to simultaneously infer the plate geometry and recover the propagation model. The resulting precision is compared with that obtained when using a predetermined propagation model that will be used as a baseline. The overall results demonstrate that this target objective is successfully achieved with our method in nominal acquisition conditions (Scenarios 1 and 2). Our approach is also tested in slightly disturbed conditions (Scenario 3) to illustrate that it remains efficient. For all the situations, the results are obtained using the full batch of measurements, which would amount to having the robot covering an important portion of the plate surface before calibrating the propagation model. As this may not be a representative scenario, the efficiency of our procedure is also assessed with simulated robot paths with acquisition points that are more sparse, and with a varying number of measurements available during the calibration.

### 7.3.1 Experimental scenarios

We test our approach on experimental acoustic data that have been acquired in three different scenarios. To generate the datasets, an emitter/receiver pair of nearly collocated contact
<figure><img src="image_75.png" /><figcaption>Figure 7.1: Picture showing the experimental setup for Scenario 3.</figcaption></figure>

piezoelectric transducers is moved to different positions on the plate surface, in each scenario.

For Scenario 1, we use the data acquired manually on our 600x450x6 mm aluminum plate. For Scenario 2, we use the dataset collected on our 1700x1000x6 mm steel plate. For Scenarios 1 and 2, the transducers were in contact with the surface during the acquisition. Coupling gel is placed at their interface to ensure a good coupling. For Scenario 3, the data are acquired on the same plate as the one used for Scenario 1. However, a thin layer of water of approximately 1 mm of thickness is put all over the plate surface, and the transducers are not placed in contact with the plate surface during the acquisition, but are kept in contact with the water only. This acquisition is performed using a customer-design five axes immersion scanner made by Inspection Technology Europe BV. It is also used to place the transducers (that are kept one beside the other) at specific positions on the measurement grid while maintaining a constant distance between the plate surface and the transducers. Due to the plastic holder, the distance between the centers of the transducers is 2.5 cm, which may not be an issue for our approach that relies on an exact pulse-echo setup, as this distance is not large comparing to the signal wavelength. For this scenario, a total of 108 measurements is collected at the same positions as those used for Scenario 1. The sampling frequency is 1.25 MHz, and the number of samples per measurement is 500. The data acquired in this setup are expected to be representative of those that would be acquired on a real structure by the robotic system, where water may also be used as coupling, and where the transducers may not be directly in contact with the structure surface to avoid damaging them by surface irregularities. Furthermore, this scenario will also allow us to highlight the potential of our approach in a slightly "disturbed scenario" comparing to the nominal acquisition conditions of Scenarios 1 and 2. A picture of the experimental setup used for Scenario 3 is available in Fig. 7.1.

### 7.3.2 Influence of the high-pass filter

First, we illustrate the benefit of applying a high-pass filter to the beamforming maps for more accurate evaluation of the focusing ability of a beamformer. We apply DAS beamforming on all the measurements acquired on the larger steel plate (Scenario 2), and retrieve a map estimate using Eq. (7.3). The propagation model parameters $d = 6$ mm, $c_L = 5880$ m/s and $c_T = 3250$ m/s are used based on prior knowledge of the plate material. The DAS beamforming results are provided in Fig. 7.2a, which depicts two visible intensity peaks. The two other peaks are much
<figure><img src="image_76.png" /><figcaption>Figure 7.2: Different beamforming maps computed using the data from Scenario 2 and the retrieved edges. a) shows the standard beamforming map, which yields correct estimates. b) depicts the map obtained with the same data sub-sampled in time. One edge is not correctly estimated. c) shows the high-pass filtered beamforming map obtained from the sub-sampled data. The correct plate geometry is recovered, while the fuzzy areas due to interference and higher-order reflections have been partially filtered out.</figcaption></figure>

less visible due to their lower intensity, and can be easily mistaken with interference and high-
order reflections, which cause areas of homogeneous intensity. Although the estimation is correct
in this case (the estimated plate dimension is 1720x980mm), it may not be robust. To illustrate
this, we sub-sample in time the data by a rate of 4 and reconstruct the beamforming map using
the same propagation model. The results can be seen in Fig. 7.2b. One may observe that the
resulting beamforming map is very similar. However, all but one line boundary are still correctly
recovered, causing the plate dimension to be underestimated. This could be detrimental to our
final objective, as the resulting map estimate, obtained by maximization of "acoustic focusing"
as we defined it, is incorrect.

Using high-pass filtered beamforming maps may increase the robustness of the geometry
estimates, and thus, may result in better definition of the focusing ability of DAS beamforming
for model calibration. The previous beamforming, after filtering, is shown in Fig. 7.2c. To
set up the high pass filter, $V(r, \alpha)$ is chosen as the rectangular area centered at $(r, \alpha)$ with a
size $6.5\text{cm} \times 24^\circ$. This window is sufficiently large to encompass the intensity peaks on the
beamforming map, and sufficiently limited to make the filtering effective. When now relying on
the filtered version of the beamforming map, the map estimation is made correct, even though
the detrimental effect of interference and higher-order reflections cannot be completely alleviated.
Hence, the introduction of the high-pass filter may provide additional robustness, which is to
be leveraged to calibrate the propagation model. It may be even more useful in the presence of
disturbances during the acquisition, as will be illustrated in the following.

7.3.3 Correlation between the loss and the mapping accuracy

We assess whether the designed loss is sufficiently correlated with the reconstruction error to
ensure that its minimization yields an accurate plate geometry estimate. We create a family
of 40 propagation models based on the approximate solutions to the Helmholtz equation, with
values of $\Omega = \{c_T, c_L, d\}$ uniformly distributed within the intervals $c_T, c_L \in [2500\text{m/s}, 6500\text{m/s}]$
and $d \in [3\text{mm}, 7\text{mm}]$. Then, we determine the filtered beamforming map using each model to
obtain a plate geometry estimate in each of the three scenarios under consideration. The results
are displayed in Fig. 7.3 as plots of the reconstruction error (both in range and orientation).
<figure><img src="image_77.png" /><figcaption>a) Scenario 1 with filtering</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>b) Scenario 1 without filtering</figcaption></figure>

<figure><img src="image_79.png" /><figcaption>c) Scenario 2 with filtering</figcaption></figure>

<figure><img src="image_80.png" /><figcaption>d) Scenario 2 without filtering</figcaption></figure>

<figure><img src="image_81.png" /><figcaption>e) Scenario 3 with filtering</figcaption></figure>

<figure><img src="image_82.png" /><figcaption>f) Scenario 3 without filtering</figcaption></figure>

Figure 7.3: 2D plots showing the correlation and loss values achieved using each of the propagation models from the predefined model set, and for the three scenarios. On the left, the loss value and geometry estimates are determined from the high-pass filtered beamforming maps. On the right, they are determined from the standard maps without filtering.
w.r.t. the loss value. We also display the same plots using the geometry estimates and loss values obtained without the filtering, for comparison. Ideally, we would obtain reconstruction errors that would monotonically decrease with the loss value.

For both Scenarios 1 and 2 in Fig. 7.3a-d, we observe that the range errors globally diminish for lower loss values. It can also be observed that models with close loss values can yield completely different reconstruction errors, or that models with high loss values can provide an accurate estimate (e.g the plot of the range error in Fig. 7.3c). Yet, these tendencies are not critical, as the major requirement for our approach to be effective is that the best estimation results are achieved for the lower losses as this is what will be minimized. The variations of the orientation error present a different aspect. There is no strict decrease in the error for lower losses. Also, the discretization of the orientation value is visible due to the resolution of the beamforming map that is limited for the sake of computational efficiency. Yet, it is still possible to achieve an accuracy in the order of one degree, which is sufficient for our application. Also, the limited resolution is not a significant problem for propagation model learning, as the best accuracies are also achieved for the lower loss values.

We notice a slightly more divergent behavior for the two losses in Scenario 3, which corresponds to our slightly disturbed scenario. Regarding Fig. 7.3f, where the beamforming maps are not filtered, the reconstruction error presents a tendency to decrease for lower loss values. However, the minimum reconstruction error is achieved for a loss value of approximately $l(\Omega_{\min}) = -86$, and it becomes higher for lower losses. This is not the desired behavior, as it would result in poorer estimation results if the loss minimization is effective. This tendency is not visible when the loss values and the geometry estimates are obtained from the high-pass filtered beamforming maps (Fig. 7.3e), where the lowest reconstruction errors are globally obtained for the lower losses. Hence, these results support the benefit of high-pass filtering the DAS beamforming maps for mapping and propagation model learning. It provides additional robustness due to its ability to lower the detrimental effect of high-order reflection and interference.

### 7.3.4 Mapping results with model learning in nominal conditions

To assess the performance of our approach for optimal beamforming for Lamb wave-based mapping, we run it 10 times for Scenario 1 and 2, with random initial parameter values for each repetition. The simulated annealing is set with the following parameters: $\Delta c = 7000\text{m/s}$, $\Delta d = 2\text{cm}$, $c_{\min} = 1500\text{m/s}$, $c_{\max} = 8000\text{m/s}$, $m = 500\text{m/s}$, $\gamma = 1$, $T = T_0/t$ where $T_0 = 200$, $d_{\min} = 3\text{mm}$ and $d_{\max} = 1\text{cm}$. Besides, the performance of our approach is compared with the results achieved when using predetermined propagation models based on prior knowledge of the plate thickness, the longitudinal and transversal velocities in nominal conditions. These predetermined models will serve as a baseline for comparison. Fig. 7.4 depicts the evolution of the minimum loss value and the corresponding reconstruction errors achieved during the optimization process. The upper and lower bounds of the filled areas represent the minimum and maximum values over the 10 runs, at each iteration step, while the solid blue lines account for the mean values. For comparison, the loss values and reconstruction errors achieved when using the predetermined propagation models for the two scenarios are also displayed.

In both scenarios, the achieved minimum loss decreases monotonically until it reaches a plateau. The final loss value is always lower than the loss value achieved with the predetermined model in Scenario 1 and significantly lower in Scenario 2. In terms of reconstruction error, the final average range error is comparable to that achieved with the predetermined model, which is 3mm for Scenario 1. It is slightly above the estimation error value obtained with the predetermined model (0.5cm) for Scenario 2. Besides, the orientation errors rapidly decrease
<figure><img src="image_83.png" /><figcaption>Figure 7.4: Evolution of the loss value and reconstruction errors for 10 repetitions of the optimization process with simulated annealing, and using the data from Scenario 1 (a) and Scenario 2 (b). The solid lines represent the mean values. The upper and lower bounds of the filled areas represent the minimum and maximum values respectively during each iteration. For comparison, the values achieved with the predetermined propagation model for each scenario are displayed as horizontal lines. The scale along the y-axis is logarithmic for the range error plots only.</figcaption></figure>

to zero for all the runs, achieving the same precision as the predetermined models in the two scenarios. Altogether, these experiments demonstrate that a propagation model can be efficiently recovered to estimate a plate geometry through optimal beamforming by relying directly on data, regardless of the plate size and material, and as long as the measurements have been acquired on a sufficient portion of the surface.

### 7.3.5 Mapping results with model learning in non-nominal conditions

We perform the same experiments using the data from Scenario 3 which have been acquired in slightly disturbed conditions due to the layer of water that is placed at the interface between the transducers and the plate surface. The evolution of the loss and the reconstruction errors for 10 repetitions of the optimization process with random initial parameter values are shown in Fig. 7.5. It can be observed that the loss rapidly decreases to a value lower than that achieved with the predetermined model. Regarding the reconstruction error, the average range error, after convergence, reaches 1cm while the error is 2cm when using the predetermined model for the geometry reconstruction. Furthermore, the orientation error rapidly decreases to a value lower than one degree.

Altogether, our approach, which automatically calibrates the propagation model based on data, yields better geometry estimates in this disturbed scenario compared to mapping with the predetermined propagation model, as illustrated by Fig. 7.5. These results are promising as the propagation in the thin layer of water has not been modeled explicitly. This further shows the potential of our approach which adapts the propagation model to achieve accurate plate reconstruction in varying conditions, as it is expected to happen for a robotic inspection
<figure><img src="image_84.png" /><figcaption>Figure 7.5: Evolution of the loss value and reconstruction errors for 10 repetitions of the optimization process and using the data from Scenario 3. The blue lines represent the mean values. The upper and lower bounds of the colored areas represent the minimum and maximum values respectively during each iteration. For comparison, the values achieved with the predetermined propagation model are displayed as horizontal red lines. The scale along the y-axis is logarithmic for the range error plot only.</figcaption></figure>

task on large structures in challenging outdoor environments. Also, with our approach, there is potentially no need to take into account slight sources of disturbance explicitly in the propagation models.

### 7.3.6 Evaluation of our approach with sparse measurements

During a real robotic inspection task, the acoustic measurements may not have been acquired over a dense grid, and over a sufficient portion of the surface during the calibration operation to fully recover the plate geometry. To evaluate the performance of our approach in more realistic scenarios, we assess how the quantity of measurements available during the optimization process affects the propagation model selection when the data are acquired along more realistic trajectories. To avoid running the optimizer every time, we consider the previous family of 40 propagation models. We evaluate, for each propagation model, the reconstruction errors obtained with and without filtering the beamforming maps, and select the error relative to the propagation model, for every set of measurements, that yields the lowest loss value. We also assess the estimation errors when the predetermined propagation model is used along with high-pass filtered beamforming maps, for comparison. We design two different trajectories (one for each plate) that are provided in Fig.7.6. The results, which are shown in Fig. 7.7, were obtained using the data acquired on the large steel plate (Scenario 2) and using the data acquired in the disturbed conditions (Scenario 3).

When only a few measurements are available, the average range errors are relatively high in
<figure><img src="image_85.png" /><figcaption>Figure 7.6: Simulated robot paths used for the experiments on the steel plate (left) and the aluminum plate (right). The path lengths (i.e. the number of acquisition positions) are respectively 36 and 26.</figcaption></figure>

<figure><img src="image_86.png" /><figcaption>Figure 7.7: Reconstruction errors w.r.t. the number of measurements considered in the loss, and using the data from Scenario 2 (a) or from Scenario 3 (b). The measurements are integrated one after the other in the loss, following the simulated paths. The error is only evaluated for the model yielding the lower loss value for the two cases (model selection with and without filtering the beamforming maps). For comparison, the errors achieved with the predetermined models are also displayed.</figcaption></figure>
<figure><img src="image_87.png" /><figcaption>Figure 7.8: Loss value computed on the high-pass filtered beamforming maps with respect to the longitudinal and transverse velocities. The data from Scenario 1 are used for the plot on the left, while the data from Scenario 2 are used for the plot on the right.</figcaption></figure>

all the cases, because the sensors need to pass by a border closely enough to detect it. The most accurate results are achieved when all the measurements are taken into account. In the middle, we can see, in Fig. 7.7a, that the minimum range error is reached faster when the beamforming maps have been filtered, and this minimum error is lower than the error achieved with the predetermined model. Furthermore, both in Fig. 7.7a and Fig. 7.7b, the reconstruction errors are seemingly lower when relying on model calibration for the reconstruction. Overall, the results show that our approach is effective for mapping, and yields similar - if not better - performance than that obtained with a predetermined propagation model, even when the measurements have been sparsely acquired by the mobile unit.

### 7.3.7 Discussions

Although we have shown that a propagation model can be recovered based on data to achieve precise localization by adapting the model parameters, our method cannot recover the real values of the physical model parameters. Indeed, different model parameters can yield equivalent loss values as illustrated by Fig. 7.8, where we represent the variation of the loss w.r.t. the two velocity parameters for a fixed thickness (d = 6mm), and using the data acquired on the two metal plates considered in this study (Scenario 1 and 2). The more specific explanation is that different model parameters may lead to similar dispersion values in the considered frequency bandwidth, after resolving the Rayleigh-Lamb equations. Yet, the fact that the actual model parameters cannot be recovered with our method is not a significant issue, as our primary objective is to achieve accurate mapping of the plate structure with a calibrated model *g*, so that the wavepackets inside the measurements can be appropriately accounted for.

The method presented here restricts the propagation models to be approximate solutions of the Helmholtz equation with dispersion relations derived from the Rayleigh-Lamb equation. This approach is not sufficient for recovering the actual model parameters, and addressing this in future work would be desirable for a complete NDT method (the knowledge of plate thickness variation would be of particular importance). Yet, the recovered propagation model $\hat{g}$ is sufficient to achieve accurate mapping results through beamforming. Whether this approach remains appropriate in more realistic conditions (i.e on a real ship hull for example) is still subject to investigation. It is expected that this parametrization would be sufficient as long as the
first-order wave-packets are sufficiently energetic within the measurements, and that the major hypotheses on the propagation model (linear propagation, homogeneous and isotropic material) approximately hold. Due to the large size of the plates on a real structure and the low reflectivity of weld joints, the signal-to-noise ratio (SNR) is expected to be low. Hence, a study evaluating the performance of our approach for various SNR conditions is needed to assess how likely it would work in practice. Also, it is to be noted that, in the present study, the sensor positions were measured accurately, whereas in practice, only estimated positions will be available.

The case of multi-modal propagation is not considered here, whereas it is likely if the frequency is not sufficiently adequate to the material in a -presumably- unknown state, or if mode conversion occurs. Having at least A0 and S0 modes propagating simultaneously is the most frequent scenario. It is believed that the algorithm could be extended by considering several hypotheses (bi-modal, A0-only, S0-only propagation, mode conversion for different paths, see [108]) to determine which one is most likely based on data. However, integrating more complex interactions such as diffraction due, for example, to complex structural features such as stiffeners, holes... would be more challenging. We expect isotropic propagation to be prevalent, as these complex wave interactions may be scarce and have a sufficiently small incidence on the signals to not affect the mapping results.

Besides, our method, as presented here, is restricted to rectangular geometries, as it facilitates edge retrieval from the beamforming map. This constraint is not a limitation for mapping storage tanks or ship hulls, as they are almost entirely made of rectangular panels. Yet, our approach could be extended, following the framework presented in Chapter 5, to make it applicable to structures with more general geometries. One may also want to adapt our approach to more conventional applications in SHM where, for example, the propagation model could be automatically calibrated by maximizing the energy focused at the estimated defect location on the imaging results, or by relying on the reflections on the sample boundaries.

## 7.4 Chapter summary

This chapter introduces a method to accurately recover the geometry of a metal plate by relying on ultrasonic measurements acquired by a mobile unit, in a pulse-echo mode, and without using a predetermined propagation model. Our approach is based on focalization in the model parameter space and DAS beamforming for localization of the plate boundaries. We restrict the propagation model to be an approximate solution of the Helmholtz equation and parameterize it with only three physical values. We introduce a loss function assessed on high-pass filtered beamforming maps to quantify the focusing ability of a candidate beamformer (i.e., candidate model parameters). Eventually, we find optimal model parameters with a simulated annealing optimization process.

We demonstrate the performance of our method on three sets of experimental acoustic data acquired on a dense grid on two metal plates of different sizes and different materials. The results illustrate the relevance of filtering the beamforming maps to reduce the detrimental effect of interference and high-order reflections in the assessment of the focusing ability of beamforming. We also show that a propagation model enabling accurate boundary localization can be recovered with simulated annealing. For the two undisturbed scenarios, the precision of the localization is found to be similar to that achieved with the model built from prior knowledge on the plate material, but it is found to be superior in non-nominal acquisition conditions. This highlights the very potential of our method for Lamb wave-based localization and mapping on a large metal structure, where the wave propagation conditions may not be known a priori. The
benefit of this approach is real for practical industrial inspection tasks, where the propagation
models could be automatically calibrated. Eventually, we assess the performance of our approach
using reduced numbers of measurements acquired on the two plates, and following more realistic
robot trajectories. The results illustrate that the proposed approach can recover the plate ge-
ometry accurately even with sparse measurements, and outperforms the mapping based on the
predetermined propagation model.

Yet, the method shall be integrated within a simultaneous localization a mapping framework,
as the sensor positions need to be estimated as well. Also, the hypothesis of rectangular plate
geometries shall be relaxed, and more complex wave phenomena such as anisotropic and/or
multi-modal propagation, diffraction, or mode conversion are to be integrated in the model. The
recovery of the real physical parameters, such as plate thickness, and the mapping of defects
shall be investigated to make possible a complete robotic NDT task. Finally, the method shall
be tested in more realistic conditions.
# Chapter 8

## Conclusion

### 8.1 Summary of the contributions

This thesis elaborates various algorithmic frameworks to evaluate the potential of omnidirectional ultrasonic guided waves for localization and mapping purposes on large metal structures, such as storage tanks or ship hulls. The underlying motivation is to contribute to the emergence of (potentially autonomous) mobile robotic systems for accurately mapping large surfaces: by processing the long-range ultrasonic measurements acquired over the robot trajectory, the presence of potential defects that may act as acoustic scatterers could eventually be revealed.

As the knowledge of the measurement positions is critical for such an application, we first address the localization-only problem on a rectangular metal panel by relying on the ultrasonic reflections of an excitation wave on the plate boundaries. We propose a simple particle filter algorithm whose observation model is obtained by correlating the ultrasonic measurements with a wave propagation model. The particularity of this approach is that, on the contrary to the mapping strategies reviewed in the literature, no explicit echo detection nor echo association with the plate boundaries are required, whereas these are challenging tasks in signal processing, especially in the presence of noise. The experimental results obtained in a laboratory environment demonstrate the potential of this approach for accurate positioning.

We subsequently address the corresponding SLAM problem to additionally recover the position of the plate boundaries. The proposed method is a FastSLAM algorithm that integrates Delay-and-Sum (DAS) beamforming to solve the online SLAM problem. The performance of this approach is compared with that of a second FastSLAM algorithm that relies on echo detection and echo-plate boundary association. The experiments demonstrate that the former method provides more robust localization and mapping results. A SLAM experiment with a magnetic crawler is also carried out to demonstrate the applicability of ultrasonic SLAM with a real system.

Despite its apparent simplicity, DAS beamforming, which can be used for mapping acoustic scatterers such as the boundaries of a metal panel, presents several limitations, namely signal interference, and the lack of uncertainty assessment on the map estimates. To address these issues, we develop a more elaborate mapping framework that can jointly recover a probabilistic grid representation along with a feature-based representation of the map. In this context, the consistency between the two representations can be leveraged to mitigate the effect of interference. Furthermore, this novel mapping framework is integrated within a frontier exploration strategy for autonomous map reconstruction, and it is also leveraged for recovering more arbitrary polygonal shapes than rectangles. The performance of the developed methods is extensively evaluated in simulation, and experimentally validated using real robotic platforms.
In the next chapter, we propose a GraphSLAM framework to solve the full SLAM problem. In our case-study, it consists in simultaneously recovering the metal plate geometry and the sensor trajectory through batch optimization, using all the available ultrasonic measurements. We first propose a novel observation model that can leverage the high-order reflections for localization. The latter model is integrated, along with DAS beamforming, into a NLS optimizer. The experiments illustrate the benefits of the proposed approach, as it can achieve more accurate SLAM results than FastSLAM in more challenging situations.

Finally, as the SLAM methods we have elaborated essentially rely on an accurate wave propagation model, the last chapter of this thesis presents an optimal beamforming approach for mapping a metal plate with little prior knowledge on the model, as it may be lacking during a practical inspection task. The experiments performed in different conditions show that, not only the proposed approach is effective, but it can lead to better mapping results, in non-nominal acquisition conditions, when compared to mapping with a predetermined propagation model. Indeed, the latter may no longer yield accurate results in the presence of variations of the acquisition conditions, for example.

## 8.2 Perspectives

Despite the promising results on ultrasonic SLAM that we were capable to obtain, it is clear that the development of a true long-range inspection robot is a challenging task which necessitates substantial work to emerge.

Although our work proposes algorithmic frameworks for mapping acoustic scatterers with a mobile robotic system, one may note that no experiments are performed with the aim to detect actual defects. The reason is that the detection of a defect signature within a multi-echo ultrasonic measurement is challenging. This is all the more true due to multi-modal propagation of the waves, and the complex scattering effects that may occur in structures with complex geometries. Hence, future works shall focus on how the signals related to geometrical structures, those related to defects, and the signals due to noise can efficiently be identified and separated for subsequent defect identification. A promising research direction in that regard would be to explore the use of more advanced and adaptive beamforming methods, instead of the simple delay-and-sum approach, such as Minimum Variance Distorsionless Response (MVDR) beamforming, or subspace-based methods.

Another remark on the present work is that different SLAM problems are addressed in this thesis, but it is not specified exactly how they should be merged together. For example, it is not yet explicit how our FastSLAM and GraphSLAM frameworks can leverage optimal beamforming to calibrate the propagation model, and whether the mapping framework with combined representations can be easily integrated as well. The actual reason is that the different methods need first to be evaluated in realistic conditions, on actual large structures, so that one may better understand their limitations for the target application. Finally, very little work in this thesis focuses on instrumentation, whereas the efficient generation and acquisition of the signals is a major stake for long-range inspection. This also ought to be further investigated.
# Appendix: Derivation of the integrals for calculation of the inverse measurement models

In this appendix, we provide calculation details on how the integrals needed to determine inverse models in Section 5.2 are determined.

## A: Derivation of the first integral

To derive the expression of the first integral in Eq. 5.2.4, we need to show that, for $0 \le r \le a$:

$$ \int h(r - \rho) \mathcal{N}_{[0,a]}(\rho|\mu, \sigma^2) d\rho = \frac{\Phi_{\mu,\sigma}(r) - \Phi_{\mu,\sigma}(0)}{\Phi_{\mu,\sigma}(a) - \Phi_{\mu,\sigma}(0)}, $$

where $\Phi_{\mu,\sigma}$ is expressed with the error function:

$$ \Phi_{\mu,\sigma}(\rho) = \erf\left(\frac{\rho - \mu}{\sigma\sqrt{2}}\right). $$

First, we make explicit the normalization factor of the truncated Gaussian distribution:

$$ \eta = \int_{0}^{a} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{(\rho - \mu)^2}{2\sigma^2}\right\} d\rho = F(a) - F(0), $$

where $F$ is the cumulative distribution of a Gaussian distribution with mean $\mu$ and standard deviation $\sigma$: $F(x) = \frac{1}{2} \left[ 1 + \erf\left(\frac{x-\mu}{\sigma\sqrt{2}}\right) \right] = \frac{1}{2} \left[ 1 + \Phi_{\mu,\sigma}(x) \right]$. Finally, we have:

$$ \int h(r - \rho) \mathcal{N}_{[0,a]}(\rho|\mu, \sigma^2) d\rho = \frac{1}{\eta} \int_{0}^{r} \frac{1}{\sqrt{2\pi\sigma^2}} \exp\left\{-\frac{(\rho - \mu)^2}{2\sigma^2}\right\} d\rho = \frac{F(r) - F(0)}{F(a) - F(0)}, $$

that yields, when substituting the expression of $F$, the desired result.

## B: Derivation of the second integral

To derive the expression of the second integral in Eq. 5.2.4, we need to show that:

$$ \int h(r - \rho) \mathcal{N}_{[0,a]}(\rho|\mu, \sigma^2) \exp\left\{-\frac{(\rho - \nu)^2}{2b^2}\right\} d\rho = \frac{\Gamma_{\mu_{\text{new}}, \sigma_{\text{new}}}(r) - \Gamma_{\mu_{\text{new}}, \sigma_{\text{new}}}(0)}{\Phi_{\mu, \sigma}(a) - \Phi_{\mu, \sigma}(0)}, $$

where we have:

$$ \Gamma_{\mu_{\text{new}}, \sigma_{\text{new}}}(x) = \frac{\sigma_{\text{new}}}{\sigma} \cdot \exp\left(-\frac{(\mu_{\text{new}} - \mu)^2}{2\sigma^2} - \frac{(\mu_{\text{new}} - \nu)^2}{2b^2}\right) \cdot \Phi_{\mu_{\text{new}}, \sigma_{\text{new}}}(x), \\ \mu_{\text{new}} = \frac{b^2\mu + \sigma^2\nu}{b^2 + \sigma^2}, \quad \sigma_{\text{new}} = \frac{\sigma b}{\sqrt{b^2 + \sigma^2}}. $$

Given the fact that the product of two Gaussian probability density functions is still Gaussian, and based on the expressions of $\mu_{new}$ and $\sigma_{new}$, we have:

$$ \exp\left\{-\frac{(\rho - \mu)^2}{2\sigma^2} - \frac{(\rho - \nu)^2}{2b^2}\right\} = \lambda \cdot \exp\left\{-\frac{(\rho - \mu_{\text{new}})^2}{2\sigma_{\text{new}}^2}\right\}, $$
where $\lambda$ is determined by assessing the two sides of the equation at $\rho = \mu_{new}$ which yields:

$$ \lambda = \exp \left( - \frac{(\mu_{new} - \mu)^2}{2\sigma^2} - \frac{(\mu_{new} - \nu)^2}{2b^2} \right) $$

Hence:

$$ \begin{aligned} & \int h(r - \rho) N_{[0,a]}(\rho|\mu, \sigma^2) \exp \left\{ -\frac{(\rho - \nu)^2}{2b^2} \right\} d\rho \\ &= \frac{\lambda}{\eta} \int_0^r \frac{1}{\sqrt{2\pi}\sigma^2} \exp \left\{ -\frac{(\rho - \mu_{new})^2}{2\sigma_{new}^2} \right\} d\rho \\ &= \frac{\lambda\sqrt{2\pi\sigma_{new}^2}}{\eta\sqrt{2\pi}\sigma^2} \int_0^r \frac{1}{\sqrt{2\pi}\sigma_{new}^2} \exp \left\{ -\frac{(\rho - \mu_{new})^2}{2\sigma_{new}^2} \right\} d\rho \\ &= \lambda \frac{\sigma_{new}}{\sigma} \frac{\Phi_{\mu_{new},\sigma_{new}}(r) - \Phi_{\mu_{new},\sigma_{new}}(0)}{\Phi_{\mu,\sigma}(a) - \Phi_{\mu,\sigma}(0)} \end{aligned} $$

which yields the desired result after substituting the products $\lambda \cdot \sigma_{new} \cdot \Phi_{\mu_{new},\sigma_{new}}/\sigma$ with $\Gamma_{\mu_{new},\sigma_{new}}$.
Bibliography

[1] Bugwright2, autonomous robotic inspection and maintenance on ship hulls and storage tanks, description of the innovative action, 2019.

[2] Jont B Allen and David A Berkley. Image method for efficiently simulating small-room acoustics. *The Journal of the Acoustical Society of America*, 65(4):943–950, 1979.

[3] David N Alleyne and Peter Cawley. The interaction of lamb waves with defects. *IEEE transactions on ultrasonics, ferroelectrics, and frequency control*, 39(3):381–397, 1992.

[4] Francesco Amigoni and Alessandro Gallo. A multi-objective exploration strategy for mobile robots. In *Proceedings of the 2005 IEEE International Conference on Robotics and Automation*, pages 3850–3855. IEEE, 2005.

[5] Franz Andert and Lukas Goormann. Combined grid and feature-based occupancy map building in large outdoor environments. In *2007 IEEE/RSJ International Conference on Intelligent Robots and Systems*, pages 2065–2070. IEEE, 2007.

[6] Steven R Anton, Daniel J Inman, and Gyuhae Park. Reference-free damage detection using instantaneous baseline measurements. *AIAA journal*, 47(8):1952–1964, 2009.

[7] Arthur B Baggeroer and William A Kuperman. Matched field processing in ocean acoustics. In *Acoustic signal processing for ocean exploration*, pages 79–114. Springer, 1993.

[8] Joan Bordoy, Christian Schindelhauer, Fabian Höflinger, and Leonhard M Reindl. Exploiting acoustic echoes for smartphone localization and microphone self-calibration. *IEEE Transactions on Instrumentation and Measurement*, 69(4):1484–1492, 2019.

[9] Jeffrey Borish. Extension of the image model to arbitrary polyhedra. *The Journal of the Acoustical Society of America*, 75(6):1827–1836, 1984.

[10] Frederic Bourgault, Alexei A Makarenko, Stefan B Williams, Ben Grocholsky, and Hugh F Durrant-Whyte. Information based adaptive robotic exploration. In *IEEE/RSJ international conference on intelligent robots and systems*, volume 1, pages 540–545. IEEE, 2002.

[11] Peter Cawley and David Alleyne. The use of lamb waves for the long range inspection of large structures. *Ultrasonics*, 34(2-5):287–290, 1996.

[12] Sonda Chaabene,mas, Mohamed Najib Ichchou, and Mohamed Haddar. Wave mode diffusion and propagation in structural wave guide under varying temperature. *Applied Acoustics*, 108:84–91, 2016.

[13] CY Chang and FG Yuan. Extraction of guided wave dispersion curve in isotropic and anisotropic materials by matrix pencil method. *Ultrasonics*, 89:143–154, 2018.