Figure 5: Utility graph with various gradient norm $C$ and noise level
$\sigma$

(a)

(b)

Figure 6: Utility graph with various gradient norm $C$ and noise level
$\sigma$ for MNIST and CIFAR$10$ data sets.

As can be seen from Figure 4 and Figure 5, for clipping constant $C=0.1$, we
can choose the maximum $\sigma$ somewhere in the range
$\sigma\in\left[18,22\right]$ if we want to guarantee there is at most about
10% accuracy loss compared to the (near)-optimal solution without noise.
Another option is $C=0.075$, where we can tolerate
$\sigma\in\left[18,30\right]$ yielding the same accuracy loss guarantee. When
the gradient bound $C$ gets smaller, our DP-SGD can tolerate bigger noise,
i.e, bigger values of $\sigma$. However, we need to increase the number $K$ of
iterations during the training process when $C$ is smaller in order to
converge and gain a specific test accuracy – this is the trade-off. For
simplicity, we intentionally choose $C=0.1$, $\sigma\leq 20$ and expected test
accuracy loss about 10% for our experiments with strongly convex and plain
convex objective functions.

The utility graph is extended to the non-convex objective function in Figure
6. To keep the test accuracy loss less or equal to 10% (of the final test
accuracy of the original model $w^{*}$), we choose $C=0.025$ and noise level
$\sigma\leq 12$ for MNIST data set (as shown in Figure 6a) and $C=0.025$ and
noise level $\sigma\leq 6.572$ for CIFAR$10$ data set (as shown in Figure 6b).
For simplicity, we use this parameter setting for our experiments with the
non-convex problem.

###  D.3 Asynchronous distributed learning with differential privacy

We consider the asynchronous DP-SGD framework with strongly convex, plain
convex and non-convex objective functions for different settings, i.e.,
different levels of privacy budget $\epsilon$ and different constant sample
size sequences.

####  D.3.1 Asynchronous DP-SGD with different constant sample size sequences

The purpose of this experiment is to investigate which is the best constant
sample size sequence $s_{i}=s$. This experiment allows us to choose a decent
sample size sequence that will be used in our subsequent experiments. To make
the analysis simple, we consider our asynchronous DP-SGD framework with
$\Upsilon(k,i)$ defined as false if and only if $k<i-1$, i.e., compute nodes
are allowed to run fast and/or have small communication latency such that
broadcast global models are at most $1$ local round in time behind (so
different clients can be asynchronous with respect to one another for 1 local
round). We also use iid data sets. The detailed parameters are in Table 2.

Table 2: Basic parameter setting for strongly and plain convex problems Parameter | Value | Note
---|---|---
$\bar{\eta}_{0}$ | $0.1$ | initial stepsize
$N_{c}$ | $10,000$ | # of data points
$K$ | $50,000$ | # of iterations
$\epsilon$ | $0.04945$ |
$\sigma$ | $19.29962$ |
$\delta$ | $0.0001$ |
$C$ | $0.1$ | clipping constant
$s$ | $\\{1,5,10,15,20,26\\}$ | constant sample size sequence
dataset | LIBSVM | iid dataset
$n$ | $5$ | # of nodes
$\Upsilon$ | $k\geq i-1$ |  $1-$asynchronous round

(a) Strong convex.

(b) Plain convex.

Figure 7: Effect of different constant sample size sequences

(a) Strong convex.

(b) Plain convex.

Figure 8: Effect of different constant sample size sequences

The results from Figure 7 to Figure 8 confirm that our asynchronous DP-SGD
framework can converge under a very small privacy budget. When the constant
sample size $s=1$, it is clear that the DP-SGD algorithm does not achieve good
accuracy compared to other constant sample sizes even though this setting has
the maximum number of communication rounds. When we choose constant sample
size $s=26$ (this meets the upper bound for constant sample sizes for our
small $N=10,000$ and small $\epsilon\approx 0.05$, see Theorem 4), our DP-SGD
framework converges to a decent test accuracy, i.e, the test accuracy loss is
expected less than or equal to 10% when compared to the original mini-batch
SGD without noise. In conclusion, this experiment demonstrates that our
asynchronous DP-SGD with diminishing step size scheme and constant sample size
sequence works well under DP setting, i.e, our asynchronous DP-SGD framework
can gain differential privacy guarantees while maintaining an acceptable
accuracy.

Table 3: Basic parameter setting for non-convex problem with MNIST data set Parameter | Value | Note
---|---|---
$\bar{\eta}_{0}$ | $0.1$ | initial stepsize
$N_{c}$ | $60,000$ | # of data points
$K$ | $360,000$ | # of iterations
$\epsilon$ | $0.15007$ |
$\sigma$ | $12.10881$ |
$\delta$ | $1.667\cdot 10^{-5}$ |
$C$ | $0.025$ | clipping constant
$s$ | $\\{10,25,50,100,200,300,370\\}$ | constant sample size sequence
dataset | MNIST | iid dataset
$n$ | $5$ | # of nodes
$\Upsilon$ | $k\geq i-1$ |  $1-$asynchronous round
Table 4: Basic parameter setting for non-convex problem for CIFAR$10$ data set Parameter | Value | Note
---|---|---
$\bar{\eta}_{0}$ | $0.1$ | initial stepsize
$N_{c}$ | $50,000$ | # of data points
$K$ | $350,000$ | # of iterations
$\epsilon$ | $0.50102$ |
$\sigma$ | $6.572$ |
$\delta$ | $2\cdot 10^{-5}$ |
$C$ | $0.025$ | clipping constant
$s$ | $\\{10,25,50,100,300,500,689\\}$ | constant sample size sequence
dataset |  CIFAR$10$ | iid dataset
$n$ | $5$ | # of nodes
$\Upsilon$ | $k\geq i-1$ |  $1-$asynchronous round
(a) Non-convex. Figure 9: Effect of different constant sample size sequences
(a) Non-convex. Figure 10: Effect of different constant sample size sequences.

We also conduct the experiment for the non-convex objective function with
MNIST and CIFAR$10$ data sets.

The detailed parameter settings can be found in Table 3 and Table 4. Here, we
again consider our asynchronous setting where each compute node is allowed to
run fast and/or has small communication latency such that broadcast global
models are at most $1$ local round in time behind. As can be seen from Figure
9 ( with MNIST data set), our proposed asynchronous DP-SGD still converges
under small privacy budget. Moreover, when we use the constant sample size
$s=370$ (this meets the upper bound for constant sample sizes for our small
$N=60,000$ and small $\epsilon\approx 0.15$, see Theorem 4), we can
significantly reduce the communication cost compared to other constant sample
sizes while keeping the test accuracy loss within 10%. The constant sample
size $s=10$ (as well as $s\leq 10$) shows a worse performance while this
setting requires more communication rounds, compared to other constant sample
sizes. We can observe the same pattern for CIFAR$10$ data set as shown in
Figure 10, where we can choose the constant sample size $s\leq 689$ with
$N=50,000$ data points and $\epsilon\approx 0.5$. While the constant sample
size $s$ satisfying $300\leq s\leq 689$, the test accuracy gets the highest
level while the constant sample size $s\leq 50$ deteriorates the performance
of accuracy significantly. This figure again confirms the effectiveness of our
asynchronous DP-SGD framework towards a strong privacy guarantee for all types
of objective function.

####  D.3.2 Asynchronous DP-SGD with different levels of privacy budget

Table 5: Different privacy budget settings for strongly and plain convex problems Privacy budget $(\epsilon,\delta)$ | $\sigma$ |  Sample size $s$
---|---|---
$(0.04945,0.0001)$ | $19.29962$ | $26$
$(0.1,0.0001)$ | $13.06742$ | $55$
$(0.25,0.0001)$ | $8.59143$ | $103$
$(0.5,0.0001)$ | $6.05868$ | $168$
$(1.0,0.0001)$ | $4.27273$ | $265$
$(2.0,0.0001)$ | $3.03241$ | $400$

We conduct the following experiments to compare the effect of our DP-SGD
framework for different levels of privacy budget $\epsilon$ including the non-
DP setting (i.e., no privacy at all, hence, no noise). The purpose of this
experiment is to show that the test accuracy degradation is at most 10% even
if we use very small $\epsilon$. The detailed constant sample sequence $s$ and
noise level $\sigma$ based on Theorem 4 are illustrated in Table 5. Other
parameter settings, such as initial stepsize $\eta_{0}$, are kept the same as
in Table 2.

(a) Strong convex.

(b) Plain convex.

Figure 11: Effect of different levels of privacy budgets $\epsilon$ and non-DP
settings

(a) Strong convex.

(b) Plain convex.

Figure 12: Effect of different levels of privacy budgets $\epsilon$ and non-DP
settings

As can be seen from Figures 11 and Figure 12, the test accuracy degradation is
about 10% for $\epsilon=0.04945$ compared to the other graphed privacy
settings and non-DP setting. Privacy budget $\epsilon=0.1$, still significant
smaller than what is reported in literature, comes very close to the maximum
attainable test accuracy of the non-DP setting.

We ran the same experiment for the non-convex objective function. The detailed
setting of different privacy budgets is shown in Table 6. Note that we also
set the asynchronous behavior to be $1$ asynchronous round, and the total of
iterations on each compute node is $K=360,000$. Other parameter settings for
the non-convex case, such as initial stepsize $\eta_{0}$, are kept the same as
in Table 3. As can be seen from Figure 13 with MNIST data set, the test
accuracy loss with $\epsilon\approx 0.15$ is less than 10% (the expected test
accuracy degradation from utility graph at Figure 6). Another pattern can be
found in Figure 14. By selecting $\epsilon\approx 0.5$ for CIFAR$10$ data set,
the test accuracy reduces less than $10\%$, compared to the non-DP setting.
Note that we use AlexNet for CIFAR$10$, which shows $\approx 0.74$ maximum
test accuracy in practice111111https://github.com/icpm/pytorch-cifar10.

Table 6: Different privacy budget settings for non-convex problem for MNIST data set Privacy budget $(\epsilon,\delta)$ | $\sigma$ |  Sample size $s$
---|---|---
$(0.15007,1.667\cdot 10^{-5})$ | $12.10881$ | $370$
$(0.2,1.667\cdot 10^{-5})$ | $10.48452$ | $460$
$(0.25,1.667\cdot 10^{-5})$ | $9.37379$ | $543$
$(0.5,1.667\cdot 10^{-5})$ | $6.63120$ | $889$
$(0.75,1.667\cdot 10^{-5})$ | $5.41887$ | $1168$
$(1.0,1.667\cdot 10^{-5})$ | $4.69244$ | $1409$
$(2.0,1.667\cdot 10^{-5})$ | $3.31648$ | $2159$
Table 7: Different privacy budget settings for non-convex problem for CIFAR$10$ data set Privacy budget $(\epsilon,\delta)$ | $\sigma$ |  Sample size $s$
---|---|---
$(0.25,2.0\cdot 10^{-5})$ | $9.29838$ | $417$
$(0.5,2.0\cdot 10^{-5})$ | $6.57192$ | $689$
$(0.75,2.0\cdot 10^{-5})$ | $5.36937$ | $909$
$(1.0,2.0\cdot 10^{-5})$ | $4.65014$ | $1099$
$(1.5,2.0\cdot 10^{-5})$ | $4.16111$ | $1267$
$(2.0,2.0\cdot 10^{-5})$ | $3.28831$ | $1690$
$(3.0,2.0\cdot 10^{-5})$ | $2.68273$ | $1994$
(a) Non-convex. Figure 13: Effect of different levels of privacy budgets
$\epsilon$ and non-DP settings (a) Non-convex. Figure 14: Effect of different
levels of privacy budgets $\epsilon$ and non-DP settings

These figures again confirm the effective performance of our DP-SGD framework,
which not only conserves strong privacy, but also keeps a decent convergence
rate to good accuracy, even for a very small privacy budget.

##  Appendix E Towards Using Proactive DP-SGD in Practice

Figure 15: The flow chart of our proactive (asynchronous) DP-SGD framework.

In order to apply our theory in practice, we introduce a flow chart of our
asynchronous DP-SGD learning from a client’s perspective in Figure 15. Each
client wants to participate in the collective learning of a global model that
achieves a sufficient target accuracy $acc_{target}$ with respect to their
test data set. That is, when a client tests the final global model against its
own private test set, then the client is satisfied if an accuracy
$acc_{target}$ is achieved. Such accuracy can only be achieved for certain
combinations of noise $\sigma$ and clipping constant $C$. In particular, the
final round introduces Gaussian noise with deviation $\sigma$ and this leads
to inherent inaccuracy of the final global model. Experiments in Section 4
have shown that reducing the clipping constant $C$ allows a larger $\sigma$
for attaining a good target accuracy $acc_{target}$ (if a sufficient larger
number $K$ of gradient computations have been executed). Here, we note that
$C$ cannot be reduced indefinitely without hurting accuracy; this is because a
reduced clipping constant into some extent plays the role of a reduced step
size (or learning rate) and we know that convergence to an accurate global
model must start with a large enough step size. Therefore, before starting any
learning, we need to understand how utility relates to $\sigma$ and the
clipping constant $C$. This is reflected in the functionality “Utility Graph”
in the flow chart of Figure 15. Based on learning a local model (and based on
a-priori information from learning models on similar data sets), “Utility
Graph” produces a set of pairs $\\{(C_{i},R_{i})\\}$ together with a total
number $K$ of gradient computations: For various clipping constants $C_{i}$, a
range $R_{i}$ of possible $\sigma$ is output. That is, if $\sigma\in R_{i}$
for clipping constant $C=C_{i}$, then there is a good indication that this
will lead to target accuracy $acc_{target}$ of the final model after the
client has contributed $K$ local gradient computations.

A second functionality in the flow chart is the “DP Calculator”. This
calculator takes the pair $(C_{i},R_{i})$ that allows the maximum possible
$\sigma$ in $R_{i}$. This $\sigma$ defines the best possible
$(\epsilon,\delta)$-DP curve through our formula
$\sqrt{2(\epsilon+\ln(1/\delta))/\epsilon}$. That is, it allows the “smallest"
possible $(\epsilon,\delta)$ pairs. The “DP Calculator” checks whether the
maximum $\sigma$ allows the target differential privacy budget of the client
defined by $(\epsilon_{target},\delta_{target})$. If not, then the client
cannot participate in the collective learning. If the target differential
privacy budget does fit, then, given $\sigma$, $\epsilon=\epsilon_{target}$,
$\delta=\delta_{target}$, $K$, and the client’s data set size $N$, the “DP
Calculator” computes the maximum possible sample size $s_{max}$ according to
the conditions in Theorem 4. This in turn results in a number of rounds
$T=K/s_{max}$. We propose to choose the maximum possible sample size as this
leads to the best accuracy/utility in our experiments. This is because
$s_{max}$ leads to the smallest number of rounds $T$, hence, the smallest
number of times noise is added and aggregated into the global model at the
server. (Also, as a secondary objective, a smaller number of rounds means less
round communication.)

As soon as the “DP Calculator” has calculated all parameters, the client
executes $T$ rounds of DP-SGD with sample (mini-batch) size $s=s_{max}$. This
is represented by the “FLwithDP” functionality in the flow chart of Figure 15.
Once the $T$ rounds are done, the client estimates the accuracy $acc$ of the
last received global model based on the client’s test data set. Next, the
client checks whether it is at least the target accuracy $acc_{target}$. If
so, then the client stops participating. That is, with each new global model
received by the server, the clients test whether the accuracy is satisfactory;
if not, then the client will want to participate again. If the test accuracy
$acc$ is not sufficient, then the “DP Calculator” will work with a new
$(C_{i},R_{i})$. The next subsection details the computations by the “DP
Calculator” and explains the feedback loop in the flow chart of Figure 15.

We remark that the client can use a complimentary differential privacy
accountant to keep track of the exact privacy budget that has been spent.

###  E.1 DP Calculator

We propose each local client to take control over its own privacy budget while
making sure the locally measured test accuracy of the final global model is
acceptable. The main idea is to start with an initial $\sigma=\sigma_{0}$ with
appropriate clipping constant $C=C_{0}$ and an estimated number $K=K_{0}$ of
local gradient computations needed for convergence to “sufficient” test
accuracy (utility). For local data set size $N$, we want to compute proper
parameter settings including the batch size $s=s_{0}$ for each round and the
total number of rounds $T=T_{0}$ (with $K=sT$). Once the $T$ rounds are
finished, local test data is used to compute the test accuracy of the final
global model. If the accuracy is not satisfactory, then $\sigma$ must be
reduced to a lower $\sigma_{1}$ (and we may re-tune to a larger clipping
constant $C_{1}$). This leads to an additional (estimated) number of local
$K_{1}$ gradient computations that need to be executed. The lower $\sigma$
corresponds to worse differential privacy since a lower $\sigma$ is directly
related to a higher $\epsilon$ for given $\delta=1/N$. The local client is in
control of what $\epsilon$ is acceptable – and if needed, the local client
simply stops participating helping the central server learn a global model.

In order to apply our theory, we pretend as if the initial $T$ rounds used the
lower $\sigma=\sigma_{1}$ – this means that our analysis provides an advantage
to the adversary as we assume less noise is used compared to what was
initially actually used. Hence, the resulting DP guarantee for $\sigma_{1}$
will hold for all $K=K_{0}+K_{1}$ local gradient computations. We use the new
$\sigma=\sigma_{1}$ and $K$ together with badge size $s_{0}$ for the first
$T_{0}$ rounds to compute a new parameter setting for the next rounds; this
includes the number $T_{1}$ of additional rounds (making $T=T_{0}+T_{1}$) and
the new batch size $s_{1}$. The new batch size implies a new average
$\bar{s}=(s_{0}T_{0}+s_{1}T_{1})/(T_{0}+T_{1})$ as well as a new variation
$\theta_{1}>1$ of the sequence of batch sizes.

Once all $T=T_{0}+T_{1}$ rounds are finished (or equivalently all
$K=K_{0}+K_{1}$ local gradient computations are finished), the local client
again computes the test accuracy of the last global model. If not acceptable
$\sigma$ is reduced again and we repeat the above process. If the test
accuracy is acceptable, then the local client stops participating, that is,
the local client stops gradient computations but continues to receive global
models from the central server. As soon as the local client measures a new
unacceptable local test accuracy, the client will continue the above process
and starts a new series of rounds based on $\sigma$.

Stopping participation and later continuing if needed best fits learning
problems over large data: Here, each local client samples its own local data
set according to the ‘client’s behavior’. The local client wants to prevent as
much leakage of its privately selected local data set as possible. Notice that
each local data set is too small for a local client to learn a global model on
its own – this is why local clients need to unite in a joint effort to learn a
global model (by using distributed SGD). Assuming all samples are iid (all
local data sets are themselves sampled from a global distribution), the final
global model is not affected by having more or less contribution from local
clients (as a result of different stopping and continuation patterns). Notice
that if local data sets would be heterogeneous, then the final global model
corresponds to a mix of all heterogeneous data sets and here it matters how
much each local client participates (as this influences the mix).

The above procedure describes a proactive method for adjusting $\sigma$ to
lower values if the locally measured test accuracy is not satisfactory. Of
course, the local client sets an a-priori upper bound $\epsilon_{target}$ on
the $\epsilon$ (with $\delta=\delta_{target}=1/N$), its privacy budget. This
privacy budget cannot be exceeded, even if the local test accuracy becomes
unsatisfactory.

We notice that our theory is general in that it can be used to analyse varying
sequences of batch sizes, which is needed for our proactive method. We now
describe in detail how to calculate parameter settings according to our
theorems:

Suppose the local client has already computed for $T_{0}+T_{1}+\ldots+T_{j-1}$
rounds with badge sizes $s_{0},s_{1},\ldots,s_{j-1}$, hence,
$K_{0}=s_{0}T_{0},K_{1}=s_{1}T_{1},\ldots,K_{j-1}=s_{j-1}T_{j-1}$. The local
client sets/fixes the total $K_{j}$ of gradient computations it wants to
compute over the next $T_{j}$ rounds. We want to compute a new $s_{j}$ and
$T_{j}$. Notice that $s_{j}=K_{j}/T_{j}$ and
$\bar{s}=\sum_{i=0}^{j}s_{i}T_{i}/T$, where $T=\sum_{i=0}^{j}T_{i}$. We want
to find a suitable $s_{j}$.

We start our calculation with $s_{j}=1$ and we rerun our calculation for
bigger batch sizes until we reach a maximum. Given a choice $s_{j}=s$, we
execute the following steps (we base the calculator on the slightly more
complex but more accurate Theorem 4 of Appendix B):

  1. 1.

Set $\delta=1/N$, compute $T_{j}=K_{j}/s_{j}$ given the input values $K_{j}$
and $s_{j}$, compute the corresponding $\bar{s}$ (see above) together with
corresponding $\theta=\max\\{s_{i}\\}/\bar{s}$. Compute
$K=\sum_{i=0}^{j}K_{i}$.

  2. 2.

Set $\gamma=2$ (because $\gamma=2+O(\bar{\alpha})$) as the initial value.

  3. 3.

According to Theorem 4, we compute $\epsilon$, $\sigma$, and $\bar{\alpha}$ as
follows:

     * •

Based on inequality (30), set $\epsilon$ as small as possible, that is,
$\epsilon=\gamma\theta^{2}\bar{s}\frac{K}{N^{2}}$.

     * •

We distinguish two cases:

$j=0$:

In case we want to determine $s_{0}$, we compute $\sigma_{0}=\sigma$ where
$\sigma$ meets (31) with equality, that is,
$\sigma=\sqrt{2(\epsilon+\log{1/\delta})/\epsilon}$.

$j>0$:

During previous computations we already selected a $\sigma_{j-1}$. As
described above, we only perform these calculations if the corresponding test-
accuracy is not satisfactory. For this reason we want a lower
$\sigma_{j}<\sigma_{j-1}$. The local client chooses a smaller $\sigma_{j}$
with possibly a larger clipping constant $C_{j}$ for which better accuracy
within $K_{j}$ local gradient computing steps is expected. We compute
$\epsilon$ as a solution of
$\sigma_{j}=\sqrt{2(\epsilon+\log{1/\delta})/\epsilon}$ and set $\epsilon$ to
the maximum of this solution and the minimal possible
$\epsilon=\gamma\theta^{2}\bar{s}\frac{K}{N^{2}}$ computed above.

     * •

Compute $\bar{\alpha}=\frac{\epsilon N}{\gamma K}$.

  4. 4.

Recompute the new
$\gamma_{new}=\frac{2}{1-\bar{\alpha}}+\frac{2^{4}\cdot\bar{\alpha}}{1-\bar{\alpha}}\left(\frac{\sigma}{(1-\sqrt{\bar{\alpha}})^{2}}+\frac{1}{\sigma(1-\bar{\alpha})-2e\sqrt{\bar{\alpha}}}\frac{e^{3}}{\sigma}\right)e^{3/\sigma^{2}}$.

  5. 5.

Repeat steps 3 and 4 with $\gamma$ replaced by $\gamma_{new}$ until
$\gamma_{new}-\gamma\leq 0.0001\gamma$, that is, $\gamma$ has converged
sufficiently.

  6. 6.

The resulting set of parameters $(\epsilon,\delta,\sigma,\gamma,\theta,K,N)$
can only be used if inequalities (28) and (29) are satisfied and
$\epsilon\leq\epsilon_{target}$.

     * •

If these conditions are satisfied, then we save the parameters
$(s,\epsilon,\sigma)$ and rerun the above calculation for bigger sample size
$s$. Otherwise, we output $(s,\epsilon,\sigma)$ of the previous run (as this
corresponds to the maximum $s$ and thus minimal number of communication
rounds) and terminate: We set $\sigma_{j}=\sigma$, $s_{j}=s$, and
$T_{j}=K_{j}/s_{j}$. Our theory proves that we satisfy
$(\epsilon,\delta=1/N)$-differential privacy.

     * •

It may be that even the minimal batch size $s_{j}=1$ does not result in valid
parameters $(s,\epsilon,\sigma)$. This means that the local client cannot
participate any more otherwise its required differential privacy guarantee
cannot be met.


# Model-based Reconstruction for Multi-Frequency Collimated Beam Ultrasound
Systems

Abdulrahman M. Alanazi1,4, Singanallur Venkatakrishnan2, Hector Santos-
Villalobos3, Gregery T. Buzzard1, and Charles Bouman1
1Purdue University-Main Campus, West Lafayette, IN 47907.
2Oak Ridge National Laboratory, One Bethel Valley Road, Oak Ridge, TN 37831.
3Amazon Prime Video, 410 Terry Ave N, Seattle 98109, WA.
4King Saud University (KSU), Riyadh, Saudi Arabia This manuscript has been
supported by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the
U.S. Department of Energy. G. Buzzard was partially supported by NSF
CCF-1763896, and C. Bouman was partially supported by the Showalter Trust. The
United States Government retains and the publisher, by accepting the article
for publication, acknowledges that the United States Government retains a non-
exclusive, paid-up, irrevocable, world-wide license to publish or reproduce
the published form of this manuscript, or allow others to do so, for United
States Government purposes. The Department of Energy will provide public
access to these results of federally sponsored research in accordance with the
DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

###### Abstract

Collimated beam ultrasound systems are a technology for imaging inside multi-
layered structures such as geothermal wells. These systems work by using a
collimated narrow-band ultrasound transmitter that can penetrate through
multiple layers of heterogeneous material. A series of measurements can then
be made at multiple transmit frequencies. However, commonly used
reconstruction algorithms such as Synthetic Aperture Focusing Technique (SAFT)
tend to produce poor quality reconstructions for these systems both because
they do not model collimated beam systems and they do not jointly reconstruct
the multiple frequencies.

In this paper, we propose a multi-frequency ultrasound model-based iterative
reconstruction (UMBIR) algorithm designed for multi-frequency collimated beam
ultrasound systems. The combined system targets reflective imaging of
heterogeneous, multi-layered structures. For each transmitted frequency band,
we introduce a physics-based forward model to accurately account for the
propagation of the collimated narrow-band ultrasonic beam through the multi-
layered media. We then show how the joint multi-frequency UMBIR reconstruction
can be computed by modeling the direct arrival signals, detector noise, and
incorporating a spatially varying image prior. Results using both simulated
and experimental data indicate that multi-frequency UMBIR reconstruction
yields much higher reconstruction quality than either single frequency UMBIR
or SAFT.

###### Index Terms:

Non-destructive testing (NDT), ultrasonic imaging, ultrasonic model-based
iterative reconstruction (UMBIR), multilayered objects, collimated beams,
multi-frequency.

##  I Introduction

Non-destructive evaluation (NDE) of multi-layered structures that can be
accessed from only a single side is important in many applications. For
example, this imaging scenario occurs when monitoring the structural integrity
of oil and geothermal wells that lie behind layers of fluid and steel casing.
While ultrasound imaging is widely used in NDE applications, multi-layered
structures present a challenge for ultrasounds systems because of the complex
propagation and reverberation of the signal through the material.

Figure LABEL:fig:MultiLayer_strucut_example illustrates an example of a
collimated beam ultrasound system [1, 2, 3, 4] that is designed to image
through multi-layered structures. The system consists of a narrow-band
collimated beam transmitter combined with an array of receivers. The system
can penetrate through heterogeneous layers because the transmitter is
collimated and the center frequency is below 100 kHz. However, since the
systems are narrow-band, they typically are operated at a few different center
frequencies, and then the data from each measurement frequency is processed
separately to image the structure.

The most popular methods to reconstruct data from ultrasound systems use a
delay-and-sum (DAS) approach because of their low computational complexity.
One such approach is the synthetic aperture focusing technique (SAFT), which
produces acceptable ultrasound images for simple objects [5]. SAFT has been
applied to single-layer [6, 7] and multi-layered structures [8, 9] but not to
collimated beam systems. Multi-layer SAFT combines DAS with techniques such as
ray-tracing [10, 11, 12] and root-mean-square velocity [13, 14] to compute the
travel time in multi-layered media. However, SAFT and its variations rely on a
simple model that often leads to artifacts such as multiple reflections and
blur. Methods to counteract these effects include [15] and [16], which use a
linear forward model for single-layer structures and a carefully constructed
sparse deconvolution approach. Finally, the conventional practice when
processing data obtained from multiple transmit frequency bands is to obtain
the reconstruction for each band separately and then visualize the results in
order to identify structures of interest.

More physically realistic inversion methods from seismology include least-
squares reverse time migration (LSRTM) [17, 18, 19] and full wave inversion
(FWI) [20]. LSRTM and FWI are iterative methods that seek the best least-
squares fit between observed and reconstructed data; a reflectivity image for
LSRTM and a velocity image for FWI. These methods have the capability to image
complex structures but rely on an iteration using a non-linear forward model,
making them computationally expensive and impractical for many imaging
applications. Finally, these methods are typically applied separately to data
from each frequency band which does not allow the end user to obtain a single
image corresponding to the structure to be imaged.

In order to reduce the reconstruction artifacts of SAFT while maintaining
computational efficiency, model-based iterative reconstruction (MBIR) can be
used with a linear propagation model. These methods combine a forward model
for the ultrasound measurement system with a prior-model/regularizer for the
unknown structure and cast the reconstruction as a maximum a-posteriori
estimation problem. In [15, 21], the forward model is designed to handle
plane-wave imaging. Recently, the MBIR approach of [22] used a propagation
model of the ultrasound through the media and combined all the data from the
source-detector pairs to jointly reconstruct a fully 3D image. However,
current MBIR approaches have limitations when imaging through heterogeneous
materials. In a single-layer structure, the time delay can be easily computed
using Snell’s law [22, 23]. However, in structures containing two layers, the
time delay no longer has a closed form that can be easily computed using
Snell’s law [24]. When the number of layers exceeds two, some techniques such
as ray-tracing or marching methods [25, 26, 27] can be used to approximate the
computations but these are computationally complex. In summary, no existing
regularized MBIR method accounts for multi-frequency collimated beam systems
used to image structures which are behind several layers.

In this paper, we propose a MBIR algorithm for multi-frequency collimated beam
ultrasound systems in order to image structures that are behind multiple
layers of materials. Our method, which we refer to as multi-frequency
ultrasound model-based iterative reconstruction (UMBIR), can be used to
accurately reconstruct heterogeneous structures behind multiple layers of
material excited at multiple frequencies. UMBIR does this by combining a novel
physics-based forward model that models multi-layered heterogeneous structures
while also coherently integrating data from of multiple frequencies to form a
single reconstructed cross section. In order to estimate the travel time of
the ultrasound signals through multiple media as a part of the forward model,
we introduce a efficient binary search-based method. Finally, the maximum a
posteriori (MAP) estimate is then computed using iterative coordinate descent
(ICD) optimization with a spatially varying prior similar to [23]. We note
that our work builds on preliminary ideas developed in the context of single
band measurements which was presented as a conference article [28]. Using
experimental and simulated data, we demonstrate that the proposed multi-layer
and multi-frequency UMBIR approach yields more accurate reconstructions with
higher spatial resolution and reduced artifacts when compared to both single
frequency UMBIR and SAFT.

##  II The UMBIR Forward Model

In order to image the structure which is separated from the source by multiple
layers of known materials, we design a ultrasound model-based iterative
reconstruction (UMBIR) approach. Assuming a linear system for simplicity, we
seek to reconstruct an image $x$ using a measurement model of the form

$y=Ax+Dg+w,$ (1)

where $y\in\mathbb{R}^{MK\times 1}$ is a vector of measurements from $K$
receivers at $M$ timepoints, $A\in\mathbb{R}^{MK\times N}$ is the system
matrix, $x\in\mathbb{R}^{N\times 1}$ is the vectorized version of the desired
image with $N$ total voxels, $D\in\mathbb{R}^{MK\times K}$ is a matrix whose
columns form a basis for the possible direct arrival signals,
$g\in\mathbb{R}^{K\times 1}$ is a scaling coefficient vector for $D$, and $w$
is a Gaussian random vector with distribution $N(0,\sigma^{2}I)$. From (1), we
will be able to formulate the reconstruction as the maximum a posteriori
estimate of $x$ and $g$ given $y$. However, to do this, we will first need to
introduce an acoustic model of propagation through the multi-layered material
that we can use to compute the matrices $A$ and $D$.

###  II-A Multi-Layer Acoustic Propagation Model

In this section, we introduce a model of the multi-layer acoustic propagation
based on an extension of the single-layer model used in [22]. Figure
LABEL:fig:MultiLay illustrates the problem. An acoustic signal is transmitted
from location $r^{t}_{o}$, reflected from the voxel, $v$, and then received by
one of $K$ possible microphones at location $r^{r}_{j}$. As the signal
propagates outward, it passes through $L$ different materials each with its
own acoustic velocity, $c_{l}$, in $\frac{m}{s}$, and attenuation coefficient,
$\alpha_{l}$, in $\frac{s}{m}$. Let $T^{t}_{\ell}(v)$, and
$T^{r}_{\ell,j}(v)$, denote that outgoing and returning propagation time of
the beam through each layer of the medium. Notice that both times will be
functions of the particular voxel, $v$, and the return time will also be a
function of the particular microphone, $j$. Then we model the frequency domain
transfer function as

$G_{j}(v,f)=\lambda(v)\prod_{\ell=1}^{L}\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}e^{-(c_{\ell}\alpha_{\ell}|f|+2\mathfrak{j}\pi
f)\left[T^{t}_{\ell}(v)+T^{r}_{\ell,j}(v)\right]}\ ,$ (2)

where $\mathfrak{j}^{2}=-1$ and

$\lambda(v)=\phi_{j}(v)\left(\prod_{\ell=2}^{L}\frac{2\zeta_{\ell}}{\zeta_{\ell-1}+\zeta_{\ell}}\right)\left(\prod_{\ell=1}^{L-1}\frac{2\zeta_{\ell}}{\zeta_{\ell}+\zeta_{\ell+1}}\right)\
,$ (3)

where $\zeta_{\ell}$ in $\frac{kg}{m^{2}s}$ is the acoustic impedance of the
$\ell^{\text{th}}$ layer and $\phi_{j}(v)$ is the scalar correction due to
beam collimation discussed in Section II-C. Note that $\lambda(v)$ in (2)
models the effects on the amplitude of the received signal by the beam
collimation and impedance mismatch between the layers, which are the formulas
defined between parentheses in (3). We can further simplify the expression,
but defining two quantities

$\displaystyle\gamma_{j}(v)$
$\displaystyle=\sum_{l=1}^{L}c_{\ell}\alpha_{\ell}\left[T^{t}_{\ell}(v)+T^{r}_{\ell,j}(v)\right]$
(4) $\displaystyle T_{j}(v)$
$\displaystyle=\sum_{l=1}^{L}\left[T^{t}_{\ell}(v)+T^{r}_{\ell,j}(v)\right]\
.$ (5)

In this case, the frequency domain transfer function has the simpler form of

$G_{j}(v,f)=\lambda(v)\,e^{-(\gamma_{j}(v)\,|f|+2\mathfrak{j}\pi
f\,T_{j}(v))}\ .$ (6)

In this form, it is clear that $T_{j}(v)$ represents the round-trip group
delay to voxel $v$, and $\gamma_{j}(v)$ represents the signal dispersion. From
this, the Fourier transform of the received signal given by

$Y_{j}(v,f)=x(v)\lambda(v)S(f)e^{-(\gamma_{j}(v)\,|f|+2\mathfrak{j}\pi
f\,T_{j}(v))}\ ,$ (7)

where $S(f)$ is the Fourier transform of the transmitted signal, and $x(v)$ is
the reflection coefficient for the voxel $v$. Then in the time-domain the
received signal is given by

$y_{j}(v,t)=x(v)\,\lambda(v)\,h(\gamma_{j}(v),t-T_{j}(v)),$ (8)

where

$h(\gamma,t)=\mathcal{F}^{-1}\left\\{S(f)e^{-\gamma|f|}\right\\}$ (9)

and $\mathcal{F}^{-1}$ is the inverse Fourier transform.

Thus, the output at each time $t$ and each receiver $j$ can be expressed as an
inner product between the input $x(v)$ and a row of the forward model’s system
matrix, $A$. This is expressed in the following formula.

$\displaystyle y_{j}(t)$
$\displaystyle=\sum_{v}\left[h(\gamma_{j}(v),t-T_{j}(v))\,\lambda(v)\right]x(v)$
(10)

So from this we see that, for the multi-layer case, we will need to first
compute, $T_{j}(v)$, the acoustic travel time of the signal through each layer
for each voxel in order to obtain a expression for the received signal. This
is the subject of the next section.

###  II-B Time Delay Computation

In single-layer structures, the computation of the time delays at the image
voxels is straightforward. However, in multilayered structures, the acoustic
speed varies with depth, which causes reflections and refractions that result
in a complex wave path. In order to compute the received signal of (10), we
need to know the time delays for each layer, $T^{t}_{\ell}(v)$ and
$T^{r}_{\ell,j}(v)$, which will be dependent on the path of the acoustic
signal through the multiple layers of the material. Figure LABEL:fig:MultiLay
illustrates the acoustic path as it passes through the layers of the medium.
Because each layer has different acoustic velocity, $c_{\ell}$, the signal
will be refracted as it passes between layers. Let $\theta^{t}_{\ell}$
represent the angle of outbound propagation through the $\ell^{th}$ material,
and let $\theta^{r}_{\ell,j}$ represent the angle of return propagation, as in
Figure LABEL:fig:MultiLay. Then by Snell’s law, we know that

$\theta^{t}_{\ell}=\text{sin}^{-1}\left(\text{sin}(\theta^{t}_{\ell-1})\frac{c_{\ell}}{c_{\ell-1}}\right)\
,$ (11)

and the return angles are given by

$\theta^{r}_{\ell-1}=\text{sin}^{-1}\left(\text{sin}(\theta^{r}_{\ell})\frac{c_{\ell}}{c_{\ell+1}}\right)\
.$ (12)

Consequently, if we know $\theta^{t}_{1}$ and $\theta^{r}_{L}$, then we can
use the recursions of (11) and (12) to compute the remaining angles. We denote
these two functions that represent the result of this calculation as

$\displaystyle\theta^{t}_{\ell}$ $\displaystyle=$ $\displaystyle
f^{t}_{\ell}\left[\theta^{t}_{1}\right]$ $\displaystyle\theta^{r}_{\ell}$
$\displaystyle=$ $\displaystyle f^{r}_{\ell}\left[\theta^{r}_{L}\right]\ .$

Let $\eta_{\ell}$ denote the thickness of the $\ell^{th}$ layer, then we can
express the vertical distance that the outbound acoustic signal travels as

$\displaystyle
Z^{t}\\!\left[\theta^{t}_{1}\right]=\sum_{\ell=1}^{L}\eta_{\ell}\text{tan}\left(f^{t}_{\ell}\left[\theta^{t}_{1}\right]\right)\
,$ (13)

and the vertical distance that the returning signal travels as

$\displaystyle
Z^{r}\\!\left[\theta^{r}_{L}\right]=\sum_{\ell=1}^{L}\eta_{\ell}\text{tan}\left(f^{r}_{\ell}\left[\theta^{r}_{L}\right]\right)\
.$ (14)

As an example, Figure 3 shows the relationship between $Z^{t}(\theta^{t}_{1})$
versus $\theta^{t}_{1}$ and $Z^{r}(\theta^{r}_{3})$ versus $\theta^{r}_{3}$ in
three layer media with parameters as shown in Table I. We note that it is
clear from this example, that $Z^{t}(\theta^{t}_{1})$ is a monotone increasing
function of $\theta^{t}_{1}$.

---
(a) (b)
Figure 3: (a) The vertical distance traveled by the outgoing acoustic signal, $Z^{t}(\theta^{t}_{1})$, as a function of $\theta^{t}_{1}$. (b) The vertical distance traveled by the returning signal, $Z^{r}(\theta^{r}_{3})$, as a function of $\theta^{r}_{3}$. Notice that both functions are monotone increasing.  TABLE I: The parameters used to show the relationship between the vertical distance and propagation angle in Figure 3. | Layer 1 | Layer 2 | Layer 3
---|---|---|---
Thickness ($\eta_{\ell}$) in m | 0.073 | 0.006 | 0.12
Acoustic velocity ($c_{\ell}$) in $\frac{m}{s}$ | 1500 | 2800 | 2620

Next, for each voxel, $v$, we must solve for the unknown angles,
$\theta^{t}_{1}(v)$, the departure angle of the outbound acoustic signal, and
$\theta^{r}_{L}(v)$, the arrival angle of the returning acoustic signal. We
can do this by solving the following two equations.

$\displaystyle Z_{v}-Z_{o}$ $\displaystyle=$ $\displaystyle
Z^{t}\left[\theta^{t}_{1}(v)\right]$ $\displaystyle Z^{r}_{j}-Z_{v}$
$\displaystyle=$ $\displaystyle Z^{r}\left[\theta^{r}_{L,j}(v)\right]\ ,$

where $Z_{v}-Z_{o}$ is the vertical distance between the transmitter,
$r^{t}_{o}$, and the voxel, $v$; and $Z^{r}_{j}-Z_{v}$ is the vertical
distance between the microphone, $r^{r}_{j}$, and the voxel, $v$. Since both
functions are monotone increasing functions of their arguments, these
equations can be easily solved using half interval search as described in
[29].

Once the values of $\theta^{t}_{1}$ and $\theta^{r}_{L}$ are determined, then
the values of the time delay can be computed using the following two
equations.

$\displaystyle T^{t}_{\ell}(v)$ $\displaystyle=$
$\displaystyle\eta_{\ell}\frac{\sqrt{1+\tan^{2}\left(f^{t}_{\ell}[\theta^{t}_{1}(v)]\right)}}{c_{\ell}}$
(15) $\displaystyle T^{r}_{\ell,j}(v)$ $\displaystyle=$
$\displaystyle\eta_{\ell}\frac{\sqrt{1+\tan^{2}\left(f^{r}_{\ell}[\theta^{r}_{L,j}(v)]\right)}}{c_{\ell}}\
.$ (16)

This provides all the values that are needed to compute the solution to (6).

###  II-C Collimated Beam Modeling

In this section, we develop a model for the collimated beam generated by the
acoustic source transducer used in our experiments. In order to accurately
model the effects of beam collimation, we use an apodization function inspired
by work of [7] as later adapted by [22]. Adapting this approach to better
model the collimated beam and defining the angle at the receiver as
$\theta^{r}_{1,j}=f_{1}^{r}[\theta^{r}_{L,j}(v)]$ our apodization function is
given by

$\phi_{j}(v)=\text{cos}^{\beta}\left(\theta^{t}_{1}(v)-\theta^{t}_{p}\right)\cos^{2}\left(\theta_{1,j}^{r}(v)\right)\
,$ (17)

where $\beta$ is a parameter that controls the beam apodization, and
$\theta^{t}_{p}$ is the pointing angle of the transmitter as shown in Figure
LABEL:fig:MultiLay. The reason behind modeling the pointing angle of the
transmitter is because in practice, the transmitter is often tilted upward
with respect to the axial direction of the microphones array to increase the
ultrasonic illumination range and ensure reflection of the transmitted signal
to the receivers [19].

###  II-D Single Frequency System Matrix Construction

In this section, we describe how the the matrices, $A$ and $D$, are
constructed from the multi-layer acoustic model for a single frequency as
described in Section II. In the following section, we will describe how these
system matrices are combined to form the full system matrix used in the multi-
frequency case.

In order to reduce computation, we window our model of the received signal in
time by replacing $h$ of (8) with

$\tilde{h}(\gamma_{j}(v),t)=h(\gamma_{j}(v),t)\,\text{rect}\left(\frac{t}{t_{0}}-\frac{1}{2}\right),$
(18)

where $t_{0}$ is a constant based on the assumption that
$h(\gamma_{j}(v),t)\approx 0$ for $t>t_{0}$ and $\text{rect}(u)=1$ for
$|u|<\frac{1}{2}$ and is $0$ otherwise. By populating $A$ and $D$ with the
windowed function, $\tilde{h}$, we ensure that the matrices are sparse so that
computation and memory usage are reduced.

Equation (10) can then be used to populate the entries of the system matrix
$A$. For each voxel, $v_{i}$, and receiver location, $r_{j}^{r}$, the
following partial column vector is formed

$\displaystyle
a^{i,j}=\left[\begin{array}[]{c}\tilde{h}(\gamma_{j}(v_{i}),t_{0})\\\
\vdots\\\ \tilde{h}(\gamma_{j}(v_{i}),t_{M-1})\end{array}\right]$ (22)

where $t_{m}=m\Delta+T_{o}$ and $\Delta$ is the time sampling period. A full
column vector of the system matrix, $A$, is then formed by concatenating the
partial columns for each receiver.

$\displaystyle A_{*,i}=\left[\begin{array}[]{c}a^{i,1}\\\ \vdots\\\
a^{i,K}\end{array}\right]$ (26)

And then the full system matrix is formed by concatenating the columns.

$\displaystyle A=\left[A_{*,1},\cdots,A_{*,N}\right]$ (27)

Note that Equation (10) can also be used to populate the entries of the matrix
$D$. However, the group delay and signal dispersion formulas introduced above
depend on the voxel $v$; which is not the case for the direct arrival signals.
Hence, we define the time delay from the transmitter, $r^{t}_{o}$, to each
receiver, $r^{r}_{j}$, and the direct arrival signal dispersion, respectively,
as

$\tau_{j}=\frac{\lVert r^{t}_{o}-r^{r}_{j}\rVert}{c_{d}}\ ,$
$\bar{\gamma_{j}}=\alpha_{d}c_{d}\tau_{j}\ ,$

where $\alpha_{d}$ and $c_{d}$ are the attenuation coefficient and acoustic
velocity of the material that the transducers are embedded in. The matrix $D$
is formed from $K$ columns, one for each detector location $r^{r}_{j}$. We
denote the columns of $D$ as $d^{j}$. Then, for the $k^{th}$ receiver, a
partial column vector, $d^{k}\in\Re^{M}$, is formed by

$\displaystyle
d^{k}=\left[\begin{array}[]{c}\tilde{h}(\bar{\gamma_{j}},t_{0})\\\ \vdots\\\
\tilde{h}(\bar{\gamma_{j}},t_{M-1})\end{array}\right]\ .$ (31)

These partial vectors can then be concatenated to form the full matrix given
by

$\displaystyle D=\left[\begin{array}[]{cccc}d^{1}&0&\cdots&0\\\
0&d^{2}&\cdots&0\\\ \vdots&\vdots&\vdots&\vdots\\\
0&0&\cdots&d^{K}\end{array}\right]$ (36)

###  II-E Multi-Frequency UMBIR

In order to do multi-frequency reconstruction, we must form a system matrix
that account for measurements at all frequencies simultaneously. Let $S$
denote the number of distinct excitation frequencies, and let $A^{s}$ and
$D^{s}$ denote the associated system matrix and direct arrival signal matrix
constructed using the methods described in Section II-D, and let $y_{s}$
denote the measurements associated measurements. Then we can form the full
measurement vector, $y$, full system matrix, $A$, and direct arrival signal
matrix, $D$, to be as follows

$\displaystyle y=\left[\begin{array}[]{c}y^{1}\\\ \vdots\\\
y^{S}\end{array}\right]\hskip 10.0ptA=\left[\begin{array}[]{c}A^{1}\\\
\vdots\\\ A^{S}\end{array}\right]\hskip
10.0ptD=\left[\begin{array}[]{cccc}D^{1}&\cdots&0\\\ \vdots&\ddots&\vdots\\\
0&\cdots&D^{S}\end{array}\right].$ (46)

##  III Prior Model of UMBIR

For the prior model, we adopt the q-generalized Gaussian Markov random field
(QGGMRF) from [30, 22]. With this design, the prior probability is

$\displaystyle p(x)=\frac{1}{z}\exp\left(-\sum_{\\{s,r\\}\in C}b_{s,r}\
\rho(x_{s}-x_{r})\right),$ (47)

where $z$ is a normalizing constant, $C$ is the set of pair-wise cliques,

$\displaystyle\rho(\Delta)=\frac{|\Delta|^{p}}{p\sigma_{s,r}^{p}}\left(\frac{|\frac{\Delta}{T\sigma_{g_{s,r}}}|^{q-p}}{1+|\frac{\Delta}{T\sigma_{s,r}}|^{q-p}}\right),$
(48) $\displaystyle\sigma_{s,r}$ $\displaystyle=$
$\displaystyle\sigma_{0}\sqrt{\nu_{s}\nu_{r}},$ (49) $\displaystyle\nu_{s}$
$\displaystyle=$ $\displaystyle
1+(\nu-1)*\left(\frac{d_{s}}{d_{max}}\right)^{a},$ (50)

where $\nu>0$, $a>0$, $d_{s}$ is the distance from the sensor assembly to
pixel $s$, and $d_{max}$ is the max of $d_{s}$ over all $s$. We use $1<p<q=2$
to insure convexity and continuity of first and second derivatives of the
prior model. The parameter $T$ is unit-less and controls the edge threshold.
The QGGMRF parameter $\nu$ is unit-less and can be adjusted to amplify
reflections at deeper regions if needed. Finally, by taking the negative log
of Equation (47), the prior model penalty function is given by

$\displaystyle-\log p(x)=\sum_{\\{s,r\\}\in C}b_{s,r}\
\rho(x_{s}-x_{r})+\text{constant}\ .$ (51)

##  IV MAP Estimation and the UMBIR Algorithm

In order to perform UMBIR reconstruction, we will need to estimate both the
image, $x$, and the direct arrival signals coefficients, $g$. Using the MAP
formulation, the multi-layer UMBIR reconstruction is then given by

$\displaystyle\left(\hat{x},\hat{g}\right)$
$\displaystyle=\arg\min_{(x,g)}\left\\{-\log p(y|x,g)-\log p(x)-\log
p_{g}(g)\right\\}\ .$ (52)

We will use an improper prior distribution for $g$ of the form

$-\log p_{g}(g)=\text{constant}\ ;$

the prior term for $x$ is given by (51) above; and the forward model term is
given by

$-\log p(y|x,g)=\frac{1}{2\sigma^{2}}\left\|y-Ax-
Dg\right\|_{2}^{2}+\text{constant}\ ,$ (53)

where the system matrix, $A$, and the direct arrival matrix, $D$, are
constructed as described in Section II above.

Putting this together, results in the

$\displaystyle(\hat{x},\hat{g})=\arg\min_{(x,g)}$
$\displaystyle\left\\{\rule{0.0pt}{20.0pt}\frac{1}{2\sigma^{2}}\left\|y-Ax-
Dg\right\|^{2}\right.$ $\displaystyle\left.+\sum_{\\{s,r\\}\in C}b_{s,r}\
\rho(x_{s}-x_{r})\right\\}.$ (54)

In order to compute the MAP estimate (54), we use Iterative Coordinate Descent
(ICD) algorithm with the majorization technique for the prior model as
described in [29].

|  |  |
---|---|---|---
(a) $\beta=1$ | (b) $\beta=4$ | (c) $\beta=8$ | (d) Measured
Figure 4:  Comparison of experimentally measured and simulated apodization functions. A simulated beam profile, $\phi(v)^{(\beta)}e^{-\alpha_{0}r(v)}$, using the apodization function of (17) for (a) $\beta=1$, (b) $\beta=4$, and (c) $\beta=8$. (d) Experimentally measured beam profile for a collimated source. Notice that the value of $\beta=8$ provides a reasonably accurate approximation to the true beam profile.  |
---|---
(a) | (b)
|
(c) | (d)
Figure 5: Illustration of experimental set up for concrete cylinder
experiment. (a) Picture of concrete cylinder that contains a central bore with
the acoustic imaging senor. Notice that a 50 mm notch in the concrete is
painted red. (b) Picture of acoustic imaging sensor removed from cylinder. (c)
Diagram of a cross section of the concrete cylinder with detailed
specification of distances. The blue rectangle indicates the region used to
generate synthetic data to evaluate the UMBIR algorithm. The red region
indicates the cross-section that is to be reconstructed. (d) Diagram of the
receiver for the acoustic imaging sensor with detailed specification of sensor
positions.

##  V Experimental Results

In this section, we present results using both synthetic and measured data
sets. The synthetic data was generated using the K-Wave simulation package
[31]. Both our synthetic and measured data experiments are designed to
evaluate the performance of the well-bore integrity inspection system designed
at the Los Alamos National Laboratory (LANL) and shown in Figure 5 and Figure
9. This system uses a collimated acoustic transmitter along with an array of
15 receiving transducers. More detailed about the designed transducer can be
found in [32].

###  V-A Methods

In order to validate our method, we performed two different experimental
studies: The concrete cylinder (CC) experiment and the granite block (GB)
experiment. In order to evaluate our method we also generate simulated data
that is similar to the CC experiment using the K-Wave software (referred to in
the results as CC-KWave). Table II provides all parameters used in our
reconstruction experiments for the three cases.

TABLE II: Parameter settings for concrete cylinder and granite block
reconstructions. Forward Model Parameters
---
Parameter | CC-KWave | CC-Exp | GB-Exp | Unit
Num. rows | 140 | 140 | 110 | -
Num. cols | 70 | 70 | 124 | -
Recon. resolution | 3 | 3 | 3 | mm
Sampling freq. | 2 | 2 | 5 | MHz
FOV height | 420 | 420 | 330 | mm
FOV depth | 210 | 210 | 370 | mm
$\beta$ | 8 | 8 | 8 | -
$\alpha_{\text{water}}$ | 2 | 2 | 2 | $1/\text{Hz m}$
$\alpha_{\text{concrete}}$ | 30 | 30 | 30 | $1/\text{Hz m}$
$\alpha_{\text{granite}}$ | - | - | 88 | $1/\text{Hz m}$
Prior Model Parameters
Parameter | CC-KWave | CC-Exp | GB-Exp | Unit
Num. iterations | 100 | 100 | 100 | -
$\sigma$ | 0.1 | 0.12 | 0.2 | Pascal
p | 1.1 | 1.1 | $1.1$ | -
q | 2.0 | 2.0 | $2.0$ | -
T | 0.01 | 1 | 0.001 | -
$\sigma_{0}$ | 2 | 2 | 2 | $m^{-2}$
$\nu$ | 10 | 10 | 5 | -
$a$ | 2 | 2 | 3 | -

Our forward model can account for the the shape of the collimated beam as
discussed in (17). Therefore, we have to pick the value of the $\beta$ in (17)
that best matches the type of profile used in the experimental data. In Figure
4, we visualize the effects of $\beta$ on beam spread and compare it with
measured data. In Figure 4(a-c), we plot
$\phi(v)^{(\beta)}e^{-\alpha_{0}r(v)}$ for $\beta$ = 1, 4, and 8, where
$\phi(v)^{(\beta)}$ is the apodization function from (17), $\alpha_{0}$ is the
attenuation coefficient in $m^{-1}$, $\theta(v)$ is the angle between beam
direction and $v$, and $r(v)$ is distance from source to $v$. Increasing
$\beta$ decreases the beam spread and makes it more collimated, with good
perceptual fit to the measured data in Fig. 4(d) when $\beta=8$.

|  |  |  |
---|---|---|---|---
(a) GT | (b) SAFT 29kHz | (c) SAFT 42kHz | (d) SAFT 58kHz | (e) GT
|  |  |  |
(f) MF-UMBIR | (g) UMBIR 29kHz | (h) UMBIR 42kHz | (i) UMBIR 58kHz | (j) MF-UMBIR
Figure 6:  Results using synthetic data generated with K-Wave for concrete
cylinder experiment. (a) ground truth without notch; (b), (c), and (d) single
frequency SAFT reconstructions at 29 kHz, 42 kHz, and 58 kHz without notch;
(e) ground truth with notch; (f) multi-frequency UMBIR reconstruction of
K-Wave data without notch; (g), (h), and (i) single frequency UMBIR
reconstructions at 29 kHz, 42 kHz, and 58 kHz without notch; (j) multi-
frequency UBMIR reconstruction of K-Wave data with notch. The red and green
dashed lines indicate the notch and backwall locations, respectively. Notice
that the two MF-UMBIR reconstructions of (f) and (j) accurately reconstruct
the location of the notch and the back wall.

###  V-B Concrete Cylinder Experiment

Figure 5 illustrates the experimental setup for the concrete cylinder
experiment that was performed at LANL [32]. The concrete cylinder is designed
to represent a concrete well bore with a central hole that contains the
acoustic imaging sensor.

Figure 5(a) shows that on one side of the concrete cylinder there is a notch
marked with red paint that is 50 mm deep and subtends an angle of
approximately $45^{\circ}$. Figure 5(c) is a detailed diagram showing the
dimensions and positions of all the components. Notice that the acoustic
imaging sensor is placed in the center of the cylinder with the receiver array
at the top and the acoustic transmitter at the bottom. The blue box in Figure
5(c) shows the region in which the K-Wave simulation is performed for the
generation of synthetic data, and the red box shows the region in which the
UMBIR reconstruction is computed.

Figure 5(b) shows the receiver array along with the collimated transmitter
hanging below. This entire assembly is positioned inside the bore hole at the
center of the concrete cylinder where it can be rotated by the computer-
controlled rotation system. Figure 5(d) is a diagram showing the position of
the transducers in the sensor array. We note that the transmitters and
receivers are immersed in water to facilitate acoustic transmission into and
out of the concrete. However, there is an vacuum barrier shown as a green line
in Figure 5(c) that is positioned between the transmitter and receivers aiming
to block the direct arrival signal.

Table III lists out the transmission parameters at the three different
excitation frequencies that were used. The data was collected at $5^{\circ}$
increments using a rotational span of $180^{\circ}$, with the sensor assembly
facing the middle of the notch at the rotational position of $90^{\circ}$.

29 kHz | 42.4 kHz | 58 kHz | MF-UMBIR
---|---|---|---
$0^{\circ}$ (no notch) |  |  |
$90^{\circ}$ (notch) |  |  |
Figure 7: Results of reconstructing experimental data collected from concrete cylinder. Column 1 to 3 of figure: Single frequency UMBIR reconstructions at excitation frequencies of 29 kHz, 42.4 kHz, 58 kHz and multi-frequency UMBIR, without and with the notch. Column 4 of figure: Multi-frequency UMBIR reconstruction without and with the notch. Notice that the multi-frequency UMBIR reconstruction results in much better localization and greater accuracy of the back wall location.  TABLE III: Transmit signal parameters used for the concrete cylinder experiment. Excitation Frequency | Duration | Pulse Shape
---|---|---
29 kHz | 200 $\mu$s | Tukey
42.4 kHz | 200 $\mu$s | Tukey
58 kHz | 50 $\mu$s | Tukey

####  V-B1 K-Wave Results

TABLE IV: K-Wave parameters used for concrete cylinder. Component | Dimensions | Units
---|---|---
Computational domain height | 855 | mm
Computational domain width | 362 | mm
Pixel pitch in x/y direction | 3 | mm
perfectly matched layer (PML) size | 20 | samples
water acoustic speed | 1.5 | km/s
Plexiglas acoustic speed | 2.82 | km/s
concrete acoustic speed | 2.62 | km/s
water density | 997 | kg/m3
Plexiglas density | 1180 | kg/m3
concrete density | 1970 | kg/m3

In order to better understand the concrete cylinder experiment, we first
perform reconstructions using synthetic data generated with the K-Wave
simulation using the parameters shown in Table IV. Table II lists the
parameter settings used in K-Wave data reconstructions under CC-KWave, and all
simulations are done using a perfectly matched layer (PML) that extends in all
directions from the outer boundary of the computational domain. The regions
beyond the top, bottom, and back wall are assumed to be air (yellow region) to
mimic the real data, and an isolating layer (void) was placed between the
source and receivers to partially block direct arrival signals.

Figure 6 shows a comparison of single frequency UMBIR, multi-frequency UMBIR,
and SAFT when reconstructing synthetically generated concrete cylinder data.
Figure 6(a) and (e) show the ground truth used to generate the data with and
without the notch, respectively. Without the notch, the back wall should be
located at a depth of 18.85 cm (shown with a red dotted line), and with the
notch it should be located at a depth of 23.85 cm (shown with a green dotted
line).

The single frequency SAFT reconstructions in Figure 6(b-d) show multiple
reflections from the back wall and other artifacts at each excitation
frequency, which leads to uncertainty in the estimated location of the back
wall. In contrast, the single frequency UMBIR reconstructions in Figure 6(g-i)
show that even single-frequency UMBIR reconstructions provide fewer artifacts
and more accurate localization of the back wall than SAFT.

The multi-frequency UMBIR reconstruction without the notch in Figure 6(f) is a
substantial improvement over the single frequency UMBIR reconstructions. The
same is true for the multi-frequency UMBIR reconstruction with the notch in
Figure 6(j). This demonstrates that joint processing of low and high
excitation frequencies results in substantially better reconstruction quality
than single-frequency reconstruction.

####  V-B2 Experimental Data Results

Figure 7 depicts a selection of cross-section reconstruction results using the
measured experimental data. The left 3 columns of Figure 7 show single
frequency UMBIR cross-section reconstructions for the excitation frequencies
29 kHz, 42.4 kHz, 58 kHz without and with the notch. Column 4 shows the multi-
frequency UMBIR reconstruction without and with the notch. We again note that
the back wall without and with the notch are located at 23.85 cm and 18.85 cm,
respectively. So from this we see that the joint reconstruction of the
multiple frequencies more accurately localizes of the back wall.

The panoramic reconstruction view of the concrete cylinder in Figure
LABEL:fig:Concrete-Cylinder-Real-Panoramic provides a more visually intuitive
interpretation of the multi-frequency UMBIR reconstruction. The panoramic
reconstruction is formed by combining the views from each measured angle (37
equi-spaced angles from $0^{\circ}$ to $180^{\circ}$) to form a 2-dimensional
horizontal cross-section at a fixed height of 27 cm. In this case, both the
back wall with and without the notch is shown as a dotted red line. Notice
that the multi-frequency reconstruction localization closely follows its true
location. There is some error in capturing the transition along the notch, but
the locations are captured accurately over a large range of angles, and there
are few reflection artifacts.

Table V gives approximate computation times for 100 iterations of single and
multi-frequency UMBIR using an optimized code running on an Intel(R) Core(TM)
i7 CPU E5-2603 0 @1.80 GHz, 32.00 GB RAM. The table lists out both the time
required to compute the system matrices, $A$ and $D$, and the time required to
reconstruct one 2D image for the concrete cylinder data set shown in Figure
6(i) and (j).

TABLE V: Computation time to obtain the UMBIR and MF-UMBIR reconstructions in Figure 6(i) and (j). Computation Time | UMBIR | MF-UMBIR
---|---|---
Time to compute system matrices | 26.6 s | 91.3 s
Time to perform reconstruction | 44.3 s | 205.2 s

###  V-C Granite Block Experiment

|  |  |
---|---|---|---
(a) | (b) | (c) | (d)
Figure 9: Set up for the Granite Block experiment. (a) Top-view photo of the
granite block used in the experiment. The two left and right defects are
approximately 21 cm away from the center of the borehole. (b) The
transmitter/receiver system used in experiment. (c) The dimensions of the
granite block and (d) A diagram for a cross-section. The red box indicates the
area reconstructed by UMBIR.

Figure 9 illustrates the experimental set up use of the granite block data set
described in this section. The granite block has a borehole with a steel
casing of 4" inner diameter and 4.5" outer diameter. The casing-cement
interface is located at 4.5" from the borehole center, and the cement-granite
interface is located at 5.25". The granite block has induced defects, as can
be seen in Figure 9(a), located at approximately 21 cm from the borehole
center.

Before imaging, the borehole was filled with water. Similar to the previous
experiment, there is one transmitter and 15 receivers. The transmitter is
aligned with the detectors and has a firing angle of $35^{\circ}$. The
receiver/transmitter system deployed to the granite block is shown in Figure
9(b).

Data was collected from 24 azimuth scans uniformly separated with an angle
$15^{\circ}$ to cover of $360^{\circ}$. Table VI lists out the transmission
parameters at the three different excitation frequencies that were used for
this experiment. Table II gives the parameter values used in reconstructions
of the granite block data, and the red box in Figure 9(c) shows the area
reconstructed by UMBIR.

TABLE VI: Transmit signal parameters used for the granite block experiment. Excitation Frequency | Duration | Pulse Shape
---|---|---
103.38 kHz | 150 $\mu$s | Tukey
162.316 kHz | 50 $\mu$s | Tukey
220.554 kHz | 50 $\mu$s | Tukey
103.38 kHz | 162.316 kHz | 220.554 kHz | MF-UMBIR |
---|---|---|---|---
$0^{\circ}$ (no defect) |  |  |  |
$90^{\circ}$ (defect) |  |  |  |
Figure 10: UMBIR and multi-frequency UMBIR reconstruction results using
experimentally measured data from the granite block (GB) experiment. Left to
right columns correspond to the excitation frequencies 103.38 kHz, 162.316
kHz, 220.554 kHz, and the multi-frequency joint reconstructions, respectively.
The first row of reconstruction corresponds to $0^{\circ}$ with no defect, and
the second row corresponds to $90^{\circ}$ with a defect. The red doted line
indicates the location of the defect. Notice that the multi-frequency UMBIR
reconstructions shows substantial improvements over single-frequency
reconstructions and that it more accurately localizes the position of the
defect and then wall.

Figure 10 shows cross-section reconstructions from the granite block
experimental data, at rotational positions $0^{\circ}$ and $90^{\circ}$, with
the defect in view at $90^{\circ}$ and not at $0^{\circ}$. Notice that the
multi-frequency UMBIR reconstruction at $0^{\circ}$ has a clear reflection
from the sidewall around 31 cm (the black dashed line), while the
reconstruction $90^{\circ}$ has a clear reflection from the defect around 21
cm (the red dashed line).

Figure LABEL:fig:LANL_II_realResults_slice shows the multi-frequency UMBIR
panoramic reconstruction at a fixed height of 25.8 cm in polar coordinates.
The two red circles indicate the defect locations, and the black dashed square
shows the wall location. Notice that the multi-frequency UMBIR reconstruction
shows clear reflections from the edges of the specimen around the wall
location along with the defects at $90^{\circ}$ and $270^{\circ}$. The
distance from center to each corner is underestimated, which we speculate is
due to multiple reflections near the corners.

##  VI Conclusion

In this paper, we proposed a multi-layer, multi-frequency collimated
ultrasound model-based iterative reconstruction (UMBIR) algorithm. To do this,
we introduced a computationally efficient method for computing an accurate
forward model system matrix for multi-layered structures that can typically
occur in practical ultrasound imaging scenarios using collimated beam
ultrasonic transducers. We also formulated the reconstruction method using MAP
reconstruction with space-varying image prior along with model of the direct
arrival signal.

We tested our method on both simulated and experimentally measured data using
two different scenarios corresponding to a concrete cylinder and a granite
block. Both scenarios were designed to represent the imaging of defects in a
well-bore.

In all cases, we found that multi-frequency UMBIR reconstructions had
substantially better quality than single-frequency UBMIR reconstructions, and
that single-frequency UMBIR reconstructions had substantially better quality
than SAFT reconstructions. Multi-frequency UMBIR reconstructions had much
better localization of image defects. In addition, the multi-frequency UMBIR
reconstructions accurately detected the location of defects and object back
walls.

## Acknowledgment

A. M. Alanazi was supported by King Saud University. C. A. Bouman was
partially supported by the Showalter Trust and by the U.S. Department of
Energy. C. A. Bouman and G.T. Buzzard were partially supported by NSF
CCF-1763896. S.V. and Hector Santos-Villalobos were supported by the U.S.
Department of Energy staff office of the Under Secretary for Science and
Energy under the Subsurface Technology and Engineering Research, Development,
and Demonstration (SubTER) Crosscut program, and the office of Nuclear Energy
under the Light Water Reactor Sustainability (LWRS) program.

## References

  * [1] V. K. Chillara, C. Pantea, and D. N. Sinha, “Low-frequency ultrasonic bessel-like collimated beam generation from radial modes of piezoelectric transducers,” _Applied Physics Letters_ , vol. 110, no. 6, p. 064101, 2017\.
  * [2] V. K. Chillara, C. Pantea, and D. Sinha, “Radial modes of laterally stiffened piezoelectric disc transducers for ultrasonic collimated beam generation,” _Wave Motion_ , vol. 76, pp. 19–27, 2018.
  * [3] V. K. Chillara, E. S. Davis, C. Pantea, and D. N. Sinha, “Collimated acoustic beams from radial modes of piezoelectric disc transducers,” in _AIP Conference Proceedings_ , vol. 2102, no. 1. AIP Publishing LLC, 2019, p. 040013.
  * [4] V. K. Chillara, J. Greenhall, and C. Pantea, “Ultrasonic waves from radial mode excitation of a piezoelectric disc on the surface of an elastic solid,” _Smart Materials and Structures_ , vol. 29, no. 8, p. 085002, 2020.
  * [5] D. Prine, “Synthetic aperture ultrasonic imaging,” in _Proceedings of the Engineering Applications of Holography Symposium, Los Angeles, CA, USA_ , vol. 1617, 1972.
  * [6] T. Stepinski, “An implementation of synthetic aperture focusing technique in frequency domain,” _ieee transactions on ultrasonics, ferroelectrics, and frequency control_ , vol. 54, no. 7, pp. 1399–1408, 2007.
  * [7] K. Hoegh and L. Khazanovich, “Extended synthetic aperture focusing technique for ultrasonic imaging of concrete,” _NDT & E International_, vol. 74, pp. 33–42, 2015.
  * [8] M. H. Skjelvareid, T. Olofsson, Y. Birkelund, and Y. Larsen, “Synthetic aperture focusing of ultrasonic data from multilayered media using an omega-k algorithm,” _IEEE transactions on ultrasonics, ferroelectrics, and frequency control_ , vol. 58, no. 5, pp. 1037–1048, 2011.
  * [9] S. Lin, S. Shams, H. Choi, and H. Azari, “Ultrasonic imaging of multi-layer concrete structures,” _NDT & E International_, vol. 98, pp. 101–109, 2018\.
  * [10] V. Cerveny, _Seismic ray theory_. Cambridge university press, 2005.
  * [11] G. F. Margrave and M. P. Lamoureux, _Numerical methods of exploration seismology: With algorithms in MATLAB®_. Cambridge University Press, 2019.
  * [12] A. Shlivinski and K. Langenberg, “Defect imaging with elastic waves in inhomogeneous–anisotropic materials with composite geometries,” _Ultrasonics_ , vol. 46, no. 1, pp. 89–104, 2007.
  * [13] J. Gazdag, “Wave equation migration with the phase-shift method,” _Geophysics_ , vol. 43, no. 7, pp. 1342–1351, 1978.
  * [14] M. H. Skjelvareid and Y. Birkelund, “Ultrasound imaging using multilayer synthetic aperture focusing,” in _Pressure Vessels and Piping Conference_ , vol. 49248, 2010, pp. 379–387.
  * [15] E. Ozkan, V. Vishnevsky, and O. Goksel, “Inverse problem of ultrasound beamforming with sparsity constraints and regularization,” _IEEE transactions on ultrasonics, ferroelectrics, and frequency control_ , vol. 65, no. 3, pp. 356–365, 2017.
  * [16] A. Tuysuzoglu, J. M. Kracht, R. O. Cleveland, M. C¸ etin, and W. C. Karl, “Sparsity driven ultrasound imaging,” _The Journal of the Acoustical Society of America_ , vol. 131, no. 2, pp. 1271–1281, 2012.
  * [17] X. Liu, Y. Liu, X. Huang, and P. Li, “Least-squares reverse-time migration with cost-effective computation and memory storage,” _Journal of Applied Geophysics_ , vol. 129, pp. 200–208, 2016.
  * [18] Y. Zhang, L. Duan, and Y. Xie, “A stable and practical implementation of least-squares reverse time migration,” _Geophysics_ , vol. 80, no. 1, pp. V23–V31, 2015.
  * [19] Y. Chen, K. Gao, E. S. Davis, D. N. Sinha, C. Pantea, and L. Huang, “Full-waveform inversion and least-squares reverse-time migration imaging of collimated ultrasonic-beam data for high-resolution wellbore integrity monitoring,” _Applied Physics Letters_ , vol. 113, no. 7, p. 071903, 2018\.
  * [20] C. Zeng, S. Dong, and B. Wang, “A guide to least-squares reverse time migration for subsalt imaging: Challenges and solutions,” _Interpretation_ , vol. 5, no. 3, pp. SN1–SN11, 2017.
  * [21] H. Wu, J. Chen, S. Wu, H. Jin, and K. Yang, “A model-based regularized inverse method for ultrasonic b-scan image reconstruction,” _Measurement Science and Technology_ , vol. 26, no. 10, p. 105401, 2015.
  * [22] H. Almansouri, S. Venkatakrishnan, C. Bouman, and H. Santos-Villalobos, “Model-based iterative reconstruction for one-sided ultrasonic nondestructive evaluation,” _IEEE Transactions on Computational Imaging_ , vol. 5, no. 1, pp. 150–164, 2018.
  * [23] H. Almansouri, S. Venkatakrishnan, D. Clayton, Y. Polsky, C. Bouman, and H. Santos-Villalobos, “Anisotropic modeling and joint-map stitching for improved ultrasound model-based iterative reconstruction of large and thick specimens,” in _AIP Conference Proceedings_ , vol. 1949 (1). AIP Publishing LLC, 2018, p. 030002.
  * [24] M. Weston, P. Mudge, C. Davis, and A. Peyton, “Time efficient auto-focussing algorithms for ultrasonic inspection of dual-layered media using full matrix capture,” _Ndt & E International_, vol. 47, pp. 43–50, 2012.
  * [25] T. Moser, “Shortest path calculation of seismic rays,” _Geophysics_ , vol. 56, no. 1, pp. 59–67, 1991.
  * [26] J. A. Sethian, “A fast marching level set method for monotonically advancing fronts.” _Proceedings of the National Academy of Sciences_ , vol. 93, no. 4, pp. 1591–1595, 1996.
  * [27] A. J. Brath and F. Simonetti, “Phased array imaging of complex-geometry composite components,” _IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control_ , vol. 64, no. 10, pp. 1573–1582, 2017\.
  * [28] A. Alanazi, S. Venkatakrishnan, H. Santos-Villalobos, G. Buzzard, and C. Bouman, “Model-based reconstruction for collimated beam ultrasound systems,” in _ICASSP 2022-2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2022, pp. 1601–1605.
  * [29] C. A. Bouman, _Foundations of Computational Imaging: A Model-Based Approach_. SIAM, 2022, vol. 180.
  * [30] J.-B. Thibault, K. D. Sauer, C. A. Bouman, and J. Hsieh, “A three-dimensional statistical approach to improved image quality for multislice helical ct,” _Medical physics_ , vol. 34, no. 11, pp. 4526–4544, 2007.
  * [31] B. E. Treeby and B. T. Cox, “k-wave: Matlab toolbox for the simulation and reconstruction of photoacoustic wave fields,” _Journal of biomedical optics_ , vol. 15, no. 2, p. 021314, 2010.
  * [32] C. Pantea, “Collimated beams for cement evaluation,” Los Alamos National Lab.(LANL), Los Alamos, NM (United States), Tech. Rep., 2019.