Bayesian Approach to Inverse Time-harmonic Acoustic Scattering with Phaseless Far-field Data

TL;DR

This paper introduces a Bayesian method for reconstructing sound-soft obstacles from phaseless far-field data in inverse acoustic scattering, utilizing MCMC techniques for improved numerical approximation.

Contribution

It develops a Bayesian framework for inverse scattering with phaseless data and demonstrates the effectiveness of pCN-MCMC algorithms for obstacle reconstruction.

Findings

Successful reconstruction of various obstacle shapes

Effective use of pCN-MCMC for convergence

Validation through numerical examples

Abstract

This paper is concerned with inverse acoustic scattering problem of inferring the position and shape of a sound-soft obstacle from phaseless far-field data. We propose the Bayesian approach to recover sound-soft disks, line cracks and kite-shaped obstacles through properly chosen incoming waves in two dimensions. Given the Gaussian prior measure, the well-posedness of the posterior measure in the Bayesian approach is discussed. The Markov Chain Monte Carlo (MCMC) method is adopted in the numerical approximation and the preconditioned Crank-Nicolson (pCN) algorithm with random proposal variance is utilized to improve the convergence rate. Numerical examples are provided to illustrate effectiveness of the proposed method.

Tables

Table 1. Table 1: Mean, standard deviation and relative error of ( x 1 , x 2 ) , r subscript 𝑥 1 subscript 𝑥 2 𝑟 (x_{1},x_{2}),r vs L 𝐿 L and M 𝑀 M .

$L$	$M$	mean	standard deviation	relative error
32	128	(0.9976, 0.2494), 0.1199	(0.0402, 0.0173), 0.0026	$(0.24\%,0.24\%),0.07\%$
32	32	(0.9192, 0.2804), 0.1203	(0.1893, 0.0847), 0.0100	$(8.08\%,12.15\%),0.28\%$
32	16	(0.9222, 0.2833), 0.1209	(0.1997, 0.0808), 0.0118	$(7.78\%,13.30\%),0.77\%$
64	128	(0.9941, 0.2514), 0.1199	(0.0278, 0.0144), 0.0016	$(0.59\%,0.56\%),0.10\%$
32	64	(0.9934, 0.2542), 0.1200	(0.0606, 0.0287), 0.0034	$(0.66\%,1.67\%),0.03\%$
16	32	(0.9023, 0.2794), 0.1201	(0.2207, 0.0841), 0.0112	$(9.77\%,11.74\%),0.11\%$
64	64	(0.9934, 0.2521), 0.1198	(0.0399, 0.0204), 0.0025	$(0.66\%,0.83\%),0.18\%$
16	64	(0.9290, 0.2660), 0.1196	(0.1857, 0.0763), 0.0104	$(7.10\%,6.41\%),0.36\%$
8	64	(0.9127, 0.2661), 0.1187	(0.1830, 0.0828), 0.0115	$(8.73\%,6.43\%),1.07\%$

Table 2. Table 2: Mean and relative error of ( x 1 , x 2 ) , r subscript 𝑥 1 subscript 𝑥 2 𝑟 (x_{1},x_{2}),r vs small L 𝐿 L and M 𝑀 M .

$L$	$M$	mean	relative error
4	2	(3.0494, 0.9982), 0.2643	(2.05%, 2.99%),  1.20%
4	8	(3.0457, 0.9949), 0.2717	(2.05%, 2.98%),  1.26%
8	4	(2.9069, 1.1143), 0.2732	(1.91%, 3.46%),  1.28%

Table 3. Table 3: The mean, standard deviation and relative error of ( x 1 , x 2 ) , r subscript 𝑥 1 subscript 𝑥 2 𝑟 (x_{1},x_{2}),r vs noise coefficient σ η subscript 𝜎 𝜂 \sigma_{\eta} .

${\sigma }_{\eta }$	mean	standard deviation	relative error
$3\%$	(0.9917, 0.2520), 0.1199	(0.0077, 0.0039), 0.0005	(0.83%, 0.79%),  0.10%
$6\%$	(0.9914, 0.2519), 0.1199	(0.0137, 0.0065), 0.0009	(0.86%, 0.77%),  0.10%
$9\%$	(0.9918, 0.2517), 0.1199	(0.0199, 0.0095), 0.0013	(0.82%, 0.69%),  0.10%

Table 4. Table 4: The mean, standard deviation and relative error of x 1 , x 2 , y 1 , y 2 subscript 𝑥 1 subscript 𝑥 2 subscript 𝑦 1 subscript 𝑦 2 x_{1},x_{2},y_{1},y_{2} vs noise coefficient σ η subscript 𝜎 𝜂 \sigma_{\eta} .

${\sigma }_{\eta }$	mean	standard deviation (${10}^{-3}$)	relative error ($\%$)
$3\%$	1.9994, 3.0055, 3.9992, 5.0029	3.5, 4.4, 3.4, 5.1	0.03, 0.18, 0.02, 0.06
$6\%$	1.9989, 3.0056, 4.0013, 5.0133	6.2, 7.7, 6.1, 8.8	0.06, 0.19, 0.03, 0.27
$9\%$	1.9936, 2.9994, 3.9968, 5.0081	8.6, 11.4, 8.9, 12.5	0.32, 0.02, 0.08, 0.16

Table 5. Table 5: Reconstruction and Hausdorff distance (HD) vs L ​ ( M = L ) 𝐿 𝑀 𝐿 L\ (M=L) .

$L$	Reconstruction	HD
5	-0.6441,  0.6287,  1.0013,  0.6460,  1.4934, -0.0051	3.6185
25	-0.6488, -3.0037,  1.0047,  0.6479,  1.5030, -0.0043	0.0021
50	-0.6499, -3.0001,  1.0000,  0.6503,  1.4998, -0.0009	$1.30\times {10}^{-4}$
100	-0.6497, -2.9995,  1.0000,  0.6508,  1.4999,  0.0006	$6.98\times {10}^{-5}$

Table 6. Table 6: Numerical solutions vs σ η subscript 𝜎 𝜂 \sigma_{\eta} with L = M = 50 𝐿 𝑀 50 L=M=50 .

${\sigma }_{\eta }$	mean of reconstructions	mean of HD	SD of HD
$3\%$	-0.6496, -2.9988, 0.9988, 0.6506, 1.4988, 0.0010	0.0039	0.0035
$6\%$	-0.6493, -2.9978, 0.9975, 0.6514, 1.4975, 0.0020	0.0062	0.0054
$9\%$	-0.6482, -2.9960, 0.9959, 0.6517, 1.4962, 0.0031	0.0088	0.0073

Equations

$$\Delta u^{in}+k^2{u}^{in}=0\text{in}{\mathbb{R}}^2,$$

$$\Delta u^{sc}+k^2{u}^{sc}=0\text{in}{\mathbb{R}}^2\backslash \bar{D},$$

$$u^{sc}=-u^{in}\text{on}\partial D,$$

$$\underset{r\to \infty }{\lim }\sqrt{r}\left(\frac{\partial u^{sc}}{\partial r}-ik{u}^{sc}\right)=0,r=|x|,$$

$$u^{sc}(x)=\frac{e^{\mathrm{i}k|x|}}{\sqrt{|x|}}\left\{u^{\infty }(\hat{\mathbf{x}})+\mathcal{O}\left(\frac{1}{\sqrt{|x|}}\right)\right\},|x|\to \infty ,$$

$$u^{\infty }(\hat{\mathbf{x}};D_z,\mathbf{d})=e^{ikz\cdot (\mathbf{d}-\hat{\mathbf{x}})}{u}^{\infty }(\hat{\mathbf{x}};D,\mathbf{d})\text{mbox}forall\hat{\mathbf{x}}\in \mathbb{S}.$$

$$|u^{\infty }(\hat{\mathbf{x}};D_z,\mathbf{d})|=|u^{\infty }(\hat{\mathbf{x}};D,\mathbf{d})|\text{mbox}forall\mathbf{d},\hat{\mathbf{x}}\in \mathbb{S}.$$

$$u_{\ell }^{in}(x):=e^{\mathrm{i}kx\cdot {\mathbf{d}}_0}+e^{\mathrm{i}kx\cdot {\mathbf{d}}_{\ell }},\ell =0,1,2,\cdots ,L,$$

$$\mathbf{Z}:=(z_1,z_2,\cdots ,z_N{)}^{\top }\in {\mathbb{R}}^N,N\in {\mathbb{N}}_0.$$

$$u^{\infty }(\hat{\mathbf{x}};\mathbf{Z},{\mathbf{d}}_0,{\mathbf{d}}_{\ell },k),\ell =1,2,\cdots ,L,\hat{\mathbf{x}}\in \mathbb{S}.$$

$$|u^{\infty }(\hat{\mathbf{x}};\mathbf{Z},{\mathbf{d}}_0,{\mathbf{d}}_{\ell },k)|,\ell =1,2,\cdots ,L,\hat{\mathbf{x}}\in \mathbb{S},$$

$${\pi }_{pr}(\mathbf{Z})$$

$${\mu }_{pr}(d\mathbf{Z})$$

$${\mu }_{pr}^n=\mathcal{N}(m_n,{\sigma }_n),n=1,2,\cdots ,N.$$

$$|u^{\infty }(\hat{\mathbf{x}})|=F(\mathbf{Z}),\hat{\mathbf{x}}\in \mathbb{S},$$

$$g_m(|u^{\infty }(\hat{\mathbf{x}})|)=|u^{\infty }({\hat{\mathbf{x}}}_m)|,m=1,2,\cdots ,M,$$

$$y_m=g_m(|u^{\infty }(\hat{\mathbf{x}})|)+{\eta }_m=|u^{\infty }({\hat{\mathbf{x}}}_m)|+{\eta }_m,$$

$$\mathbf{Y}=\mathcal{G}(\mathbf{Z})+\eta ,Y,\eta \in {\mathbb{R}}^M,\mathbf{Z}\in {\mathbb{R}}^N.$$

$$\Phi (\mathbf{Z};\mathbf{Y})=\frac{1}{2}|\mathbf{Y}-\mathcal{G}(\mathbf{Z}){|}_{{\Sigma }_{\eta }}^2,$$

$$\rho\big{(}\mathbf{Y}-\mathcal{G}(\mathbf{Z})\big{)}=\frac{1}{\big{(}(2\pi)^{M}\;\mbox{det}(\Sigma_{\eta})\big{)}^{1/2}}\;e^{-\Phi(\mathbf{Z};\mathbf{Y})}.$$

$$\frac{d{\mu }_{post}}{d{\mu }_{pr}}(\mathbf{Z})\propto e^{-\Phi (\mathbf{Z};\mathbf{Y})}.$$

$$|\mathcal{G}(\mathbf{Z}){|}_{{\Sigma }_{\eta }}\le e^{\epsilon \Vert \mathbf{Z}{\Vert }_2^2+\hat{M}},$$

$$|\mathcal{G}({\mathbf{Z}}_1)-\mathcal{G}({\mathbf{Z}}_2){|}_{{\Sigma }_{\eta }}\le K\Vert {\mathbf{Z}}_1-{\mathbf{Z}}_2{\Vert }_2.$$

$$d_{\text{Hell}}\left({\mu }_1,{\mu }_2\right):=\sqrt{\frac{1}{2}\int {\left(\sqrt{\frac{d{\mu }_1}{d{\mu }_0}}-\sqrt{\frac{d{\mu }_2}{d{\mu }_0}}\right)}^{2}d{\mu }_0},$$

$$d_{\text{Hell}}({\mu }_{post}^1,{\mu }_{post}^2)\le C\Vert {\mathbf{Y}}_1-{\mathbf{Y}}_2{\Vert }_2,$$

$$\mathbf{X}={\color[rgb]{0.000,0.000,0.000}\mathbf{m}_{pr}+}(1-\beta^{2})^{1/2}{\color[rgb]{0.000,0.000,0.000}(\mathbf{Z}-\mathbf{m}_{pr})}+\beta\omega,$$

$$\mathbf{X}={\color[rgb]{0.000,0.000,0.000}\mathbf{m}_{pr}+(1-\beta_{j}^{2})^{1/2}(\mathbf{Z}_{j}-\mathbf{m}_{pr})}+\beta_{j}\omega,\hskip 8.5359pt\omega\sim\mathcal{N}(\mathbf{0},\Sigma_{pcn});$$

$$\alpha\big{(}\mathbf{X},\mathbf{Z}_{j}\big{)}=min\{1,e^{\Phi(\mathbf{Z}_{j};\mathbf{Y})-\Phi(\mathbf{X};\mathbf{Y})}\};$$

$$\mathbf{Z}_{j+1}=\left\{\begin{array}[]{ll}\mathbf{X},&\hbox{if\; $U\leq\alpha\big{(}\mathbf{X},\mathbf{Z}_{j}\big{)}$,}\\ \mathbf{Z}_{j},&\hbox{if\, otherwise;}\end{array}\right.$$

$${\beta }_{new}=(1-{\gamma }^2{)}^{1/2}{\beta }_j+\gamma ({\omega }_{\beta }-0.5),{\omega }_{\beta }\sim \mathcal{U}(0,1),$$

$$\beta_{j+1}=\left\{\begin{array}[]{ll}\beta_{new},&\mbox{if}\quad\beta_{new}\in[0,1],\\ -\beta_{new},&\mbox{if}\quad\beta_{new}<0,\\ \beta_{new}-1,&\mbox{if}\quad\beta_{new}>1.\end{array}\right.$$

Full text

1 Introduction

Time-harmonic inverse scattering problems have attracted extensive attention due to their numerous applications in many areas such as radar and sonar detection, geophysical prospection, medical imaging, nondestructive testing and so on. In this paper, we are interested in the inverse problem of reconstructing the location and shape of an acoustically sound-soft obstacle using phaseless far-field data.

The propagation of a time-harmonic incident field $u^{in}$ in a homogeneous and isotropic medium is governed by the Helmholtz equation

[TABLE]

where $k>0$ is the wavenumber. Let $D\subset\mathbb{R}^{2}$ be a sound-soft scatterer, which occupies a bounded subset with $C^{2}$-smooth boundary $\partial D$ such that the exterior $\mathbb{R}^{2}\backslash\bar{D}$ of $D$ is connected. In this paper $D$ maybe a domain or a curve, which represents an extended obstacle or a crack in acoustics. The forward scattering problem is to find the scattered (perturbed) field $u^{sc}$ to the Helmholtz equation

[TABLE]

which satisfies the Dirichlet boundary condition

[TABLE]

and the Sommerfeld radiation condition

[TABLE]

uniformly in all directions $\hat{\mathbf{x}}=x/|x|\in\mathbb{S}:=\{x:|x|=1\}$, $x\in\mathbb{R}^{2}\backslash\bar{D}$, and $i=\sqrt{-1}$ is the imaginary unit. The total field $u$ is defined as $u=u^{in}+u^{sc}$ in $\mathbb{R}^{2}\backslash\bar{D}$. The Sommerfeld radiating solution $u^{sc}$ has an asymptotic behavior of the form

[TABLE]

where $u^{\infty}(\hat{\mathbf{x}})$ is called the far-field pattern at the observation direction $\hat{\mathbf{x}}\in\mathbb{S}$. Note that $u^{\infty}:\mathbb{S}\rightarrow\mathbb{C}$ is an analytic function with phase information. The above model also appears in the TE polarization of time-harmonic electromagnetic scattering from infinitely long and perfectly conducting cylinders.

The uniqueness, stability and inversion algorithms for recovering $\partial D$ from phased far-field patterns have been extensively studied with one or many incoming plane and point source waves. We refer to the monographs [6, 8, 24, 23, 37] for historical remarks, an overview of recent progresses and the comparison between different approaches. In many practical applications, the phase information of the far-field pattern cannot be measured accurately compared with its modulus or intensity. For instance, in optics it is not trivial to measure the phase of electromagnetic waves incited at high frequencies. One of the essential difficulties in using phaseless far-field data lies in the translation invariance property for plane wave incidence, which we state as follows. Let $u^{\infty}(\hat{\mathbf{x}};D,\mathbf{d})$ be the far-field pattern corresponding to the incident plane wave $e^{ikx\cdot\mathbf{d}}$ ($\mathbf{d}\in\mathbb{S}$ is the incident direction) and the sound-soft obstacle $D$. For the shifted obstacle $D_{z}:=\{x+z:x\in D\}$, the corresponding far-field pattern is given by (see [29])

[TABLE]

Hence, we get

[TABLE]

This implies that it is impossible to recover the location of $D$ from the phaseless far-field pattern of a plane wave.

There has been tremendous interest in inverse scattering with phaseless data or in phase retrieval problems in optics and other physical and engineering areas (see, e.g. [1, 2, 12, 17, 20, 25, 26, 34, 35] and the references therein). In a deterministic setting where randomness are not taken into account, Kress & Rundell and Ivanyshyn & Kress proposed a Newton-type iterative approach to reconstruct the shape of sound-soft obstacles from only the modulus of the far-field pattern in [17, 20, 29]. The approach of [17, 20] was based on a pair of nonlinear and ill-posed integral equations motivated by an inverse boundary value problem for the Laplace equation with phase information [18, 30]; see also [12, 19, 28] for inverse scattering from sound-soft cracks using a single far-field pattern with phase or phaseless information. Klibanov proved unique determination of a compactly supported potential of the stationary three-dimensional Schr$\ddot{\text{o}}$dinger equation from the phaseless near-field data incited by an interval of frequencies [25]. This was later extended in [26] to the reconstruction of a smooth wave speed in the three-dimensional Helmholtz equation. To broke the translation invariance property, it was recently prosed in [41] that, phaseless far-field patterns generated by infinitely many sets of superpositions of two plane waves with different directions can be used to uniquely determine a penetrable or impenetrable scatterer; see [42] for a fast imaging algorithm based on this idea. Similar uniqueness results were derived by Zhang & Guo [43] where the superposition of a fixed plane wave and some point sources was taken as incident waves and a reference ball technique was proposed. Uniqueness and direct sampling algorithms using the superposition of plane waves and fixed source location point sources were considered in [21]. In this paper we propose to generate the phaseless data using the following superposition of two plane waves (see [41])

[TABLE]

and then to recover a sound-soft disk, a line crack or a kite-shaped obstacle through the Bayesian approach. In (7), we fix $\mathbf{d}_{0}\in\mathbb{S}$ and change $\mathbf{d}_{\ell}\in\mathbb{S}$ as incident directions, due to the a priori information of the obstacle; see Theorem 4 for a uniqueness proof for sound-soft disks.

In recent years, the Bayesian method has received increasing attention for inverse problems [10, 16, 22, 31, 39], which also has been applied to the inverse scattering problems [3, 5, 14, 32, 33, 40] with phase far-field data. In particular, the authors of [5] adopt the Bayesian framework of [39] to shape identification problems in inverse scattering and establish a framework for proving well-posedness of the Bayesian formulation using a suitable shape parametrization and the regularity of shape derivatives. The aim of this paper is to propose the Bayesian method using a single far-field pattern without phase information. The Bayesian method provides us a new perspective to view the inverse scattering problem in the form of statistical inferences. In this statistical approach, all parameters are random variables and the key issue is to estimate the posterior distribution of the unknown quantities based on the Bayes’ formula [39] and the known prior distribution. The Bayesian method could be an alternative method to overcome the challenges in deterministic inverse problems, although it usually leads to expensive computational cost. The advantageous over deterministic inversion schemes in inverse scattering (for example, optimization-based iterative schemes and non-iterative sampling methods) are summarized as follows. (i) Instead of deterministic reconstructions, the Bayesian approach gives rise to statistic characteristics of the posterior distribution of unknown parameters and provides a quantification of the uncertainties arising from the corresponding model predictions. (ii) The Bayesian method could lead to all possible solutions of the inverse problem. For example, using one plan wave rather than the superposition of two plane waves, the Bayesian method could reconstruct the shape of an unknown obstacle which is located at every possible position due to the translation invariance (6). (iii) For inverse scattering problems, the theoretical analysis and numerical methods in the Bayesian framework are only based on the deterministic forward model. Hence, it is easy to perform theoretical analysis and numerical examples. (iv) The Bayesian method needs less measurement data without phase information and does not require a good initial guess. In this paper, the Markov chain Monte Carlo (MCMC) method [4, 11, 13] is proposed to accomplish the characterization of the posterior distribution, while the preconditioned Crank-Nicolson (pCN) algorithm [9] is adopted to improve the convergence rate in the iteration of MCMC method. Since the MCMC method and the pCN algorithm are adopt to calculate the numerical approximation, the numerical method in this paper is insensitive to the initial guess of the obstacle shape. In our numerical examples, we exhibit that accurate reconstructions can be achieved when the number of incident waves and observation directions is small.

This paper is organized as follows. In section 2, we adapt the Bayesian framework to inverse scattering problems with phaseless data. In section 3, we exhibit numerical results for recovering a disk, a line crack and a kite-shaped obstacle. Conclusions are given in section 4.

2 Bayesian Framework

In this paper we want to recover an unknown sound-soft obstacle from phaseless far-field patterns corresponding to a set of superposition of two plane waves. We propose the Bayesian approach to solve this inverse scattering problem. First of all, we set a suitable parameterization of the position and the shape of an obstacle. Then, in the Bayesian framework, we estimate the posterior distribution of the unknown obstacle parameters. By the Bayes’ theorem [31, 39], the posterior distribution of these parameters can be obtained from the prior distribution and the likelihood function to be specified in this section. Numerically, we will adopt the Markov chain Monte Carlo method (MCMC) to get an approximation of the posterior distribution.

2.1 Parameterization of the obstacle

Since the boundary of the underlying obstacle is a $C^{2}$-smooth curve, we can represent or approximate its geometrical shape by a finite set $\mathbf{Z}$ of variables

[TABLE]

For example, we can use four parameters $\mathbf{Z}:=(z_{1},z_{2},z_{3},z_{4})^{\top}$ to represent a line segment, where $(z_{1},z_{2})^{\top}$ and $(z_{3},z_{4})^{\top}$ denote respectively the two ending points, or we can use $\mathbf{Z}:=(a_{1}$, $b_{1}$, $a_{2}$, $b_{2}$, $\cdots$, $a_{N}$, $b_{N})^{\top}$ to approximate a star-shaped closed curve where $\{(a_{j},b_{j}):j=1,\cdots,N\}$ stand for the Fourier coefficients in the truncated Fourier expansion.

Recalling the incident waves $u^{in}_{\ell}(x),\ell=1,2,\cdots,L$ in the form of a set of superpositions of two plane waves (7), we express the (phased) far-field patterns of the scattering model (2)-(4) by

[TABLE]

Correspondingly, the phaseless far-field pattern are denoted by

[TABLE]

where $|\cdot|$ is the modulus of a complex number.

2.2 Prior distribution

By (8), the prior distribution of the obstacle parameters $\mathbf{Z}$ depends on the distribution of $z_{n},n=1,2,\cdots,N$. Let $\{z_{n}\}_{n=1}^{N}$ be independent variables with the prior density $\pi_{pr}^{n}$ and prior measure $\mu_{pr}^{n}$. Then the prior density $\pi_{pr}$ and prior measure $\mu_{pr}$ of $\mathbf{Z}$ are respectively given by

[TABLE]

In this paper we assume that $z_{n}$ are random variables with the Gaussian distribution, that is,

[TABLE]

For simplicity, we assume that $\sigma_{1}=\cdots=\sigma_{N}=\sigma_{pr}$, implying that $\mu_{pr}=\mathcal{N}(\mathbf{m}_{pr},\sigma_{pr}\mathbf{I})$, where $\mathbf{m}_{pr}=(m_{1},m_{2},\cdots,m_{N})^{\top}$ and $\mathbf{I}\in\mathbb{R}^{N\times N}$ is the identity matrix.

2.3 Observation of far-field pattern

To bridge the parameterization (8) of the obstacle and the associated phaseless far-field data (10), we define an operator F:\mathbb{R}^{N}\rightarrow{\color[rgb]{0.000,0.000,0.000}C^{\infty}(\mathbb{S})} as

[TABLE]

which can be regarded an abstract map from the space of obstacle parameters to the space of observation data in the continuous sense. From the well-posedness of forward scattering, $F$ is continuous but highly non-linear.

Let G=(g_{1},g_{2},\cdots,g_{M})^{\top}:{\color[rgb]{0.000,0.000,0.000}C^{\infty}(\mathbb{S})}\rightarrow\mathbb{R}^{M} be a bounded linear observation operator with g_{m}:{\color[rgb]{0.000,0.000,0.000}C^{\infty}(\mathbb{S})}\rightarrow\mathbb{R}_{+} given by

[TABLE]

where $\{\hat{\mathbf{x}}_{m}\}_{m=1}^{M}\subset\mathbb{S}$ is the set of discrete observation directions. Then the observation at the observation direction $\hat{\mathbf{x}}_{m}$ can be rephrased as

[TABLE]

where $\eta_{m}$ represents the noise polluting the observation data at the direction $\hat{\mathbf{x}}_{m}$.

Set $\mathbf{Y}=(y_{1},y_{2},\cdots,y_{M})^{\top}\in\mathbb{R}^{M}$ and $\eta=(\eta_{1},\eta_{2},\cdots,\eta_{M})^{\top}\in\mathbb{R}^{M}$. Denote by $\mathcal{G}=G\circ F$ the map from the obstacle parameter space $\mathbb{R}^{N}$ to observation space $\mathbb{R}^{M}$, that is,

[TABLE]

Our inverse problem in this paper is to determine the obstacle parameters $\mathbf{Z}\in\mathbb{R}^{N}$ from the observation data $\mathbf{Y}\in\mathbb{R}^{M}$ with the noise pollution $\eta\in\mathbb{R}^{M}$.

2.4 Likelihood

We assume the observation pollution $\eta$ is independent of $u^{\infty}$ and drawn from the Gaussian distribution $\mathcal{N}(\mathbf{0},\Sigma_{\eta})$ with the density $\rho$, where $\Sigma_{\eta}\in\mathbb{R}^{M\times M}$ is a self-adjoint positive matrix. By the observation of the phaseless data with noise (15), we can get the relationship $\mathbf{Y}|Z\sim\mathcal{N}(\mathcal{G}(\mathbf{Z}),\Sigma_{\eta})$. Define the model-data misfit function $\Phi(\mathbf{Z};\mathbf{Y}):\mathbb{R}^{N}\times\mathbb{R}^{M}\rightarrow\mathbb{R}$ as

[TABLE]

where $|\cdot|_{\Sigma_{\eta}}=|\Sigma_{\eta}^{-\frac{1}{2}}\cdot|$. Hence, the likelihood function is given by

[TABLE]

Furthermore, the posterior density $\pi_{post}$ and the posterior measure $\mu_{post}$ are connected to the prior measure $\mu_{pr}$ through the Radon-Nikodym derivative [38], given by

[TABLE]

2.5 Well-posedness of Bayesian framework

The well-posedness arguments of [5, 39] can be applied to deal with our inverse scattering problem with the Bayesian approach. In our phaseless case, we are required to justify the following Assumption 1, relying on regularity properties of the forward operator $\mathcal{G}$.

Assumption 1**.**

The map $\mathcal{G}:\mathbb{R}^{N}\rightarrow\mathbb{R}^{M}$ satisfies

• (i)

For every $\varepsilon>0$, there is an $\hat{M}=\hat{M}(\varepsilon)\in\mathbb{R}$ such that, for all $\mathbf{Z}\in\mathbb{R}^{N}$,

[TABLE]

where $\|\cdot\|_{2}$ is the Euclidean norm.

• (ii)

*For every $r>0$, there is a $K=K(r)>0$ such that, for all $\mathbf{Z}_{1},\mathbf{Z}_{2}\in\mathbb{R}^{N}$ with *

\max\big{\{}\|\mathbf{Z}_{1}\|_{2},\|\mathbf{Z}_{2}\|_{2}\big{\}}<r, it holds that

[TABLE]

We remark that there is no essential difference in proving Assumption 1 (i) between the phased and phaseless inverse scattering problems. The Assumption 1 (ii) follows directly from the triangle inequality \big{|}|a|-|b|\big{|}\leq|a-b| for complex numbers $a,b\in\mathbb{C}$ and the corresponding assumption for phased inverse scattering problems. Hence, when $D$ is sound-soft scatterer, Assumption 1 can be proved following the phased arguments of [5]; see also [32, 40] for the proofs in the case of limited aperture data and for interior scattering problems. If $D$ is sound-soft crack, the same results can be verified by applying the Fr$\acute{\text{e}}$chet differentiability with respect to the boundary of the far field operator ([27, 28]). The above assumptions together with the choice of the Gaussian prior measure (which satisfies $\mu_{pr}(\mathbb{R}^{N})=1$) lead to well-posedness of the Bayesian inverse problem, which is a result of application of Lemma 2.8, Theorem 4.1, Theorem 4.2 and Theorem 6.31 in [39]. Before stating the well-posedness (see Theorem 2 below), we recall the Hellinger distance defined by

[TABLE]

where $\mu_{1},\mu_{2}$ are two measures that are absolutely continuous with respect to $\mu_{0}$.

Theorem 2**.**

If the operator $\mathcal{G}$ satisfies the Assumption 1 and the prior measure $\mu_{pr}$ satisfies $\mu_{pr}(\mathbb{R}^{N})=1$, then the posterior measure $\mu_{post}$ is a well-defined probability measure on $\mathbb{R}^{N}$ and absolutely continuous with respect to prior measure $\mu_{pr}$. What’s more, the posterior measure $\mu_{post}$ is Lipschitz in the data $\mathbf{Y}$, with respect to the Hellinger distance: if $\mu_{post}^{1}$ and $\mu_{post}^{2}$ are two posterior measures corresponding to data $\mathbf{Y}_{1}$ and $\mathbf{Y}_{2}$, then there exists $C=C(r)>0$ such that,

[TABLE]

for all $\mathbf{Y}_{1}$, $\mathbf{Y}_{2}$ with \max\big{\{}\|\mathbf{Y}_{1}\|_{2},\|\mathbf{Y}_{2}\|_{2}\big{\}}<r.

2.6 Preconditioned Crank-Nicolson (pCN) algorithm with random proposal variance

This subsection is devoted to the numerical approximation of the posterior distribution. We adopt the Markov chain Monte Carlo method (MCMC) [4, 11, 13] to generate a large number of samples subject to the posterior distribution. The numerical approximation of the posterior distribution of unknown obstacle parameters can be obtained by statistic analysis on these samples. The Metropolis-Hastings [15, 36] algorithm will be used to construct MCMC samples. To improve the convergence rate of the MCMC method, we apply the preconditioned Crank-Nicolson algorithm [9].

According to the pCN algorithm, the new obstacle parameter $\mathbf{X}$ can be iteratively updated by the old parameter (initial guess) $\mathbf{Z}$ through the formula

[TABLE]

where $\beta\in[0,1]$ is the proposal variance coefficient and $\omega\sim\mathcal{N}(\mathbf{0},\Sigma_{pcn})$ is a zero-mean normal random vector with covariance matrix $\Sigma_{pcn}\in\mathbb{R}^{N\times N}$. We remark that it is important and very tricky to select a suitable $\beta$, because the value of $\beta$ dominates the proposal variance in the pCN algorithm. If $\beta\ll 1$ is small, the parameter $\mathbf{Z}$ will be updated slightly in the MCMC sequence, leading to a time-consuming iteration process to get the ergodic in the space of obstacle parameters. If $\beta$ is big, the parameter $\mathbf{Z}$ may stay at one state for quite a long time with a huge number of iterations in the MCMC method. Consequently, one cannot get enough number of samples to approximate the posterior distribution, due to the computational cost prohibition. To over come this difficulty, we recommend the pCN algorithm with a random proposal variance [9] to obtain good MCMC sequences (see below for the description).

Algorithm 3**.**

pCN Algorithm with Random Proposal Variance

• •

Initialize $\mathbf{Z}_{0}\in\mathbb{R}^{N}$ and $\beta_{0}\in[0,1]$.

• •

Repeat

Draw new obstacle parameter $\mathbf{X}$ from the old state $\mathbf{Z}_{j}$ by the pCN algorithm (19) with the proposal variance coefficient $\beta_{j}$ as:

[TABLE] 2. 2.

Compute Hasting ratio $\alpha(\cdot,\cdot):\mathbb{R}^{N}\times\mathbb{R}^{N}\rightarrow[1,\infty)$ as:

[TABLE] 3. 3.

Accept or reject $\mathbf{X}$: draw $U\sim\mathcal{U}(0,1)$ and then update $\mathbf{Z}_{j}$ by the criterion

[TABLE] 4. 4.

Generate new proposal variance coefficient $\beta_{j+1}$ from $\beta_{j}$. First we set

[TABLE]

with $\gamma\in[0,1]$. In our case we choose $\gamma=0.1$. Then $\beta_{j+1}$ can be updated by

[TABLE]

• •

Select $\mathbf{Z}_{\hat{j}},\ \hat{j}=J_{1}+(\tilde{j}-1)J_{2},\ \tilde{j}=1,2,\cdots,J_{3}$, $J_{1},J_{2},J_{3}\in\mathbb{N}_{0}$.**

In the Algorithm 3, the randomness of the proposal variance has the potential advantage of including the possibility of large and small proposal variance coefficients $\beta$. The large proposal variance coefficient $\beta$ helps the pCN algorithm explore the state space efficiently, and the small proposal variance coefficient $\beta$ protects the iterations from dropping into a fixed state. Then the random proposal variance gives rise to ergodic Markov chains.

The sequence $\mathbf{Z}_{\hat{j}}$ in the Algorithm 3 is selected to approximate the posterior distribution. Here, $J_{1}$ is the number of initial states; we take every $J_{2}$ sates to guarantee the selected sates are independent; $J_{3}$ is the number of the total selected states.

3 Numerical Examples

In this section we exhibit numerical examples to demonstrate the effectiveness of the method described in the previous section. To save computational costs, we consider three types of acoustically sound-soft scatterers in two dimensions:

• •

Sound-soft disks with unknown centers and radii;

• •

Line cracks with unknown starting and ending points;

• •

Kite-shaped obstacles with unknown position and shape.

There are totally three unknown parameters for disks, four parameters for cracks and six parameters for kite-shaped obstacle in 2D, implying that the unknown parameters always lie in a finite space with low dimensions.

With the definition of the observation (15), we construct two types of observations. Since the observation noise $\eta$ in (15) is assumed to be a Gaussian random vector, the zero vector is a special sample of the observation noise. Then in the first type observations, we consider an ideal model where the phaseless far-field data corresponding to the exact obstacle are polluted by this special sample ($\eta=\mathbf{0}$) of observation noise. In the second type observations, we consider a practical model with noise-polluted far-field data, which will be used to discuss the robustness of the numerical method in practical applications. The Hausdorf difference between the reconstructed and exact scatterers indicates the accuracy of our numerical method.

3.1 Disk

In this subsection, we assume the underlying obstacle is a sound-soft disk. The inverse problem of recovering disks arises from, for instance, the polarization model of time-harmonic electromagnetic scattering from perfectly conducting cylinders whose cross-section is a disk. The parameterization of a disk is given by

[TABLE]

where $(x_{1},x_{2})^{\top}$ is the center and $r$ is the radius of the disk. Since $r>0$, we assume that $r$ is a lognormal random variable, i.e., $z_{3}=\log{r}$ is a Gaussian random variable.

Let the incident wave be given by the sum of two plane waves of the form (7). If the disk is located at the origin, it is well-known that the corresponding far-field pattern incited by the plane wave $u^{in}(x)=e^{ikx\cdot\mathbf{d}_{\ell}}$ is given by the convergent series

[TABLE]

Here, $\theta_{\ell}=\angle(\hat{\mathbf{x}},\mathbf{d}_{\ell})$ denotes the angle between the observation direction $\hat{\mathbf{x}}$ and the incident direction $\mathbf{d}_{\ell}$, $J_{n}(\cdot)$ is the Bessel function of order $n$ and $H_{n}^{(1)}(\cdot)$ is the Hankel function of the first kind of order $n$. If the disk is located at $(x_{1},x_{2})^{\top}\in\mathbb{R}^{2}$, by the translational formula (6) and the linear superposition principle, the exact far-field pattern with phase information of the scattered waves can be expressed as

[TABLE]

Note that the first line on the right hand side of (27) denotes the far-field pattern corresponding to the incoming plane wave $e^{ikx\cdot\mathbf{d}_{0}}$, while the second line corresponding to $e^{ikx\cdot\mathbf{d}_{\ell}}$.

The following theorem states that our phaseless data set is sufficient to uniquely identify a sound-soft disk.

Theorem 4**.**

Let $k>0$ and $\mathbf{d}_{0}\in\mathbb{S}$ be fixed. Then the data $\{|u^{\infty}(\hat{\mathbf{x}};\mathbf{Z},\mathbf{d}_{0},\mathbf{d},k)|:\mathbf{d}\in\mathbb{S}\}$ uniquely determine a sound-soft disk (that is, the center and radius of a disk).

Proof.

Suppose that $D_{j}:=\{z=(z_{1},z_{2})^{\top}\in\mathbb{R}^{2}:\|z-x^{(j)}\|_{2}<r_{j}\}$ ($j=1,2$) are two sound-soft disks centered at $x^{(j)}=(x^{(j)}_{1},x^{(j)}_{2})^{\top}\in\mathbb{R}^{2}$ with the radius $r_{j}>0$ ($j=1,2$). Set $\mathbf{Z}^{(j)}=(x^{(j)}_{1},x^{(j)}_{2},r_{j})^{\top}$ and denote by $u^{\infty}(\hat{\mathbf{x}};\mathbf{Z}^{(j)},\mathbf{d}_{0},\mathbf{d},k)$ the far-field data corresponding to $D_{j}$ and the incident wave (7). Suppose that the phaseless far-field pattern are identical, i.e.,

[TABLE]

In particular, choosing $\mathbf{d}=\mathbf{d}_{0}$ in the previous relation yields

[TABLE]

for all $\hat{\mathbf{x}}\in\mathbb{S}$, where $u^{\infty}(\hat{\mathbf{x}};\mathbf{Z}^{(j)},\mathbf{d}_{0},k)$ stands for the far-field pattern corresponding to the plane wave $e^{ikx\cdot\mathbf{d}_{0}}$ incident onto $D_{j}$. This implies that

[TABLE]

Next, we shift the center of the disk $D_{j}$ to the origin and set $\mathbf{Z}^{(j)}_{0}:=(0,0,r_{j})^{\top}$. Recalling the translational formula (6), we obtain

[TABLE]

Since the shifted disks with the parameters $\mathbf{Z}_{0}^{(j)}$ are rotationally invariant, the far-field pattern $u^{\infty}(\hat{\mathbf{x}};\mathbf{Z}_{0}^{(j)},\mathbf{d},k)$ only depends on the angle between the incident direction $\mathbf{d}$ and the observation direction $\hat{\mathbf{x}}$. For any $\mathbf{d}\in\mathbb{S}$, there exist a orthogonal matrix $Q$ such that $\mathbf{d}=Q\mathbf{d}_{0}$. It then follows that (see e.g., [8, Chapter 5.1])

[TABLE]

Combining the previous two identities yields

[TABLE]

which together with the translational formula implies

[TABLE]

As a consequence of [41, Theorem 2.2], the relations (28) and (29) lead to the coincidence of $D_{1}$ and $D_{2}$, which proves Theorem 4. ∎

To apply the pCN algorithm (19) or (20), we need to set the key parameters $\beta$ and $\Sigma_{pcn}$. For sound-soft disks, we set $\beta_{j}=\beta=0.1$ and let $\Sigma_{pcn}=\mathbf{I}$ be the $N$-by-$N$ identity matrix. Then the proposal is given by

[TABLE]

We describe the settings of our computational performance as follows:

• •

Unless otherwise specified, the wave number is always taken as $k=1$;

• •

The incident directions are

[TABLE]

• •

The observation directions are

[TABLE]

• •

To compute the far-field pattern (27), we truncate the infinite series of (27) by using the Bessel and first-kind Hankel functions of order $n=0,1,2,\cdots,100$;

• •

In the Algorithm 3, we choose $J_{1}=9000$, $J_{2}=5$, $J_{3}=201$;

• •

In the setting of the prior distribution $\pi_{pr}$, we assume $\sigma_{pr}=1$;

• •

For the observation $\mathbf{Y}_{\ell}$ corresponding to the incident wave $u^{in}_{\ell}(x)$ in (7), $\ell=1,2,\cdots,L$, we assume the observation pollution $\eta^{\ell}=(\eta^{\ell}_{1},\eta^{\ell}_{2},\cdots,\eta^{\ell}_{M})^{\top}$ is a M-dimensional Gaussian variable, given by

[TABLE]

where $\omega^{\ell}_{m}\sim\mathcal{N}(0,1)$, $m=1,2,\cdots,M$ and $\sigma_{\eta}$ is the noise coefficient. In our numerical tests we choose $\sigma_{\eta}=3\%,6\%,9\%$. In other words, we take $\eta^{\ell}\sim\mathcal{N}(\mathbf{0},\mathbf{\Sigma_{\eta^{\ell}}})$ and the diagonal matrix $\mathbf{\Sigma_{\eta^{\ell}}}=\mbox{diag}(\sigma^{\ell}_{1},\sigma^{\ell}_{2},\cdots,\sigma^{\ell}_{M})$ with \sigma^{\ell}_{m}=\big{(}\sigma_{\eta}\times\left|u^{\infty}(\hat{\mathbf{x}}_{m};\mathbf{Z},\mathbf{d}_{0},\mathbf{d}_{\ell},k)\right|\big{)}^{2},m=1,2,\cdots,M;

• •

Unless otherwise specified, the accurate obstacle parameters are set as $\hat{\mathbf{Z}}=(\hat{x}_{1},\hat{x}_{2},\hat{r})^{\top}=(1,0.25,0.12)^{\top}$, that is, a disk centered at $(1,0.25)^{\top}$ with the radius $0.12$;

• •

We assume the initial guess is a disk centered at origin $(0,0)^{\top}$ with the radius $0.05$. Then the mean of the prior distribution is $\mathbf{m}_{pr}=(0,0,\log{0.05})^{\top}$.

In the first part, we adopt the ideal setting. Noting that in the formula (33) $\omega^{\ell}_{m}\sim\mathcal{N}(0,1)$, we can gain a special sample of the observation noise with $\omega^{\ell}_{m}=0$ and $\sigma_{\eta}=3\%$, $\ell=1,2,\cdots,L$, $m=1,2,\cdots,M$. This would help us to investigate the accuracy of the numerical method and to verify the above uniqueness result in Theorem 4.

At first, we discuss the accuracy of the numerical solutions for different choice of $L$ and $M$. Recall that the parameter $L$ denotes the number of incident waves and $M$ the number of observation directions. In Table 1, we show the mean, standard deviation and the relative error of the reconstructed parameters with different choice of $L,M$. The numerical solutions of the recovered centers and radii are shown in Figures 1 and 2. The histograms of numerical solutions of the centers and radii are shown in Figures 4, 3 and 5, respectively.

Based on these results, we find that the reconstructed parameters are getting more accurate as the number of incident or observation directions become larger. The numerical solutions with $(L,M)=(32,32)$, $(32,16)$, $(16,32)$, $(16,64)$, $(8,64)$ are relatively inaccurate, since the resulted relative errors are larger than $5\%$ in these cases. In contrast, the relative error is less $1\%$ if we choose $L$ and $M$ large enough such as $(L,M)=(64,128)$, $(64,64)$, $(32,128)$. On the other hand, it can be observed from Table 1 that the standard deviation decreases as $L$ or $M$ increases. Table 2 and Figure 6 illustrate that small $L$ and $M$ may lead to unreliable reconstructions.

In the following we suppose that the location of the center $(\hat{x}_{1},\hat{x}_{2})^{\top}=(2,2)^{\top}$ is known and the knowledge of the radius needs to be recovered. Since only one parameter of the scatter remains unknown, we make use of minimal number of incident and observation directions by setting $L=M=1$. In our tests we set incident directions $\mathbf{d}_{0}=(0,-1)^{\top},\mathbf{d}_{1}=(0,1)^{\top}$, observation direction $\hat{\mathbf{x}}_{1}=(0,-1)^{\top}$ and accurate radius $\hat{r}=1$. The numerical approximations of radius $r$ vs different wave numbers $k$ are exhibited in Figure 7. For each fixed $k$, we plot the phaseless far-field pattern $|u^{\infty}(\hat{\mathbf{x}}_{1})|$ against the radius $r$ in Figure 8.

From the numerical results we conclude that an accurate approximation of the radius can be obtained if the wave number $k$ is less than a threshold. It is seen from Figure 8 that the function $r\rightarrow|u^{\infty}(\hat{\mathbf{x}}_{1})|$ is monotonically increasing in $(0,R(k))$ where $R(k)\rightarrow 0^{+}$ as $k\rightarrow+\infty$. This suggests that for large $k$ such as $k=97,200,2000$, there are more than one radii corresponding to the measured phaseless far-field pattern at $\hat{\mathbf{x}}_{1}$. Hence, the reconstructed radii are inaccurate. These findings are consistent with the uniqueness result of [34], which states that a sound-soft disk can be uniquely determined from the phaseless far-field pattern at one observation direction, provided the radius is sufficiently small for a fixed wave number. The monotonicity property of the backscattered phaseless data with respect to the radius was rigorously justified in [34].

Having verified the accuracy of our inversion scheme at the ideal setting with a special sample of the observation noise, we now consider the inverse problem with a general sample of the observation noise at the noise coefficient $\sigma_{\eta}=3\%$. In the second part, we estimate the obstacle parameters by setting $k=1$, $L=32$ and $M=64$. We generate one sample of the observation noise, which is a matrix with $M\times L$ elements constructed by the formula (33). The numerical approximations from the polluted observation data are exhibited in Figures 9 and 10. The mean of the numerical center and radius are $(1.0089,0.2527)^{\top}$ and $0.1211$, respectively. The standard deviations of these parameters are $(0.0601,0.0259)^{\top},0.0039$ and the relative errors are $(0.89\%,1.08\%)^{\top},0.92\%$.

To demonstrate the robustness of the numerical scheme, we generate 1000 samples of the observation noise. For each sample of the observation noise, one can gain a corresponding reconstruction of the parameters $x_{1},x_{2},r$. Hence, we can perform statistical analysis over totally 1000 reconstructions of $x_{1},x_{2},r$. In our tests, we pollute the phaseless data at different levels $\sigma_{\eta}=3\%,6\%,9\%$ and exhibit the numerics in Table 3, Figure 11 and Figure 12. From these reconstructed parameters we conclude that the mean and relative error are robust against the noise pollution, but the standard deviation is very sensitive to the noisy level. Further, the phaseless data with less noise give rise to a more reliable reconstruction result.

3.2 Line cracks

A crack or an open arc can be used to model the defects inside elastic and solid bodies such as bridge structures, aircraft engines and wings etc. Detection of such scatterers is important in safety and health assessment and is one of the fundamental topics in ultrasonic non-destructive testing. In this subsection, we want to recover a sound-soft crack of line-segment-type with the starting point at $x=(x_{1},x_{2})^{\top}\in\mathbb{R}^{2}$ and the ending point $y=(y_{1},y_{2})^{\top}\in\mathbb{R}^{2}$. Hence, such line cracks can be characterized by $N=4$ parameters:

[TABLE]

Unlike the scattering from disks, we do not have an analytical expression of the far-field pattern corresponding to a line crack. Below we describe the integral equation method to solve the forward scattering problem, following the numerical scheme of [28] for general cracks. Denote by $\Gamma\subset\mathbb{R}^{2}$ an open arc of class $C^{3}$ in 2D, which can be parameterized as

[TABLE]

Using the integral equation method, the solution $u^{sc}$ to the Helmholtz equation (2) in $\mathbb{R}^{2}\backslash\Gamma$ can be expressed as a single-layer potential ([8])

[TABLE]

where $\Phi(x,y)$ is the fundamental solution to the Helmholtz equation in two dimensions given by

[TABLE]

Due to the Dirichlet boundary condition (3) on $\Gamma$, the unknown density function $\varphi$ is sought as a solution to the integral equation

[TABLE]

Once $\varphi$ is calculated from (38), the far-field pattern could be expressed in the form

[TABLE]

To describe the numerical scheme of [28], we first introduce two functions defined on $\mathbb{R}\times\mathbb{R}$ as follows

[TABLE]

and

[TABLE]

Here, $C\approx 0.577216$ is the Euler’s constant. Then the integral equation (38) can be rephrased as

[TABLE]

for $t\in[0,2\pi]$. Here,

[TABLE]

The quadrature method [7] can be employed to discretize the integral equation (42), based on the trigonometric interpolation with $2n$ equidistant nodal points $t_{j}:=\frac{j\pi}{n}$, $j=0,1,\cdots,2n-1$. Then the unknown solution $\psi$ to the integral equation (42) can be approximated by the $2n$ discrete nodal values $\left\{\psi_{j}=\psi(t_{j})\right\}_{j=0}^{2n-1}$. Since $\psi_{k}=\psi_{2n-k},k=1,2,\ldots,n-1$ with $\psi_{j}=|\sin(t_{j})|\left|z^{\prime}(\cos(t_{j}))\right|\varphi(z(\cos(t_{j})))$, it suffices to compute the $n+1$ discrete nodal values $\left\{\psi_{j}\right\}_{j=0}^{n}$ from the following $(n+1)\times(n+1)$ algebraic system

[TABLE]

with

[TABLE]

Note that there are totally $n+1$ unknown discrete nodal values $\left\{\psi_{j}\right\}_{j=0}^{n}$ in (48), because $\psi_{j}=\psi_{n-|n-j|}$ for all $j=0,1,\cdots,2n-1$. Now the far-field pattern can be approximated by

[TABLE]

If the right hand side of (48) (or (38)) is given by the incident wave (7), we obtain the far-field pattern $u^{\infty}(\hat{\mathbf{x}};\mathbf{Z},\mathbf{d}_{0},\mathbf{d}_{\ell},k)$ where $\mathbf{Z}$ denotes the crack parameter (34).

To set the parameters $\beta$ and $\Sigma_{pcn}$, we let $\Sigma_{pcn}=\mathbf{I}$ be the identity matrix, which is the same as the case of sound-soft disks. However, in this section the proposal variance coefficient $\beta$ is not a fixed number, but is taken as a random variable. This suggests that a random proposal variance is adopted to reconstruct line cracks. Then the proposal takes the form

[TABLE]

and the proposal variance coefficients $\beta_{j}$ need to be updated by the formulas (23) and (24). It should be noted that, the MCMC method with a fixed proposal variance coefficient converges slowly or even does not converge after a large number of iterations, which is in contrast to the efficient MCMC method for recovering disks. This could partly be due to the number of reconstructed parameters, which is four in the line crack case while three for a disk. The trace of the iterations of MCMC (shown in Figure 15) verifies the efficiency of the random proposal variance. In the first 10000 iteration steps, the trace converges fast but always drops into some fixed states, when the proposal variance coefficients $\beta_{j}$ are not appropriately updated. Numerics show that the trace can converge to and oscillate around the accurate state only after a large number of iterations.

As in the previous subsection, we set some computational parameters as follows:

• •

The wave number $k$, the incident directions $\mathbf{d}_{\ell},\ell=0,1,\cdots,L$, the observation directions $\hat{\mathbf{x}}_{m},m=1,2,\cdots,M$, $\sigma_{pr}$ of the prior distribution and the observation pollution $\eta^{\ell},\ell=1,\cdots,L$ are given as same as those for recovering disks;

• •

We choose $L=M=40$;

• •

In the Algorithm 3, we choose $J_{1}=18000$, $J_{2}=5$, $J_{3}=401$;

• •

The accurate obstacle is $\hat{\mathbf{Z}}=(\hat{x}_{1},\hat{x}_{2},\hat{y}_{1},\hat{y}_{2})^{\top}=(2,3,4,5)^{\top}$, that is, a line segment with the starting point $(2,3)^{\top}$ and the ending point $(4,5)^{\top}$;

• •

We assume the initial guess is a line segment with the starting point $(0,0)^{\top}$ and the ending point $(1,1)^{\top}$. Then the mean of the prior distribution is $\mathbf{m}_{pr}=(0,0,1,1)^{\top}$.

Unfortunately we do not have the uniqueness result analogous to Theorem 4 for recovering cracks. A local uniqueness result for general cracks was proved in [27] using a single far-field pattern with information. In the idea setting ($\omega^{\ell}_{m}=0$ and $\sigma_{\eta}=3\%$, $\ell=1,2,\cdots,L$, $m=1,2,\cdots,M$), the numerical approximations of the crack parameters are exhibited in Figure 13. The mean solutions of the starting and ending points are $(2.0012,3.0005)$, $(3.9990,4.9992)$, the standard deviations are $(0.0139,0.0164)$, $(0.0129,0.0170)$ and the relative errors are $(0.06\%,0.02\%)$, $(0.02\%,0.02\%)$.

Setting the noise coefficient $\sigma_{\eta}=3\%$, we generate one general sample of the observation noise by the formula (33), which takes the form of an $M\times L$ matrix. The numerical approximations from the polluted observations is exhibited in Figure 14. The mean solutions of the numerical starting and ending points are $(2.0052,3.0031)$, $(3.9958,4.9986)$, the standard deviations are $(0.0135,0.0155)$, $(0.0130,0.0184)$ and the relative errors are $(0.26\%,0.10\%)$, $(0.10\%,0.03\%)$. Besides, the trace of the iterations in MCMC are shown in Figure 15.

As done for recovering disks, we also demonstrate the robustness of the numerical scheme with 1000 samples of the observation noise at different levels $\sigma_{\eta}=3\%,6\%,9\%$. For each sample of the observation noise, one can gain a corresponding reconstruction of the parameters $x_{1},x_{2},y_{1},y_{2}$. Hence, we can perform statistical analysis over totally 1000 reconstructions of $x_{1},x_{2},y_{1},y_{2}$. The corresponding results are exhibited in Table 4, Figure 16 and Figure 17. From these reconstructed parameters, we can draw almost the same conclusions as those for determining a sound-soft disk. The mean and relative error are robust against the noise pollution, but the standard deviation is very sensitive to the noisy level. It follows that the phaseless data with less noise give rise to a more reliable reconstruction result.

3.3 Kite-shaped obstacle

In this subsection, we consider the following kite-shaped sound-soft obstacle $D_{1}$,

[TABLE]

which is a benchmark acoustically impenetrable scatterer in inverse scattering problems. Obviously, the boundary of $D_{1}$ can be parameterized by six parameters

[TABLE]

Suppose that the exact parameters are given by $\hat{\mathbf{Z}}$ $=(\hat{z}_{1}$, $\hat{z}_{2}$, $\cdots$, $\hat{z}_{6})^{\top}=(-0.65,-3,1,0.65,1.5,0)^{\top}$. To calculate numerical solutions of the forward problem, we adopt the MATLAB code given by [37, Chapter 8]. Since there are six unknown parameters, the random proposal variance is adopted as same as in recovering line cracks. As in the previous subsections, we set the computational settings as follows.

• •

The wave number is $k=2$;

• •

The directions of plane incident waves are

[TABLE]

• •

The observation directions are

[TABLE]

• •

Unless otherwise specified, we choose $L=M=50$;

• •

In the Algorithm 3, we choose $J_{1}=180000$, $J_{2}=20$, $J_{3}=1001$;

• •

Since the initial guess is assumed to be a unit circle centered at the origin, the mean of the prior distribution $\pi_{pr}$ is $\mathbf{m}_{pr}=(0,1,0,0,1,0)^{\top}$. In this example, we also assume $\sigma_{pr}=1$;

• •

The observation pollution $\eta^{\ell},\ell=1,\cdots,L$, are given as same as those for recovering disks.

Since the sixth exact obstacle parameter $\hat{z}_{6}$ is [math], we can not use the relative error to evaluate the accuracy of numerical results in this example. Instead, the Hausdorff distance (HD) is chosen to compute the distance between reconstructed and exact boundaries. Recall that the Hausdorff distance between two obstacles $\partial D_{2}$ and $\partial D_{3}$ is defined by

[TABLE]

In the first part, we consider the ideal setting ($\omega^{\ell}_{m}=0$ and $\sigma_{\eta}=3\%$, $\ell=1,2,\cdots,L$, $m=1,2,\cdots,M$) of observations to investigate the accuracy of our numerical method. As done for recovering disks, we discuss the accuracy of numerical solutions for different choice of $L\ (M=L)$. In Table 5 and Figure 18, we exhibit the numerical reconstructions and the Hausdorff distances (HD) between the numerical reconstructions and the exact boundaries. We find that the reconstructed parameters are getting more accurate as the number of incident waves and observation directions becomes larger. The Hausdorff distance is less than $0.001$ if we choose $L=M=50,100$. The numerical solution with $L=M=5$ is unreliable as illustrated in the Figure 18.

To show that this method is not sensitive to initial guess, we exhibit the initial guess, the exact boundary and numerical reconstruction in Figure 19, where we set $L=M=50$. We can obtain an accurate numerical solution with the Hausdorff distance (HD) being $1.30\times 10^{-4}$, even if the initial guess of the obstacle is separated from the exact one.

We draw the histogram of the selected 1001 sates with $L=M=50$ in the Figure 20, which are used to construct the posterior density. From Figure 20, we conclude that the posterior density is a Gaussian distribution and the mean of the posterior density approaches the exact obstacle parameters.

Then we consider the practical setting with noise-polluted phaseless far-field data. As done for recovering disks and line cracks, we also demonstrate robustness of the numerical scheme with 1000 samples of the observation noise at different noise levels $\sigma_{\eta}=3\%,6\%,9\%$ with $L=M=50$. For every sample of the observation noise, one can gain a corresponding reconstruction of the parameters $z_{1},z_{2},\cdots,z_{6}$, and then gain the corresponding Hausdorff distance (HD) between the numerical reconstruction and the exact boundary. Then we can use standard statistical tools to analyze these 1000 reconstructions to discuss robustness of our numerical scheme. In Table 6, we exhibit the mean of reconstructions, the mean and standard deviation (SD) of Hausdorff distances (HD) at different noise levels. In Figure 21, we describe the 1000 reconstructions at different noise levels $\sigma_{\eta}$. In Figure 22, we show the histogram of these 1000 reconstructions, which correspond to the 1000 samples of the observation noise. We can find that the mean and standard deviation (SD) of Hausdorff distances (HD) become larger as the noise level $\sigma_{\eta}$ is getting bigger. However, our inversion scheme is still robust against the noise pollution, since the mean and standard deviation (SD) of Hausdorff distances (HD) are very small as shown in Table 6. Further, the phaseless data with less noise give rise to a more reliable reconstruction result.

4 Conclusion

In this paper, we propose the Bayesian approach to inverse acoustic scattering from sound-soft disks, line cracks and kite-shaped obstacles with phaseless far-field data. Motivated by [41], the incoming waves are properly chosen in order to break the translational invariance of the far-field patten. Uniqueness of the inverse solution is proven for recovering a disk. When the Gaussian prior measure is given, we discuss well posedness of the posterior measure based on regularity properties of the deterministic direct scattering problem. Our numerics verify the efficiency of the preconditioned Crank-Nicolson algorithm with the random proposal variance. Further, increasing the number of incident and observation directions would lead to more accurate and reliable reconstructions. It is shown that the Bayesian method is robust for phaseless inverse scattering problems with respect to the observation noise. In this paper the obstacle boundary can be easily parameterized in a finite dimensional space. Our future efforts will be devoted to recovering the shape and physical properties of more general acoustic obstacles with a large number of unknown parameters from phaseless far-field patterns. Noe that in this paper the number of unknown obstacle parameters is not larger than six, which has reduced the computational cost. For more complex scatterers, the increased computational cost of the Markov chain Monte Carlo method needs to be improved, for example, by combining the Gibbs sampling method and the stochastic surrogate model of the forward solver. Besides the idea of using superposition of two plane wave, one can also make use of a single spherical incident wave within the Bayesian framework. Research outcomes along these directions will be reported in our forthcoming publications.