Inverse scheme for sound source localization

The inverse problem as presented in [^1] is to reconstruct acoustic sources from microphone measurements. The wave equation in the frequency domain (Helmholtz equation) is solved with the corresponding given boundary conditions using the finite element method. The optimization based inverse scheme is based on minimizing a Tikhonov functional matching measured microphone signals with simulated ones. This method identifies the amplitude and phase information of the acoustic sources such that the prevailing sound field can be reconstructed with high level of accuracy.

Governing equations

Inverse Problem

We assume that the original geometry of the setup and Fourier-transformed acoustic pressure signals $p_i^\mathrm{ms}(\omega)$ ($i=1,...,M$) are given. We define a domain $\Omega$ that consists of a domain $\Omega_\mathrm{acou} \subset \Omega$ and $\Omega_\mathrm{sc} \subseteq \Omega_\mathrm{acou}$. We assume the sound sources to be located in $\Omega_\mathrm{sc}$ and on $\partial\Omega_\mathrm{sc}$ and the microphones to be at positions $\boldsymbol{x}_i \in \Omega_\mathrm{acou}$,.

At the outer boundary $\partial \Omega$ we can either apply sound hard boundary conditions that lead to full reflection or we apply boundary conditions that lead to $-$ in practice more realistic $-$ partial reflections. This can be accomplished e.g. by applying impedance boundary conditions or by modeling acoustic absorbers as a separate region as an equivalent fluid with complex material properties.

The forward problem in its strong form is defined in frequency domain as

\begin{align}\label{eqn:Helmholtz} \nabla\cdot\nabla p_\mathrm{a} + k^2 p_\mathrm{a} = \sigma^\mathrm{}\quad\text{in } \Omega, \end{align}

also known as the Helmholtz equation. The reconstruction of sound sources is performed for a fixed angular sound frequency $\omega$. In \eqref{eqn:Helmholtz} $p_\mathrm{a}$ denotes the acoustic pressure, $k = \omega/c_0$ the wave number whereby $c_0$ is the speed of sound. Since we solve the forward problem via the finite element method, we derive the weak form of \eqref{eqn:Helmholtz} by multiplying with an arbitrary (complex valued) test function $p'$ and integrating by parts \begin{align} \int\limits_{\Omega}\left(\nabla p_\mathrm{a} \cdot \nabla \bar{p}' - k^2 p_\mathrm{a} \bar{p}'\right) \,\mathrm{d}\boldsymbol{x} = \ -\int\limits_{\Omega} \sigma \bar{p}' \,\mathrm{d}\boldsymbol{x} + \int\limits_{\partial\Omega}\nabla p_\mathrm{a}\bar{p}'\cdot\mathrm{d} \boldsymbol{s}\,. \label{eqn:weakForm} \end{align}

Sound sources $\sigma$ are modeled as point sources as

$$\sigma = \sum_{n=1}^{N} a_n\mathrm{e}^{\,\mathrm{j}\varphi_n}\delta_{\boldsymbol{x}_n}\,,$$

with $N$ being the number of possible sources which corresponds to the number of nodes in $\Omega_\mathrm{sc}$, the amplitudes $a_n$ and the phases $\varphi_n$.

Optimization problem

Defining

$$A(p_\mathrm{a},p') := \Re\left(\int\limits_{\Omega}\left(\nabla p_\mathrm{a} \cdot \nabla \bar{p}' - k^2 p_\mathrm{a} \bar{p}'\right) \,\mathrm{d}\boldsymbol{x}\right)\,,$$

fitting of the parameters $a_j$ and $\varphi_j$ by means of Tikhonov regularization amounts to solving the constrained optimization problem

\begin{align} \min \ J(p,a,\varphi) \quad &\text{s.t. } \\\nonumber \quad A(p_\mathrm{a},p')=-&\Re\left( \sum\limits_{n=1}^{N^\mathrm{in}} a_n \mathrm{e}^{\,\mathrm{j}\varphi_n} \bar{p}'(\boldsymbol{x}_n) \right) \,. \end{align}

The cost functional $J$ is defined as

\begin{align} \label{eqn:functionalJ} \nonumber J(p,a,&\varphi)= \frac{\varPhi}{2} \sum\limits_{i=1}^M |p_\mathrm{a}(\boldsymbol x_i) - p^{\rm ms}_i |^2 \\ &+ \alpha \sum\limits_{n=1}^N | a_n|^q + \beta \sum\limits_{n=1}^N \varphi_n^2 \\ \nonumber &- \varrho \sum\limits_{n=1}^N \Bigl(\ln\left(\frac{\pi}{2}+\varphi_n\right)+\ln\left(\frac{\pi}{2}-\varphi_n\right)\Bigr), \end{align}

with appropriately chosen regularization parameters $\alpha$, $\beta$ and $\varrho$. The exponent $q\in (1,2]$ leads to sparse source reconstruction when chosen close to 1. The scaling factor $\varPhi$ should be chosen appropriately to ensure convergence and avoiding numerical rounding errors. If $\varPhi \neq 1$ the identified amplitudes are scaled back after the identification in order to receive the correct values for the unscaled problem.

The regularization parameters are defined in the xml-Scheme within <analysis> as follows:

<analysis>
<inverseSource>
...
<alphaPar>   alpha          </alphaPar>
<betaPar>    beta           </betaPar>
<rhoPar>     rho            </rhoPar>
<expPar>     q              </expPar>
<scale2Val>  Phi/max(p_ms)  </scale2Val>
...
</inverseSource>
</analysis>


Note that the scaling value $\Phi$ is not defined directliy but one has scale2Val$~= \frac{\Phi}{\max(p^{\rm ms}_i)}$.

The minimization is realized by a gradient method with Armijo line search. The maximum amount of reducing the step size of the line search is limited by an interior tolerance or a defined maximum amount $N_\mathrm{i,max}$ <maxNumLineSearch>. The gradient is computed using an adjoint approach. After reaching the desired tolerance $tol$ of the norm of the gradient or reaching the specified maximal amount of gradient steps $N_\mathrm{g,max}$ <maxGradSteps>, the regularization parameters $\alpha$, $\beta$ and $\varrho$ are multiplyed by $0.5$. These outer iteration steps are limited by the defined error $\varepsilon$ <resStopCritRel> which is defined by the difference between measured and simulated acoustic pressures or again a maximum number of steps $N_\mathrm{o,max}$ <maxReduceParSteps>.

The above mentioned stopping criteria are also defined within <analysis>:

<analysis>
<inverseSource>
...
<resStopCritRel>     epsilon  </resStopCritRel>
<maxReduceParSteps>  N_o,max  </maxReduceParSteps>
<maxNumLineSearch>   N_i,max  </maxNumLineSearch>
...
</inverseSource>
</analysis>


For further details on the optimization problem as well as a pseudo code of the implementation in openCFS see [^1]. For apllications in real measurement environments see [^2].

The pressure measurements are provided in a text document that includes the node number, real and imaginary part of the Fourier transformed acoustic pressure.

Node No. $\Re(p_{\mathrm{a},i})$ $\Im(p_{\mathrm{a},i})$
... ... ...

The path to this text document is again defined within <analysis>:

<analysis>
<inverseSource>
...
<measDataFilename>  pathToTxt  </measDataFilename>
...
</inverseSource>
</analysis>


Refernces

[^1]: M. Kaltenbacher, B. Kaltenbacher, and Stefan Gombots. Inverse scheme for acoustic source localization using microphone measurements and finite element simulations. Acta Acustica United With Acustica, 104:647–656, 2018. [^2]: Stefan Gombots. Acoustic source localization at low frequencies using microphone arrays. PhD thesis, TU Wien, Vienna, Austria, 2018.

[^1]: M. Kaltenbacher, B. Kaltenbacher, and Stefan Gombots. Inverse scheme for acoustic source localization using microphone measurements and finite element simulations. Acta Acustica United With Acustica, 104:647–656, 2018. [^2]: Stefan Gombots. Acoustic source localization at low frequencies using microphone arrays. PhD thesis, TU Wien, Vienna, Austria, 2018.