Locally Optimal Subsampling Strategies for Full Matrix Capture Measurements in Pipe Inspection

Abstract

In ultrasonic non-destructive testing, array and matrix transducers are being employed for applications that require in-field steerability or which benefit from a higher number of insonification angles. Having many transmit channels, on the other hand, increases the measurement time and renders the use of array transducers unfeasible for many applications. In the literature, methods for reducing the number of required channels compared to the full matrix capture scheme have been proposed. Conventionally, these are based on choosing the aperture that is as wide as possible. In this publication, we investigate a scenario from the field of pipe inspection, where cracks have to be detected in specific areas near the weld. Consequently, the width of the aperture has to be chosen according to the region of interest at hand. On the basis of ray-tracing simulations which incorporate a model of the transducer directivity and beam spread at the interface, we derive application specific measures of the energy distribution over the array configuration for given regions of interest. These are used to determine feasible subsampling schemes. For the given scenario, the validity/quality of the derived subsampling schemes are compared on the basis of reconstructions using the conventional total focusing method as well as sparsity driven-reconstructions using the Fast Iterative Shrinkage-Thresholding Algorithm. The results can be used to effectively improve the measurement time for the given application without notable loss in defect detectability.

1. Introduction

Multi-channel ultrasound measurements allow a high level of flexibility in the data acquisition of imaging and defect recognition. A common acquisition scheme is the full matrix capture (FMC), which denotes consecutively transmitting with one transducer element after another, while all elements receive on each cycle. This acquisition scheme is then commonly paired with the Total Focusing Method (TFM) [1], which is an imaging algorithm applied in post-processing.

When being applied for inline inspection, where speed is crucial, FMC measurements suffer from a higher acquisition time compared to, e.g., a transmit-focused beam, as each measurement cycle requires sending with each array element. The purpose of this study is to derive a simple heuristic for choosing only a subset of the elements for transmission to improve the acquisition speed while preserving the achievable imaging quality. The algorithm takes a priori knowledge of the problem geometry, i.e., a forward model, into account. Furthermore, only a specific region of interest is considered, limiting the width of the usable aperture.

In the literature, reducing the number of transmit or receive channels of a FMC acquisition has been studied both in the light of reducing acquisition as well as processing times and hardware cost. The design of such sparse ultrasound arrays is extensively being studied for medical applications, where the most common optimization criteria are side lobe levels and point spread functions over the full aperture. In non-destructive testing, transmit subsampling is being investigated in [2], where the authors show that using two or four elements near the end of the array for excitation is to be favored over transmitting with all elements as one. In [3,4], equispaced element choices as well as using only the outer elements are being investigated, with the conclusion that using fewer elements can be an application-dependent choice to improve measurement speed. In [5], sparse arrays are designed using a combination of almost different sets and a genetic algorithm using attributes of the point spread function as the fitness measure. In [6], using fewer receivers is investigated for near surface defects. In [7], particle swarm optimization is performed to reduce the computation time of FMC by using fewer elements. In [8], a Compressive Sensing (CS)-based approach is used for reducing the data both in receive as well as in frequency based on Cramer Rao Bounds.

All aforementioned publications deal with the optimization of the patterns in terms of the full aperture width. In contrast, we consider a scenario where, due to the limited extent of the region of interest, only a subset of the full width of the aperture is of relevance. Similar investigations were carried out in [9] for acoustic source localization using microphones, prioritizing certain regions in the image. In the same spirit, we derive a simpler heuristic. In contrast, our approach is based on a more sophisticated forward model, including effects such as refraction and beam patterns, for an application in non-destructive testing.

We investigate a scenario where a bowed transducer array is positioned over a moving pipe in production. While the transducer encloses a section of the pipe to ensure recording echoes from a large variety of angles, as well as to enable in-field adjustability, the region where defects are expected is limited to the weld. This also limits the usable aperture, posing the question of how to choose the active array elements.

Consequently, we derive a scheme where the transmit elements are chosen on the basis of the a priori knowledge of the problem geometry and area of interest. The reconstruction as well as the array element selection are performed on the basis of a linearized forward model, which, calculated by a raycasting scheme, incorporates effects as transducer directivity and refraction at the pipe boundary. The reconstruction is performed by TFM. We additionally also compare to the Fast Iterative Shrinkage-Thresholding Algorithm (FISTA) [10] for reconstruction and tailor a forward model-based weighting scheme for the scenario at hand that benefits both algorithms.

2. Materials and Methods

2.1. Forward Model

The post-processing algorithms considered in this publication combine data from multiple observations for imaging. The bowed linear array transducer that is used for the measurements is assumed to be stationary during each cycle of transmit and receive events. The data acquisition is performed such that one transducer element shoots at a time and multiple transducer elements record the echoes [1]. The data from all transmit and receive events at one location is then fused in post-processing.

We consider each observation to be a superposition of multiple time shifted echoes as

$$g_{m_{tx},m_{rx}}\left(t\right)=\Re \left\{{\displaystyle \sum }_d{\displaystyle \sum }_p{\alpha }_{m_{tx},m_{rx},p,d}{h}^{\left(a\right)}\left(t-{\tau }_{m_{tx},m_{rx},p,d}\right)\right\}+n\left(t\right),$$

(1)

where $g\left(t\right)$ is an A-scan and $m_{tx}$ and $m_{rx}$ index a transmitter and receiver, respectively. The analytic signal of the transmit pulse is given by $h^{\left(a\right)}\left(t\right)$, $\Re \{.\}$ takes the real part. Each A-scan consists of a superposition of echoes from multiple defects $d$. The ultrasonic wave can have taken different propagation paths and modes from transmitter to defect to receiver, which are indexed by $p$ (see also [11,12]). The complex factor $\alpha$ models the amplitude and phase change on the path, constructed from the transducer directivity, the interfaces, and the interaction with the defect itself. The time-shift $\tau$ indicates the propagation time on the path. The variable $n\left(t\right)$ indicates zero-mean additive white Gaussian noise.

2.2. Total Focusing Method

The Total Focusing Method (TFM) allows the estimation of the reflectivity of the specimen at arbitrary points of interest $r\in {ℝ}^2$ with $r={[x,y]}^T$. Starting from the model given in Equation (1), we assume a perfect point scatterer to be located at $r$. The reflectivity is then given by correlating the measurement data with the estimated point scatterer response. Repeating this procedure for a set of points of interest yields an image of the reflectivity of the medium.

For the calculation of the reconstruction, we limit ourselves to longitudinal propagation, as responses including transversal propagation can be omitted from the data. Furthermore, multi-view data fusion is an open research topic [11] and out of the scope of this publication. Amplitude changes that are due to interfaces have been shown to be neglectable for the investigated scenario; furthermore, the influence of attenuation is considered to be small enough to be neglected as well, which is in accordance with the literature [11].

The TFM reconstruction at a point $r$ is then given by [1]

$$I\left(r\right)=\left|{\displaystyle \sum }_{m_{tx}}{\displaystyle \sum }_{m_{rx}}{\alpha }_{m_{tx},m_{rx}}\left(r\right)g_{m_{tx},m_{rx}}^{\left(a\right)}\left({\tau }_{m_{tx},m_{rx}}\left(r\right)\right)\right|,$$

(2)

where ${\alpha }_{m_{tx},m_{rx}}\left(r\right)={\alpha }_{m_{tx}}\left(r\right){\alpha }_{m_{rx}}\left(r\right)$ is a real weighting factor constructed from the transducer directivities ${\alpha }_{m_{tx}}\left(r\right)$ and ${\alpha }_{m_{rx}}\left(r\right)$ at the transmitter element and the receiver element, respectively. We consider directivities to follow a Gaussian shape as $\alpha \left(r\right)=\exp \left(\frac{-\varphi \left(r\right)}{2\cdot {(0.5{\varphi }_{max})}^2}\right)$ here, where $\varphi \left(r\right)$ is the angle of the according propagation path at the transducer element, which can be calculated once the path from the transducer element to $r$ has been determined.

2.2.1. Matrix Formulation

For brevity of notation we now limit ourselves to reconstructions which are performed on an equispaced grid, i.e., to $r=r\left(n_z,n_x\right)={[n_z{d}_z+z_0,n_x{d}_x+x_0]}^T$ where $n_z$ and $n_x$ are integers in the half-open intervals $\left[0,N_z\right]$ and $\left[0,N_x\right]$, respectively, $d_z$ and $d_x$ denote the sampling intervals and $z_0$ and $x_0$ a start point in the according directions. Furthermore, we denote the number of transducer elements as $N_{el}$ and, in the case of using fewer transmit cycles, the number of transmit cycles $N_{tx}$. The A-scans are considered to be time-discrete equivalents of the continuous definition used before, holding $N_t$ samples acquired at a sampling frequency $f_S$.

Equation (2) can be rewritten as matrix-vector product [8,13,14]

$$i_{tfm}=\left|Mg\right|,$$

(3)

where $g={\left[g_{0,0}^{{\left(a\right)}^T},g_{0,1}^{{\left(a\right)}^T},\dots ,{g_{0,N_{el}-1}^{\left(a\right)}}^T,g_{1,0}^{{\left(a\right)}^T},\dots ,{g_{N_{tx}-1,N_{el}-1}^{\left(a\right)}}^T\right]}^T\in {ℂ}^{N_{tx}{N}_{el}{N}_t}$ are the stacked analytic representations of the measurement data with $g_{m_{tx},m_{rx}}^{\left(a\right)}\in {ℝ}^{N_t}$ being a single A-scan, $i_{tfm}=[I(\left[z_0,x_0{]}^T\right),I(\left[z_0+d_z,x_0{]}^T\right),\dots ,I{(\left[z_0,x_0+d_x{]}^T\right),\dots ]}^T\in {ℝ}^{N_z{N}_x}$ is the vectorized image on the equispaced grid, $M\in {ℝ}^{N_z{N}_x\times N_{tx}{N}_{el}{N}_t}$ is the model matrix, and $|.|$ denotes the elementwise absolute value.

The matrix is constructed from blocks $M_{\left(n_z,n_x\right),\left(m_{tx},m_{rx}\right)}\in {ℝ}^{1,N_t}$ as

$$M=\left[\begin{array}{ccccc}{M}_{\left(0,0\right),\left(0,0\right)}& M_{\left(0,0\right),\left(0,1\right)}& \dots & M_{\left(0,0\right),\left(1,0\right)}& \dots M_{\left(0,0\right),\left(N_{tx}-1,N_{el}-1\right)}\\ M_{\left(1,0\right),\left(0,0\right)}& \dots & & & \\ \dots & & & & \\ M_{\left(0,1\right),\left(0,0\right)}& & & & \\ \dots & & & & \\ M_{\left(N_z-1,N_x-1\right),\left(0,0\right)}& & & & \end{array}\right],$$

(4)

with each block being defined as

$${\left[M_{\left(n_z,n_x\right),\left(m_{tx},m_{rx}\right)}\right]}_{m_t}=\{\begin{array}{ll}{\alpha }_{m_{tx},m_{rx}}\left(r\left(n_z,n_x\right)\right)\hfill & \mathrm{if}\lfloor f_S{\tau }_{m_{tx},m_{rx}}\left(r\left(n_z,n_x\right)\right)\rfloor =m_t\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array},$$

(5)

where $m_t\in \left[0,N_t\right]$ denotes an integer index, ${[.]}_m$ denotes access to the particular element of the column vector, and $\lfloor .\rfloor$ is the floor operation. In contrast to the state of the art TFM formulations e.g., in [1,11] the times of flights are rounded to the temporal sampling grid here. Formulating the problem this way allows for an efficient precalculation of the matrix, reducing the reconstruction to a simple matrix vector multiplication [14] and opening the door for the application of algorithms like FISTA.

The model matrix $M$ can be calculated for different transmit schemes. In the remainder of the paper, we assume that $M$ was calculated according to the given transmit pattern, furthermore, as a special case, we denote $M_{fmc}$ with $N_{tx}=N_{el}$ as the matrix that corresponds to the reconstruction of a FMC dataset. Reducing the size of the matrix also consequently reduces the computation time of the reconstruction, which is a simple matrix vector product here.

2.2.2. Weighting

In the given scenario, the number of A-scans contributing to each reconstruction voxel varies strongly. This is due to the limited insonification angle of each transducer element, the influence of the specimen geometry and refraction, as well as the reconstruction being performed on a cartesian grid, whereas the measurement setup is laid out in a circular pattern. In [15], a weighting to counter the effects of transducer directivity and beamspread was applied. In a similar spirit, we introduce a weighting that corrects for the influence of the refracting boundary in combination with the directivity of our transducer elements.

A measure of the visibility of a reconstruction voxel from each active transducer is given by the number of nonzero elements in its corresponding row in $M$. Summing over these elements further includes the effect of the directivity of the transducers. The vector $e\in {ℝ}^{N_z{N}_x}$ captures this measure for all voxels and can be computed by:

$$e=M_{fmc}\hspace{0.17em}\mathbf{1}.$$

(6)

The weight matrix to scale the reconstruction can then be constructed as the inverse of a diagonal matrix constructed from $e$

$$W={\left(\mathrm{diag}\left(e\right)\right)}^{-1},$$

(7)

yielding the weighted TFM reconstruction as

$$i_{wtfm}=\left|W\hspace{0.17em}M\hspace{0.17em}g\right|.$$

(8)

2.3. FISTA

The task of recovering the point scatterer locations from the measurements using a linear model of the reflectivity of the medium can also be considered as a sparse recovery problem, when the number of point scatterers is small. The reconstruction can therefore also be obtained by approximating a solution to

$$\underset{i_{FISTA}}{\min }||g-M^H{W}^H{i}_{FISTA}|{|}_2^2+\lambda ||i_{FISTA}|{|}_1,$$

(9)

where $i_{FISTA}\in {ℂ}^{N_z{N}_x}$ is the sparse representation, which is equivalent to an image of the reflectivity of the medium, $||.|{|}_n$ denotes the $n$-norm, $\lambda$ is a regularization parameter controlling the sparsity of the result, and ${.}^H$ denotes the Hermitian transpose. Such problems can be approximated by the FISTA algorithm [10], yielding a regularized image of the reflectivity according to $M$ for our usecase. The matrix $W$ serves as a normalization to the dictionary $M$. Note that this is common practice in the Simultaneous Algebraic Reconstruction Technique (SART) [16].

We use the implementation of the algorithm provided in the fastmat package for python [17]. Following [18], the regularization parameter is chosen as

$$\lambda =\beta |\left|WMg\right|{|}_{\infty }$$

(10)

with $0\le \beta \le 1$.

2.4. Raycasting

All the above formulations require the ray paths to be known for the calculation of $M$. In the literature, two approaches to determine the propagation paths between transmitters, receivers (which are interchangeable due to acoustic reciprocity), and $r$ are prevalent. In many publications, iterative schemes for finding a propagation path that fulfils Fermats principle given one or multiple surfaces are used [11,12,19,20]. On the other hand, e.g., in [21], the Eikonal equations are solved directly to determine the propagation times.

In contrast, we use a scheme where a set of rays is shot for each transmit and receive event. The rays are propagated through the medium and refracted and reflected at each surface on their path. The propagation times and amplitudes are extrapolated to the grid of points of interest on each step using an adaption of the Bresenham line tracing algorithm [22], expanded by a linear interpolation of the amplitudes and propagation times along the line. A sketch that outlines the procedure is given in Figure 1. The rays are traced on a map that represents the time of flights and amplitudes from which the model is built.

2.5. Transmit Subsampling Strategy

By using fewer transmit cycles, the measurement time can be improved. Our goal is to use the information derived from the raycasting procedure and forward modelling to get an optimal choice of transmit elements for imaging a region of interest. The steps of the algorithm are also summarized in Figure 2. The model in Equation (3) can be used to derive the linear forward model. The following equation

$$\widehat{g}=M_{fmc}^H\widehat{i}$$

(11)

yields a synthetic dataset $\widehat{g}$ given by a superposition of all point scatterer responses modeled by $M_{fmc}^H\widehat{i}$. We now choose $\widehat{i}$ such that all entries corresponding to voxels in our region of interest have value 1 and 0 elsewhere. In consequence, $\widehat{g}$ now holds data which, though following a simplified model, reflects a scenario where the full region of interest is reflecting.

We can now sum these theoretical responses to get a measure of how much the excitation of each transducer element contributes to imaging the region of interest as

$$[c{]}_{m_{tx}}={\displaystyle \sum }_{m_{rx}}{\displaystyle \sum }_{m_t}|[{\widehat{g}}_{m_{tx},m_{rx}}{]}_{m_t}|$$

(12)

with $c\in {ℝ}^{N_{el}}$. We assume here that the elements indexed by $m_{tx}$ are ordered such that for adjacent transducer elements the contributions are given by ${[c]}_i$ and ${[c]}_{i+1}$. In order to get the best trade off between collecting the maximum energy from the region of interest (choosing only the center elements) and using an as wide as possible aperture (using only the outer elements) we use the following heuristic. We divide $c$ into $N_{tx}$ sections where each section has the same area under curve. We then use the center transducer element corresponding to the center index of each section as transmitter. This can be expressed by the following computation steps. We first compute the cumulative sum $c^{\left(c\right)}\in {ℝ}^{N_{el}}$ of the elements of $c$,

$${\left[c^{\left(c\right)}\right]}_{m_{tx}}=\frac{1}{{\mathbf{1}}^{T}c}{\displaystyle \sum }_{i=0}^{m_{tx}}{[c]}_i$$

(13)

normalized such that the maximum is 1. We then find our $k$-th transmitter element index ${[p]}_k$ with $p\in {\mathbb{N}}^{N_{tx}}$

$$p_k=\underset{\xi }{\min } \xi \mathrm{s}.\mathrm{t}. {\left[c^{\left(c\right)}\right]}_{\xi }\ge \frac{1}{N_{tx}}\left(\frac{1}{2}+k\right)$$

(14)

which is the smallest index for which $c^{\left(c\right)}$ exceeds the given threshold.

3. Results

3.1. Simulation Scenario

Simulations are performed in CIVA 2020 from EXTENDE [23]. The scenario is depicted in Figure 3. A summary of the simulation parameters is given in

Table 1. Simulation scenario parameters.

Parameter	Value
transducer radius transducer element count transducer element width pitch outer pipe radius inner pipe radius pipe material coupling material	35 mm 128 0.3 mm 0.05 mm 16.5 mm 11.7 mm Steel Water

. The simulation mimics an inline inspection application where a circular array is used to monitor a pipe in production. Water is used as an acoustic couplant between transducer and pipe. The weld region is configured to behave like the surrounding material, except for defects being placed in that region, as the main focus of the publication is the investigation of the transmit pattern choice rather than thoroughly modelling a weld.

For each simulation, one inclusion is modelled by a $0.5 \mathrm{mm}$diameter circular defect. Simulations with the defect being placed at radius $13 \mathrm{mm}$ and $15 \mathrm{mm}$ are performed. The defect is placed centered with respect to the transducer in the angular domain.

Furthermore, $0 \mathrm{dB}$ zero-mean white Gaussian noise is added to the simulated measurement data to account for imperfections in the electronics and measurement. The term $0 \mathrm{dB}$ refers here to the noise variance being equal to the energy of the first echo from the defect at the center transducer element in pulse echo mode.

3.2. Transmit Subsampling

3.2.1. Impact of the Region of Interest Selection

In order to optimize our subsampling, we choose a region of interest which overlays the weld. The variation of the visibility of the region of interest over the transducer elements is displayed in Figure 4. When choosing a map that reflects a region near the front wall, as given by (b), only the center transducer elements contribute to the imaging. Vice versa, a region near the inner pipe radius (c) can be imaged using a much wider aperture.

Consequently, since imaging with the widest possible aperture is beneficial, as defects like cracks can show very direction-specific reflection patterns, the width of the subsampled aperture needs to be chosen carefully. For the remainder of the paper we use the region of interest depicted in (c), which overlays the full weld as our best a priori guess to the defect positions to be observed in the data. It spans an area of $3.56\hspace{0.17em}\mathrm{mm}\times 1.02 \mathrm{mm}$. From Figure 4, it is apparent that this region results in an aperture choice that is a mix of (a) and (b).

The 16, 8, and 4 transmit elements chosen on the basis of the region of interest in our scenario are displayed in Figure 5. The center of the array is being sampled more densely, as it contributes more energy to our region of interest. Still, a wide aperture is used which can help with the imaging of very specular reflecting defects.

3.2.2. Verification of the Transmit Element Selection Scheme

To verify the approach to subsampling we perform reconstructions of the simulated measurement data with wider or narrower transmit patterns. To generate these patterns, the vector $c$ is linearly interpolated or extrapolated to reflect a wider or narrower energy spread, the patterns are then chosen from that deviated distribution. This is done in order to verify that the chosen approach of taking the distribution from the model matrix is valid.

In Figure 6, the interpolated distributions $c$ are related to the corresponding scale parameters. For larger values of the scaling factor, the curve becomes wider, yielding a more uniform transmitter distribution over the full aperture. In contrast, a small scale parameter leads to a narrower distribution that results in only picking the center elements of the array as transducers.

We evaluate the transmit element selection scheme in Figure 7 by choosing the transmit pattern from a scaled $c$ (the scaled $c$ are displayed in Figure 6 and performing a reconstruction both with FISTA as well as TFM on the given data. The reconstructions are then compared by the ratio of the energy in a square region around the defect to the energy in the whole image. The evaluation has been run both for a defect at $z=13 \mathrm{mm}$ in (a) as well as for a defect at $z=15\hspace{0.17em}\mathrm{mm}$ averaging over 10 noise realizations for each. All plots have been scaled to maximum 1 for better comparison.

The curves show that the TFM reconstruction performs best for narrow transmit patterns, except for a slight increase for 8 and 16 elements from scaling factor 0.1 to 0.25. The performance decreases for wider aperture choices. The original distribution at scaling factor 1.0 is approximately at half of the peak to minimum performance value. For FISTA, the reconstruction of the defect at $z=13\hspace{0.17em}\mathrm{mm}$ benefits from using a pattern with scaling factor 0.5 to 1, even failing to reconstruct when picking only the center elements. For the other defect, the algorithm performs in a similar fashion as TFM. The results for FISTA show relatively large jumps for small changes of the scaling factor. This illustrates that the algorithm reacts strongly to small changes in the presence of multiple scattering events, leading to variations of the proportion of the front wall energy to defect energy in the reconstruction. For future investigations, the effect could be decreased by separating the front wall and back wall from the remainder of the data.

The measures indicate that, except for FISTA on the $z=13 \mathrm{mm}$ defect, choosing a narrower distribution yields more energy for the defect in the reconstruction, compared to the overall energy. When choosing a wider aperture, the performance slowly starts to decrease, as the ultrasonic waves from the outer transmitters that are reflected by the defect have less amplitude due to the transducer directivity and boundary. On the other hand, choosing a wider distribution is beneficial if the defect reflects stronger to the side than back to the center, as might be the case for a crack. This mismatch is due to the model being derived for perfect point scatterers, which have, unlike the side-drilled holes used in the simulations, an omnidirectional reflection pattern.

Consequently, as the measures in Figure 7 are derived for side-drilled holes as reflectors, the results indicate the aperture has been chosen slightly too large. Still, the choice seems suitable, as the algorithm performance is good, and will be more robust to directional reflectors or off-center than a more narrow aperture. The given measures reflect the performance of the algorithms for a side-drilled hole that is located perfectly in the center of the aperture. It is also indicated that the aperture choice should not be too narrow by the FISTA reconstruction in Figure 7a and the slight increase of most curves from scaling factor 0.1 to 0.25. We therefore deem the transmit pattern selection scheme suitable for the application at hand.

3.2.3. FISTA and the Dictionary Normalization

FISTA requires tuning the regularization parameter as well as number of steps for the scenario at hand. We empirically determined that $\beta =0.01$ as well as running 20 steps of FISTA is suitable for our application and kept these values constant for all investigations. To justify the dictionary normalization introduced by the matrix $W$ in Equations (9) and (10) the performance of the algorithm is compared with and without the weighting. The comparison is carried out using four transmit elements that were determined by the proposed selection scheme. The resulting B-scans are displayed in Figure 8.

It can be seen that without the weighting, i.e., a uniform $W=I$, the level of the noise compared to the maximum amplitude of the defect in Figure 8a is higher. This is due to the rows of $M$ corresponding to scatterers in that region having a larger norm, as the region is visible over a wider aperture (compare also to Figure 4). Consequently, FISTA yields a less sparse solution for scatterers being predominantly in that region, as a slight change to values of $i_{FISTA}$ associated with that region introduces a large change to the data fidelity term $||.|{|}_2^2$ but a small change to the regularization term. For the defect at $z=15 \mathrm{mm}$ in Figure 8b, the weighted reconstruction has a better defect to front wall separation. A likely cause for this is the uniformly weighted reconstruction favoring the backwall reflection much more than the weighted reconstruction, leading to a more faint defect image.

Concluding, without the introduction of $W$ as given in Equation (7), the FISTA shows more noise, preferably near the inner pipe radius. This is due to the non-uniform effect of the L1-norm regularization term, which is meant to promote sparsity, for non-normalized dictionaries. In praxis, the proposed weighting scheme is capable of countering that effect.

3.2.4. Reconstruction of the Simulated Data

The reconstructions for 4, 8, and 16 transmit elements chosen according to the proposed selection scheme as well as of the full FMC data are given in Figure 9. For all defect locations, transmit element counts, as well as both algorithms, the defects are clearly visible. The TFM reconstruction shows more background noise, predominantly for the reconstructions of the defect at $z=13 \mathrm{mm}$. The difference of the perceived noise level of the TFM reconstruction for both defects is due to the colormap being adjusted to the maximum amplitude in the region of the defect of the according image. As the defect at $z=13 \mathrm{mm}$ has a smaller amplitude, the noise is more predominant here.

The FISTA reconstructions show a better separation of defect and frontwall/backwall reflection as well as less noise than the TFM reconstructions, which is due to the algorithm favoring sparser images. Consequently, the algorithm shows a better behavior with regards to the noise in the measurement data. Still, this requires a suitably chosen $\mathit{\beta }$ as well as number of iterations.

Despite the strong subsampling, imaging the according defect is still possible for all settings. Using fewer transmit cycles increases the noise level, as there is less data to be averaged over. For (a) TFM using four transmit cycles the noise floor becomes more pronounced, for FISTA, the defect echo is a bit more faint than in the other images. Consequently, when using four elements, imaging the defects is still possible for the given setup. Using 8 or 16 transmit cycles improves the result so far that no apparent difference to the reconstruction of the FMC data is apparent in the region of interest.

3.2.5. Performance in the Presence of Structural Noise

In addition to the investigations in the last subsection, which display the performance when the data is subject to $0 \mathrm{dB}$ additive white Gaussian noise, in order to emulate a highly scattering medium we also added more realistic structural noise in CIVA. The scatterers were placed with a density of $100 \mathrm{pts}/{\mathrm{mm}}^3$ and amplitude 1. The reconstructions of data subject to both noise sources are displayed in Figure 10.

Using fewer transmit cycles clearly degrades the imaging quality in the presence of structural noise, with the effect being more pronounced than for the dataset that only contained additive white Gaussian noise in Figure 9. This is due to the fact that the additional small scatterers in the medium that emulate the structural noise show up over multiple measurements according to the forward model and therefore cannot simply be averaged out like the zero-mean Gaussian noise. In the presence of such small-scale scattering events, FISTA achieves a visually better defect to noise separation due to the sparsity constraint.

In consequence, the simulation shows that for a real world application of the subsampling algorithm to a highly scattering specimen the noise level from e.g., the grain needs to be taken into account. Based on FMC measurement data, a tuning of the number of transmit channels should be performed.

4. Discussion

We show that a significant reduction of transmit cycles is possible for the simulation scenario at hand. Using only 8 instead of 128 elements in transmit, the TFM and FISTA reconstructions still yield a clear image of the defects. Consequently, the measurement time for the acquisition of one dataset can be reduced by factor 16 with virtually no impairment of the imaging quality for the given region of interest. Furthermore, the computation time of the reconstruction algorithms is scaled down by approximately the same factor. For highly scattering media, a tuning of the number of transmit channels according to the noise level is advised.

The transmit elements are chosen according to a deterministic scheme which, for a given region of interest, selects an as wide as suitable aperture while emphasizing those transmit elements that contribute more energy. The scheme is verified by evaluating the performance of the reconstruction algorithms for wider and narrower choices of the aperture. From these observations, it becomes apparent that using a slightly narrower aperture would benefit the imaging of the side-drilled holes used in the simulations. Still, as less well-behaving and more directional reflections are expected when imaging actual cracks or inclusions, the given aperture choice seems reasonable and in a good performance range over the full region of interest.

Inspired by SART, we introduce a normalization of the model matrix, which improves the performance of FISTA in the given scenario. The algorithm does not perform optimally for problems that give varying weight to contributions in the imaging region. Though, this is natural for ultrasound imaging scenarios, as different areas in the imaging region are in most cases being observed by a varying number of transducers. We show empirically that the proposed normalization improves the results in the given scenario.

The reconstructions, normalization for FISTA, as well as transmit element selection are based on a linear forward model which is derived using a ray casting scheme. It includes effects such as the directivity of the transducer elements as well as the refraction at the pipe boundary. The algorithms are expressed in terms of linear algebra, allowing easy and fast parallel computation and precomputation.

As next step, the findings will be applied and verified on real measurement data. Furthermore, including other propagation modes than only longitudinal propagation seems beneficial, further influencing the aperture choice as well as the possibilities for imaging from a larger aperture. Including multiple propagation modes is inherently possible given the way the algorithms are expressed here, yet the fusion of the multi-mode images will require further investigation [11]. Furthermore, in order to reduce the processing time further, including subsampling on the receive side might be beneficial.