Inverse Random Source Problem for the Helmholtz Equation from Statistical Phaseless Data

TL;DR

This paper develops a novel method to reconstruct and analyze a random source for the Helmholtz equation using statistical phaseless data, overcoming non-uniqueness with phase retrieval and Bayesian techniques.

Contribution

It introduces a reference source technique, phase retrieval formulas, and Bayesian reconstruction for the inverse random source problem, providing theoretical analysis and numerical validation.

Findings

Successfully reconstructs the expectation and variance of the random source.

Proves stability and uniqueness of the inverse problem solutions.

Demonstrates effectiveness through numerical experiments.

Abstract

This paper investigates the problem of reconstructing a random source from statistical phaseless data for the two-dimensional Helmholtz equation. The major challenge of this problem is non-uniqueness, which we overcome through a reference source technique. Firstly, we introduce some artificially added point sources into the inverse random source system and derive phase retrieval (PR) formulas for the expectation and variance of the radiated fields. This paper rigorously analyze the uniqueness and stability of the recovered statistics of the radiated fields. Afterwards, since the direct problem has a unique mild solution, by examining the expectation and variance of this solution and combined with the phase retrieval formulas, we derive the Fredholm integral equations to solve the inverse random source problem (IRSP). We prove the stability of the corresponding integral equations. To quantify the uncertainty of the random source, we utilize the Bayesian method to reconstruct the random source and establish the well-posedness of the posterior distribution. Finally, numerical experiments demonstrate the effectiveness of the proposed method and validate the theoretical results.

Tables

Table 1. Table 1: The relative L 2 L^{2} errors between 𝐄 ​ ( u ​ ( ⋅ , k ) ) | ∂ B R \mathbf{E}(u(\cdot,k))|_{\partial B_{R}} and 𝐄 ϵ ​ ( u ​ ( ⋅ , k ) ) | ∂ B R \mathbf{E}^{\epsilon}(u(\cdot,k))|_{\partial B_{R}} for different wavenumbers and noise levels

	$k^{\ast }=\frac{\pi }{30}$	$k=\pi$	$k=60.5\pi$	$k=84.5\pi$
$ϵ=0.01$	0.2034	0.1272	0.6575	0.8143
$ϵ=0.005$	0.1014	0.0686	0.3585	0.4411
$ϵ=0.001$	0.0247	0.0313	0.1872	0.2287

Table 2. Table 2: The relative L 2 L^{2} errors between 𝐕𝐚𝐫 ​ ( u ​ ( ⋅ , k ) ) | ∂ B R \mathbf{Var}(u(\cdot,k))|_{\partial B_{R}} and 𝐕𝐚𝐫 ϵ ​ ( u ​ ( ⋅ , k ) ) | ∂ B R \mathbf{Var}^{\epsilon}(u(\cdot,k))|_{\partial B_{R}} for different wavenumbers and noise levels

	$k^{\ast }=\frac{\pi }{30}$	$k=\pi$	$k=60.5\pi$	$k=84.5\pi$
$ϵ=0.01$	0.0265	0.0122	0.0140	0.0141
$ϵ=0.005$	0.0135	0.0092	0.0113	0.0113
$ϵ=0.001$	0.0028	0.0082	0.0105	0.0107

Table 3. Table 3: The relative L 2 L^{2} errors between 𝐂𝐨𝐯 ​ ( Re ⁡ ( u j , k ) , Im ⁡ ( u j , k ) ) | ∂ B R \mathbf{Cov}(\operatorname{Re}(u_{j,k}),\operatorname{Im}(u_{j,k}))|_{\partial B_{R}} and 𝐂𝐨𝐯 ϵ ​ ( Re ⁡ ( u j , k ) , Im ⁡ ( u j , k ) ) | ∂ B R \mathbf{Cov}^{\epsilon}(\operatorname{Re}(u_{j,k}),\operatorname{Im}(u_{j,k}))|_{\partial B_{R}} for different wavenumbers and noise levels

	$k^{\ast }=\frac{\pi }{30}$	$k=\pi$	$k=60.5\pi$	$k=84.5\pi$
$ϵ=0.01$	0.0484	0.0620	0.3256	0.2438
$ϵ=0.005$	0.0247	0.0482	0.2352	0.1696
$ϵ=0.001$	0.0049	0.0429	0.1870	0.1363

Table 4. Table 4: The relative L 2 L^{2} errors with various noise levels ϵ \epsilon .

		$ϵ=5\%$	$ϵ=1\%$	$ϵ=0.1\%$
$g_{1_{err}}$	$f_{1,k}^{ϵ}$	0.5954	0.1499	0.1211
	$f_{2,k}^{ϵ}$	0.5773	0.1475	0.1151
${\sigma }_{1_{err}}$	$F_{k,1,1}^{ϵ}$	0.3950	0.1102	0.1030
	$F_{k,2,2}^{ϵ}$	0.3541	0.1183	0.1077
	$F_{k,1,2}^{ϵ}$	0.1761	0.0832	0.0830

Equations

$$\{\begin{array}{rl}\Delta u+k^{2}u& =f\text{in}{\mathbb{R}}^2,\\ \frac{\partial u}{\partial r}-\mathrm{i}ku& =o(r^{-1/2}),r=|x|\to \infty ,\end{array}$$

$$f(x)=g(x)+\sigma (x){\dot{W}}_x,$$

$$u(x,k)={\int }_{{\Omega }_0}{G}_k(x,y)g(y)\mathrm{d}y+{\int }_{{\Omega }_0}{G}_k(x,y)\sigma (y)\mathrm{d}{W}_y,$$

$$G_k(x,y):=-\frac{\mathrm{i}}{4}{H}_0^{(1)}(k|x-y|),$$

$${\Omega }_0=(-a,a{)}^2,a>0.$$

$${\mathbb{K}}_N:=\left\{\frac{\pi }{a}x,x⩾1,x\in \mathbb{R}\right\}\bigcup \left\{k^{\ast }\right\}.$$

$$B_j=\left\{r(\cos \theta ,\sin \theta ):0⩽r⩽R,{\nu }_{2j-2}⩽\theta ⩽{\nu }_{2j}\right\},j=1,\cdots ,m,$$

$$\partial B_j=\left\{R(\cos \theta ,\sin \theta ):{\nu }_{2j-2}⩽\theta ⩽{\nu }_{2j}\right\},j=1,\cdots ,m.$$

$$B_R=\bigcup _{j=1}^m{B}_j\text{ and }\partial B_R=\bigcup _{j=1}^m\partial B_j.$$

$$\Delta G_{j,\ell }+k^2{G}_{j,\ell }={\delta }_{j,\ell }\mathrm{in}{\mathbb{R}}^2.$$

$$\mathrm{Re}\left(G_{j,\ell }(x,k)\right)=\frac{1}{4}{Y}_0\left(k{r}_{j,\ell }\right),\mathrm{Im}\left(G_{j,\ell }(x,k)\right)=-\frac{1}{4}{J}_0\left(k{r}_{j,\ell }\right),$$

$$\{\begin{array}{ll}\Delta v_{j,\ell }+k^2{v}_{j,\ell }=f+{\delta }_{j,\ell },& \text{in }{\mathbb{R}}^2,\\ \frac{\partial v_{j,\ell }}{\partial r}-\text{i}k{v}_{j,\ell }=o\left(r^{-1/2}\right),& r=|x|\to \infty .\end{array}$$

$$u_{j,k}(\cdot ):={u(\cdot ,k)|}_{\partial B_j,}{v}_{j,\ell ,k}(\cdot ):={v_{j,\ell }(\cdot ,k)|}_{\partial B_j,}\ell =1,2,$$

$$\mathrm{Re}\left(v_{j,\ell ,k}\right)=\mathrm{Re}(u_{j,k})+\frac{1}{4}{Y}_0\left(k{r}_{j,\ell }\right),\mathrm{Im}\left(v_{j,\ell ,k}\right)=\mathrm{Im}(u_{j,k})-\frac{1}{4}{J}_0\left(k{r}_{j,\ell }\right).$$

$$|v_{j,\ell ,k}{|}^2$$

$$=(\mathrm{Re}(u_{j,k}){)}^2+(\mathrm{Im}(u_{j,k}){)}^2+\frac{1}{2}{Y}_0(k{r}_{j,\ell })\cdot \mathrm{Re}(u_{j,k})-\frac{1}{2}{J}_0(k{r}_{j,\ell })\cdot \mathrm{Im}(u_{j,k})$$

$$+\frac{1}{16}{Y}_0^2(k{r}_{j,\ell })+\frac{1}{16}{J}_0^2(k{r}_{j,\ell })$$

$$=|u_{j,k}{|}^2+\frac{1}{16}{\left|H_0^{(1)}(k{r}_{j,\ell })\right|}^2+\frac{1}{2}{Y}_0(k{r}_{j,\ell })\cdot \mathrm{Re}(u_{j,k})-\frac{1}{2}{J}_0(k{r}_{j,\ell })\cdot \mathrm{Im}(u_{j,k}).$$

$$\mathbf{E}(|v_{j,\ell ,k}{|}^2)=\mathbf{E}(|u_{j,k}{|}^2)+\frac{1}{16}{\left|H_0^{(1)}(k{r}_{j,\ell })\right|}^2+\frac{1}{2}{Y}_0(k{r}_{j,\ell })\cdot \mathbf{E}(\mathrm{Re}(u_{j,k}))-\frac{1}{2}{J}_0(k{r}_{j,\ell })\cdot \mathbf{E}(\mathrm{Im}(u_{j,k})).$$

$$Y_0(k{r}_{j,1})\cdot \mathbf{E}(\mathrm{Re}(u_{j,k}))-J_0(k{r}_{j,1})\cdot \mathbf{E}(\mathrm{Im}(u_{j,k}))=f_{j,1,k},$$

$$Y_0(k{r}_{j,2})\cdot \mathbf{E}(\mathrm{Re}(u_{j,k}))-J_0(k{r}_{j,2})\cdot \mathbf{E}(\mathrm{Im}(u_{j,k}))=f_{j,2,k},$$

$$f_{j,\ell ,k}=2\left(\mathbf{E}(|v_{j,\ell ,k}{|}^2)-\mathbf{E}(|u_{j,k}{|}^2)\right)-\frac{1}{8}{\left|H_0^{(1)}(k{r}_{j,\ell })\right|}^2,\ell =1,2.$$

$$\mathbf{E}(\mathrm{Re}(u_{j,k}))=\frac{\det (A_{j,k}^R)}{\det (A_{j,k})},\mathbf{E}(\mathrm{Im}(u_{j,k}))=\frac{\det (A_{j,k}^I)}{\det (A_{j,k})},$$

$$A_{j,k}=\left(\begin{array}{cc}{Y}_0(k{r}_{j,1})& -J_0(k{r}_{j,1})\\ Y_0(k{r}_{j,2})& -J_0(k{r}_{j,2})\end{array}\right),A_{j,k}^{\mathrm{R}}=\left(\begin{array}{cc}{f}_{j,1,k}& -J_0(k{r}_{j,1})\\ f_{j,2,k}& -J_0(k{r}_{j,2})\end{array}\right),$$

$$A_{j,k}^{\mathrm{I}}=\left(\begin{array}{cc}{Y}_0(k{r}_{j,1})& f_{j,1,k}\\ Y_0(k{r}_{j,2})& f_{j,2,k}\end{array}\right).$$

$$\mathbf{Cov}(|u_{j,k}{|}^2,|v_{j,\ell ,k}{|}^2)$$

$$=$$

$$4\mathbf{Var}(|v_{j,\ell ,k}{|}^2)=$$

$$=$$

$$+4Y_0(k{r}_{j,\ell })\cdot \mathbf{Cov}(|u_{j,k}{|}^2,\mathrm{Re}(u_{j,k}))-4J_0(k{r}_{j,\ell })\cdot \mathbf{Cov}(|u_{j,k}{|}^2,\mathrm{Im}(u_{j,k}))$$

$$-2Y_0(k{r}_{j,\ell })\cdot J_0(k{r}_{j,\ell })\cdot \mathbf{Cov}(\mathrm{Re}(u_{j,k}),\mathrm{Im}(u_{j,k}))$$

$$=$$

$$-4\mathbf{Var}(|u_{j,k}{|}^2)-2Y_0(k{r}_{j,\ell })\cdot J_0(k{r}_{j,\ell })\cdot \mathbf{Cov}(\mathrm{Re}(u_{j,k}),\mathrm{Im}(u_{j,k}))$$

$$4\mathbf{Cov}(|v_{j,1,k}{|}^2,|v_{j,2,k}{|}^2)=$$

$$|u_{j,k}{|}^2+\frac{1}{2}{Y}_0(k{r}_{j,2})\cdot \mathrm{Re}(u_{j,k})-\frac{1}{2}{J}_0(k{r}_{j,2})\cdot \mathrm{Im}(u_{j,k}))$$

$$=$$

$$+Y_0(k{r}_{j,1})Y_0(k{r}_{j,2})\cdot \mathbf{Var}(\mathrm{Re}(u_{j,k}))+J_0(k{r}_{j,1})J_0(k{r}_{j,2})\cdot \mathbf{Var}(\mathrm{Im}(u_{j,k}))$$

Full text

Inverse Random Source Problem for the Helmholtz Equation from Statistical Phaseless Data

Qiao-Ping Chen School of Mathematics, Hunan University, Changsha 410082, China. Email: s230600375@hnu.edu.cn   

Hongyu Liu Department of Mathematics, City University of Hong Kong, Kowloon, Hong Kong, China.  Email: hongyu.liuip@gmail.com; hongyliu@cityu.edu.hk   

Zejun Sun School of Mathematics, Hunan University, Changsha 410082, China. Email: sunzejun@hnu.edu.cn   

Li-Li Wang School of Mathematics, Hunan University, Changsha 410082, China. Email: lilywang@hnu.edu.cn   

Guang-Hui Zheng School of Mathematics, Hunan University, Changsha 410082, China. Email: zhenggh2012@hnu.edu.cn; zhgh1980@163.com

ABSTRACT

This paper investigates the problem of reconstructing a random source from statistical phaseless data for the two-dimensional Helmholtz equation. The major challenge of this problem is non-uniqueness, which we overcome through a reference source technique. Firstly, we introduce some artificially added point sources into the inverse random source system and derive phase retrieval (PR) formulas for the expectation and variance of the radiated fields. This paper rigorously analyze the uniqueness and stability of the recovered statistics of the radiated fields. Afterwards, since the direct problem has a unique mild solution, by examining the expectation and variance of this solution and combined with the phase retrieval formulas, we derive the Fredholm integral equations to solve the inverse random source problem (IRSP). We prove the stability of the corresponding integral equations. To quantify the uncertainty of the random source, we utilize the Bayesian method to reconstruct the random source and establish the well-posedness of the posterior distribution. Finally, numerical experiments demonstrate the effectiveness of the proposed method and validate the theoretical results.

keywords: Stochastic Helmholtz equations; Inverse random source problem; Phase retrieval; Bayesian inference

1 Introduction

The inverse source problem, which involves reconstructing an unknown source from measurements of the radiated fields, has broad applications in scientific and engineering areas, such as acoustic tomography [26, 24], seismic exploration [32, 6], medical imaging [1, 12], and antenna design [29, 33]. Existing studies of the inverse source problem mostly consider that the source is a deterministic function [13, 5, 35]. However, in many practical situations, the source is better modeled as a random field. Compared to deterministic inverse source problems, inverse random source problems (IRSP) involve inherent randomness and uncertainties, which make the problem more difficult. In this case, the solution of the forward problem is a random function. The statistics of the solution, such as expectation, variance and higher-order moments, play a more significant role. Recent studies have proposed various theoretical and numerical methods for addressing IRSP. In particular, the IRSP associated with the Helmholtz equation or time-harmonic acoustic model is addressed in this paper, which refers to [3, 4, 7, 10, 20, 21, 22, 23, 25].

In existing research, most studies on reconstructing sources rely on complete datasets, including both intensity and phase information. However, acquiring the full data is often unfeasible or even impossible. Typically, only the intensity or modulus of the field can be measured, espectially in high-frequency regime. This limitation has motivated the research for inverse source problems with phaseless data [8, 37, 15]. Xu et al. [34] is concerned with uniqueness in inverse acoustic scattering with phaseless far-field data at a fixed frequency. O. Ivanyshyn and R. Kress [14] propose a modification of the Newton-type algorithm to recover a real-valued surface impedance from phase-less far field data. Inverse scattering problems using phaseless data have been widely examined by many researchers [2, 17, 18, 19, 38, 39].

Considering IRSP when only phaseless data is available, the problem becomes even more challenging. The reconstruction process is more complex, and the theoretical analysis of uniqueness and stability encounters greater difficulties. It is precisely for this reason that related results are nearly absent.

This paper is concerned with the problem of reconstructing a random acoustic source in the two-dimensional Helmholtz equation from multi-frequency phaseless data. More precisely, consider the scattering problem of the Helmholtz equation in a homogeneous medium

[TABLE]

where the wavenumber $k>0$ is a constant, $f$ is the electric current density, and $u$ is the radiated random wave field that satisfies the Sommerfeld radiation condition. Here, $f$ is assumed to be a random function driven by an additive noise, given by

[TABLE]

where $g$ and $\sigma\geq 0$ are two deterministic real functions with compact supports contained in the rectangular domain $\Omega_{0}\subset\mathbb{R}^{2}$, and $\dot{W}_{x}$ represents the white noise. In this random source model, the mean, standard deviation, and variance of $f$ are $g$, $\sigma$, and $\sigma^{2}$, respectively. We assume that $g\in L^{2}(\Omega_{0})$ and $\sigma\in L^{4}(\Omega_{0})$. Then, there exists a unique continuous stochastic process $u$ satisfying

[TABLE]

which is referred to as the mild solution of the stochastic scattering problem (1.1) [4, Theorem 2.7]. Here $G_{k}$ is Green’s function of the Helmholtz equation. Specifically, we have

[TABLE]

where $H_{0}^{(1)}$ is the Hankel function of the first kind of order zero.

In the following, the notation $u(x,k)$ is used to indicate the explicit dependence of the radiated field $u(x)$ on the wavenumber $k$. Suppose that the radiated fields are measured on $\partial B_{R}:=\{x\in\mathbb{R}^{2}:|x|=R\}$, with $\Omega_{0}\subset\subset B_{R}:=\{x\in\mathbb{R}^{2}:|x|<R\}$. Then the IRSP considered in this paper is described as follows:

(IRSP)  Reconstruct the source $f(x)$ from the multi-frequency phaseless statistics data $\{\mathbf{E}(|u(x,k)|^{2});\\ \mathbf{Var}(|u(x,k)|^{2})\}:x\in\partial B_{R},\,k\in\mathbb{K}_{N}\}$, where $N\in\mathbb{N}_{+}$ and $\mathbb{K}_{N}$ is an admissible set consisting of a finite number of wavenumbers.

This problem is challenging due to its inherent non-uniqueness. Specifically, the source term $f(x)$ and its negative counterpart $-f(x)$ generate identical phaseless data, leading to ambiguity in source recovery.

To address this issue, motivated by reference scatterer techniques that are applied to inverse scattering problems [24, 27, 28, 36, 37], we develop a reference source technique by using some artificially added point sources as additional reference sources to retrieve the phase information. By adding these point sources, we derive two systems of equations, one system for recovering the expectation and the other for recovering the variance and the covariance. The point sources are placed at suitable locations, which guarantees that the absolute value of the determinant of the coefficient matrix for each of the two systems of equations has a strictly positive lower bound. Thereby, the two systems are ensured to be uniquely solvable, leading to the stability results for the phase information of both the mean and variance of the radiated fields on $\partial B_{R}$.

Once the phase retrieval formulas are obtained, the problem becomes a standard IRSP. The inverse problem is to reconstruct the random source by utilizing the measured phaseless statistical information together with the phase retrieval formulas. Given the random source, the direct problem is to determine the random wave field. The authors have demonstrated that a unique mild solution exists for the stochastic direct scattering problem in [4]. By analyzing the statistics of the mild solution and incorporating the phase retrieval formulas, we derive Fredholm integral equations to solve the IRSP. We have established the stability results for the phaseless IRSP. To quantify the uncertainty of the random source, we employ the Bayesian method to address the IRSP and have shown the well-posedness of the posterior distribution.

The major contributions of this work can be summarized as follows:

• •

We investigate a class of inverse problems that utilize phaseless statistical data to reconstruct two-dimensional random wave sources. To date, there has been relatively little research compared to the inverse problems involving phase data. Specifically, by introducing reference point sources, we ingeniously integrate phase retrieval techniques, integral equation theory, and Bayesian inference methods to address the challenges posed by inherent ill-posedness and randomness of inverse problems. Consequently, we develop effective numerical algorithms for IRSP.

• •

We prove the stability estimates for the phase retrieval formulas and the corresponding integral equations, and we establish the well-posedness of the posterior measure within the Bayesian inference framework. This provides a theoretical foundation and guidance for the subsequent design of effective reconstruction algorithms.

• •

Based on phase retrieval techniques, integral equation theory, and Bayesian inference methods, we have constructed several reconstruction formulas for random wave sources within the Bayesian framework, relying on different artificial point sources. Finally, we employed a dimension-independent Preconditioned Crank-Nicolson (PCN) sampling algorithm to validate the effectiveness and robustness of these reconstruction formulas against measurement errors.

The outline of this paper is as follows. In Section 2, we introduce the reference source technique and derive phase retrieval formulas for the mean and variance of the radiated fields. We also prove the unique solvability of the phase retrieval formulas and the stability of the recovered statistics. Section 3 is dedicated to the stochastic inverse problem, where we derive Fredholm integral equations. The stability of integral equations are established. In Section 4, we apply the Bayesian method to solve the inverse problem and demonstrate the well-posedness of the posterior distribution. Section 5 verifies the phase retrieval formulas through numerical experiments, and evaluates the performance of the Bayesian method in reconstructing the random source. Finally, conclusions are provided in Section 6.

2 Phase retrieval formulas

In this section, we mainly introduce how to derive two systems of equations for the phase retrieval formulas, one for the expectation and the other for the variance and the covariance, by incorporating artificially added point sources.

Compared to the literature [37], which uses the reference source technique to invert a deterministic source from phaseless data, our study discusses the IRSP with phaseless data. Some necessary notation and definitions are introduced. Without loss of generality, let

[TABLE]

We introduce the definition of an admissible set, which is slightly modified from [37].

Definition 2.1**.**

Let $N\in\mathbb{N}_{+}$ and $k^{*}\in\mathbb{R}_{+}$ be a small wavenumber such that $0<k^{*}R<1.$ Then the admissible set of wavenumbers is given by

[TABLE]

Let $m$ be a positive integer and $\nu_{s}=s\pi/m$ for $s\in\mathbb{R}$. Denote

[TABLE]

For each $j\ (j=1,\cdots,m)$, take two points $z_{j,\ell}:=\lambda_{j,\ell}R\left(\cos\nu_{2j-1},\sin\nu_{2j-1}\right)$, where $\lambda_{j,\ell}\in\mathbb{R},\sqrt{2}a/R\leqslant\left|\lambda_{j,\ell}\right|<1$ and $\ell=1,2$. Then, we have

[TABLE]

The geometrical configuration is illustrated in Figure 2.1.

Let $\delta_{j,\ell}(j=1,\cdots,m,\ell=1,2)$ denote the Dirac distributions at the points $z_{j,\ell}$ and define $\|\cdot\|_{j,\infty}:=\|\cdot\|_{L^{\infty}(\partial B_{j})}.$ $G_{k}(x,z_{j,\ell})$ satisfies the following inhomogeneous Helmholtz equation

[TABLE]

Since $H_{0}^{(1)}(z)=J_{0}(z)+\mathrm{i}Y_{0}(z)$ for $z\in\mathbb{C}\backslash\{0\}$, we have

[TABLE]

where $r_{j,\ell}=|x-z_{j,\ell}|$, $J_{0}$ and $Y_{0}$ are the Bessel functions of the first and second kind with order zero, respectively. Further, $v_{j,\ell}:=u+G_{j,\ell}$ is the unique solution to the problem

[TABLE]

For simplicity, we denote

[TABLE]

for each fixed $j=1,\cdots,m$ and $k\in\mathbb{K}_{\mathrm{N}}$. With these preparations, we consider the following phase retrieval problem:

(PRP)  Given $N\in\mathbb{N}_{+},k\in\mathbb{K}_{N}$ and the phaseless statistical data

$\{\mathbf{E}^{\epsilon}(|u(x,k)|^{2})$, $\mathbf{Var}^{\epsilon}(|u(x,k)|^{2})\}$, $x\in\partial B_{R}$;

$\{\mathbf{E}^{\epsilon}(|v_{j,1,k}|^{2})$, $\mathbf{E}^{\epsilon}(|v_{j,2,k}|^{2})$, $\mathbf{Cov}^{\epsilon}(|v_{j,1,k}|^{2},|v_{j,2,k}|^{2})$, $\mathbf{Cov}^{\epsilon}(|u_{j,k}|^{2},|v_{j,1,k}|^{2})$, $\mathbf{Cov}^{\epsilon}(|u_{j,k}|^{2},|v_{j,2,k}|^{2})$, $\mathbf{Var}^{\epsilon}(|v_{j,1,k}|^{2})$, $\mathbf{Var}^{\epsilon}(|v_{j,2,k}|^{2})\}$ for $j=1,\cdots,m$;

Recover the mean of the radiated fields $\mathbf{E}(u_{j,k})$, the covariance $\mathbf{Cov}(\operatorname{Re}(u_{j,k}),\operatorname{Im}(u_{j,k}))$, and the variance: $\mathbf{Var}(\operatorname{Re}(u_{j,k}))$, $\mathbf{Var}(\operatorname{Im}(u_{j,k}))$ for $j=1,\cdots,m$

Now, we turn to the derivation of the phase retrieval formulas. Notice

[TABLE]

Therefore, we have

[TABLE]

We have established a relationship between the phaseless and phase data, enabling us to derive the phase retrieval formulas.

2.1 Recovering expectation

In this section, we derive the phase retrieval formula for the expectation. Taking the expectation on both sides of (2.1), we have

[TABLE]

Hence,

[TABLE]

where

[TABLE]

Thus, by a direct calculation, the phase retrieval formula for the expectation on $\partial B_{j}$ is obtained:

[TABLE]

where the function matrices $A_{j,k},A_{j,k}^{R},A_{j,k}^{I}$ are defined as follows

[TABLE]

Therefore, the mean of the radiated fields $\mathbf{E}(u_{j,k})$ can be recovered from $\mathbf{E}(u_{j,k})=\mathbf{E}(\operatorname{Re}(u_{j,k}))+i\mathbf{E}(\operatorname{Im}(u_{j,k}))$ for $j=1,\cdots,m.$

2.2 Recovering variance and covariance

In this section, we will derive the phase retrieval formula for the variance.

Based on (2.1), for each fixed $j=1,\cdots,m$ and $k\in\mathbb{K}_{\mathrm{N}}$, we have

[TABLE]

From (2.1) and (2.5), we obtain

[TABLE]

and

[TABLE]

Using equations (2.6) and (2.7), we derive the variance retrieval formula

[TABLE]

where

[TABLE]

Thus, by a simple calculation, the variance and covariance of real and imaginary parts on $\partial B_{j}$ is given by

[TABLE]

where the function matrices $D_{j,k},D_{j,k}^{(1)},D_{j,k}^{(2)},D_{j,k}^{(3)}$ are defined as follows

[TABLE]

Therefore, the variance of the radiated fields $\mathbf{Var}(u_{j,k})$ can be recovered from $\mathbf{Var}(u_{j,k})=\mathbf{Var}(\operatorname{Re}(u_{j,k}))+\mathbf{Var}(\operatorname{Im}(u_{j,k}))$ for $j=1,\cdots,m.$

2.3 Phase retrieval algorithm

Similar to [37], the parameters $R,m$ and $\lambda_{j,\ell}$ for $\ell=1,2$ are taken as:

[TABLE]

The selection of parameters (2.11) ensures the equations (2.2) and (2.8) are uniquely solvable. It will be discussed in the next section. However, in most cases, the measured data contain noise, which means that formulas (2.4) and (2.10) cannot be directly applied in computations. Therefore, we present an algorithm with perturbed data. Algorithm 1 describes the PR algorithm process.

2.4 Stability of phase retrivel

In this section, we will prove that the norm of the determinant of the coefficient matrix for each system of equations (2.2) and (2.8) has a strictly positive lower bound. This ensures the uniqueness of phase information of the radiated fields. Additionally, we will also further analyze stability results of the recovered statistics.

Theorem 2.1**.**

Under (2.11), we have the following estimate

[TABLE]

and

[TABLE]

where $M=\frac{20\tau-7}{20\tau}$ and $M^{*}=\frac{4}{9}.$

Proof.

Regarding (2.12), the proof can be found in [37]. Now, we proof (2.13). Set $s_{\ell}=kr_{j,\ell}$, we have

[TABLE]

Since

[TABLE]

it follows that

[TABLE]

∎

Remark 2.1*.*

We would like to emphasize that, aside from the parameter settings selected in (2.11), there are alternative strategies available that can ensure the inequalities (2.12) and (2.13) hold, thereby guaranteeing the unique solvability of the two systems of equations (2.2) and (2.8).

Now, we proceed to investigate the stability of the phase retrieval formula. For a fixed $k$ and $j$, we consider the following perturbed linear system with the unknowns $\mathbf{E}^{\epsilon}(\mathrm{Re}(u_{j,k}))$ and $\mathbf{E}^{\epsilon}(\mathrm{Im}(u_{j,k}))$:

[TABLE]

where

[TABLE]

It is easy to observe that the solutions to the perturbed equations (2.16) can also be expressed using formula (2.4), with $f_{j,\ell,k}^{\epsilon}$ replacing $f_{j,\ell,k}$ ($\ell=1,2$). Let

[TABLE]

and

[TABLE]

Before analyzing stability, we state a boundedness lemma for $J_{0}$ and $Y_{0}$, whose estimates are derived in the proof of Theorem 3.2 of [37].

Lemma 2.2**.**

Let $t_{\ell}=kr_{j,\ell}$ for $\ell=1,2$.

(i) If $k\in\mathbb{K}_{N}\setminus\{k^{*}\}$, then

[TABLE]

(ii) If $k=k^{*}=\dfrac{\pi}{30a}$, then

[TABLE]

We first establish stability for the expectations.

Theorem 2.3**.**

Under (2.11), it follows that

[TABLE]

where $C>0$ is a constant depends on $k,R,M,M^{*},\lambda_{j,2}$.

Proof.

From (2.2), the linear systems are given by

[TABLE]

By inverting $A_{j,k}$, we obtain

[TABLE]

which implies

[TABLE]

Using the explicit form

[TABLE]

and applying Lemma 2.2 and (2.12), we deduce that

[TABLE]

Finally, we have

[TABLE]

∎

We now turn to analyze the stability of the variance. Define

[TABLE]

[TABLE]

Theorem 2.4**.**

Under (2.11), we have

[TABLE]

where $C>0$ is a constant depends on $k,R,M,M^{*},\lambda_{j,2}$.

Proof.

First, consider the real part.

Case 1: $k\in\mathbb{K}_{N}\backslash\{k^{*}\}$. From (2.10), (2.13) and (2.18), we obtain

[TABLE]

with

[TABLE]

Case 2: $k=k^{*}$. A uniform bound holds:

[TABLE]

with

[TABLE]

Analogous estimates hold for the imaginary part and for the covariance term. Hence, we obtain the $2$-norm bound:

[TABLE]

where

[TABLE]

This completes the proof. ∎

3 Stochastic inverse problem

In this section, we derive the Fredholm integral equations to reconstruct random source by examining the expectation and variance of the mild solution. Furthermore, the stability of the IRSP has been proved, and thus the stability of the phaseless IRSP is also established.

3.1 Integral equations

By taking the expectation on both sides of (1.3) and applying

[TABLE]

it deduces

[TABLE]

This equation can be employed to reconstruct $g$. The reconstruction formula for $g$ is similar to a deterministic inverse problem. The difference is that the known boundary data is obtained from the expectation of the radiation wave field. To simplify the solution process, complex-valued quantities are separated into their real and imaginary parts. More explicitly, we have

[TABLE]

and

[TABLE]

Taking the expectation on both sides of (3.3) and (3.4), we obtain

[TABLE]

[TABLE]

Substituting the real and imaginary parts of the two-dimensional Green function into the above equations yields

[TABLE]

[TABLE]

We can combine two equations to perform the inversion of $g$.

From Lemma 7.1 of Appendix 7.1, it finds

[TABLE]

[TABLE]

and

[TABLE]

By taking the variance on both sides of (3.3) and (3.4), it leads to the following equations

[TABLE]

[TABLE]

Furthermore, substituting the real and imaginary parts of Green function into (3.12) and (3.13), we obtain

[TABLE]

[TABLE]

Then we can recover the variance from above Fredholm integral equations. Using (3.1), (3.3)-(3.6) and (3.11), it follows that

[TABLE]

Hence, we obtain

[TABLE]

that is,

[TABLE]

3.1.1 Stability of integral equations

In this section, we establish the stability result of the integral equations (3.2) and (3.13).

Theorem 3.1**.**

Assume that $g\in L^{2}(\Omega_{0})$ and $\|g\|_{H_{p}(\Omega_{0})}\leqslant M$. Then, for any $k\in\mathbb{K}_{N}$, the following estimate holds

[TABLE]

where

[TABLE]

Proof.

Extend $g(x)$ by zero from $\Omega_{0}$ to $\mathbb{R}^{2}$, defining it as follows:

[TABLE]

Setting

[TABLE]

From (3.2), it follows that

[TABLE]

where $G_{k}(x,y)$ is the Green’s function and $*$ denotes convolution. Applying the Fourier transform to both sides, and using the convolution theorem, we obtain

[TABLE]

where $\widehat{G_{k}}(\omega)$ is the Fourier transform of $G_{k}(x)$, and $\hat{g}(\omega)$ is the Fourier transform of $g(x)$.

Using Parseval’s theorem, Hölder inequality and $\widehat{G_{k}}(\omega)=\frac{1}{\sqrt{2\pi}}\frac{1}{k^{2}-|\omega|^{2}}$, we have

[TABLE]

Combining the above estimates, we conclude that

[TABLE]

where

[TABLE]

∎

It shows that the norm of $g(x)$ can be controlled by the norm of $\mathbf{E}(u(x,k))$ in $L^{2}(\mathbb{R}^{2})$. We now proceed to prove that $\sigma^{2}(x)$ can similarly be controlled.

Lemma 3.2**.**

Let $f(x)\in L^{2}(\mathbb{R}^{2})$ with $\|f(x)\|^{2}_{L^{2}(\mathbb{R}^{2})}=A>0$. Then, for any $\epsilon>0$, there exists $N^{\prime}(\epsilon)>0$ such that for all $N\geqslant N^{\prime}$,

[TABLE]

Proof.

By Parseval’s theorem, we have

[TABLE]

Taking the limit as $N\to\infty$, we obtain

[TABLE]

Thus, for any $\epsilon>0$, there exists $N_{1}(\epsilon)>0$ such that for all $N>N_{1}$,

[TABLE]

Moreover, we also have

[TABLE]

Hence, for any $\epsilon>0$, there exists $N_{2}(\epsilon)>0$ such that for all $N>N_{2}$,

[TABLE]

which implies

[TABLE]

Setting $\epsilon=\frac{A}{2}$ and choosing $N_{3}=\max\{N_{1},N_{2}\}$, we obtain that for all $N>N_{3}$,

[TABLE]

∎

Remark 3.1*.*

It demonstrate that the energy of $f(x)$ is concentrated in the low-frequency region, meaning that the low-frequency energy can effectively control the high-frequency energy.

Theorem 3.3**.**

Assume that $\sigma^{2}\in L^{2}(\Omega_{0})$. Then, for any $k\in\mathbb{K}_{N}$, the following estimate holds

[TABLE]

where

[TABLE]

Proof.

Extend $\sigma^{2}(x)$ by zero from $\Omega_{0}$ to $\mathbb{R}^{2}$ as

[TABLE]

From (3.15), it follows that

[TABLE]

Taking Fourier transform on both sides and using the convolution theorem, we get

[TABLE]

Recall that

[TABLE]

Using Parseval’s theorem and Lemma 3.2, we have

[TABLE]

By Parseval’s identity, this implies

[TABLE]

where

[TABLE]

∎

Remark 3.2*.*

It shows that the norm of $\sigma^{2}(x)$ can be controlled by the norm of $\mathbf{Var}(\operatorname{Im}(u(x,k))$ in $L^{2}(\mathbb{R}^{2})$. However, the same conclusion does not directly apply to $\mathbf{Var}(\operatorname{Re}(u(x,k))$ and $\mathbf{Cov}(\operatorname{Re}(u(x,k),\operatorname{Im}(u(x,k))$. In addition, it should be noted that the right-hand sides of stability estimates (3.19) and (3.20) are based on statistical data of phase-scattered fields in the whole space, rather than data from measurement curve (circle). Therefore, we cannot derive the stability for phaseless IRSP from the stability estimates of the phase retrieval formulas. Unless the scattered field exhibits exceptional regularity, such as analyticity, we cannot control the whole space data using the data from low-dimensional measurement curve. However, the scattered fields induced by random sources often lack such regularity. Thus, we need to develop new methods to obtain stability estimates for IRSP with phaseless statistical data, which will be an important focus for us in the future.

4 Bayesian inference

Bayesian inference has become a powerful framework for solving inverse problems, as it allows for the quantification of uncertainty in the inferred parameters [31, 16]. Consider the inverse problem of finding $\phi$ from $\psi$, where the relationship between $\phi$ and $\psi$ is given by the model:

[TABLE]

here, $\mathcal{T}:X\to\mathbb{R}^{m}$ is a measurable mapping, referred to as the forward operator. $X$ is a separable Hilbert space with inner product $\langle\cdot,\cdot\rangle_{X}$. $\phi\in X$ is the unknown function. The observed data $\psi\in\mathbb{R}^{m}$ is finite-dimensional. The noise $\eta$ follows an $m$-dimensional Gaussian noise with zero mean and covariance operator as $\Sigma$, i.e., $\eta\thicksim N(0,\Sigma)$.

The Bayesian approach relies on deriving the posterior distribution of the unknown parameter $\phi$ based on a prior distribution and the observed data, with the forward model connecting the parameter space to the data space. Let $\mu_{pr}$ be the prior measure on $\phi$. The posterior distribution of $\phi$, conditioned on the observed data $\psi$, is given by the measure $\mu^{\psi}$, which satisfies the following relation:

[TABLE]

where the left hand side of (4.2) is the Radon-Nikodym derivative of the posterior distribution $\mu^{\psi}$ with respect to the prior $\mu_{pr}$ and

[TABLE]

which is commonly referred to as the data fidelity term in deterministic inverse problems.

It is well known that Markov Chain Monte Carlo (MCMC) methods are widely used in Bayesian inference to explore the posterior distribution. In our work, we implement the Preconditioned Crank-Nicolson (pCN) algorithm to generate samples from the posterior distribution $\mu^{\Psi}$ given by equation (4.2). This method, originally developed in [11, 9, 30], is particularly advantageous due to its dimension-independent properties.

4.1 Well-posedness

In this section, we discuss the well-posedness of the posterior distribution arising from the Bayesian inference with Gaussian prior. The regularity properties of Green’s function play an important role in the subsequent analysis.

Lemma 4.1**.**

Let $\Omega_{0}\subset\mathbb{R}^{2}$ be a bounded domain. It holds that for any $p>0$, $G_{k}(x,y)\in L^{p}(\Omega_{0})$ for any $y\in\Omega_{0}$.

Proof.

Set $\rho=\sup_{x,y\in\Omega_{0}}|x-y|$, it implies $\bar{\Omega}_{0}\subset B_{\rho}(y)$. Notice that

[TABLE]

here $V$ is a continuous function.

For any $p>0$, it follows that

[TABLE]

This completes the proof. ∎

We now derive the forward operator equation. By combining equations (2.2), (3.5), and (3.6), we can obtain

[TABLE]

where

[TABLE]

with

[TABLE]

and

[TABLE]

Hence the noisy measured data $f_{\ell,k}$ are generated by

[TABLE]

where $\eta_{\ell,k}\thicksim\mathcal{N}(0,\Sigma_{\ell,k})$ .

The forward operator $T_{\ell,k}:L^{2}(\Omega_{0})\rightarrow\mathbb{R}^{m}$ for the inversion of $g$ needs to satisfy the following lemma as in [30]:

Lemma 4.2**.**

(i) For every $\delta>0$ there exists $M_{\ell,k}=M_{\ell,k}(\delta)\in\mathbb{R}$ such that, for all $g\in L^{2}(\Omega_{0})$,

[TABLE]

(ii) For every $r>0$, there exists $K_{\ell,k}=K_{\ell,k}(r)>0$ such that, for all $g_{1},g_{2}\in L^{2}(\Omega_{0})$ with

[TABLE]

we have

[TABLE]

Proof.

Throughout the proof, the constant $C_{\ell,k}$ changes from occurrence to occurrence. For any $g\in L^{2}(\Omega_{0})$, according to Triangle inequality, Cauchy-Schwarz inequality, Lemma 2.2 and Lemma 4.1, we estimate

[TABLE]

From Young’s inequality, for every $\delta>0$, it follows that

[TABLE]

The combination of (4.6) and (4.7) yield

[TABLE]

Using the linearity of $T_{\ell,k}$, we can easily obtain

[TABLE]

This completes the proof. ∎

Definition 4.1** ([30]).**

Assume that two probability measures $\mu$ and $\mu^{\prime}$ both absolutely continuous with respect to the same reference measure $\nu$, then Hellinger metric with respect to measure $\mu$ and $\mu^{\prime}$ is defined by

[TABLE]

The Lemma 4.2 about $T_{\ell,k}$ can derive the bounds and Lipschitz properties of the corresponding data fidelity term $\Psi_{\ell,k}$ as Assumptions 2.6 in [30]. As a result, the probability measure $\mu^{f^{\epsilon}_{\ell,k}}$ defined by equation (4.2) is well defined on $L^{2}(\Omega_{0})$ and it is Lipschitz continuous with respect to the data $f^{\epsilon}_{\ell,k}$ under the Hellinger metric, as stated in the following theorem.

Theorem 4.3**.**

*Let forward operator $T_{\ell,k}:L^{2}(\Omega_{0})\rightarrow\mathbb{R}^{m}$ satisfies Lemma 4.2. For a given $f^{\epsilon}_{\ell,k}\in\mathbb{R}^{m}$, $\mu^{f^{\epsilon}_{\ell,k}}$ is given by equation (4.2). Then we have the following:

$(i)$ $\mu^{f^{\epsilon}_{\ell,k}}$ defined as equation (4.2) is well defined on $L^{2}(\Omega_{0})$.

$(ii)$ $\mu^{f^{\epsilon}_{\ell,k}}$ is Lipschitz in the data $f^{\epsilon}_{\ell,k}$ with respect to the Hellinger metric. Specifically, if $\mu^{f^{\epsilon}_{\ell,k}}$ and $\mu^{{f_{\ell,k}^{\prime,\epsilon}}}$ are two measures corresponding to data $f^{\epsilon}_{\ell,k}$ and ${f_{\ell,k}^{\prime,\epsilon}}$ respectively, then for every $r>0$, there exists $C_{\ell,k}=C_{\ell,k}(r)>0$ such that for all $f^{\epsilon}_{\ell,k}$, ${f_{\ell,k}^{\prime,\epsilon}}\in\mathbb{R}^{m}$ with $\max\left\{\|f^{\epsilon}_{\ell,k}\|_{L^{2}(\Omega_{0})},\|{f_{\ell,k}^{\prime,\epsilon}}\|_{L^{2}(\Omega_{0})}\right\}<r$, we have*

[TABLE]

The above theorem is a direct consequence of the fact that $\Psi_{\ell,k}$ satisfies the Assumptions 2.6 in [30], and therefore we omit the proof here.

Combining equations (2.8), (3.12), (3.13) and (3.17), we obtain the following representation

[TABLE]

where

[TABLE]

with

[TABLE]

and

[TABLE]

Hence the noisy measured data $F_{k,\ell_{1},\ell_{2}}$ are generated by

[TABLE]

where $\eta_{k,\ell_{1},\ell_{2}}\thicksim\mathcal{N}(0,\Sigma_{k,\ell_{1},\ell_{2}})$. The forward operator $T_{k,\ell_{1},\ell_{2}}:L^{2}(\Omega_{0})\rightarrow\mathbb{R}^{m}$ for the inversion of $\sigma^{2}$ needs to satisfy the following lemma as in [30]. The proof follows a similar procedure as in Lemma 4.2, thus we omit it here.

Lemma 4.4**.**

(i) For every $\delta>0$ there exists $M_{k,\ell_{1},\ell_{2}}=M_{k,\ell_{1},\ell_{2}}(\delta)\in\mathbb{R}$ such that, for all $\sigma^{2}\in L^{2}(\Omega_{0})$,

[TABLE]

(ii)For every $r>0$, there exists $K_{k,\ell_{1},\ell_{2}}=K_{k,\ell_{1},\ell_{2}}(r)>0$ such that, for all $\sigma^{2}_{1},\sigma^{2}_{2}\in L^{2}(\Omega_{0})$ with

[TABLE]

we have

[TABLE]

Lemma 4.4 imply that the forward operator $T_{k,\ell_{1},\ell_{2}}$ is bounded and Lipschitz continuous, similar as Assumptions 2.6 in [30]. Consequently, the probability measure $\mu^{F^{\epsilon}_{k,\ell_{1},\ell_{2}}}$ defined by equation (4.2) is well-defined on $L^{2}(\Omega_{0})$ and is Lipschitz in the data $F^{\epsilon}_{k,\ell_{1},\ell_{2}}$ with respect to the Hellinger metric, as shown in the following theorem.

Theorem 4.5**.**

*Let forward operator $T_{k,\ell_{1},\ell_{2}}:L^{2}(\Omega_{0})\rightarrow\mathbb{R}^{m}$ satisfies Lemma 4.4. For a given $F^{\epsilon}_{k,\ell_{1},\ell_{2}}\in\mathbb{R}^{m}$, $\mu^{F^{\epsilon}_{k,\ell_{1},\ell_{2}}}$ is given by equation (4.2). Then we have the following:

$(i)$ $\mu^{F^{\epsilon}_{k,\ell_{1},\ell_{2}}}$ defined as equation (4.2) is well defined on $L^{2}(\Omega_{0})$.

$(ii)$ $\mu^{F^{\epsilon}_{k,\ell_{1},\ell_{2}}}$ is Lipschitz in the data $F^{\epsilon}_{k,\ell_{1},\ell_{2}}$ with respect to the Hellinger metric. Specifically, if $\mu^{F^{\epsilon}_{k,\ell_{1},\ell_{2}}}$ and $\mu^{{F_{k,\ell_{1},\ell_{2}}^{\prime,\epsilon}}}$ are two measures corresponding to data $F^{\epsilon}_{k,\ell_{1},\ell_{2}}$ and $F_{k,\ell_{1},\ell_{2}}^{\prime}$ respectively, then for every $r>0$, there exists $C_{k,\ell_{1},\ell_{2}}=C_{k,\ell_{1},\ell_{2}}(r)>0$ such that $\max\left\{\|F^{\epsilon}_{k,\ell_{1},\ell_{2}}\|_{L^{2}(\Omega_{0})},\|{F_{k,\ell_{1},\ell_{2}}^{\prime,\epsilon}}\|_{L^{2}(\Omega_{0})}\right\}<r$ for all $F^{\epsilon}_{k,\ell_{1},\ell_{2}}$, ${F_{k,\ell_{1},\ell_{2}}^{\prime,\epsilon}}\in\mathbb{R}^{m}$, we have*

[TABLE]

Remark 4.1*.*

Theorem 4.3 and Theorem 4.5 ensure that the posterior measure (or posterior density) in the corresponding Bayesian formula is well-defined. In particular, the stability of the posterior measure with respect to the measurement data provides a theoretical foundation for the subsequent design of numerical algorithms.

5 Numerical examples

In the section, we present several numerical results of a two-dimensional example to evaluate the effectiveness of the proposed method. Next, we solve the IRSP using Bayesian inference methods.

First, the numerical approximation to the forward problem is obtained by discretizing the integral equations. We utilize Trapezoidal rules to approximate integrals, i.e.,

[TABLE]

[TABLE]

where $y_{1}^{i}=a+i\Delta y_{1},i=0,1,\cdots,N_{y_{1}}$ and $y_{2}^{j}=c+j\Delta y_{2},j=0,1,\cdots,N_{y_{2}}$ are the grid points in the $y_{1}$ and $y_{2}$ directions. $I_{1}$ is a deterministic integral and $I_{2}$ is a stochastic integral. We producing $v_{j,\ell}$ by merely adding $u_{j,\ell}$ and $\Phi_{j,\ell}$ together. To evaluate the stability of the PR algorithm, we consider not only the discretization error of the forward solver but also add Gaussian distributed random noises to the measurements. The noisy phaseless data is generated by the following formulas:

[TABLE]

where $r$ is a standard Gaussian distributed random number, and $\epsilon>0$ is the noise level.

In our experiment, the test source function is chosen as the following

[TABLE]

and

[TABLE]

and reconstruct the mean $g_{1}$ and the variance $\sigma_{1}$ given by

[TABLE]

where $x_{1},x_{2}\in\Omega_{0}=[-1,1]\times[-1,1]$. See Figure 5.2 for the surface plot of the exact $g_{1}$ and $\sigma_{1}$ in the domain $\Omega_{0}$.

5.1 Phase retrieval

In this section, a numerical validation of the phase retrieval technique will be shown. First we specify the parameters used in (2.11): $a=1,\tau=6,R=\tau a=6,m=10,\lambda_{1,1}=0.5,\nu=0.1\pi$ and thus $k^{*}=\pi/30$. The other wavenumbers are taken as $k=\pi$, $k=60.5\pi$ and $k=84.5\pi$. We conducted experiments to confirm the accuracy and stability of the phase retrieval formulas. Table 1 and 2 show the relative errors of phase retrieval for $\mathbf{E}(u)$ and $\mathbf{Var}(u)$ with different wavenumbers and noise levels. The discrete relative errors are calculated as follows:

[TABLE]

[TABLE]

[TABLE]

here, $x_{n}\in\partial B_{j}$ $(j=1,\dots,m)$, $N_{j}$ denotes the number of measurements on $\partial B_{j}$, $\mathbf{E}(u)$ and $\mathbf{E}^{\epsilon}(u)$ represent the exact and retrieved means of the radiated fields, $\mathbf{Var}(u)$ and $\mathbf{Var}^{\epsilon}(u)$ represent the exact and retrieved variances, $\mathbf{Cov}(\operatorname{Re}(u),\operatorname{Im}(u))$ and $\mathbf{Cov}^{\epsilon}(\operatorname{Re}(u),\operatorname{Im}(u))$ represent the exact and retrieved covariance. The specifications are as follows: $N_{j}=120,j=1,\cdots,m.$ $N_{y_{1}}=60,N_{y_{2}}=60,\Delta t=0.0022,M=5\times 10^{4}.$ Here, $M$ denotes the total number of samples in the pCN algorithm.

From Table 1, 2 and 3, the PR algorithm can accurately recover the phase information at different wavenumbers in small noise level. As we expect, the relative errors gradually decrease as the noise level decreases when the wavenumber is fixed. The conclusions of Theorems 2.3 and 2.4 can be verified by these results. This demonstrates the effectiveness and feasibility of PR algorithm.

5.2 Inverse random source

In this section, we present the numerical results of reconstructing the random source using the Bayesian method.

We can reconstruct the mean $g_{1}$ via the equations (4.5) for a set of wavenumbers $k\in\{k_{1},k_{2},\cdots,k_{N}\}$. Similarly, the variance $\sigma_{1}$ can be reconstructed using the equations (4.9) for the same set of wavenumbers. In the following examples, we take $k_{j}=(j+82.5)\pi,j=0,1,2,3$. We can simply discretize equation (4.5) and (4.9) to obtain a discrete linear system by using Trapezoidal rules on a uniform grid.

Specifically, we choose the covariance operator $\mathcal{C}_{0}$ of the Gaussian prior $\mathcal{N}(0,\mathcal{C}_{0})$ is given by

[TABLE]

where $\gamma=0.001$ and $d=0.2$ in the subsequent numerical experiment of reconstructing $g_{1}$. Additionally, we fix $\beta=0.05$ and draw $10^{5}$ samples from the posterior distribution using pCN algorithm, with a burn-in period of $5\times 10^{3}$. The error of the numerical experiments is measured using relative $L^{2}$ error

[TABLE]

where $g_{1}$ is the true solution and $\hat{g}_{1}$ is the reconstructed result obtained using our algorithm. When solving inverse problems, we employ a coarse grid to avoid the “inverse crimes” , i.e., $N_{j}=60,\ j=1,\cdots,m$. $N_{y_{1}}=30$, $N_{y_{2}}=30.$

The numerical results for reconstructing the source term $g_{1}$ are shown in Figure 5.3. For $N=4$, the inversion attains a satisfactory outcome. The inversion results using $f^{\epsilon}_{1,k}$ or $f^{\epsilon}_{2,k}$ are similar.

The inversion process for $\sigma_{1}$ is the same as that for $g_{1}$. When we use $F^{\epsilon}_{k,1,1}$ or $F^{\epsilon}_{k,2,2}$ to reconstruct $\sigma_{1}$, we take $\gamma=0.01,d=0.2,\beta=0.01$ and draw $10^{5}$ samples from pCN algorithm with $5\times 10^{3}$ in burn-in period. When we use $F^{\epsilon}_{k,1,2}$, we change paramater $\gamma$ to 0.05. From Table 4 and Figure 5.4, we can accurately recover the true solution in all cases.

In Table 4, the reconstruction of $g_{1}$ becomes more accurate as the noise level decreases. Note that even with low levels of noise, the reconstruction remains quite precise and stable, regardless of whether $f^{\epsilon}_{1,k}$ or $f^{\epsilon}_{2,k}$ is used. However, when the noise level is high, the inversion results deteriorate significantly. When reconstructing $\sigma_{1}$, similar results are obtained.

6 Conclusion

In this paper, we have investigated methods to determine acoustic random sources from multi-frequency phaseless data. The reconstruction process consists of two stages. On the first stage, phase retrieval formulas for the mean and variance of the scattered field are derived using the reference source technique, and the uniqueness and stability of the retrieved statistics are proved. On the second stage, based on the mild solution of the Helmholtz equation and phase retrieval formulas, we derive Fredholm integral equations for the inverse scattering problem to reconstruct the mean and variance of the random source. The Bayesian method is employed to solve the inverse problem for the reconstruction of random sources. The stability of the integral equations problems are proved, thereby leading to the stability results for the phaseless IRSP. A two-dimensional numerical example is provided to demonstrate the effectiveness and validity of the proposed method. The proposed framework and methodology can be directly applied to solving a wide range of IRSP with phaseless data, such as the heat and wave equations. Moreover, the methodology can be extended to random source problems in inhomogeneous media.

Acknowledgements

The research of H. Liu was supported by NSFC/RGC Joint Research Scheme, N CityU101/21, ANR/RGC Joint Research Scheme, A-CityU203/19, and the Hong Kong RGC General Research Funds (projects 11311122, 11303125 and 11300821). The research of G. Zheng was supported by the NSF of China (12271151).

7 Appendix

7.1 Appendix A

Lemma 7.1**.**

Let $f(x)$ be a test function with a compact support in $\mathbb{R}^{d}.$ It holds that

[TABLE]

If $g(x)$ is also a fucntion with a compact support in $\mathbb{R}^{d}$, it holds that

[TABLE]

Proof.

The proof of (7.1) can be found in [4]. Now, we prove the third equation. According to [4], we have

[TABLE]

for any $x=(x_{1},\ldots,x_{d})\in\mathbb{R}_{+}^{d}\mathrm{~and~}y=(y_{1},\ldots,y_{d})\in\mathbb{R}_{+}^{d},\mathrm{~where~}x_{j}\wedge y_{j}=\min\{x_{j},y_{j}\}.$ We also have

[TABLE]

and

[TABLE]

From (7.3), (7.4) and (7.5), it follows that

[TABLE]

First, apply integration by parts to $y_{1}$. Since $g(y)$ has a compact support in $\mathbb{R}^{d}$, we have

[TABLE]

Then for $y_{2}$, we also obtain

[TABLE]

Repeat the process for $y_{3},\cdots,y_{d}$, hence

[TABLE]

which completes the proof. ∎