vibro-acoustography]On the inverse problem of vibro-acoustography

[1]Barbara<EMAIL_ADDRESS>
The aim of this paper is to put the problem of vibroacoustic imaging into the mathematical framework of inverse problems (more precisely, coefficient identification in PDEs) and regularization. We present a model in frequency domain, prove uniqueness of recovery of the spatially varying nonlinearity parameter from measurements of the acoustic pressure at multiple frequencies, and derive Newton as well as gradient based reconstruction methods.

[MSC Classification]35R30,65J20

§ INTRODUCTION

Vibro-acoustography by means of ultrasound  was developed [1, 2] to achieve
the high resolution by high frequency waves while avoiding the drawbacks of scattering from small inclusions and of stronger attenuation at higher frequencies.

schematic of the experimental setup

The experiment for image acquisition is illustrated in Figure <ref>:

Two ultrasound beams of high and slightly different frequencies $\omega_1$ and $\omega_2$ are excited at two parts $\Sigma_1$, $\Sigma_2$ of an array of piezoelectric  transducers (emitters). They interact nonlinearly at a focus and this interaction excites a wave that basically propagates at the difference frequency $\omega_1-\omega_2$ and is eventually measured by a receiver array $\Gamma$ (in experiments consisting of hydrophones, in imaging this would be a piezoelectric transducer array as well).
After each measurement, the focal region is shifted to scan the overall region of interest.
Inhomogeneity of the medium leads to spatial dependence of two coeffients in the governing models: The speed of sound $c=c(x)$ and the nonlinearity parameter $\gamma=\gamma(x)$.
Both parameters are susceptible to local variations in the acoustic medium (e.g., human tissue in medical applications) and thus their reconstruction yields a spatial image of the region of interest.
In case of reconstructing the individual coefficients $c=c(x)$ or $\gamma=\gamma(x)$, this is related to ultrasound tomography and nonlinearity parameter imaging, respectively, cf. e.g., [3, 4, 5, 6] and the citing literature.

A modeling and simulation framework for this methodology has been devised in
[7, 8].
In this paper we put an emphasis on the inverse problem of reconstructing $c=c(x)$ and $\gamma=\gamma(x)$.

§ MODEL

Two ultrasound beams with acoustic velocity potentials $\phi_1$, $\phi_2$ are excited by transducers

and their interaction in turn excites a wave field with velocity potential $\psi$. This is discribed by a system of PDEs with inhomogeneous Neumann conditions

\begin{eqnarray}
&&\partial_t^2 \phi_k-c^2\Delta \phi_k=0 \mbox{ in }\Omega,\nonumber\\
&&\partial_\nu \phi_k=g_k \mbox{ on } \Sigma_k\,,\ k\in\{1,2\}\label{eq:phik}\\
&&\partial_t^2 \psi-c^2\Delta \psi= \tilde{f}(\phi_1,\phi_2,\gamma,c) \mbox{ in }\Omega
\label{eq:psi}\\
&&\partial_t\left(\left\vert\nabla(\phi_1+\phi_2)\right\vert^2 + \frac{\gamma-1}{2c^2}
\left\vert\partial_t(\phi_1+\phi_2)\right\vert^2\right)
\nonumber
\end{eqnarray}

see [9] for the derivation of the nonlinear forcing $f$.
In here, $c=c(x)$ and $\gamma=\gamma(x)$ are the spatially varying sound speed and nonlinearity parameter, respectively, and the manifold $\Sigma_k$ represents the emitting transducer array with given time harmonic excitation $g_k(x,t)=\hat{g}_k(x)\, e^{\imath\omega_k t}$.

The system (<ref>), (<ref>) is not fully nonlinear but the task of its solution can be decoupled into two linear subproblems: First compute $\phi_1,\phi_2$ from (<ref>), then insert them into the right hand side of (<ref>), and finally solve (<ref>) for $\psi$.

Transformation into frequency domain
Linearity of the subproblems allows to easily transfer the time domain formulation (<ref>), (<ref>) into frequency domain.
With the time harmonic ansatz $\phi_k(x,t)=\hat{\phi}_k(x)\, e^{\imath\omega_k t}$, $\psi(x,t)=\Re\Bigl(\hat{\psi}(x)\, e^{\imath(\omega_1-\omega_2)t}\Bigr)$, where the latter is induced by real-valuedness of the right hand side of (<ref>)

\begin{eqnarray*}
%&&\partial_t\left[\left\vert\nabla\left(\hat{\phi}_1\, e^{\imath\omega_1t}+\hat{\phi}_2\, e^{\imath\omega_2t}\right)\right\vert^2 + \frac{\gamma(x)-1}{2c(x)^2}
%\left\vert\partial_t \left(\hat{\phi}_1\, e^{\imath\omega_1t}+\hat{\phi}_2\, e^{\imath\omega_2t}\right)\right\vert^2\right]\\
+2\Re\left(\nabla\hat{\phi}_1\cdot\nabla\overline{\hat{\phi}_2}\,  e^{\imath(\omega_1-\omega_2)t}\right)\\
&&\qquad\qquad+2\omega_1\omega_2\Re\left(\hat{\phi}_1\overline{\hat{\phi}_2}\, e^{\imath(\omega_1-\omega_2)t}\right)
\Bigr)\Bigr]\\
&&\ \cdot\Re\Bigl(\imath\left(
\nabla\hat{\phi}_1\cdot\overline{\nabla\hat{\phi}_2}
+\omega_1\omega_2 \frac{\gamma(x)-1}{2c(x)^2}\hat{\phi}_1\overline{\hat{\phi}_2}\right) e^{\imath(\omega_1-\omega_2)t}\Bigr),
\end{eqnarray*}

we get

\begin{eqnarray}
&&-\frac{\omega_k^2}{c(x)^2} \hat{\phi}_k-\Delta \hat{\phi}_k=0 \mbox{ in }\Omega \nonumber\\
&&\partial_\nu \hat{\phi}_k=\hat{g}_k \mbox{ on } \Sigma_k\,,\quad k\in\{1,2\}
\label{eq:phikhat0}\\
-\frac{(\omega_1-\omega_2)^2}{c(x)^2} \hat{\psi}-\Delta \hat{\psi}=
f(\hat{\phi}_1,\hat{\phi}_2,\gamma,c)\mbox{ in }\Omega\quad\label{eq:psihat0}\\
\nabla\hat{\phi}_1\cdot\overline{\nabla\hat{\phi}_2}
+\omega_1\omega_2 \frac{\gamma(x)-1}{2c(x)^2}\hat{\phi}_1\overline{\hat{\phi}_2}\Bigr)\nonumber
\end{eqnarray}

which nicely illustrates the physical fact that the propagating wave described by $\psi$ is concentrated at the difference frequency $\omega_1-\omega_2$.

We mention in passing that in fact also in the harmonic ansatz for $\phi_1$, $\phi_2$ taking the real part would be demanded by physics. This would lead to certain (actually higher frequency) correction terms, that we neglect here, though, as they are not relevant for reconstructions.

Boundary conditions
We consider a bounded computational domain $\Omega$, where the excitation surfaces $\Sigma_k$ are part of the boundary $\Sigma_k\subseteq\partial\Omega$ and the rest of $\partial\Omega$ is subject to impedance
boundary conditions in order to damp reflected waves

\[
%c\partial_\nu \hat{\phi}_k=-\imath\omega_k \hat{\phi}_k  \mbox{ on }\partial\Omega\setminus\Sigma_k\,,
%\quad c\partial_\nu \hat{\psi}=-\imath(\omega_1-\omega_2) \hat{\psi} \mbox{ on }\partial\Omega\,.
\partial_\nu \hat{\phi}_k=-\imath\sigma_k \hat{\phi}_k  \mbox{ on }\partial\Omega\setminus\Sigma_k\,,
\quad \partial_\nu \hat{\psi}=-\imath\sigma \hat{\psi} \mbox{ on }\partial\Omega\,.
\]

with nonnegative $L^\infty$ impedance coefficients $\sigma,\sigma_k$ that are bounded away from zero on an open subset of $\partial\Omega$ or $\partial\Omega\setminus\Sigma_k$, respectively.
Note that with the choice
these would be first order absorbing boundary conditions; however, later on in the definition of the operator $\mathcal{A}_c$ we wish to avoid explicit frequency dependence.
An alternative scenario that allows to work on a bounded domain $\Omega$ as well is to use a perfectly matched layer PML, (see. e.g., <cit.> and the references therein,) making the replacements
\Delta\, \leftrightarrow \, \nabla\cdot(D(x)\nabla)$,
$\frac{\omega_k^2}{c(x)^2} \, \leftrightarrow \, d(x) \frac{\omega_k^2}{c(x)^2}
in the above Helmholtz equations on an augmented domain
$\Omega=\Omega_{\rm acou}\cup\Omega_{\rm PML}$ with space dependent coefficients $D$ (matrix valued) and $d$;

for details see, e.g., [10].
Note that the real part of $D$  and $d$ is close to unity also in the PML region.
The boundary condition on the outer boundary can then be set to homogeneous Neumann.

The pressure data taken at the receiver array can, via the identity

\[
\varrho\partial_t\psi=p
\]

be expressed by an observation operator

\begin{equation}\label{eq:C}
C:(\hat{\phi}_1,\hat{\phi}_1,\hat{\psi})\mapsto \imath(\omega_1-\omega_2)\mbox{tr}_\Gamma\hat{\psi}
\end{equation}

where $\Gamma$ is a manifold representing the receiver array and lying inside the acoustic domain $\Omega$.

Inverse problem
The inverse problem of vibro-acoustography consists of determining the spatially varying coefficients $c$ and $\gamma$ from observations (<ref>) of the low frequency wave field.
We assume that $c$ is known on the outer boundary and needs to be reconstructed only in a subdomain (region of interest) $\widetilde{\Omega}\subseteq\Omega$ of the computational domain. With a slight abuse of notation we write

\[\frac{1}{c^2}=:\kappatil=\kappa_0+\chi_{\widetilde{\Omega}}\kappa\,, \quad \frac{\gamma-1}{2c^4}=:\gammatil=\gamma_0+\chi_{\widetilde{\Omega}}\gamma\]

where the background $\kappa_0$, $\gamma_0\in L^\infty(\widetilde{\Omega})$, $\kappa_0(x)\geq\underline{\kappa}>0$, and the subdomain $\tilde{\Omega}$ are known and $\chi_{\widetilde{\Omega}}$ is the extension by zero operator from $\widetilde{\Omega}$ to $\Omega$, defined by $(\chi_{\widetilde{\Omega}}\kappa)(x)=\kappa(x)$ for $x\in\widetilde{\Omega}$ and zero else.

Therewith our aim is to recover $\kappa,\gamma\in L^2(\widetilde{\Omega})$, in the weak form of (<ref>), (<ref>)

\begin{equation} \label{eq:A0}
\begin{aligned}
&0=\langle A(\hat{\phi}_1,\hat{\phi}_2,\hat{\psi},\kappa,\gamma),(v_1,v_2,w)\rangle:=\\
&\Re\Bigl((1-\imath)\Bigl(\sum_{k=1}^2\int_{\Omega} (-\omega_k^2\kappatil\hat{\phi}_k\overline{v}_k+\nabla\hat{\phi}_k\cdot\nabla \overline{v}_k)\, dx\\
&+\int_{\Omega} (-(\omega_1-\omega_2)^2\kappatil\hat{\psi}\overline{w}+\nabla\hat{\psi}\cdot\nabla \overline{w})\, dx\\
&+\imath\int_{\partial\Omega\setminus\Sigma}\sigma_k\hat{\phi}_k \overline{v}_k \, ds
-\int_\Sigma \hat{g}_k \overline{v}_k\, ds
+\imath\int_{\partial\Omega}\sigma\hat{\psi} \overline{w} \, ds\\
&-2(\omega_1-\omega_2) \imath \\
&\quad\cdot\int_\Omega \left(
\kappatil\nabla\hat{\phi}_1\cdot\overline{\nabla\hat{\phi}_2}
+\omega_1\omega_2 \gammatil \hat{\phi}_1\overline{\hat{\phi}_2}\right)\,\overline{w}\, dx \Bigr)\Bigr)
\\&\mbox{for all } v_1,v_2,w\in H^1(\Omega;\mathbb{C})
\end{aligned}
\end{equation}

(it suffices to take the real part here since $v_1,v_2,w$ vary over complex valued functions).
This is the weak form of

\begin{eqnarray}
&-\omega_k^2\kappatil \hat{\phi}_k-\Delta \hat{\phi}_k=0 \mbox{ in }\Omega
\nonumber\\
&\partial_\nu \hat{\phi}_k=\hat{g}_k \mbox{ on } \Sigma_k\,,\quad k\in\{1,2\}
\label{phihatk}\\
&-(\omega_1-\omega_2)^2\kappatil \hat{\psi}-\Delta \hat{\psi}=f(\kappatil,\gammatil,\hat{\phi}_1,\hat{\phi}_2)\mbox{ in }\Omega
\label{psihat}\\
& \quad f(\kappatil,\gammatil,\hat{\phi}_1,\hat{\phi}_2)=2\imath(\omega_1-\omega_2)
\nonumber\\
\nabla\hat{\phi}_1\cdot\overline{\nabla\hat{\phi}_2}
+\omega_1\omega_2 \gammatil\hat{\phi}_1\overline{\hat{\phi}_2}\Bigr)
\nonumber
\end{eqnarray}

with homogeneous impedance boundary conditions on (the rest of) $\partial\Omega$, which we also tacitly assume to hold in the following.

In Section <ref> we will prove that for every $\kappa,\gamma\in L^2(\widetilde{\Omega})$, there exists a unique solution $\hat{\phi}_1,\hat{\phi}_2,\hat{\psi}$ of the operator equation $A(\hat{\phi}_1,\hat{\phi}_2,\hat{\psi},\kappa,\gamma)=0$ in appropriate function spaces, such that also the observation operator $C$ according to (<ref>) can be applied and yields an element of $L^2(\Gamma)$.

This justifies the use of the function spaces

\begin{equation}\label{defXY}
X=L^2(\widetilde{\Omega})\times L^2(\widetilde{\Omega})\,, \quad  Y=L^2(\Gamma),
\end{equation}

to define the forward operator

\begin{equation}\label{defF}
\begin{aligned}
&F:X\to Y, \quad \kappagamma\mapsto F(\kappagamma)=\imath(\omega_1-\omega_2)\mbox{tr}_\Gamma\hat{\psi}\\ &\mbox{ where $\phi_1,\phi_2,\psi$ solve \eqref{eq:A0} }
\end{aligned}
\end{equation}

and write the inverse problem in reduced form as

\begin{equation}\label{Fxy}
\end{equation}

where $y\in L^2(\Gamma)$ is the pressure distribution measured at the receiver array.
Concerning the choice of spaces (<ref>), working in $L^2$ spaces makes definition of methods most convenient. This is on one hand due to their Hilbert space structure, on the other hand due to the fact that no derivatives are involved, which avoids having to solve additional PDEs for evaluating the adjoint operator.

Alternatively, using the model and observation operators $A$ and $C$ defined in (<ref>), (<ref>), we may write the inverse problem as an all-at-once system for the parameters $\kappagamma$ and the states $u:=(\hat{\phi}_1,\hat{\phi}_2,\hat{\psi})$ as

\begin{equation}\label{AuxCu}
\begin{aligned}
\end{aligned}
\end{equation}

The two formulations are related via the identity $F=C\circ S$, where the parameter-to-state map $S:\kappagamma\to u:=(\hat{\phi}_1,\hat{\phi}_2,\hat{\psi})$ is implicitly defined by the identity

\begin{equation}\label{eq:S}
\end{equation}

Thus, using the forward operator $F$ requires an analysis of the operator $S$.

§ FORWARD PROBLEM AND FUNCTION SPACE SETTING

In the following, function spaces such as $L^2(\Omega;\mathbb{C})$ or $H^1(\Omega;\mathbb{C})$ will be regarded as spaces of functions with values in $\mathbb{C}$, but treated as real Hilbert spaces with a real valued inner product, e.g. $(v,w)_{L^2(\Omega)}=\Re(\int_\Omega v\overline{w}\, dx)$.
The $L^2$ space of real valued functions will simply denoted by $L^2(\Omega)$.

Consider the Laplace operator equipped with impedance boundary conditions, defined in its weak form by

\[
\begin{aligned}
&\langle D_\sigma(\hat{\psi},w\rangle:=B_\sigma(\hat{\psi},w)\\
\Re\Bigl((1-\imath)\Bigl(\int_{\Omega} \nabla\hat{\psi}\cdot\nabla \overline{w})\, dx+\imath\int_{\partial\Omega}\sigma\hat{\psi} \overline{w} \, ds\Bigr)\Bigr), \\
&\forall w\in H^1(\Omega)\,.
\end{aligned}
\]

Here $B_\sigma$ is a symmetric, bounded and coercive bilinear form on on $H^1(\Omega;\mathbb{C})$ by the identity

\[
\int_{\Omega} \vert\nabla\hat{\psi}\vert^2\, dx+\int_{\partial\Omega}\sigma\vert\hat{\psi}\vert^2 \, ds
\]

and Poincaré's inequality. Thus, by the Lax-Milgram Lemma, $D_\sigma:H^1(\Omega;\mathbb{C})\to H^1(\Omega;\mathbb{C})^*$ is boundedly invertible and its inverse is compact as an operator from $L^2_{\kappatil}(\Omega;\mathbb{C})$ into itself, where $L^2_{\kappatil}(\Omega;\mathbb{C})$ is the weighted $L^2$ space with weight function $\kappatil\in L^2(\Omega)$, $\kappatil\geq0$ almost everywhere. Thus, by spectral theory for compact operators, $D_\sigma$ has a countable sequence of positive real eigenvalues tending to infinity, which we will denote by $\{\lambda^\sigma_n\, : n\in \mathbb{N}\}$.
Likewise, the eigenvalues of the operators $D_{\sigma_k}$ defined by the Laplacian on $\Omega$ with impedance boundary conditions (coefficient $\sigma_k$) on $\partial\Omega\setminus\Sigma_k$ are given by the countable set $\{\lambda^{\sigma_k}_n\, : n\in \mathbb{N}\}$, $k\in\{1,2\}$.
Thus (<ref>) is uniquely solvable provided $\hat{g}_k\in H^{-\frac12}(\Sigma_k)$ and $\omega_k\notin \{\lambda^{\sigma_k}_n\, : n\in \mathbb{N}\}$, $k\in\{1,2\}$, $\omega_1-\omega_2\notin\{\lambda^\sigma_n\, : n\in \mathbb{N}\}$.

Higher regularity (actually only higher summability) can be achieved
under the additional assumption $\hat{g}_k\in (W^{1-\frac{1}{q},\frac{q}{q-1}}(\partial\Omega\setminus\Sigma_k))^*\subseteq(\mbox{tr}_{\Sigma_k}(W^{1,\frac{q}{q-1}}(\Omega)))^*$, meaning that the linear map $v\mapsto \int_\Sigma\hat{g}_k v\, ds$ lies in $(W^{1,\frac{q}{q-1}}(\Omega)))^*$.
Therefore according to elliptic regularity (e.g. <cit.>), (<ref>)
admits weak solutions $\hat{\phi}_k\in W^{1,q}(\Omega;\mathbb{C})$, $k\in\{1,2\}$.
Thus, the right hand side of (<ref>) has the following regularity.
From $\kappatil\in L^2(\Omega)$ and $\nabla\hat{\phi}_1$, $\overline{\nabla\hat{\phi}_2} \in L^{q}(\Omega)$ we conclude by Hölder's inequality

\begin{equation}\label{Hoelderabc}
\begin{aligned}
&\|a\,b\,\overline{c}\|_{L^r}\leq \|a\|_{L^2}\|b\|_{L^{\frac{4r}{2-r}}}\|c\|_{L^{\frac{4r}{2-r}}}\\
&\mbox{ for any }a\in L^2(\Omega), \ b,\,c\in L^{\frac{4r}{2-r}}(\Omega,\mathbb{C})
%&\leq \|\kappatil\|_{L^2}\|\nabla\hat{\phi}_1\cdot\overline{\nabla\hat{\phi}_2}\|_{L^{\frac{2r}{2-r}}}\\
%&\leq \|\kappatil\|_{L^2}\|\nabla\hat{\phi}_1\|_{L^{\frac{4r}{2-r}}} \|\nabla\hat{\phi}_2\|_{L^{\frac{4r}{2-r}}}
\end{aligned}
\end{equation}

that $\kappatil\nabla\hat{\phi}_1\cdot\overline{\nabla\hat{\phi}_2} \in L^r(\Omega)\subseteq W^{-1,p}(\Omega)$, provided

\begin{equation}\label{eq:qrp}
r\leq\min\left\{2,\frac{q}{2}\right\}\mbox{ and }
\frac{2r}{2-r}\leq\frac{q}{2} \mbox{ and }
1-\frac{d}{p^*}\geq -\frac{d}{r^*}
\end{equation}

where $p^*=\frac{p}{p-1}$ denotes the dual index.
This regularity (and even more) also holds true for the second quadratic term $\gammatil\hat{\phi}_1\overline{\hat{\phi}_2}$, with $\gammatil\in L^2(\Omega)$. Thus we conclude $\psi\in W^{1,p}(\Omega)$ (cf. <cit.>) and hence, by the Trace Theorem,  $\mbox{tr}_\Gamma\psi\in W^{1-\frac{1}{p},p}(\Gamma)\subseteq L^2(\Gamma)$ provided

\begin{equation}\label{eq:p}
1-\frac{1}{p}-\frac{d-1}{p}\geq -\frac{d-1}{2}\,.
\end{equation}

It is readily checked that conditions (<ref>), (<ref>) can be satisfied for $\Omega\subseteq\mathbb{R}^d$,

by choosing,

\[
\begin{aligned}
%q\geq \frac{4d}{3}\,, \quad p=\frac{2d}{d+1}\,, \quad p=\frac{2d}{d+3}\,,
&\frac{2d}{d+2}\leq p^*\leq\frac{2d}{d-1}, \\
&2\leq r^*\leq\frac{d p^*}{d-p^*} (\mbox{or }p^*\geq d), \quad
q\geq \frac{4r}{2-r},
\end{aligned}
\]

that is, in the physically relevant case $d\leq 3$, e.g.
$p=\frac32$, $r=1$, $q=4$.

Thus we have proven

Let $\hat{g}_k\in H^{-\frac12}(\Sigma_k)$ and $\omega_k\notin \{\lambda^{\sigma_k}_n\, : n\in \mathbb{N}\}$, $k\in\{1,2\}$, $\omega_1-\omega_2\notin\{\lambda^\sigma_n\, : n\in \mathbb{N}\}$ the sets of eigenvalues of the Laplacians $D_{\sigma_k}$, $D_\sigma$ with impedance boundary conditions.

Then the parameter-to-state map $S:\mathcal{D}(F)\to W^{1,q}(\Omega;\mathbb{C})^2\times W^{1,p}(\Omega;\mathbb{C})$, and the forward operator $F:\mathcal{D}(F)\to L^2(\Gamma)$ are well-defined by (<ref>), (<ref>), (<ref>), (<ref>) on
$\mathcal{D}(F)=\{\kappa\in L^2(\widetilde{\Omega})\, : \, \kappatil\geq0\mbox{ a.e. }\}\times L^2(\widetilde{\Omega})$.

The domain $\mathcal{D}(F)$ has empty interior with respect to the $L^2$ topology and this prevents applicability of convergence results for the Newton and gradient methods to be discussed below.
To avoid this, we restrict $F$ to an open ball around a strictly positive $L^\infty$ function $\kappa_0>0$ (e.g., the background)

\begin{equation}\label{DF}
\begin{aligned}
\tilde{\mathcal{D}}(F)=&\{\kappa\in L^2(\widetilde{\Omega})\, : \, \|\kappatil-\kappa_0\|_{L^2}\leq\rho\}\\
&\times L^2(\widetilde{\Omega})
\end{aligned}
\end{equation}

for $\rho$ sufficiently small, and apply a fixed point argument to obtain well-definedness of $F$ on $\tilde{\mathcal{D}}(F)$, see, e.g., [12].

For use in Newton and gradient type methods we also need differentiability of $F$. It sufficies to prove that the parameter-to-state map $S$ is differentiable, since $F=C\circ S$ with $C$ being a bounded linear operator.
It is straightforward to see that for
the difference $(d\hat{\phi}_1,d\hat{\phi}_2,d\hat{\psi}):=S(\kappagamma+\delta\kappagamma)-S(\kappagamma)$ satisfies the weak form of

\begin{equation}\label{phihatkpsihat_diff}
\begin{aligned}
&-\omega_k^2\kappatil \, d\hat{\phi}_k-\Delta \, d\hat{\phi}_k=\omega_k^2\,\underline{\widetilde{\delta\kappa}}\, \hat{\phi}_k^+ \mbox{ in }\Omega \\
&\partial_\nu \, d\hat{\phi}_k=0 \mbox{ on } \Sigma_k\,,\quad k\in\{1,2\}
\\
&-(\omega_1-\omega_2)^2\kappatil \, d\hat{\psi}-\Delta \, d\hat{\psi}=f_d
\mbox{ in }\Omega\\
& \quad f_d= (\omega_1-\omega_2)^2\kappatil \, \hat{\psi}^+ +2\imath(\omega_1-\omega_2)\\
+\kappatil(\nabla\, d\hat{\phi}_1\cdot\overline{\nabla\hat{\phi}_2^+}
+\nabla\hat{\phi}_1\cdot\overline{\nabla\, d\hat{\phi}_2})\\
&\qquad + \omega_1\omega_2 \left(
\underline{\widetilde{\delta\gamma}}\,\hat{\phi}_1^+\overline{\hat{\phi}_2^+}
+\gammatil( d\hat{\phi}_1\overline{\hat{\phi}_2^+}
\right)\Bigr)
\end{aligned}
\end{equation}

and therefore, formally

\begin{equation}\label{phihatkpsihat_deriv}
\begin{aligned}
&-\omega_k^2\kappatil \, \delta\hat{\phi}_k-\Delta \, \delta\hat{\phi}_k=\omega_k^2\,\underline{\widetilde{\delta\kappa}}\, \hat{\phi}_k \mbox{ in }\Omega \\
&\partial_\nu \, \delta\hat{\phi}_k=0 \mbox{ on } \Sigma_k\,,\quad k\in\{1,2\}
\\
&-(\omega_1-\omega_2)^2\kappatil \, \delta\hat{\psi}-\Delta \, \delta\hat{\psi}=
f_\delta\mbox{ in }\Omega\\
& \quad f_\delta=(\omega_1-\omega_2)^2\kappatil \, \hat{\psi} +2\imath(\omega_1-\omega_2)\\
+\kappatil(\nabla\, \delta\hat{\phi}_1\cdot\overline{\nabla\hat{\phi}_2}
+\nabla\hat{\phi}_1\cdot\overline{\nabla\, \delta\hat{\phi}_2})\\
&\qquad + \omega_1\omega_2 \left(
\underline{\widetilde{\delta\gamma}}\,\hat{\phi}_1\overline{\hat{\phi}_2}
+\gammatil( \delta\hat{\phi}_1\overline{\hat{\phi}_2}
\right)\Bigr)\,.
\end{aligned}
\end{equation}

Hence the first order Taylor remainder

\begin{equation}\label{phihatkpsihat_remainder}
\begin{aligned}
&-\omega_k^2\kappatil \, \phiTay_k -\Delta \, \phiTay_k=\omega_k^2\,\underline{\widetilde{\delta\kappa}}\, d\hat{\phi}_k \mbox{ in }\Omega \\
&\partial_\nu \, \phiTay_k=0 \mbox{ on } \Sigma_k\,,\quad k\in\{1,2\}\\
&-(\omega_1-\omega_2)^2\kappatil \, \psiTay-\Delta \, \psiTay=f_{rest}\mbox{ in }\Omega\\
& \quad f_ {rest}=
(\omega_1-\omega_2)^2\kappatil \, d\hat{\psi} +2\imath(\omega_1-\omega_2)\\
(\nabla\, d\hat{\phi}_1\cdot\overline{\nabla\hat{\phi}_2^+}
+\nabla\hat{\phi}_1\cdot\, \overline{d\nabla\hat{\phi}_2})\\
&\ +\kappatil(\nabla\, \phiTay_1\cdot\overline{\nabla\hat{\phi}_2}
+\nabla\, d\hat{\phi}_1\cdot\overline{\nabla\, d\hat{\phi}_2}
+\nabla\hat{\phi}_1\cdot\overline{\nabla \phiTay_2})\\
&\ + \omega_1\omega_2 \Bigl(
\underline{\widetilde{\delta\gamma}}\,
+\hat{\phi}_1\, \overline{d\hat{\phi}_2})\\
\Bigr)\Bigr).
\end{aligned}
\end{equation}

Here we have used the identities

\[
\begin{aligned}
&(a+\, \delta a) (b+db) (c+\,dc) - abc\\
&\quad= \delta a\, (b+db)(c+\,dc)\, + \, a\, db \, (c+dc)\, +\, a\,b\,dc\\
&(a+\, \delta a) (b+db) (c+\,dc) - abc\\
&\quad-(\delta a\, b\, c\, + \, a\, \delta b\, c\, +\, a\,b\,\delta c)\\
&\quad=\delta a\,(db\, (c+dc)\, + \, b\, dc) \\
&\qquad+ a[(db-\delta b)c+db\,dc+b(dc-\delta c)].
\end{aligned}
\]

Regularity arguments as in the proof of Theorem <ref> lead to estimates of the form

\[
\begin{aligned}
&\|d\hat{\phi}_k\|_{W^{1,q}}\leq C \|\underline{\widetilde{\delta\kappa}}\, \hat{\phi}_k^+\|_{(W^{1,q^*})^*}\\
&\|d\hat{\psi}\|_{W^{1,p}}\leq C \|f_d\|_{(W^{1,p^*})^*}\\
&\|\delta\hat{\phi}_k\|_{W^{1,q}}\leq C \|\underline{\widetilde{\delta\kappa}}\, \hat{\phi}_k\|_{(W^{1,q^*})^*}\\
&\|\delta\hat{\psi}\|_{W^{1,p}}\leq C \|f_\delta\|_{(W^{1,p^*})^*}\\
&\|\phiTay_k\|_{W^{1,q}}\leq C \|\underline{\widetilde{\delta\kappa}}\, d\hat{\phi}_k\|_{(W^{1,q^*})^*}\\
&\|\psiTay\|_{W^{1,p}}\leq C \|f_{rest}\|_{(W^{1,p^*})^*}\,,
\end{aligned}
\]

where $f_d$, $f_\delta$, $f_{rest}$ can be estimated by the same Hölder inequalities and Sobolev embeddings as those used for the proof of Theorem <ref>.

This proves Fréchet differentiability.

Under assumptions of Theorem <ref>, the parameter-to-state map $S$ and the forward operator $F$ are Fréchet differentiable on $\tilde{\mathcal{D}}(F)$ as defined in (<ref>) with respect to the $L^2$ topology in preimage space, as mappings to $W^{1,q}(\Omega;\mathbb{C})^2\times W^{1,p}(\Omega;\mathbb{C})$ and $L^2(\Gamma)$, respectively.

Concerning further convergence conditions for Newton and gradient type methods, cf. e.g. [13], we briefly comment on the tangential cone condition

\begin{equation}\label{tcc}
\begin{aligned}
&\quad\leq c_{tc}\|F(\kappagamma+\delta\kappagamma)-F(\kappagamma)\|.
\end{aligned}
\end{equation}

In case of known speed of sound $c$, when we seek to identify $\gamma=\gamma(x)$ only,
the inverse problem becomes an inverse source problem, see (<ref>) below, and is therefore affinely linear, thus trivially satisfying (<ref>) with $c_{tc}=0$.
Conversely, if $c=c(x)$ is to be determined, the inverse probems is closely related to the well-known and well-investigated model problem of recovering the potential $c$ in the Schrödinger equation $-\Delta u + c u=0$. This is known to satisfy the tangential cone condition only in case of complete observations of $u$ on all of $\Omega$ [12]. Thus (<ref>) cannot be expected to be verifiable in our boundary observation setting.

In the definition of gradient type methods (and also in the implementation of Newton type methods) we will need the adjoint of $F'(\kappagamma)$, which we therefore derive here.
First of all, note that by $F=C\circ S$ with $S$ defined by (<ref>) and the Implicit Function Theorem we can write

\begin{equation}\label{eq:SpFpKL}
S'(\kappagamma)=-K^{-1} L\,, \quad
F'(\kappagamma)=-CK^{-1} L\,,
\end{equation}

\begin{equation}\label{eq:KL}
K=\frac{\partial A}{\partial u}(u^{(n)},\kappagamma^{(n)})\,, \quad
L=\frac{\partial A}{\partial \kappagamma}(u^{(n)},\kappagamma^{(n)})
\end{equation}

are the linearizations of the operator $A$ from (<ref>) with respect to the states and the parameters, respectively. They are given by

\begin{equation}\label{eq:Kmod}
\begin{aligned}
&\langle K\, (\delta\hat{\phi}_1,\delta\hat{\phi}_2,\delta\hat{\psi}),(v_1,v_2,w)\rangle=\\
&\Re\Bigl((1-\imath)\Bigl(\sum_{k=1}^2\int_{\Omega} (-\omega_k^2\kappatil\, \delta\hat{\phi}_k\overline{v}_k+\nabla\ \delta\hat{\phi}_k\cdot\nabla \overline{v}_k)\, dx\\
&+\int_{\Omega} (-(\omega_1-\omega_2)^2\kappatil\, \delta\hat{\psi}\overline{w}+\nabla\, \delta\hat{\psi}\cdot\nabla \overline{w})\, dx\\
&+\imath\int_{\partial\Omega\setminus\Sigma}\sigma_k\, \delta\hat{\phi}_k \overline{v}_k \, ds
+\imath\int_{\partial\Omega}\sigma\, \delta\hat{\psi} \overline{w} \, ds\\
&-\int_\Omega 2(\omega_1-\omega_2) \imath \Bigl(
\kappatil
\Bigl(\nabla\, \delta\hat{\phi}_1\cdot\overline{\nabla\hat{\phi}_2}
+\nabla\hat{\phi}_1\cdot\overline{\nabla\, \delta\hat{\phi}_2}\Bigr)\\
&\quad+\omega_1\omega_2 \gammatil \Bigl(
\, \delta\hat{\phi}_1\,\overline{\hat{\phi}_2}
+\hat{\phi}_1\overline{\, \delta\hat{\phi}_2}
\Bigr)
\Bigr)\Bigr)\,\overline{w}\, dx \Bigr)
\end{aligned}
\end{equation}

\begin{equation}\label{eq:Lmod}
\begin{aligned}
&\langle L\, (\delta\kappa,\delta\gamma),(v_1,v_2,w)\rangle=\\
&\Re\Bigl((1-\imath)\Bigl(\sum_{k=1}^2\int_{\Omega} -\omega_k^2 \underline{\widetilde{\delta\kappa}} \,\hat{\phi}_k\overline{v}_k\, dx\\
& +\int_{\Omega} -(\omega_1-\omega_2)^2 \underline{\widetilde{\delta\kappa}} \,\hat{\psi}\overline{w}\, dx\\
&-2(\omega_1-\omega_2) \imath \\
&\cdot\int_\Omega \left(
\underline{\widetilde{\delta\kappa}} \,\nabla\hat{\phi}_1\cdot\overline{\nabla\hat{\phi}_2}
+\omega_1\omega_2 \, \underline{\widetilde{\delta\gamma}}\, \hat{\phi}_1\overline{\hat{\phi}_2}\right)\,\overline{w}\, dx \Bigr)\Bigr)
\end{aligned}
\end{equation}

for any $v_1,v_2,w\in H^1(\Omega;\mathbb{C})$.
The identity (<ref>) with (<ref>) can also be used to determine the adjoint operator $F'(\kappagamma)^*=-(CK^{-1}L)^*$ as a Hilbert space adjoint in $L^2$.
To this end, for given $r\in L^2(\Gamma)$ we want to find $\xizeta:=(CK^{-1}L)^*r$ such that

\[
\begin{aligned}
&\langle CK^{-1}L\,\delta\kappagamma,r\rangle_{L^2(\Gamma)} = \langle \delta\kappa,\xi\rangle_{L^2(\widetilde{\Omega})}
+\langle \delta\gamma,\zeta\rangle_{L^2(\Omega)} \\
&\mbox{ for all }\delta\kappagamma\in L^2(\widetilde{\Omega})\times L^2(\Omega)\,.
\end{aligned}
\]

We introduce the auxiliary variables $(\widetilde{\phi}_1,\widetilde{\phi}_2,\widetilde{\psi}):=K^{-1}L\,\delta\kappagamma$, which allows us to use the identity

\begin{equation}\label{eq:stateeq}
\begin{aligned}
&\langle K(\widetilde{\phi}_1,\widetilde{\phi}_2,\widetilde{\psi}),(v_1,v_2,w)\rangle
=\langle L\delta\kappagamma,(v_1,v_2,w)\rangle\\
&\mbox{ for all }(v_1,v_2,w)
\end{aligned}
\end{equation}

and define $(p_1,p_2,q)$ as the solution to the adjoint equation

\begin{equation}\label{eq:adjeq}
\begin{aligned}
&\langle K(\delta\hat{\phi}_1,\delta\hat{\phi}_2,\delta\hat{\psi}),(p_1,p_2,q)\rangle\\
&\qquad=\langle C(\delta\hat{\phi}_1,\delta\hat{\phi}_2,\delta\hat{\psi}),r\rangle_{L^2(\Gamma)}
\\&
\mbox{ for all } (\delta\hat{\phi}_1,\delta\hat{\phi}_2,\delta\hat{\psi})\,.
\end{aligned}
\end{equation}

Using (<ref>), and (<ref>) together with (<ref>), (<ref>), we get

\[
\begin{aligned}
&\langle C(\widetilde{\phi}_1,\widetilde{\phi}_2,\widetilde{\psi}),r\rangle_{L^2(\Gamma)}
=\langle L\,\delta\kappagamma,(p_1,p_2,q)\rangle\\
&=-\Re\Bigl(\sum_{k=1}^2\int_{\widetilde{\Omega}} \omega_k^2 \,\delta\kappa \,\hat{\phi}_k^{(n)}\overline{p}_k\, dx\\
&\qquad+\int_{\widetilde{\Omega}} (\omega_1-\omega_2)^2 \,\delta\kappa \,\hat{\psi}^{(n)}\overline{q}\, dx\\
&\qquad+2(\omega_1-\omega_2) \imath \\
&\hspace*{-0.5cm}\cdot\int_{\widetilde{\Omega}} \left(
\delta\kappa \,\nabla\hat{\phi}_1^{(n)}\cdot\overline{\nabla\hat{\phi}_2^{(n)}}
+\omega_1\omega_2 \, \delta\gamma\, \hat{\phi}_1^{(n)}\overline{\hat{\phi}_2^{(n)}}\right)\,\overline{q}\, dx \Bigr)\\
&=\langle \delta\kappa,\xi\rangle_{L^2(\widetilde{\Omega})}
+\langle \delta\gamma,\zeta\rangle_{L^2(\widetilde{\Omega})}
\end{aligned}
\]

\[
\begin{aligned}
&\xi=-\Re\Bigl(\sum_{k=1}^2 \omega_k^2 \hat{\phi}_k^{(n)}\overline{p}_k\\
&+(\omega_1-\omega_2)^2 \hat{\psi}^{(n)}\overline{q}
+2(\omega_1-\omega_2) \imath \nabla\hat{\phi}_1^{(n)}\cdot\overline{\nabla\hat{\phi}_2^{(n)}} \overline{q} \Bigr)\vert_{\widetilde{\Omega}}
\\
&\zeta=-\Re\Bigl(2(\omega_1-\omega_2) \imath \omega_1\omega_2 \hat{\phi}_1^{(n)}\overline{\hat{\phi}_2^{(n)}}\overline{q}\Bigr)\vert_{\widetilde{\Omega}}.
\end{aligned}
\]

Thus we end up with an explicit expression for $\xizeta:=F'(\kappagamma)^*r$.
For this purpose, the adjoint states have to be computed as solutions to (<ref>), that is, the weak form of

\[
\begin{aligned}
&-(\omega_1-\omega_2)^2\kappatil q-\Delta q=0 \mbox{ in }\Omega\setminus\Gamma\,, \\
&\partial_\nu q=-\imath\sigma q \mbox{ on } \partial\Omega\,, \quad
\left[\partial_\nu q\right]=-\imath(\omega_1-\omega_2) r \mbox{ on }\Gamma
\end{aligned}
\]

where $\left[\partial_\nu q\right]$ denotes the jump of the normal derivative over the interface $\Gamma$, as well as

\[
\begin{aligned}
&-\omega_k^2\kappatil p_k-\Delta p_k
=f_k \mbox{ in }\Omega\,, \\
&\partial_\nu p_k=-\imath\sigma_k p_k \mbox{ on } \partial\Omega\setminus\Sigma_k\,, \quad
\partial_\nu p_k=0 \mbox{ on } \Sigma_k\\
\nabla \bigl( \kappatil q\nabla \phi_{\not k}\bigr)
+\omega_1\omega_2 \gammatil q\phi_{\not{k}}\Bigr)
\end{aligned}
\]

with $\not{\!k}=\begin{cases}2\mbox{ for }k=1\\1\mbox{ for }k=2\end{cases}$.

§ UNIQUENESS

In case the speed of sound $c$ is known, reconstruction of $\gamma=\gamma(x)$ in the time domain (<ref>), (<ref>) or the frequency domain formulation (<ref>), (<ref>) amounts to an inverse source problem.
Indeed, setting

\[
\begin{aligned}
h(x,\omega_d;\omega_2)&=2\imath\omega_d\omega_2 \nabla\hat{\phi}(x,\omega_2+\omega_d)\cdot\overline{\nabla\hat{\phi}(x,\omega_2)}\\
\multip(x,\omega_d;\omega_2)&=2\imath
\hat{\phi}(x,\omega_2+\omega_d)\overline{\hat{\phi}(x,\omega_2)}\\
\gamma'(x)&=\frac{\gamma(x)-1}{2c(x)^2}
\end{aligned}
\]

and multiplying with $c^2$, we can write (<ref>) as

\begin{equation}\label{eq:psihat_hf}
-\omega_d^2 \hat{\psi}+\mathcal{A}_c \hat{\psi}= h(\omega_d;\omega_2)+ \multip(\omega_d;\omega_2) \gamma'
\end{equation}

Here we denote the difference frequency $\omega_1-\omega_2$ by $\omega_d$ and the solution of (<ref>) with boundary excitation $\hat{g}_k=\hat{g}(\omega)$ by $\hat{\phi}_k(\omega)$.
Note that the functions $\multip$ and $h$ are known from the known excitations $\hat{g}$.
Moroever, we denote by $\mathcal{A}_c$
either of
the elliptic differential operator
$-c^2\Delta$ with homogeneous impedance boundary conditions;
or $-\frac{c^2}{d}\nabla\cdot(D\nabla)$ with homogeneous Neumann boundary conditions. In both cases
$\mathcal{A}_c$ is a selfadjoint nonnegative definite operator with respect to the weighted $L^2$ inner product with weight function $w=\frac{1}{c^2}$

or $w=\frac{d}{c^2}$, respectively.

By $\{(\lambda_k,(\varphi_j^k)_{j\in \indset_k}\, : \, k\in\mathbb{N}\}$ we denote the corresponding eigensystem, where in case of multiple eigenvalues we collect the eigenfunctions corresponding to $\lambda_k$ in the set $\{\varphi_j^k\, : \, j\in \indset_k\}$ with some finite index set $\indset_k$. (Note that in one space dimension, the eigenfunctions are single and so $\indset_k=\{1\}$.)
The requirements on $c$ for this purpose are

\begin{equation}\label{eq:cond-c}
c,\tfrac{1}{c}\in L^\infty(\Omega)\,.
\end{equation}

Since the eigenfunctions form an orthonormal basis of $L^2_{w}$, we can expand $\hat{\psi}$ with respect to this basis

\[
\hat{\psi}(x,\omega_d;\omega_2)
=\sum_{k=1}^\infty\sum_{j\in \indset_k}\langle \hat{\psi}(\omega_d;\omega_2),\varphi_j^k
\rangle_{L^2_{w}} \varphi_j^k(x)\,.
\]

This allows us to express the observations according to (<ref>) by

\begin{equation}\label{eq:obs_eig}
\begin{aligned}
&y(x_0,\omega_d;\omega_2)=\hat{\psi}(x_0,\omega_d;\omega_2) \\
&=\sum_{k=1}^\infty\sum_{j\in \indset_k}\langle \hat{\psi}(\omega_d;\omega_2),\varphi_j^k
\rangle_{L^2_{w}} \varphi_j^k(x_0) \\
&x_0\in\Gamma, \quad \omega_d\in \freqset\,,
\end{aligned}
\end{equation}

where we assume that we can take observations for all difference frequencies in some set $\freqset $, while $\omega_2$ is fixed.
On the other hand, taking inner products of (<ref>) with $\varphi_j^k$ and using the eigenvalue equation $\mathcal{A}_c\varphi_j^k=\lambda_j\varphi_j^k$ we obtain the identity

\begin{equation*}%\label{eq:psihat_j}
\begin{aligned}
&\langle \hat{\psi}(\omega_d;\omega_2),\varphi_j^k\rangle_{L^2_{w}}\\
&= \frac{1}{-\omega_d^2+\lambda_j}
\langle h(\omega_d;\omega_2)+ \multip(\omega_d;\omega_2) \gamma',\varphi_j^k\rangle_{L^2_{w}}.
\end{aligned}
\end{equation*}

Combining this with (<ref>)  we get

\begin{equation}\label{eq:obs_eig_sumj0}
\begin{aligned}
=\sum_{k=1}^\infty\frac{1}{-\omega_d^2+\lambda_k}\sum_{j\in \indset_k}
\langle \multip(\omega_d;\omega_2) \gamma',\varphi_j^k\rangle_{L^2_{w}}\\
&x_0\in\Gamma, \quad \omega_d\in \freqset \,,
\end{aligned}
\end{equation}

where $\tilde{y}$ is the modified observation function

\[
\begin{aligned}
&=y(x_0,\omega_d;\omega_2)- (-\omega_d^2+\mathcal{A}_c)^{-1} h(\omega_d;\omega_2)\,,
\end{aligned}
\]

thus a known quantity.
In order to obtain from this the desired information on $\gamma'$, we assume that $\hat{g}(\omega)$ has been chosen such that $\multip$ factorizes into a frequency dependent and a space dependent part

\begin{equation}\label{eq:fab}
\multip(x,\omega_d;\omega_2)=a(\omega_d;\omega_2) b(x)
\end{equation}

so that (<ref>) becomes

\begin{equation}\label{eq:obs_eig_sumj}
\begin{aligned}
=\sum_{k=1}^\infty\frac{a(\omega_d;\omega_2)}{-\omega_d^2+\lambda_k}\sum_{j\in \indset_k}
\langle b \gamma',\varphi_j^k\rangle_{L^2_{w}}
\\& x_0\in\Gamma, \quad \omega_d\in \freqset \,.
\end{aligned}
\end{equation}

Both sides of this equality have sigularities at $\omega_d=\pm\sqrt{\lambda_\ell}$.
Thus, these poles provide the location of the eigenvalues of $\mathcal{A}_c$ and therewith some information on $c$ (see Remark <ref> below).
Moreover, multiplying with $(\omega_d-\sqrt{\lambda_\ell})$ and taking the limit $\omega_d\to\sqrt{\lambda_\ell}$, we can extract the contribution due to the $\ell$th eigenfunction

\begin{equation}\label{eq:obs_eig_ell}
\begin{aligned}
&\lim_{\omega_d\to\sqrt{\lambda_\ell}} (\omega_d-\sqrt{\lambda_\ell})
\tilde{y}(x_0,\omega_d;\omega_2)\\
&= -\frac{a(\sqrt{\lambda_\ell};\omega_2)}{2\sqrt{\lambda_\ell}}\sum_{j\in \indset_k}
\langle b \gamma',\varphi_j^\ell\rangle_{L^2_{w}}
\quad x_0\in\Gamma
\end{aligned}
\end{equation}

For this to work out, we need to assume that

\begin{equation}\label{eq:cond_I}
\sqrt{\lambda_\ell}\mbox{ is an interior point of }\freqset \mbox{ for all }\ell\in\mathbb{N}.
\end{equation}

Finally (<ref>) allows to uniquely determine the coefficients $\langle b \gamma',\varphi_j^\ell\rangle_{L^2_{w}}$ in

\begin{equation}\label{gamma_recon}
\gamma'(x)=\frac{1}{b(x)}
\sum_{\ell=1}^\infty\sum_{j\in \indset_\ell} \langle b \gamma',\varphi_j^\ell\rangle_{L^2_{w}}
\end{equation}

provided $b$ vanishes nowhere and

\begin{equation}\label{eq:cond_linindep}
\{\varphi_j^\ell\vert_\Gamma\,:\, j\in \indset_\ell\} \mbox{ is linearly independent}
\end{equation}

Thus we have proven the following uniqueness result on recovery of $\gamma(x)$.

Assume that $c$ is known and satisfies (<ref>),
that $\freqset $ and $\Gamma$ are chosen such that (<ref>), (<ref>)
hold, and that $\hat{g}(\omega_2+\omega_d)$, is chosen such that (<ref>)
holds for all $\omega_d\in \freqset $ with $b\in L^\infty$.

Then $\gamma\in L^2(\Omega)$ is uniquely determined on the set $\{x\in\widetilde{\Omega}\, : \, b(x)\not=0\}$ by the observations
$y(x_0,\omega_d;\omega_2)=\hat{\psi}(x_0,\omega_d;\omega_2)$, $x_0\in\Gamma$, $\omega_d\in \freqset $.

Obviously, if (<ref>), (<ref>) only hold with $\mathbb{N}$ replaced by $\{1,\ldots,N\}$, we can recover the first $N$ coefficients of $b\gamma'$.

Note that no regularity assumptions with respect to $\omega_d$ need to be imposed here.

Condition (<ref>) has been discussed in detail in <cit.>. It is trivially satisfied with $\Gamma$ containing a single point $\{x_0\}$ in one space dimension, since the eigenvalues of $\mathcal{A}_c$ are single then.
Moroever, it can be extended to higher space dimensions and geometric settings in which the eigenfunctions allow for separation of variables. A simple 2-d example is a disc with radius $r$, where using polar coordinates, the eigenfunctions can be written in terms of Bessel functions. A circle with almost any radius $r_*\in(0,r]$ can then be used as observation manifold $\Gamma$, as shown in  <cit.>.

To achieve the separability (<ref>) of $\multip$ we supplement the boundary excitation $\hat{g}_k(\omega)$ by an interior one $f_g(\omega)$, which we view as an approximation of a source $\tilde{g}(\omega)\,\delta_\Sigma$ concentrated on $\Sigma$, cf., e.g., [15].

The resulting equation for $\hat{\phi}_k(\omega)$

\[
\begin{aligned}
&-\omega^2\hat{\phi}_k -\frac{c^2}{d}\Delta \hat{\phi}_k = f_g(\omega)\mbox{ in }\Omega\\
&\partial_\nu \hat{\phi}_k=-\sigma_k\hat{\phi}_k\mbox{ on }\partial\Omega\setminus\Sigma_k\,,\
\partial_\nu \hat{\phi}_k=\hat{g}_k(\Omega)\mbox{ on }\Sigma_k
\end{aligned}
\]

then has a solution of the form $\hat{\phi}_k(x,\omega)=\tilde{a}(\omega)\tilde{b}(x)$
if, e.g., we choose $\tilde{b}$ such that $\Delta \tilde{b}=0$ in $\Omega$, $\partial_\nu \tilde{b}=-\sigma_k\tilde{b}$ on $\partial\Omega\setminus\Sigma$, and set
$\hat{g}_k(\omega):=\tilde{a}(\omega)\partial_\nu \tilde{b}\vert_\Sigma$,

In case of constant sound speed $c$, uniqueness for the above inverse source problem for $\gamma(x)$ in the time domain formulation (<ref>), (<ref>) from boundary observations under a space-time separability assumption (similar to the space-frequency one (<ref>)) can be conclused from <cit.>, provided $\Gamma$, $c$ and $T$ satisfy <cit.>, which is basically a condition on sufficient size of $\Gamma$ and $T$, depending on the speed $c$ of sound propagation.

Other related uniqueness results for $\gamma$ have been found recently in the context of nonlinearity imaging in [17, 18].

Note that from the poles on both sides of (<ref>) we also obtain the eigenvalues of $\mathcal{A}_c$. According to Sturm-Liouville theory, applied as in <cit.>, this uniquely determines $c(x)$ in one space dimension, provided we can take measurements at two different impedance values $\sigma$, $\tilde{\sigma}$. Note however, that we need the eigenfunctions of $\mathcal{A}_c$ for reconstructing $\gamma'(x)$ according to (<ref>), so this only gives a uniqueness result for $c$ alone and no simultaneous uniqueness of $c$ and $\gamma$. Also, its restriction to the 1-d setting limits applicability to our experimental setting.

For uniqueness of $c=c(x)$ in higher space dimensions from boundary measurements,

results on uniqueness of the space-dependent index of refraction $n(x)=\frac{c_0^2}{c(x)^2}$ in inverse scattering, e.g.,  <cit.> or of the potential in the Schrödinger equation <cit.> are relevant.
Note however, that $c$ appears not only in the equation for the observed quantity $\psi$ but also governs the two excitation wave fields $\phi_1$, $\phi_2$ that enter the $\psi$ equation through a source term. This makes the uniqueness question for $c$ more involved than in the mentioned references.

A proof of unique recovery of both $c$ and $\gamma$ is widely open and subject of future research. We will nevertheless in the remainder of this paper discuss some simultaneous numerical reconstruction techniques.

§ ITERATIVE RECONSTRUCTION METHODS

We return to the general case in which both $c$ and $\gamma$ are unknown

and consider the Helmholtz model with absorbing boundary conditions. The PML setting could be treated analogously.

Iteratively regularized Gauss-Newton method IRGNM
A regularized Gauss-Newton step for solving (<ref>) defines $\kappagamma^{(n+1)}$ as a minimizer of

\[
\begin{aligned}
\\&+\alpha^{(n)}\|(\kappagamma-\kappagamma^{(0)})\|_{L^2}^2
\end{aligned}
\]

where $T^{(n)}=F'(\kappagamma^{(n)})$
thus, with $^*$ denoting the Hilbert space adjoint in $L^2$, the Newton step reads as

\[
\begin{aligned}
&\kappagamma^{(n+1)} = \kappagamma^{(n)}+({T^{(n)}}^* T^{(n)}+\alpha^{(n)} I)^{-1}
\\&\hspace*{0.5cm}\Bigl({T^{(n)}}^* (y-F(\kappagamma^{(n)})+\alpha^{(n)}(\kappagamma^{(0)}-\kappagamma^{(n)})\Bigr)
\end{aligned}
\]

or as $\kappagamma^{(n+1)}=\kappagamma^{(n)}+\delta \kappagamma$ where $\delta \kappagamma$ solves the variational equation

\begin{equation}\label{eq:redIRGNM0}
\begin{aligned}
&\langle T^{(n)} \delta\kappagamma, T^{(n)} \xizeta\rangle + \alpha^{(n)} \langle \delta\kappagamma,\xizeta \rangle \\
&= \langle y-F(\kappagamma^{(n)}, T^{(n)} \xizeta\rangle+ \alpha^{(n)} \langle \kappagamma^{(0)}\!\!\!-\kappagamma^{(n)}\!\!,\xizeta \rangle \\
&\mbox{for all } \xizeta\in L^2(\widetilde{\Omega})\times L^2(\widetilde{\Omega})\,,
\end{aligned}
\end{equation}

where $\langle\cdot,\cdot\rangle$ denotes the $L^2$ inner products

\[
\begin{aligned}
&\langle \kappagamma,\xizeta\rangle=\langle \kappa,\xi\rangle_{L^2(\widetilde{\Omega})}+\langle \xi,\zeta\rangle_{L^2(\widetilde{\Omega})},\\
&\langle y,z\rangle=\langle y,z\rangle_{L^2(\Gamma)}\,.
\end{aligned}
\]

Levenberg-Marquardt method
A slightly different version of Newton's method is the Levenberg-Marquardt method defining the new iterate as a minimizer of

\[
\begin{aligned}
\\&+\alpha^{(n)}\|(\kappagamma-\kappagamma^{(n)})\|_{L^2}^2
\end{aligned}
\]

and thus reads as

\[
\begin{aligned}
&\kappagamma^{(n+1)} = \kappagamma^{(n)}+({T^{(n)}}^* T^{(n)}+\alpha^{(n)} I)^{-1} \\&\hspace*{4cm} {T^{(n)}}^* (y-F(\kappagamma^{(n)})
\end{aligned}
\]

i.e., as $\kappagamma^{(n+1)}=\kappagamma^{(n)}+\delta \kappagamma$ with

\begin{equation}\label{eq:redLM0}
\begin{aligned}
&\langle T^{(n)} \delta\kappagamma, T^{(n)} \xizeta\rangle + \alpha^{(n)} \langle \delta\kappagamma,\xizeta \rangle \\
&= \langle y-F(\kappagamma^{(n)}, T^{(n)} \xizeta\rangle\\
&\mbox{for all } \xizeta\in L^2(\widetilde{\Omega})\times L^2(\Omega)\,.
\end{aligned}
\end{equation}

All-at-once Newton type methods
Alternatively, we may apply (regularized versions of) Newton's method directly to (<ref>), i.e., solve the linear system

\begin{equation}\label{eq:AClin}
\begin{aligned}
&A(u^{(n)},\kappagamma^{(n)})+K\,\delta u + L\,\delta\kappagamma=0\\
&Cu^{(n)}+C\delta u=y
\end{aligned}
\end{equation}

for $(\delta u,\delta \kappagamma)=(\delta \hat{\phi}_1,\delta \hat{\phi}_2, \delta \hat{\psi},\delta \kappa, \delta\gamma)$ and set

\begin{equation}\label{eq:step}
u^{(n+1)}=u^{(n)}+\delta u,\quad  \kappagamma^{(n+1)}=\kappagamma^{(n)}+\delta \kappagamma\,.
\end{equation}

\begin{equation}\label{eq:unp12}
\end{equation}

for $u^{(n+\frac12)}=u^{(n)}-K^{-1}A(u^{(n)},\kappagamma^{(n)})$
and eliminating $\delta u=-K^{-1}(A(u^{(n)},\kappagamma^{(n)}+L\, \delta\kappagamma)=u^{(n+\frac12)}-u^{(n)}-K^{-1}L\, \delta\kappagamma$ in (<ref>) yields an equation for $\delta\kappagamma$ only

\begin{equation}\label{eq:elim}
CK^{-1}L \, \delta \kappagamma = Cu^{(n+\frac12)}-y
\end{equation}

from which, after solving for $\delta\kappagamma$, the new state can be computed as

\begin{equation}\label{eq:unp1}
u^{(n+1)}=u^{(n+\frac12)}-K^{-1}L\, \delta\kappagamma
\end{equation}

A regularized version of (<ref>) is (after computing $u^{(n+\frac12)}$ from (<ref>)) to determine $\delta\kappagamma$ from the variational equation

\begin{equation}\label{eq:aaoIRGNM}
\begin{aligned}
&\langle CK^{-1}L\, \delta\kappagamma, CK^{-1}L\, \xizeta\rangle + \alpha^{(n)} \langle \delta\kappagamma,\xizeta \rangle \\
&= \langle Cu^{(n+\frac12)}-y, CK^{-1}L\, \xizeta\rangle\\
\left[+ \alpha^{(n)} \langle \kappagamma^{(0)}-\kappagamma^{(n)},\xizeta \rangle \right]\\
&\mbox{for all } \xizeta\in L^2(\widetilde{\Omega})\times L^2(\Omega)\,,
\end{aligned}
\end{equation}

which is followed by computation of $u^{(n+1)}$ from (<ref>) and of $\kappagamma^{(n+1)}=\kappagamma^{(n)}+\delta\kappagamma$.

In (<ref>), the term in large brackets may be skipped to obtain a Levenberg-Marquardt type version of the iteration.

Also the reduced versions (<ref>), (<ref>) can be rewritten in terms of the operators $A$, $C$, $L$, $K$, applying the Implicit Function Theorem to (<ref>) which yields $S'(\kappagamma^{(n)})=-K^{-1} L$, hence $T^{(n)}=-CK^{-1} L$ with $K,L$ according to (<ref>).
This leads to almost the same formulation as in (<ref>)

\begin{equation}\label{eq:redIRGNM}
\begin{aligned}
&\langle CK^{-1}L\, \delta\kappagamma, CK^{-1}L\, \xizeta\rangle + \alpha^{(n)} \langle \delta\kappagamma,\xizeta \rangle \\
&= \langle Cu^{(n)}-y, CK^{-1}L\, \xizeta\rangle\\
&\qquad\left[+ \alpha^{(n)} \langle \kappagamma^{(0)}-\kappagamma^{(n)},\xizeta \rangle \right]\\
&\mbox{for all } \xizeta\in L^2(\widetilde{\Omega})\times L^2(\Omega)\,,
\end{aligned}
\end{equation}

(note that here $u^{(n)}$ is defined as the solution to $A(u^{(n)},\kappagamma^{(n)})=0$ and therefore $u^{(n+\frac12)}=u^{(n)}$).
The key difference between (<ref>) and (<ref>) lies in the fact that in (<ref>), the state $u^{(n)}$ has to be precomputed as a solution to the nonlinear state equation $A(u^{(n)},\kappagamma^{(n)})=0$ for given $\kappagamma^{(n)}$, whereas in (<ref>), $u^{(n)}$ just comes from the previous iterate according to (<ref>) (with $n$ replaced by $n-1$). Instead, in (<ref>), $u^{(n+\frac12)}$ has to be precomputed from a linear state equation (<ref>). This would make a considerable difference between the two methods in case of a fully nonlinear model. However here, due to the fact that the model decouples into linear subproblems, the difference in computational effort is insignificant.

Gradient type methods
A Landweber step for solving (<ref>) is defined by a gradient descent step for the least squares functional

\[
\|F(\kappagamma^{(n)})-y\|_{L^2}^2
\]

i.e., by

\[
\begin{aligned}
\kappagamma^{(n+1)} &= \kappagamma^{(n)}+\mu {T^{(n)}}^* (y-F(\kappagamma^{(n)})
\end{aligned}
\]

with an appropriately chosen step size $\mu$.

We just mention that Landweber iteration could as well be applied to the all-at-once version (<ref>) and would completely avoid PDE solutions, thus itself act as an iterative Helmholtz solver (however, probably a very slow one).

Some remarks on the implementation are in order. For details we refer to, e.g., [13].

In the Newton type methods described in Tables <ref>, <ref>, we use a discretization $\kappa^{(n)}(x)=\sum_{i=1}^I a_i d^i(x)$, $\gamma^{(n)}(x)=\sum_{j=1}^J b_j e^j(x)$ on subspaces of $L^2(\widetilde{\Omega})$, by bases $\{d^1,\ldots,d^I\}$, $\{e^1,\ldots,e^J\}$.

Choice of $\alpha^{(n)}$
The regularization parameter in the Newton type methods may be simply chosen along a geometric sequence $\alpha^{(n)}=c \rho^n$ for some $c>0$, $\rho\in(0,1)$ in case of the IRGNM versions (both reduced and all-at-once).
For the Levenberg-Marquardt method, the choice is somewhat more complicated, namely it has to balance nonlinear and linearized residual in the sense of an inexact Newton method such that

\[
\underline{\theta}\|F(\kappagamma^{(n)}-y\|\leq \mbox{res}(\alpha^{(n)}) \leq\overline{\theta}\|F(\kappagamma^{(n)}-y\|
\]

for some constants $0<\underline{\theta}<\overline{\theta}<1$, where

$\mbox{res}(\alpha)=\|F(\kappagamma^{(n)})+T^{(n)}(\kappagamma^{(n+1)}(\alpha)-\kappagamma^{(n)})-y\|$, cf. [22].

Stopping rule
To avoid unbounded propagation of the measurement noise through the iterations, the methods defined above

have to be stopped at an appropriate index $n$. A widely used and well-investigated method for this is the discrepancy principle, which for a given noise level $\delta$ and a safety factor $\tau>1$ defines $n$ as the first index such that

\[
\|F(\kappagamma^{(n)})-y\|\leq\tau\delta\,.
\]

Multiple observations
As we have seen in Section <ref>, unique recovery of even just one of the two coefficients $c$ and $\gamma$ requires boundary measurements for several frequencies - a fact that is evident from a simple dimension count.
Also the fact that the focal point where the high frequency beams interact is moved through the region of interest should be taken into account by incorporating multiple excitations.
This corresponds to using Neumann conditions $g_k^\ell$ at transducer locations $\Sigma_k^\ell$ for $\ell\in\{1,\ldots,L\}$.
Finally, several receiver array locations $\Gamma^m$, $m\in\{1,\ldots,M\}$, might be used to recover a single pair of $\kappa$ and $\gamma$.
Thus, we actually deal with a set of several model and observation operators $A^\ell$, $\ell\in\{1,\ldots,L\}$, $C^m$, $m\in\{1,\ldots,M\}$ respectively. Labelling the resulting forward operators $F_p=C^m \circ S^\ell$ and data $y_p$ for $p=(m-1)L+\ell$, we can write the inverse problem of reconstructing $\kappagamma$ as a system of operator equations

\begin{equation}\label{Fxy_sys}
F_p(\kappagamma)=y_p \, \quad p\in\{1,\ldots,P=L\cdot M\}
\end{equation}

and apply Kaczmarz type methods as follows:

(a) parallelly apply one step of an iterative reconstruction method to each of the equations in (<ref>) and then combine the resulting reconstructions $\kappagamma_p$ in a proper way, e.g.,

\begin{equation}\label{Kacz_par}
\begin{aligned}
\kappagamma_p^{(n+1)}&=\kappagamma^{(n)}\!\!\!+G_p(\kappagamma^{(n)}\!\!\!,F_p,y_p),\ p=1,\ldots P, \\
\kappagamma^{(n+1)}&=\frac{1}{P}\sum_{p=1}^P \kappagamma_p^{(n+1)} %\,, \quad n=1,2,3,\ldots
\end{aligned}\end{equation}

(b) sequentially perform one step of an iterative reconstruction method in a
cyclically repeated manner

\begin{equation}\label{Kacz_seq}
\kappagamma^{(n+1)}=\kappagamma^{(n)}+G_p(\kappagamma^{(n)},F_p,y_p)
\end{equation}

where $p=\text{mod}(n-1,P)+1$,
(the order in which the indices $p$ are addressed could as well be randomized)

In here, $G_p(\kappagamma^{(n)},F_p,y_p)$ is defined by one of the Newton or gradient steps

from Tables <ref>, <ref>, <ref> below.

For a pseudocode description of the methods discussed above,

see Tables <ref>, <ref>, <ref>

§ OUTLOOK

In this paper we have made some first steps towards putting the problem of vibroacoustic imaging into the mathematical framework of inverse problems and regularization. We have presented a model in frequency domain, proven uniqueness of recovery of the spatially varying nonlinearity parameter $\gamma(x)$ from pressure measurements at multiple frequencies, and derived Newton as well as gradient based reconstruction methods.

Natural next steps are on one hand to refine and implement the devised numerical methods and on the other hand to answer important analytical questions. Among the latter, there is uniqueness of
simultaneous reconstruction of $c(x)$ and $\gamma(x)$.
To this end, the use of multiple excitation locations (instead of or in addition to multiple frequencies), corresponding to shifting the focus of the interacting high-frequency beams around the region of interest, needs to be further investigated.
Moreover, a priori information should be taken into account. Indeed, an important special case is the one of piecewise constant coefficients, in which only the shapes of finitely many subdomains and finitely many values of $c$ and $\gamma$ are to be found: Here one would expect uniqueness even from  boundary data at just a few frequencies, resulting from appropriately chosen excitations.

A computational framework for the reconstruction of piecewise constant coefficients could be based on the by now standard approach of alternatingly recovering the support and the value of inclusions in a homogeneous background. For a simultaneous recovery of both support and value, the known advantages of total variation regularization can be made use of. In case of known parameter values, also regularization by bound constraints (using the known values as bounds) is a promising approach [23].

Concerning forward simulation, we point to the fact that the high frequency waves $\phi_1$, $\phi_2$ have a strongly preferred direction of propagation, which can justify the use of a parabolic approximation, cf., e.g., [24].
Indeed, for efficient numerical simulation a decomposition approach has been devised in [7, 8] that splits the forward problem into a three components: (a) directed high frequency propagation of the two beams described by $\phi_1$, $\phi_2$,
(b) nonlinear interaction of these at the focal point, and
(c) undirected low frequency propagation to the measurement array via $\psi$.
This could also be implemented in our framework;

the adjoint equations for Landweber iteration would have to be re-derived for this purpose.

Also the model itself might have to be modified. Besides the use of a parabolic approximation in phase (a), also fractional damping e.g., [25, 19] is relevant in ultrasonics.

§ ACKNOWLEDGMENT

The author wishes to thank Alison Malcolm, University of New Foundland for fruitful discussions that have in fact been crucial in setting up the model considered in this paper.
The work was supported by the Austrian Science Fund fwf
under the grants P30054 and DOC78.

#1ISBN #1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1<https://doi.org/#1>et al.#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1#1<><#>1#1#1#1#1#1#1#1#1#1#1#1#1#1PreBibitemsHook

[1]

Fatemi, M.,
Greenleaf, J.F.:
Ultrasound-stimulated vibro-acoustic spectrography.

[2]

Fatemi, M.,
Greenleaf, J.F.:
Vibro-acoustography: An imaging modality based on
ultrasound-stimulated acoustic emission.
Proceedings of the National Academy of Sciences

[3]

Bjørnø, L.:
Characterization of biological media by means of their non-linearity.

[4]

Cain, C.A.:
Ultrasonic reflection mode imaging of the nonlinear parameter B/A:
I. a theoretical basis.
The Journal of the Acoustical Society of America

[5]

Ichida, N.,
Sato, T.,
Linzer, M.:
Imaging the nonlinear ultrasonic parameter of a medium.
Ultrasonic Imaging
PMID: 6686896

[6]

Byram, B.C.,
Ruiter, N.V. (eds.):
Medical Imaging 2021: Ultrasonic Imaging and Tomography.
Proceedings of SPIE,
vol. 11602.
Bellingham WA

[7]

Malcolm, A.E.,
Reitich, F.,
Yang, J.,
Greenleaf, J.F.,
Fatemi, M.:
Numerical modeling for assessment and design of ultrasound
vibro-acoustography systems.
In: Biomedical Applications of Vibration and Acoustics for Imaging and
ASME Press,
New York

[8]

Malcolm, A.E.,
Reitich, F.,
Yang, J.,
Greenleaf, J.F.,
Fatemi, M.:
A combined parabolic-integral equation approach to the acoustic
simulation of vibro-acoustic imaging.

[9]

Westervelt, P.J.:
Parametric acoustic array.
The Journal of the Acoustic Society of America

[10]

Kaltenbacher, M.:
Numerical Simulation of Mechatronic Sensors and Actuators: Finite
Elements for Computational Multiphysics,
3rd edn.

[11]

Tröltzsch, F.:
Optimal Control of Partial Differential Equations: Theory, Methods,
and Applications.
Graduate studies in mathematics.
American Mathematical Society,
Providence, RI

[12]

Hanke, M.,
Neubauer, A.,
Scherzer, O.:
A convergence analysis of the Landweber iteration for nonlinear
ill-posed problems.
Numer. Math.

[13]

Kaltenbacher, B.,
Neubauer, A.,
Scherzer, O.:
Iterative Regularization Methods for Nonlinear Problems.
de Gruyter,
Berlin, New York
Radon Series on Computational and Applied Mathematics

[14]

Kaltenbacher, B.,
Rundell, W.:
On an inverse problem of nonlinear imaging with fractional damping.
Mathematics of Computation
to appear; see also arXiv:2103.08965 [math.AP]

[15]

Kaltenbacher, B.:
Periodic solutions and multiharmonic expansions for the Westervelt
Evolution Equations and Control Theory EECT
19 pages

[16]

Isakov, V.:
Inverse Source Problems.
Mathematical Surveys and Monographs,
vol. 34,
p. 193.
American Mathematical Society,
Providence, RI

[17]

Kaltenbacher, B.,
Rundell, W.:
On the identification of the nonlinearity parameter in the Westervelt
equation from boundary measurements.
Inverse Problems & Imaging

[18]

Yamamoto, M.,
Kaltenbacher, B.:
An inverse source problem related to acoustic nonlinearity parameter
In: Kaltenbacher, B.,
Wald, A.,
Schuster, T. (eds.)
Time-dependent Problems in Imaging and Parameter Identification.
New York

[19]

Kaltenbacher, B.,
Rundell, W.:
Some inverse problems for wave equations with fractional derivative
Inverse Problems

[20]

Kirsch, A.:
An Introduction to the Mathematical Theory of Inverse Problems.
New York

[21]

Isakov, V.:
Inverse Problems for Partial Differential Equations.
New York
2nd Edition

[22]

Hanke, M.:
A regularization Levenberg–Marquardt scheme, with applications to
inverse groundwater filtration problems.
Inverse Problems

[23]

Hungerländer, P.,
Kaltenbacher, B.,
Rendl, F.:
Regularization of inverse problems via box constrained minimization.
Inverse Problems and Imaging

[24]

Tappert, F.D.:
In: Keller, J.B.,
Papadakis, J.S. (eds.)
The parabolic approximation method,
pp. 224–287.
Berlin, Heidelberg

[25]

Cai, W.,
Chen, W.,
Fang, J.,
Holm, S.:
A Survey on Fractional Derivative Modeling of Power-Law Frequency-Dependent
Viscous Dissipative and Scattering Attenuation in Acoustic Wave Propagation.
Applied Mechanics Reviews


Eungbeom Kim1⋆, Yunkee Chae1⋆, Jaeheon Sim1, Kyogu Lee1,2††thanks: ⋆ Equal
contribution. This work was partly supported by Institute of Information &
communications Technology Planning & Evaluation (IITP) grant funded by the
Korea government(MSIT) [No. 2022-0- 00641, XVoice: Multi-Modal Voice Meta
Learning, 90%] and [NO.2021-0-01343, Artificial Intelligence Graduate School
Program (Seoul National University), 10%].

# Debiased Automatic Speech Recognition for Dysarthric Speech
via Sample Reweighting with Sample Affinity Test

###### Abstract

Automatic speech recognition systems based on deep learning are mainly trained
under empirical risk minimization (ERM). Since ERM utilizes the averaged
performance on the data samples regardless of a group such as healthy or
dysarthric speakers, ASR systems are unaware of the performance disparities
across the groups. This results in biased ASR systems whose performance
differences among groups are severe. In this study, we aim to improve the ASR
system in terms of group robustness for dysarthric speakers. To achieve our
goal, we present a novel approach, sample reweighting with sample affinity
test (Re-SAT). Re-SAT systematically measures the debiasing helpfulness of the
given data sample and then mitigates the bias by debiasing helpfulness-based
sample reweighting. Experimental results demonstrate that Re-SAT contributes
to improved ASR performance on dysarthric speech without performance
degradation on healthy speech.

Index Terms: speech recognition, debiasing, dysarthric speech

##  1 Introduction

Automatic speech recognition (ASR) performance across groups should be fair
regardless of disorder, race, age, and dialect for a trustworthy and inclusive
system. Although ASR has been improved along with the success of deep
learning, ASR systems based on deep learning tend to be vulnerable to biases
[1, 2]. Consequently, the bias in ASR systems interferes with the ASR system's
trustworthiness by causing poor worst group performance. Applying ASR on
dysarthria, which is a type of motor speech disorder, also suffers from the
lack of group robustness between healthy and dysarthric speakers. In other
words, a dysarthric speaker group has a lower performance than a healthy
speaker group on ASR. This is because dysarthric speech has low
intelligibility, whereas deep learning has a tendency to be easily fitted to
shortcuts [3, 4]. However, most of the previous dysarthric speech recognition
studies mainly focus on the sole performance of the ASR system on dysarthric
speech. In this study, we aim to improve the performance of ASR on dysarthric
speech from a debiasing point of view. To the best of our knowledge, this is
the first work to endeavor to make a debiased ASR system for dysarthric
speech.

In this regime, we focus on a debiasing method based on sample reweighting.
Recently, sample reweighting such as [5, 6] is proposed as a promising
solution to handle the bias problem. Sample reweighting methods mainly consist
of 1) estimating the helpfulness of a sample for debiasing and 2) reweighting
the sample based on helpfulness. To estimate the debiasing helpfulness, [5]
defines bias-conflicting sample and bias-aligned sample which denote an
incorrectly classified sample and a correctly classified sample from an
unintended decision rule of a biased model, respectively. The strategies
utilizing those definitions are based on the assumption that upweighting the
data samples with a large loss from the biased model, i.e. bias-conflicting
samples, leads to improved generalization performance of underfitted groups.
This assumption, however, hardly holds in real-world problems. For example,
although an outlier might have a large loss value, it is an undesirable
strategy to upweight the outlier.

Therefore, we define a new taxonomy to directly categorize the samples by
debiasing helpfulness; bias-blocking sample denotes a sample that mitigates
the model's bias, and bias-accelerating sample denotes a sample that
accelerates the model's bias when the model is trained on each sample. To
estimate whether a given sample is bias-blocking sample, we propose a sample
affinity test (SAT) which systematically measures debiasing helpfulness. SAT
is based on the novel metric dubbed sample affinity, which denotes a training
effect of the given sample to the other samples' loss. This is inspired by
task affinity [7] which estimates inter-task affinity using the loss shift of
each task after lookahead training of a given task. Intuitively, SAT measures
debiasing helpfulness using sample affinity of the given sample on the bias-
conflicting sample set, based on a one-step lookahead training of a given
sample. Unlike the loss-based methods, SAT can filter unhelpful bias-
accelerating samples even if the samples are bias-conflicting.

Based on SAT, we propose a novel sample reweighting method, Re-SAT, to
challenge the debiasing problem of dysarthric speech recognition. Re-SAT
consists of four sequential components for implementation. The first component
estimates bias-conflicting samples based on loss, following [6]. Secondly, the
SAT component is activated based on the estimated bias-conflicting samples to
accurately measure the debiasing effect of each sample. In the third
component, the SAT result is normalized through sorting and mapped to the
weights. Finally, Re-SAT trains the models with the reweighted samples.

In summary, we observe the biased performance of ASR on disordered speech and
analyze the ASR system through the lens of debiasing. Furthermore, we present
sample affinity test (SAT) to directly measure the debiasing helpfulness of
samples. We also propose Re-SAT, a novel method for sample reweighting that is
available under the group label-free environment and mitigates the margin
between the ASR performances of healthy speakers and dysarthric speakers.

##  2 Related work

Dysarthric speech recognition ASR on dysarthric speech still remains a
challenging problem. Since the disordered speech dataset has limited
scalability due to the difficulty of data collection, data augmentation on
dysarthric speech has been widely studied [8, 9, 10, 11, 12]. Other
investigations adopt pre-trained self-supervised learning models [13, 14, 15,
16] for improved disordered speech recognition. On the other hand, [17, 18,
19] utilize some prior knowledge that a dysarthric speaker has distinct
characteristics compared to healthy or other dysarthric speakers by leveraging
speaker adaptation using the spectro-temporal level. In this study, we
approach ASR on dysarthric speech from a different point of view, debiasing,
which has never been used to the best of our knowledge.

Debiasing Debiasing is a challenging but essential area for fair and inclusive
deep learning. Group robustness, which denotes an ability to perform
satisfactorily across different groups in terms of the given task, is one of
the key factors of a debiased model. We focus on group robustness without
group annotations, for a wider real-world application. Sample reweighting is
one of the promising methods for debiasing. While empirical risk minimization
(ERM) uniformly averages all of the performances across the given samples,
sample reweighting methods such as Learning from Failure (LfF) [5] and Just
Train Twice (JTT) [6] upweight the samples that are expected to belong to the
poor performance group called bias-conflicting samples. These methods are
based on the intuition that upweighting the bias-conflicting samples with
large losses leads to a debiased model, which is not always true. In
particular, we test ASR systems on dysarthric speakers, each of who can be
regarded as a respective group due to the variety of dysarthria. In this
complex real-world environment, reweighting the samples using a fine-grained
metric is an important issue. For this reason, we aim to propose direct and
accurate criteria for reweighting beyond the estimation of bias-conflicting
samples based on loss.

Task affinity Given a set of tasks in a multi-task learning setup, task
grouping aims to find the clusters of tasks that should be trained
simultaneously for improved performance. Task affinity [7] is proposed to
address the task grouping problems. To examine the task affinity of task $i$
on task $j$, they leverage a lookahead update on task $i$. Then, they compute
the loss shift after the lookahead update to explore the effect of task $i$ on
$j$. Inspired by task affinity, we propose sample affinity to investigate the
inter-sample effect for debiasing.

Figure 1: Illustration of Re-SAT for debiased dysarthric speech recognition.
Re-SAT consists of 1) bias-conflicting sample estimation, 2) sample affinity
test, 3) normalized reweighting, and 4) training blocks.

##  3 Method

In this section, we introduce our sample reweighting with sample affinity test
(Re-SAT), which is proposed to challenge the ASR on dysarthric speech. Re-SAT
consists of four components: 1) bias-conflicting sample estimation, 2) sample
affinity test, 3) normalized reweighting and 4) training, as shown in Figure
1. The data samples are reweighted through the first three components, and
then Re-SAT trains the model with the reweighted data samples in the 4)
training block. In summary, Re-SAT upweights the bias-blocking samples and
downweights the bias-accelerating samples using sample affinity test for
debiasing. Details are as follows.

###  3.1 Bias-conflicting sample estimation

At the first step, Re-SAT estimates the bias-conflicting samples using the
losses of samples in each batch. Re-SAT regards the samples with the largest
$K$ losses in a batch as the bias-conflicting samples. Unlike the previous
research [6] which fixes the estimation results, Re-SAT keeps estimating the
bias-conflicting samples during training to accurately reflect the current
state of the model. That is, the estimation result for each sample can be
modified throughout training.

Although Re-SAT estimates the bias-conflicting samples, Re-SAT does not map
the estimation results for reweighting because the large loss does not
guarantee the debiasing ability, as we introduced in the above sections. We
address this issue by designing a novel criterion that directly estimates the
bias-blocking samples, for debiasing ability.

###  3.2 Sample affinity test

We propose a novel sample affinity test (SAT) to evaluate the debiasing
ability of samples to filter the bias-accelerating samples. Given the model
$f_{\theta}$, the samples in the mini-batch $x_{1},...,x_{N}$, and the
estimated bias-conflicting samples $\hat{b}_{1},...,\hat{b}_{K}$, SAT utilizes
single step lookahead updating with respect to $x_{i}$ as

$\theta^{LA}_{x_{i}}=\theta-\eta\nabla L(x_{i};f_{\theta}),$ (1)

for $i=1,...,N$ where $\eta$ is learning rate, $N$ is batch size, and $L$ is a
loss function. After the lookahead updating process, SAT compares the loss of
1) the lookahead model $\theta^{LA}_{x_{i}}$ and 2) the original model
$\theta$ on the bias-conflicting samples to calculate the averaged sample
affinity:

$\text{SA}(x_{i}\rightarrow\\{\hat{b}_{1},...,\hat{b}_{K}\\};f_{\theta})=\frac{1}{K}\sum_{\forall
k}\left(1-\frac{L(\hat{b}_{k};f_{\theta^{LA}_{x_{i}}})}{L(\hat{b}_{k};f_{\theta})}\right).$
(2)

We present sample affinity inspired by task affinity [7] for debiasing. Note
that SAT computes the sample affinity every step to take into account the
current state of the model, while task affinity is averaged across the whole
training step.

###  3.3 Normalized reweighting

Sample affinity on the bias-conflicting samples approximates the bias-blocking
effect of a given sample. However, sample affinity is not an absolute but
relative score which depends on the bias-conflicting sample set, the current
state of the model, and the maturity of the model. Therefore, Re-SAT does not
directly map the sample affinity to weights. Instead, Re-SAT normalizes the
sample affinity by extracting the rank of the sample affinity with respect to
the descending order within the batch. This normalization stabilizes unwanted
shifts in the learning rate.

For sorted samples $x_{1},...,x_{r},...,x_{N}$ with respect to sample affinity
in descending order, Re-SAT finally reweights the samples using the function
$w:\mathbb{N}\rightarrow\mathbb{R}^{+}$, which is defined as follows:

$w(r)=\frac{\exp(s(N-r)/(N-1))}{\sum_{r=1}^{N}\exp(s(N-r)/(N-1))}$ (3)

for $r=1,...,N$ where $r$ is the rank of the sample $x_{r}$ and $s$ is the
constant hyperparameter. For our experiments, we set batch size $N=32$ and
$s=4$.

###  3.4 Training

In this component, Re-SAT trains the model $f_{\theta}$ with the reweighted
data samples. Unlike the vanilla training that averages the loss of each data
sample, Re-SAT leverages weighted average for total loss as:

$\theta^{\prime}=\theta-\eta\nabla\dfrac{1}{N}\sum_{r=1}^{N}w(r)L(x_{r};f_{\theta}).$
(4)

for sorted samples $x_{1},...,x_{N}$ with respect to sample affinity
$\text{SA}(x_{i})$ in descending order.

##  4 Experiments

###  4.1 Dataset

UASpeech corpus [20] is used for our experiments. UASpeech corpus is one of
the largest English dysarthric speech datasets, composed of speech from 15
dysarthric speakers and 13 healthy control speakers. Each dysarthric speaker
is classified into very low, low, mid, and high levels of intelligibility. The
dataset consists of three blocks B1, B2, and B3. Each block contains 155
common words and 100 uncommon words spoken by dysarthric and healthy control
speakers. We used block 1 and block 3 of all healthy and dysarthric speech
recorded on microphone M5 as a training set, and all data in block 2 as the
test set. Note that we used not only dysarthric speech but also healthy speech
as the test set. Basically, UASpeech provides the speech data denoised by
noisereduce [21, 22]; therefore we used this version for our experiments. We
excluded the uncommon words from both the training and test set and no
additional data augmentation was conducted.

###  4.2 Model

We fine-tuned Whisper [23], the recent state-of-the-art ASR model, trained on
a large amount of labeled audio-transcription data. It employs a simple
encoder-decoder Transformer architecture. We used the Whisper-tiny, pre-
trained on English only with 39M parameters, which is available at HuggingFace
transformers repository [24].

###  4.3 Experimental Setups

We trained the model using AdamW optimizer [25] with a learning rate of
$10^{-5}$, weight decay of 0.1, and batch size of 32 for 30 epochs. For the
proposed method, Re-SAT, we investigated the effect of the number of bias-
conflicting samples in batch as $K\in\\{2,4,8,16\\}$. For comparison, we
tested the JTT method [6] with upweighting value 25 and bias-conflicting
sample identification epoch 3. For the ablation studies, we tested the model
named Re-Loss, which utilizes the loss-based sample reweighting instead of
sample affinity and follows Re-SAT in other settings, to figure out the impact
of the sample affinity test.

##  5 Results

|  | WER (%)
---|---|---
|  | Speaker
---
(Intelligibility %)
ERM | JTT | Re-Loss |  | Re-SAT
---
(Ours)
VL | M04 (2%) | 84.39 | 79.87 | 82.45 | 80.26
F03 (6%) | 44.33 | 42.21 | 40.28 | 38.25
M12 (7%) | 44.41 | 42.04 | 41.72 | 35.38
M01 (17%) | 50.97 | 51.77 | 47.10 | 44.19
Avg. | 54.66 | 52.46 | 51.50 | 48.09
L | M07 (28%) | 16.22 | 23.78 | 16.96 | 16.04
F02 (29%) | 19.72 | 22.86 | 16.59 | 15.94
M16 (43%) | 22.04 | 23.55 | 20.32 | 19.57
Avg. | 19.19 | 23.39 | 17.84 | 17.06
M | M05 (58%) | 12.63 | 13.55 | 12.90 | 12.72
M11 (62%) | 14.30 | 12.80 | 12.04 | 10.11
F04 (62%) | 13.88 | 13.88 | 13.22 | 14.64
Avg. | 13.57 | 13.44 | 12.75 | 12.59
H | M09 (86%) | 9.86 | 11.15 | 8.66 | 8.76
M14 (90%) | 12.17 | 13.73 | 11.98 | 9.86
M10 (93%) | 3.50 | 3.32 | 2.40 | 2.76
M08 (95%) | 4.61 | 5.53 | 4.70 | 3.78
F05 (95%) | 2.03 | 4.33 | 2.03 | 2.76
Avg. | 6.43 | 7.61 | 5.95 | 5.59
| Avg.
---
(Dysarthric)
21.49 | 22.25 | 20.15 | 19.93
| Avg.
---
(Healthy)
3.81 | 4.21 | 3.75 | 3.40
| Avg.
---
(Overall)
12.93 | 13.51 | 12.20 | 12.08
Table 1: ASR results in terms of Word Error Rate (WER) on UASpeech corpus. VL/L/M/H refer to very low/low/mid/high intelligibility groups, respectively. The percentage of intelligibility is measured based on how accurately naive human listeners can transcribe isolated words produced by speakers, according to [20]. $K$ of Re-SAT is set to 4.  | Intelligibility group |
---|---|---
K | VL | L | M | H |  | Avg.
---
(D)
| Avg.
---
(H)
Avg.
ERM | 54.66 | 19.19 | 13.57 | 6.43 | 21.49 | 3.81 | 12.93
2 | 50.12 | 17.49 | 14.15 | 5.75 | 19.93 | 3.72 | 12.08
4 | 48.09 | 17.06 | 12.59 | 5.59 | 19.05 | 3.40 | 11.47
8 | 49.44 | 17.61 | 12.62 | 6.36 | 19.75 | 3.84 | 12.05
16 | 55.31 | 21.10 | 13.60 | 6.91 | 22.21 | 3.96 | 13.36
Table 2: Comparison on hyperparameter $K$ with word error rate (WER). D and H
refer to dysarthric and healthy speech, resepctively.

We investigated with the effect of the proposed method, Re-SAT, by comparing
it to empirical risk minimization (ERM) and JTT for debiasing. In contrast to
most of the dysarthric speech recognition studies that focus only on
dysarthric speech, we investigated the effect of Re-SAT to ASR system on
healthy speech for a fair comparison. This is important for debiasing because
the performance disparity on dysarthric speech and healthy speech determines
the group robustness of the systems.

Results in Table 1 show word error rate (WER) of the models on each speaker
and intelligibility group. It is shown that Re-SAT contributes to performance
improvements for each group compared to ERM by relatively decreasing 12.02%
(54.66 $\rightarrow$ 48.09), 11.10% (19.19 $\rightarrow$ 17.06), 7.22% (13.57
$\rightarrow$ 12.59), and 13.06% (6.43 $\rightarrow$ 5.59) in terms of the
averaged WER for the very low, low, mid, and high intelligibility group,
respectively. Surprisingly, Re-SAT also surpasses ERM on the healthy speaker
group. From a speaker-wise view, it is observed that the performances of 12
out of 15 speakers show enhanced results in terms of WER on the Re-SAT-based
ASR model. These results demonstrate the robust performance gain of Re-SAT
even under a diverse group environment. On the contrary, the other debiasing
method, JTT, increases the averaged WER on healthy speech, low and high
intelligibility speech although JTT reduces the very low and mid
intelligibility group's WER, compared to those of ERM.

Figure 2: Illustrations of the averaged loss ranks sorted by descending order
in the mini-batch for each intelligibility group. The x-axis denotes training
epochs and the y-axis denotes the loss rank in mini-batch.

Figure 2 shows the rank sorted in descending order of training loss within the
batch with respect to each intelligibility group. For very low, low, and mid-
intelligibility groups that are vulnerable to bias, the increasing and large
rank value is desired because it means that the loss for the biased group is
decreased, even though it is approximated by training loss. Interestingly,
debiasing of JTT works only on the single worst (very low) intelligibility
group and even deteriorates the low and mid groups' ranks. This problem can be
caused by the structure of JTT's 1) binary classification of bias-conflicting
samples and 2) which are fixed at an early stage. On the other hand, Re-SAT
leverages 1) fine-grained reweighting and 2) updated bias-conflicting set
through training and shows desirable results across diverse groups. Figure 2
also demonstrates that ERM is inappropriate for debiasing. Although the ranks
on very low, low, and mid intelligibility groups of ERM increase at the early
stage for all training schemes, the ranks decrease at the late stage of ERM.
Even worse, the healthy speaker's rank shows the highest rank compared to
other methods, which means the biased performance of ERM. We also explore the
Re-Loss, which substitutes the sample affinity-based reweighting of Re-SAT
with loss-based reweighting, to investigate the importance of the sample
affinity test. Although Re-Loss outperforms ERM and JTT by a large margin, Re-
SAT still shows a surpassing performance. This demonstrates the precise bias-
blocking sample estimation ability of Re-SAT, as we expected.

Intuitively, selecting decent bias-conflicting samples takes an important role
in Re-SAT. For this reason, we investigate the size of the selected bias-
conflicting sample set $K$ in Re-SAT as shown in Table 2. Re-SAT shows stable
improvements across the intelligibility groups where $K=2,4,8$ compared to ERM
for all intelligibility groups. This supports our design choice of Re-SAT
which re-estimates the bias-conflicting samples through training. Even though
Re-SAT selects small bias-conflicting samples, Re-SAT can re-estimate the new
underfitted data sample as a bias-conflicting sample through training. On the
contrary, Re-SAT with $K=16$ even degrades the ASR results, since setting a
large $K$ as 16 interferes with the accurate estimation of the bias-
conflicting samples.

##  6 Conclusions

In this study, we address a debiasing problem of ASR on dysarthric speech
towards a fair and inclusive ASR system. In contrast to the previous research
that focuses only on the ASR performance for dysarthric speakers, we explore
the fairer validation system by analyzing the performance of healthy and
dysarthric speech at the same time. To achieve our goal, we propose a novel
debiasing method based on sample reweighting, Re-SAT. The ASR system using Re-
SAT surpasses the other baselines across the diverse dysarthric speakers and
shows a robust performance gain over various hyperparameters. For future work,
it is worth exploring the integrated system of Re-SAT and other dysarthric
speech recognition systems by leveraging the versatile structure of Re-SAT.

## References

  * [1] J. L. Martin and K. E. Wright, ``Bias in automatic speech recognition: The case of african american language,'' _Applied Linguistics_ , 2022.
  * [2] S. Feng, O. Kudina, B. M. Halpern, and O. Scharenborg, ``Quantifying bias in automatic speech recognition,'' _arXiv preprint arXiv:2103.15122_ , 2021.
  * [3] S. Beery, G. Van Horn, and P. Perona, ``Recognition in terra incognita,'' in _Proceedings of the European conference on computer vision (ECCV)_ , 2018, pp. 456–473.
  * [4] S. Asgari, A. Khani, F. Khani, A. Gholami, L. Tran, A. Mahdavi-Amiri, and G. Hamarneh, ``Masktune: Mitigating spurious correlations by forcing to explore,'' in _Advances in Neural Information Processing Systems_ , A. H. Oh, A. Agarwal, D. Belgrave, and K. Cho, Eds., 2022.
  * [5] J. Nam, H. Cha, S. Ahn, J. Lee, and J. Shin, ``Learning from failure: De-biasing classifier from biased classifier,'' _Advances in Neural Information Processing Systems_ , vol. 33, pp. 20 673–20 684, 2020.
  * [6] E. Z. Liu, B. Haghgoo, A. S. Chen, A. Raghunathan, P. W. Koh, S. Sagawa, P. Liang, and C. Finn, ``Just train twice: Improving group robustness without training group information,'' in _International Conference on Machine Learning_. PMLR, 2021, pp. 6781–6792.
  * [7] C. Fifty, E. Amid, Z. Zhao, T. Yu, R. Anil, and C. Finn, ``Efficiently identifying task groupings for multi-task learning,'' _Advances in Neural Information Processing Systems_ , vol. 34, pp. 27 503–27 516, 2021.
  * [8] M. Geng, X. Xie, S. Liu, J. Yu, S. Hu, X. Liu, and H. Meng, ``Investigation of data augmentation techniques for disordered speech recognition,'' _Proc. Interspeech_ , 2022.
  * [9] C. Bhat, A. Panda, and H. Strik, ``Improved asr performance for dysarthric speech using two-stage dataaugmentation,'' _Proc. Interspeech_ , pp. 46–50, 2022.
  * [10] Z. Jin, M. Geng, X. Xie, J. Yu, S. Liu, X. Liu, and H. Meng, ``Adversarial data augmentation for disordered speech recognition,'' _Proc. Interspeech_ , 2021\.
  * [11] Z. Jin, X. Xie, M. Geng, T. Wang, S. Hu, J. Deng, G. Li, and X. Liu, ``Adversarial data augmentation using vae-gan for disordered speech recognition,'' _arXiv preprint arXiv:2211.01646_ , 2022.
  * [12] B. Vachhani, C. Bhat, and S. K. Kopparapu, ``Data augmentation using healthy speech for dysarthric speech recognition.'' _Proc. Interspeech_ , 2018.
  * [13] M. K. Baskar, T. Herzig, D. Nguyen, M. Diez, T. Polzehl, L. Burget, J. Černockỳ _et al._ , ``Speaker adaptation for wav2vec2 based dysarthric asr,'' _Proc. Interspeech_ , 2022.
  * [14] S. Hu, X. Xie, Z. Jin, M. Geng, Y. Wang, M. Cui, J. Deng, X. Liu, and H. Meng, ``Exploring self-supervised pre-trained asr models for dysarthric and elderly speech recognition,'' _arXiv preprint arXiv:2302.14564_ , 2023.
  * [15] L. P. Violeta, W.-C. Huang, and T. Toda, ``Investigating self-supervised pretraining frameworks for pathological speech recognition,'' _Proc. Interspeech_ , 2022.
  * [16] A. Hernandez, P. A. Pérez-Toro, E. Nöth, J. R. Orozco-Arroyave, A. Maier, and S. H. Yang, ``Cross-lingual self-supervised speech representations for improved dysarthric speech recognition,'' _Proc. Interspeech_ , 2022.
  * [17] M. Geng, X. Xie, R. Su, J. Yu, Z. Ye, X. Liu, and H. Meng, ``On-the-fly feature based speaker adaptation for dysarthric and elderly speech recognition,'' _arXiv preprint arXiv:2203.14593_ , 2022.
  * [18] M. Geng, S. Liu, J. Yu, X. Xie, S. Hu, Z. Ye, Z. Jin, X. Liu, and H. Meng, ``Spectro-temporal deep features for disordered speech assessment and recognition,'' _Proc. Interspeech_ , 2021.
  * [19] M. Geng, X. Xie, Z. Ye, T. Wang, G. Li, S. Hu, X. Liu, and H. Meng, ``Speaker adaptation using spectro-temporal deep features for dysarthric and elderly speech recognition,'' _IEEE/ACM Transactions on Audio, Speech, and Language Processing_ , vol. 30, pp. 2597–2611, 2022.
  * [20] H. Kim, M. Hasegawa-Johnson, A. Perlman, J. Gunderson, T. S. Huang, K. Watkin, and S. Frame, ``Dysarthric speech database for universal access research,'' in _Ninth Annual Conference of the International Speech Communication Association_ , 2008.
  * [21] T. Sainburg, M. Thielk, and T. Q. Gentner, ``Finding, visualizing, and quantifying latent structure across diverse animal vocal repertoires,'' _PLoS computational biology_ , vol. 16, no. 10, p. e1008228, 2020.
  * [22] T. Sainburg, ``timsainb/noisereduce: v1.0,'' Jun. 2019. [Online]. Available: https://doi.org/10.5281/zenodo.3243139
  * [23] A. Radford, J. W. Kim, T. Xu, G. Brockman, C. McLeavey, and I. Sutskever, ``Robust speech recognition via large-scale weak supervision,'' _arXiv preprint arXiv:2212.04356_ , 2022.
  * [24] T. Wolf, L. Debut, V. Sanh, J. Chaumond, C. Delangue, A. Moi, P. Cistac, T. Rault, R. Louf, M. Funtowicz _et al._ , ``Huggingface's transformers: State-of-the-art natural language processing,'' _arXiv preprint arXiv:1910.03771_ , 2019.
  * [25] I. Loshchilov and F. Hutter, ``Decoupled weight decay regularization,'' in _International Conference on Learning Representations_ , 2019.