1 Introduction↩︎

Passive imaging with ambient noise aims at extracting information about unknown media or reflectors by analyzing wave fields generated by uncontrolled noise sources. A central idea of this approach is that, assuming the ambient noise is randomly distributed in space and statistically stationary in time, the cross correlation of recorded signals converges to the Green’s function of the underlying wave operator. This principle forms the foundation of classical passive imaging and has led to successful applications in helioseismology, seismology, and related fields; see, for example, [1]–[4].

Inverse source problem concerns determining the radiating source from measurements of the wave fields away from its support, which is an important research subject in inverse scattering theory [5]–[8]. In passive imaging, the stationary uncontrolled noise source is usually assumed to be mean-zero with covariance of the form \(F(t-s)K(x)\delta(x-y)\). Using the passive imaging techniques, the inverse problem of determining the spatial support of the noise source from the boundary cross correlation was addressed by [5]. The strategy is to build imaging functionals which utilize the empirical cross correlation to qualitatively locate the noise source. The passive imaging methodology relies on the assumption that the source is stationary in time. Under this assumption, the temporal stationarity leads to a diagonalization of the noise covariance operator in the frequency domain, allowing the time-domain wave equation to be decomposed into independent single-frequency Helmholtz equations. Based on this decoupling, Helmholtz–Kirchhoff identities can be employed to develop backpropagation-based imaging methods. For quantitative numerical methods for reconstructing the white noise random sources of time-harmonic acoustic waves, we refer the reader to [6], [8].

In some practical scenarios, the noise sources may be active only during finite time intervals or exhibit time-dependent intensities. Hence, the source is not stationary. In this work, we consider the acoustic wave equation in \(\mathbb{R}^d\), \(d=2,3\), \[\begin{cases} \partial_t^2 U(x,t) - \Delta U(x,t) = S(x,t), & (x,t)\in \mathbb{R}^d\times (0,\infty),\\ U(x,0)=\partial_t U(x,0)=0, & x\in\mathbb{R}^d, \end{cases}\] where \(S(x,t)\) denotes a random noise source. We assume that \(S(x,t)\) is a mean-zero stochastic process with covariance function \[\label{eq:cov1} C(S(x,t)S(y,s)) = F(t-s)\, K(x,t)\,K(x,s)\,\delta(x-y),\tag{1}\] where the temporal correlation function \(F(t-s)\) depends only on the time lag. \(K(x,t)\) is a deterministic function describing the time-space support and intensity of the source. We assume that \(K(x,t)\) has a compact support \(D\times[t_1,t_2]\), where \(D\subset\mathbb{R}^d\) is a bounded Lipschitz domain such that \(\mathbb{R}^d\backslash D\) is connected and \(0<t_1<t_2<\infty\). The time interval \([t_1,t_2]\) represents the radiating period of the time-dependent source and \(D\) characterizes the time-invariant shape and location.

The goal of this paper is to reconstruct the spatial support of the random source from the correlation data. It should be noticed that, although the temporal correlation function \(F(t-s)\) is stationary, the time dependence of \(K(x, t)\) implies that the random source is, in general, not stationary in time. As a consequence, the classical passive imaging theory based on temporal stationarity and Helmholtz–Kirchhoff identities may not be directly applicable. On the other hand, in many practical applications, it is required that the measurement points are limited. For instance, only a sparse network of seismometers are available to detect the location of the uncontrolled earthquake [4]. Such a requirement naturally leads to inverse problems with limited measurement positions. To the best of our knowledge, inverse source problems driven by temporally non-stationary random sources remain largely unexplored in the existing literature.

To address the inverse source problem and fulfill the requirement of using only limited measurement points, we employ the multi-frequency based factorization method. The multi-frequency data can be achieved by Fourier transforming the time-domain recordings. Specifically, letting \[u(x,k)=\int_{0}^{\infty} U(x,t)e^{\mathrm{i}kt}\,\mathrm{d}t, \qquad f(x,k)=\int_{t_1}^{t_2} S(x,t)e^{\mathrm{i}kt}\,\mathrm{d}t.\] We have that \(u(\cdot,k)\) satisfies the inhomogeneous Helmholtz equation \[\label{eq:helm} -\Delta u(x,k) - k^2 u(x,k) = f(x,k), \quad x\in\mathbb{R}^d.\tag{2}\] Notice that the resulting random source is frequency-dependent, and the non-stationary temporal behavior of the time-domain source is encoded in the frequency dependence of its second-order statistics.

Let \(G_d(x, t)\) be the fundamental solution to the wave equation with initial condition \(G_d(x, t)=0\) for all \(t<0\), where \(G_2(x, t) = \frac{H(t-|x|)}{2\pi\sqrt{t^2 - |x|^2}}\) with \(H\) denoting the Heaviside function, and \(G_3(x, t) = \frac{\delta(t-|x|)}{4\pi|x|}\). The Fourier transform of the time-dependent Green function is given by \[\hat{G}_d(x, k) = \int_0^\infty G_d(x, t)e^{\mathrm{i}kt}{\rm d}t,\] which has the explicit form \(\hat{G}_2(x, k) = \frac{i}{4}H_0^{(1)}(k|x|)\) and \(\hat{G}_3(x, k) = \frac{e^{\mathrm{i}k|x|}}{4\pi|x|}\). Here \(H_0^{(1)}\) is the Hankel function of the first kind of order zero. Therefore, the solution \(u\) to the Helmholtz equation can be represented by \[u(x,k) = \int_{\mathbb{R}^d}\hat{G}_d(x-y, k)f(y, k){\rm d}y.\] Denoting \(\mathbb{S}^{d-1}:=\{x\in\mathbb{R}^d: |x|=1\}\), due to the asymptotic behavior of the Green function, the solution satisfies the Sommerfeld radiation condition \[\begin{align} \label{src} \lim_{r\to\infty}r^{\frac{d-1}{2}}(\partial_r u - iku) = 0, \quad r=|x|, \end{align}\tag{3}\] which holds uniformly in all directions \(\hat{x} = x/|x|\in\mathbb{S}^{d-1}\). Further, \(u\) admits the asymptotic behavior at infinity \[\begin{align} u(x,k) = \frac{e^{\mathrm{i}k|x|}}{|x|^{\frac{d-1}{2}}}\left\{u^\infty(\hat{x},k) + \mathcal{O}(r^{-\frac{d+1}{2}})\right\}, \quad \text{as}\,\,|x|\to\infty. \end{align}\] Here \(u^\infty(\cdot,k)\in C^\infty(\mathbb{S}^{d-1})\) is the far-field pattern associated with the scattering problem 2 –3 with explicit form \[\label{eq:far} u^\infty(\hat{x},k) =\int_{\mathbb{R}^d} e^{-\mathrm{i}k\hat{x}\cdot y} f(y,k)\,\mathrm{d}y=\int_{D}e^{-{\rm i} k\hat{x}\cdot y}\int_{t_1}^{t_2} S(y,t)e^{{\rm i}kt} {\rm d}t {\rm d}y.\tag{4}\] Under the setting in the frequency domain, the inverse problem is to determine the spatial support of \(K(x,t)\), or equivalently the support of the random source \(S(x,t)\), from multi-frequency far-field random data \[\{u^\infty(\hat{x}_j,k): k\in(k_{\min},k_{\max}),\; j=1,\dots,J\},\] where \(\hat{x}_j\in\mathbb{S}^{d-1}\) are finitely many, possibly sparse, observation directions.

The factorization method was originally introduced by Kirsch in the context of deterministic time-harmonic inverse scattering problems [9], [10]. Its multi-static formulation, based on far-field operators constructed from multi-incident and multi-receiver data, provides a powerful non-iterative tool for reconstructing geometric information about unknown scatterers or sources. A multi-frequency factorization method was proposed by Griesmaier [11], which allows one to recover the \(\Theta-\)convex polygon of the support from the multi-frequency far-field data over sparse observation directions. However, in most existing literature, it is required that the scatterers are frequency-independent. Recently, a novel factorization method has been developed in [12] which is able to deal with inverse source problems in the time domain. Motivated by [12], to locate the random source we aim to develop a multi-frequency factorization method without knowing the temporal correlation function \(F\). Compared with its deterministic counterpart, the current inverse random source problem addresses a fundamentally different inverse problem with space-time random sources and exploits the expectation of the correlated far-field pattern. The reconstruction is therefore performed at the covariance level rather than on single realizations, which requires a new operator formulation and a nontrivial extension of the deterministic factorization framework.

To implement the factorization method, a far-field operator is constructed from the expectation of correlated far-field patterns, which encodes the second-order statistical information of the random source. Then we show that this correlation-based far-field operator admits a symmetric factorization and satisfies a suitable range identity, thereby enabling the reconstruction of the spatial support of the source. A notable advantage of the associated multi-frequency factorization framework is that it requires measurements from only finitely many observation directions, in contrast to the passive imaging method used in [5] which relies on full-aperture boundary measurements. Indeed, for a fixed observation direction, we derive a criterion for characterizing the smallest strip containing the source support and perpendicular to the observation direction. Therefore, by combining data from several observation directions, a convex polygonal approximation of the support can be recovered. Our proposed imaging algorithm is non-iterative, computationally efficient and well-suited for sparse measurement configurations.

The remainder of this paper is organized as follows. Section 2 introduces the far-field operator, presents its decomposition, and constructs the corresponding indicator function. Section 3 validates the theoretical findings through numerical experiments in both two and three dimensions. The factorization method using near-field random data is discussed in the Appendix.

2 Factorization for far-field case↩︎

We propose the factorization method based on the random multi-frequency far-field pattern. As the noise source \(S(x,t)\) is a mean-zero stochastic process, the far-field pattern \(\eqref{eq:far}\) is also a mean-zero random field. Thus, for each fixed observation direction \(\hat{x}\in \mathbb{S}^{d-1}\), we introduce the covariance function \[\label{eq:cov} {C}(\hat{x},\tau, k)={C}(u^{\infty}(\hat{x}, \tau+k), {u^{\infty}(\hat{x}, k)}):=\mathbb{E}(u^{\infty}(\hat{x}, \tau+k)\overline{u^{\infty}(\hat{x}, k)}),\tag{5}\] which captures the accessible statistical information contained in the random multi-frequency far-field data. Substituting the expression of far-field pattern 4 into 5 yields \[\begin{align} \label{eq:cov-g} &&{C}(\hat{x},\tau, k)=\mathbb{E}(u^{\infty}(\hat{x}, \tau+k)\overline{u^{\infty}(\hat{x}, k)})\nonumber\\ &=&\int_{D} \int_{t_1}^{t_2}\int_{t_1}^{t_2} F(t-s)K(y,t)K(y,s)e^{{\rm i}(k+\tau)t}e^{-{\rm i} ks} \,{\rm d} t \,{\rm d}s \, e^{-{\rm i}\tau\hat{x}\cdot y}{\rm d}y\nonumber\\ &=&\int_{D} g(y,\tau,k)\, e^{-{\rm i}\tau\hat{x}\cdot y}{\rm d}y, \end{align}\tag{6}\] where the function \(g(\cdot,\tau,k)\) is given by \[g(y,\tau,k)=\int_{t_1}^{t_2}\int_{t_1}^{t_2} F(t-s)K(y,t)K(y,s)e^{{\rm i}(k+\tau)t}e^{-{\rm i} ks} {\rm d}t {\rm d}s .\] The identity 6 provides the basis for defining a far-field correlation operator and developing a factorization framework for the inverse problem. We are concerned with identifying the position and shape of the spatial support \(D\) of \(K(x,t)\) from the data set of the covariance of the multi-frequency far-field data at sparse observation directions, which is \[\{C(\hat{x}_j, \tau, k): some fixed \;k>0, \tau\in [\tau_{\min},\tau_{\max}], j=1,2,...,J\}\] for some \(J\in \mathbb{N}\).

2.1 Definition of far-field correlation operator and factorization↩︎

Based on the covariance representation 6 and the multi-frequency far-field operator proposed by [11], [12], we introduce the central frequency \(\tau_c\) and half of the bandwidth of the given data as \[\tau_c:=\frac{\tau_{\min}+\tau_{\max}}{2},\quad \tau_w=\frac{\tau_{\max}-\tau_{\min}}{2}.\] For a given reference frequency \(k=k_0>0\), we define a far-field correlation operator \[\begin{align} \label{FO} \mathcal{F}: L^2(0,\tau_w)\rightarrow L^2(0,\tau_w), \end{align}\tag{7}\] which is given by \[\begin{align} &&(\mathcal{F}\phi)(\tau):=(\mathcal{F}^{(\hat{x},k_0)}\phi)(\tau) \\ &=&\int_{0}^{\tau_w} C(\hat{x}, \tau_c+\tau-\gamma,k_0)\,\phi(\gamma)\,\mathrm{d}\gamma, \quad {\tau \in (0,\tau_w)}, \\ &=&\int_{0}^{\tau_w} \int_D g(y,\tau_c+\tau-\gamma,k_0 ) e^{-{\rm i}(\tau_c+\tau-\gamma)\hat{x}\cdot y}\,\phi(\gamma){\rm d}y\,\mathrm{d}\gamma, \\ &=& \int_{0}^{\tau_w } \int_{D} \int_{t_1}^{t_2}\int_{t_1}^{t_2} F(t-s)K(y,t)K(y,s)e^{{\rm i}(k_0+\tau_c+\tau-\gamma)t}e^{-{\rm i} k_0s} \, {\rm d} t {\rm d}s\, e^{-{\rm i}(\tau_c+\tau-\gamma)\hat{x}\cdot y}{\rm d}y\,\phi(\gamma)\,{\rm d}\gamma. \end{align}\]

We proceed to investigate the symmetry properties of the far-field operator \(\mathcal{F}\). Under the covariance-based formulation introduced above, the operator can be decomposed in a symmetric manner, which reveals its intrinsic self-adjoint structure and enables a rigorous factorization analysis. We state and prove a factorization result for the far-field operator in the following proposition.

It is noteworthy that the operator \(\mathcal{L}\) maps a time-spatial-dependent source function supported on \(D\times[t_1,t_2]\) to the multiple-frequency far-field correlation data at the far-field observation direction \(\hat{x}\). It plays the same role as the data-to-pattern operator for multi-static inverse scattering problems at a fixed energy [10]. To establish a rigorous connection between the range of the far-field operator \(\mathcal{F}\) and that of \(\mathcal{L}\), we introduce the self-adjoint and positive-definite operator \(\mathcal{F}_{\#} := |\operatorname{Re}\mathcal{F}| + |\operatorname{Im}\mathcal{F}|\). Here, the real and imaginary parts are defined as \(\operatorname{Re}\mathcal{F} = (\mathcal{F} + \mathcal{F}^*)/2\) and \(\operatorname{Im}\mathcal{F} = (\mathcal{F} - \mathcal{F}^*)/2\), respectively, ensuring that the square root \(\mathcal{F}_{\#}^{1/2}\) is well-defined.

To facilitate the subsequent analysis of the operator’s properties, we introduce the following hypotheses:

• Hypothesis (H1): The temporal correlation function \(F(t)\) is continuous and satisfies \(F(t)>0\) for all \(t\geq0\).

• Hypothesis (H2): \(|K(x,t)|>0\) a.e. on \(D\times[t_1, t_2]\) and the nodal set of \(K\) has the Lebesgue measure of zero. Furthermore, for almost every given \(x \in D\), the mapping \(t \mapsto K(x,t)\) does not change sign on the interval \([t_1, t_2]\).

• Hypothesis (H3): The wave number \(k_0 \ge 0\) and the time window length \(t_2 - t_1\) satisfy the following inequality:\[k_0 (t_2 - t_1) \le \pi - \delta_0,\]where \(\delta_0 \in (0, \pi]\) is a positive constant independent of time and space.

The following proposition from [12] is useful in the subsequent analysis of the factorization.

We are now in a position to establish a fundamental range identity between the far-field operator \(\mathcal{F}\) and the data-to-pattern operator \(\mathcal{L}\), which forms the analytical basis of the indicator functions introduced in the subsequent section.

To translate the abstract range characterization into geometric information, we introduce, for a given observation direction \(\hat{x} \in \mathbb{S}^{d-1}\), the strip \[K_D^{(\hat{x})} := \bigl\{ y \in \mathbb{R}^d : \inf_{z\in D} (\hat{x} \cdot z) < \hat{x} \cdot y < \sup_{z\in D} (\hat{x} \cdot z) \bigr\}.\] Geometrically, \(K_D^{(\hat{x})}\) is the smallest strip containing the support \(D\) and orthogonal to the observation direction \(\hat{x}\). For a given sampling point \(y \in \mathbb{R}^d\), we further define the test function \[\label{eq:test} \phi_{y,\epsilon}^{(\hat{x})}(\tau) = \frac{1}{(t_2-t_1)\,|B_\epsilon(y)|} \int_{t_1}^{t_2} \int_{B_\epsilon(y)} e^{-\mathrm{i}\tau\hat{x}\cdot z} \,e^{\mathrm{i}\tau t} \,\mathrm{d}z\,\mathrm{d}t, \qquad \tau \in (0,\tau_w),\tag{11}\] where \(B_\epsilon(y)\) denotes the ball of radius \(\epsilon\) centered at \(y\). This test function will be used to probe the range of \(\mathcal{L}\) and to characterize the location of the source support.

2.2 Indicator functions and reconstruction of the \(\Theta\)-convex hull↩︎

By Picard’s Theorem and Theorem 4, we introduce an indicator function based on the spectral decomposition of the multi-frequency far-field operator. Since the underlying range characterization follows directly from the established factorization framework, we state the corresponding characterization theorem.

For a fixed observation direction \(\hat{x} \in \mathbb{S}^{d-1}\), we define \[\label{eq:In} \mathcal{I}^{(\hat{x})}(y) := \left[ \sum_{n=1}^{\infty} \frac{ \bigl|\langle \phi^{(\hat{x})}_{y,\epsilon}, \psi^{(\hat{x})}_n\rangle_{L^2(0,\tau_w)}\bigr|^2 }{ |\lambda_n^{(\hat{x})}| } \right]^{-1}, \qquad y \in \mathbb{R}^d,\tag{12}\] where \(\{(\lambda_n^{(\hat{x})},\psi^{(\hat{x})}_n)\}_{n\ge1}\) denotes the eigensystem of the self-adjoint operator \(\mathcal{F}^{(\hat{x})}_\#\), and \(\phi^{(\hat{x})}_{y,\epsilon}\) is the test function centered at the sampling point \(y\).

The indicator function \(\mathcal{I}^{(\hat{x})}\) provides a computational criterion for determining whether a sampling point belongs to \(K_D^{(\hat{x})}\). This criterion is based on whether the associated test function \(\phi_{y,\epsilon}^{(\hat{x})}\) belongs to the range of the data-to-pattern operator \(\mathcal{L}^{(\hat{x})}\), which encodes the geometric information of the source support. Specifically, we have the following result.

Furthermore, we have the following computational criterion:

Theorem 5 shows that, for a single observation direction, the indicator function \(\mathcal{I}^{(\hat{x})}\) yields an exact characterization of the strip \(K_D^{(\hat{x})}\) determined by the multi-frequency far-field correlation data.

In practical applications, the wave fields are measured at only finitely many observation directions \(\{\hat{x}_j\}_{j=1}^J\). To incorporate all available data, we define the combined indicator function \[\mathcal{I}(y) := \left[ \sum_{j=1}^{J} \frac{1}{\mathcal{I}^{(\hat{x}_j)}(y)} \right]^{-1} = \left[ \sum_{j=1}^{J} \sum_{n=1}^{\infty} \frac{ \bigl|\langle \phi^{(\hat{x}_j)}_{y,\epsilon}, \psi^{(\hat{x}_j)}_n\rangle_{L^2(0,\tau_w)}\bigr|^2 }{ |\lambda_n^{(\hat{x}_j)}| } \right]^{-1}, \qquad y \in \mathbb{R}^d .\] Correspondingly, we define the \(\Theta\)-convex hull of \(D\) associated with the observation directions \(\{\hat{x}_j\}_{j=1}^J\) by \[\Theta_D := \bigcap_{j=1}^{J} K_D^{(\hat{x}_j)} .\]

Geometrically, \(\Theta_D\) is the smallest convex set obtained as the intersection of all strips determined by the given observation directions. By Theorem 5, each strip \(K_D^{(\hat{x}_j)}\) can be uniquely determined from the multi-frequency far-field correlation data, and hence their intersection \(\Theta_D\) is uniquely determined as well. In summary, we have the following result.

To conclude this section, we point out that the factorization with multi-frequency far-field data can be carried over to the case of using near-field data. We discuss this extension in the Appendix.

3 Numerical examples↩︎

In this section, we present numerical experiments in both two and three dimensions to demonstrate the effectiveness of the proposed factorization-based method for identifying the support of the noise source from multi-frequency far-field correlation data. All numerical examples are shown to be consistent with the theoretical findings in previous sections.

For a fixed observation direction \(\hat{x}\) and a reference wavenumber \(k_0\), we consider the continuous far-field operator introduced in the theoretical analysis, which acts on the frequency variable through an integral formulation. In the numerical implementation, the frequency interval \((0,\tau_w)\) is uniformly discretized as \[\omega_n = (n-0.5)\Delta\tau, \quad n = 1,\ldots,N,\] with step size \(\Delta\tau = \tau_w/N\). Applying a quadrature approximation to the continuous integral operator, we obtain a finite-dimensional matrix representation \(\mathcal{F}^{(\hat{x},k_0)} \in \mathbb{C}^{N\times N}\), and the far-field correlation operator can be discretized by \[(\mathcal{F}^{(\hat{x},k_0)}\phi)(\tau_n) \approx \Delta\tau\,\sum_{m=1}^{N} C\!\left(\hat{x},\, \tau_c + \tau_n - \gamma_m,\, k_0\right) \phi(\gamma_m),\] where \(\tau_n:=n\Delta \tau\) and \(\gamma_m:=(m+0.5)\Delta\tau\), \(m,n=1,\ldots,N\). Equivalently, the resulting data matrix can be written in the explicit form \[\mathcal{F}^{(\hat{x},k_0)} := \Delta \tau \begin{pmatrix} C(\hat{x},\tau_c-\omega_1,k_0) & C(\hat{x},\tau_c-\omega_2,k_0) & \cdots & C(\hat{x},\tau_c-\omega_N,k_0) \\ C(\hat{x},\tau_c+\omega_1,k_0) & C(\hat{x},\tau_c-\omega_1,k_0) & \cdots & C(\hat{x},\tau_c-\omega_{N-1},k_0) \\ \vdots & \vdots & \ddots & \vdots \\ C(\hat{x},\tau_c+\omega_{N-1},k_0) & C(\hat{x},\tau_c+\omega_{N-2},k_0) & \cdots & C(\hat{x},\tau_c-\omega_1,k_0) \end{pmatrix}.\]

Let \[\bigl\{(\tilde{\lambda}_n^{(\hat{x},k_0)}, \phi_n^{(\hat{x},k_0)}) : n = 1,\ldots,N \bigr\}\] denote an eigensystem of the discrete far-field data matrix \(\mathcal{F}^{(\hat{x},k_0)}\). Based on this spectral decomposition, the associated factorization operator admits the eigensystem \[\bigl\{(\lambda_n^{(\hat{x},k_0)}, \phi_n^{(\hat{x},k_0)}) : n = 1,\ldots,N \bigr\},\] with eigenvalues given explicitly by \[\lambda_n^{(\hat{x},k_0)} = \bigl| \operatorname{Re}(\tilde{\lambda}_n^{(\hat{x},k_0)}) \bigr| + \bigl| \operatorname{Im}(\tilde{\lambda}_n^{(\hat{x},k_0)}) \bigr|.\] For any sampling point \(y \in \mathbb{R}^{d-1}\), we define the discrete test function vector \(\tilde{\phi}_y^{(\hat{x})} \in \mathbb{C}^N\) by discretizing the continuous test function 11 with respect to the frequency variable. The \(n\)-th component of the test vector is given by as \(\epsilon\rightarrow0\) \[\bigl(\phi_y^{(\hat{x})}\bigr)_n = \frac{1}{t_2-t_1} \int_{t_1}^{t_2} e^{-\mathrm{i}\tau_n\hat{x}\cdot y } \, e^{\mathrm{i}\tau_n t}\,\mathrm{d}t, \qquad n = 1,2,\ldots,N.\] Equivalently, the test vector can be written in the vector form \[\tilde{\phi}_y^{(\hat{x})} = \bigl( {\phi}_y^{(\hat{x})}(\tau_1), {\phi}_y^{(\hat{x})}(\tau_2), \ldots, {\phi}_y^{(\hat{x})}(\tau_N) \bigr)^\top.\] Then, substituting the eigensystem of \(\mathcal{F}^{(\hat{x}, k_0)}_\#\) and the discretized test function into the indicator function 12 , we can compute the truncated indicator function numerically.

Noting that the time-dependent source \(S(x,t)\) is real valued, we have \(u^\infty (\hat{x}, - k) = \overline{u^\infty (\hat{x}, k)}\) for all \(k > 0\). We suppose \(\tau_{\min}=0\), then the bandwidth of frequency can be extended from \((0,\tau_{\max})\) to \((-\tau_{\max},\tau_{\max})\) by \(u^{\infty}(\hat{x}, -k)=\overline{u^{\infty}(\hat{x}, k)}\), thus one deduces from these new measurement data with \(\tau_{\min}=-\tau_{\max}\) that \(\tau_c=0\) and \(\tau_w=\tau_{\max}\). Therefore, we obtain \[\begin{align} &&C(\hat{x}, -\tau, k)\\ &=&\int_D e^{\mathrm{i}\tau \hat{x}\cdot y}\int_{t_1}^{t_2} \int_{t_1}^{t_2} F(t-s)K(y,t)K(y,s) e^{{\rm i}(k-\tau)t} e^{-{\rm i} ks}\mathrm{d} s\,\mathrm{d} t\,\mathrm{d} y \nonumber\\ &=&\overline{\int_D \overline{e^{\mathrm{i}\tau \hat{x}\cdot y}}\int_{t_1}^{t_2} \int_{t_1}^{t_2} F(t-s)K(y,t)K(y,s) \overline{e^{{\rm i}(k-\tau)t}} \,\overline{e^{-{\rm i} ks}}\mathrm{d}s\,\mathrm{d}t\,\mathrm{d}y } \nonumber\\ &=&\overline{\int_D {e^{-\mathrm{i}\tau \hat{x}\cdot y}}\int_{t_1}^{t_2} \int_{t_1}^{t_2} F(t-s)K(y,t)K(y,s) {e^{{\rm i}(\tau-k)t}} e^{-{\rm i} (-k)s}\mathrm{d}s\,\mathrm{d}t\,\mathrm{d}y } \nonumber\\ &=& \overline{\mathbb{E} [u^{\infty}(\hat{x},\tau-k) \overline{u^{\infty}(\hat{x}, -k)}]}\\ &=& C(\hat{x}, \tau,-k)= \overline{\mathbb{E} [u^{\infty}(\hat{x},\tau-k) u^{\infty}(\hat{x}, k)]} \end{align}\]

In the following numerical experiments, the random source is modeled as \(S(x,t) = \mathbb{G}(t)K(x,t)\dot{W}(x)\). Here we assume that \(\mathbb{G}(t)\) is a mean-zero Gaussian process with covariance \(F(t-s)=e^{-\frac{(t-s)^2}{0.1^2}}\), and \(\dot{W}(x)\) is spatial white noise which can be interpreted in the distributional sense as the derivative of a Brownian sheet. It is straightforward to verify that the covariance function of \(S(x,t)\) is consistent with the theoretical setting (1 ). For the discretization of the random source, we introduce a regular grid of points \(\{x_i\}_{i=1}^{N_x}\), \(\{t_j\}_{j=1}^{N_t}\) with grid step \(h_x\), \(h_t\) covering the spatial and temporal support of \(K\), respectively. Then Gaussian process \(\mathbb{G}(t)\) can be discretized by a joint Gaussian distribution and spatial white noise \(\dot{W}(x)\) can be approximated by a piecewise function with Gaussian variables \[\dot{W}(x)= \sum_{i=1}^{N_x} |K_i|^{-1}\chi_i(x)\int_{K_i}{\rm d} W(x)={h_x}^{-\frac{d}{2}}\sum_{i=1}^{N_x}\chi_i(x)Z_{i},\] where \(K_i\) represents a cell centered at \(x_i\) with side length \(h_x\), \(\chi_i\) denotes the indicator function of \(K_i\), and \(Z_{i}\) are independent Gaussian random variables with mean 0 and variance 1. We employ the Monte Carlo method to evaluate the expectation of correlation far-field data \(C(\hat{x}, \tau, k)\), and we use (4 ) to compute the far-field pattern for each sample. The number of Monte Carlo samples is set to 10000.

3.1 Numerical examples in \(\mathbb{R}^2\)↩︎

In this subsection, we present several two-dimensional numerical experiments. The frequency bandwidth is chosen as \[\tau \in (0,2\pi), \qquad \tau_n = \frac{n\pi}{10}, \quad n=1,2,\ldots,N_\tau,\quad N_\tau=20,\] with the reference wavenumber \(k_0=\pi/20\) and frequency step \(\Delta\tau=\pi/10\). The observation directions are uniformly distributed on the unit circle and given by \[\hat{x}_n = (\cos\theta_n,\sin\theta_n), \qquad \theta_n = \frac{(n-1)\pi}{N_{\hat{x}}}, \quad n=1,2,\ldots,N_{\hat{x}},\] where \(N_{\hat{x}}=12\). For each observation direction, multi-frequency far-field correlation data are generated and assembled to form the discrete far-field operator used in the factorization algorithm. The sampling region is taken as a square domain containing the support of the source, which is uniformly discretized into a Cartesian grid. At each sampling point, the indicator function associated with the factorization method is evaluated. In all examples below, the reconstructed support is visualized by plotting the indicator function and the boundary of the source support is highlighted by pink solid line, where higher values indicate a higher likelihood that the sampling point belongs to the source support.

Example 1 (Disk-shaped support). In the first example, the noise source \(S(x,t)\) is supported in a disk, \[D=\{(x_1,x_2)\in\mathbb{R}^2:\;x_1^2+x_2^2\le 1\}.\] Figure 1 shows the reconstruction results using far-field data from one or 12 observation directions for a disk-shaped source support with \(K(x,t)=t+1\) and different radiating time period \([0, t_2]\). The smallest strips \(K_D^{(\hat{x})}\) and circular shape of the source support are clearly identified, and the boundary of the disk is well resolved. This example demonstrates that the proposed method is capable of accurately recovering simple convex geometries from a limited number of observation directions and frequencies. We also observe that, as the radiation period of the noise signal increases, the reconstruction quality slightly deteriorates, leading to a mild loss of accuracy in both the shape and the location of the identified support. Nevertheless, the overall geometry of the disk remains clearly recognizable.

Example 2 (Kite-shaped support). In the second example, the source is supported in a kite-shaped domain defined by the parametric curve \[x(\theta) = (x_1,x_2) + \bigl( r\cos\theta + a r (\cos 2\theta - 1), \; b r \sin\theta \bigr), \qquad \theta\in[0,2\pi],\] with parameters \(a=0.65\), \(b=1.5\), and \(r=1\).

The reconstructed images using far-field data from 12 observation directions for a kite-shaped source support with different \(K(x,t)\) and radiating time period \([0, 0.1]\) are shown in Figure 2. For subfigures (a), (c), and (e), the center of the kite-shaped support is located at the origin, whereas in subfigures (b), (d), and (k) the support is shifted away from the origin. It can be observed from subfigures (a), (d), (e), and (k) that the proposed factorization-based method is able to accurately recover the location of the kite-shaped support as well as its convex hull, even in the presence of strong non-convexity and sharp geometric features. In contrast, for subfigures (b) and (c), only partial information on the support boundary is obtained. In these cases, the imaging functional mainly highlights segments of the support edges, while the interior region is not fully reconstructed. This phenomenon can be explained by the relatively weak effective signal intensity associated with these configurations. From an operator-theoretic point of view, the corresponding multi-frequency far-field correlation operator exhibits faster spectral decay, which reduces the contribution of informative singular components and consequently limits the reconstruction capability of the factorization method. Overall, this example demonstrates that the proposed approach can reliably recover both the position and the convex hull of a kite-shaped support in most configurations, while also illustrating the intrinsic dependence of the reconstruction quality on the effective strength of the measured data.

Example 3 (Rounded square support). The third example considers a rounded square-shaped support given by \[x(\theta) = (x_1,x_2) + r(\cos^3\theta+\cos\theta,\;\sin^3\theta+\sin\theta), \qquad \theta\in[0,2\pi],\] with \(r=0.8\). The reconstruction result is shown in Figure 3. Although the boundary contains flat regions and rounded corners, the proposed method still yields a stable and accurate identification of the overall shape. However, it can be observed that the reconstruction near the four rounded corners is less accurate. This phenomenon can be attributed to the limited angular coverage and the relatively small number of observation directions.

Example 4 (Elliptic support). In the final example, the noise source is supported in an ellipse described by \[x(\theta) = (x_1,x_2) + (a\cos\theta,\;b\sin\theta), \qquad \theta\in[0,2\pi],\] where \(a=1.5\), and \(b=0.8\). Figure 4 presents the corresponding reconstruction. The elongation and orientation of the ellipse are clearly recovered, further confirming the capability of the method to distinguish different aspect ratios of source supports. It is worth noting that the reconstruction quality in this example is very similar to that observed for the kite-shaped support in Example 2. This is consistent with the fact that the signal intensity functions are constructed in a comparable manner in both cases, leading to multi-frequency data of similar effective strength. As a result, the proposed method exhibits nearly identical reconstruction performance for these two geometrically different source supports.

3.2 Numerical examples in \(\mathbb{R}^3\)↩︎

In this subsection, we present three-dimensional numerical experiments to further validate the effectiveness of the proposed multi-frequency factorization method. Compared with the two-dimensional case, the three-dimensional setting poses additional challenges due to the increased complexity of the geometry and the limited-view nature of practical observations.

The frequency bandwidth is chosen in the same way as that used in the two-dimensional setting. In each numerical experiment, multi-frequency data are computed for the following three sets of observation directions:

• The first set consists of three canonical directions, \[\hat{x}=\{(1,0,0),(0,1,0),(0,0,1)\}.\]

• The second set contains \(10\) observation directions uniformly distributed on the upper hemisphere.

• The third set contains \(20\) observation directions uniformly distributed on the upper hemisphere.

These three configurations are used to illustrate the influence of the number of observation directions on the reconstruction quality. The sampling region is chosen as a cubic domain containing the support of the source and discretized by a uniform three-dimensional grid. For each sampling point, the indicator function associated with the factorization method is evaluated for visualization. The reconstructed and original source supports is further projected onto the coordinate planes \(y_1 y_2\), \(y_1 y_3\), and \(y_2 y_3\), where the corresponding projection boundaries are delineated by red curves to facilitate comparison.

Example 1 (Cube-shaped support). In the first example, the noise source \(S(x,t)\) is supported in a cube (see Figure 6 (a)), \[D=\{(x_1,x_2,x_3)\in\mathbb{R}^3:\;|x_1|\le 1,\;|x_2|\le 1,\;|x_3|\le 1\}.\] The reconstruction results for different sets of observation directions are shown in Figure 5 and  6. Figure 6 illustrates the reconstructions of \(K_D^{(\hat{x})}\) using far-field data under different observation directions. When data from a single observation direction \(\hat{x}\) are employed, the reconstructed indicator clearly exhibits a slab-like structure, which is consistent with the projection of the true source support onto the coordinate planes. In Figure 6 (b), three carefully chosen observation directions are used, and the overall shape of the source support is well reconstructed. In Figures 6 (c) and 6 (d), far-field data from \(10\) and \(20\) observation directions, respectively, uniformly distributed on the upper hemisphere are adopted. In these cases, the location of the source support is accurately identified; however, the reconstructed shape appears smoother than the true cube. This smoothing effect is mainly caused by measurement noise and data discretization, which tend to suppress high-frequency components and thus round off the sharp corners of the cube.

Figures 7 and 8 present the reconstruction results for the same cube-shaped source support as in Figures 5 and 6, except that the cube is no longer centered at the origin. Figure 7shows the reconstruction obtained from a single observation direction, while Figure 8 corresponds to multiple observation directions. Compared with Figures 5 and 6, the spatial location of the translated source support is still accurately identified; however, for multiple uniformly distributed observation directions, the reconstruction of the shape becomes less accurate, with the sharp corners appearing smoother. This effect is mainly due to the loss of geometric symmetry and the presence of noise, which reduce the sensitivity to high-frequency shape features.

It is observed that, after translating the cube away from the origin, the reconstruction obtained from multiple uniformly distributed observation directions exhibits a degradation in shape recovery, while the location of the source support remains accurately identified. This phenomenon can be attributed to the loss of geometric symmetry and the limited observation aperture. For an off-centered target, the uniformly distributed observation directions on the upper hemisphere no longer provide balanced angular coverage of the source support, which leads to a reduced sensitivity to high-frequency geometric features, such as sharp edges and corners. Moreover, in the presence of noise and finite bandwidth, the multi-directional data tend to suppress high-frequency components, resulting in a smoothing effect on the reconstructed shape.

Example 2 (Ball-shaped support). In the second example, the noise source is supported in a ball, \[D=\{(x_1,x_2,x_3)\in\mathbb{R}^3:\;x_1^2+x_2^2+x_3^2\le 1\}.\] Figures 9 show the reconstruction results for a ball-shaped source support, where the true support is displayed in Figure 9 (a). Figure 9 (b)-(d) present the reconstructions obtained from far-field data corresponding to \(3\), \(10\), and \(20\) observation directions, respectively. Figure 9 considers a ball centered at the origin, while Figure 10 illustrates the case with an off-origin center. In all these configurations, the reconstructed supports are in good agreement with the true ones, and increasing the number of observation directions leads to more stable reconstructions. Compared with the cube-shaped support, the ball-shaped support tends to exhibit slightly improved reconstruction quality, which may be attributed to its smooth boundary and the absence of sharp edges or corners, leading to a more favorable behavior of the inversion process.

Example 3 (Ellipsoidal support). The third example considers an ellipsoidal support given by \[D=\{(x_1,x_2,x_3)\in\mathbb{R}^3:\;x_1^2+9x_2^2+4x_3^2\le 1\}.\] The reconstruction results are shown in Figure 11. Despite the anisotropic geometry of the ellipsoid, a satisfactory reconstruction can already be achieved using data from a small number of observation directions. As the number of observation directions increases, the elongation and orientation of the ellipsoid are recovered with improved accuracy. This example demonstrates that the proposed method is capable of resolving anisotropic source supports in three dimensions, even with limited angular information. Figures 11 and 12 show the reconstruction results for an ellipsoidal source support, corresponding to the same experimental settings as in Figures 9 and 10. The reconstructions exhibit similar behavior and comparable accuracy to the ball-shaped case.

The three-dimensional numerical experiments confirm that the proposed multi-frequency factorization method remains effective in higher dimensions. While limited observation data may lead to reduced resolution, incorporating a larger number of observation directions substantially improves the accuracy and stability of the reconstruction. These results are in full agreement with the theoretical analysis and further demonstrate the applicability of the method to practical three-dimensional inverse source problems.

4 Conclusion↩︎

In this paper, we consider the inverse problem of determining the uncontrolled noise source for the acoustic wave equation in two and three dimensions. The noise source is a random process with compact spatial support and will vanish after a finite period of time. Therefore, it is not stationary in time. By transforming the wave equation into the Helmholtz equation via Fourier transform, we propose an efficient factorization method to locate the support of the random source based on the multi-frequency far-field correlation data. The proposed method only requires sparse observation directions and has low computational cost. Numerical examples in both two and three dimensions are presented to demonstrate the validity and effectiveness of the proposed method. Although this paper concerns the inverse random source scattering problem for the acoustic wave equation, we believe that the proposed framework and methodology can be directly applied to solve inverse random source problems for other wave equations, such as the electromagnetic and elastic wave equations. It is a more challenging direction to consider the inverse random medium scattering problem where the medium should be modeled as a random function. We hope to be able to report the progress on these problems in the future.

5 Factorization for the near-field case in three dimensions↩︎

In this section, we will carry the factorization with multi-frequency far-field correlation data to the case with near-field measurements in three dimensions. More precisely, assuming \(D\subset B_R\), the proposed factorization method can be slightly modified to get an image of the annular region \[\tilde{K}^{(x)}_D:=\{ y\in \mathbb{R}^3: \inf_{z\in D}|x-z|<|x-y|<\sup_{z\in D}|x-z|\}\subset \mathbb{R}^3,\] for every fixed measurement position \(x\in \partial B_R\). Define the near-field operator \(\mathcal{N}^{(x)}: L^2(0,\tau_w)\rightarrow L^2(0,\tau_w)\) by \[\begin{align} &&(\mathcal{N}\phi)(\tau) =(\mathcal{N}^{(x,k_0)}\phi)(\tau) \\ &=&\int_{0}^{\tau_w} \tilde{C}(x, \tau_c+\tau-\gamma,k_0)\,\phi(\gamma)\,{\rm d}\gamma, \quad {\tau \in (0,\tau_w)}, \\ &=& \int_{0}^{\tau_w} \int_{D} \int_{t_1}^{t_2}\int_{t_1}^{t_2} F(t-s)K(y,t)K(y,s)e^{{\rm i}(k_0+\tau_c+\tau-\gamma)t}e^{-{\rm i} k_0 s} \, {\rm d} t {\rm d}s\, \frac{e^{{\rm i}(\tau_c+\tau-\gamma)|x-y|}}{4\pi|x-y|}{\rm d}y\,\phi(\gamma)\,{\rm d}\gamma \end{align}\] where \(\tilde{C}(x, \tau, k):= {\mathbb{E}}[u(x, k+\tau)\overline{u(x, k)}]\). Following the proof of Proposition 1, we obtain a factorization of the near-filed operator: \[\mathcal{N}^{(x)}=\tilde{\mathcal{L}} \tilde{\mathcal{T}} \tilde{\mathcal{L}}^*\] where \(\mathcal{\tilde{L}}=\mathcal{\tilde{L}}_D^{(x,k_0)}: L^2(D\times[t_1,t_2])\rightarrow L^2(0,\tau_w)\) is defined by \[\begin{align} (\mathcal{\tilde{L}}u)(\tau)=\int_{t_1}^{t_2}\int_D e^{-{\rm i} \tau (|x-y|+t)} u(y,t)\,{\rm d}y {\rm d}t,\qquad \tau\in (0, \tau_w), \end{align}\] for all \(u\in L^2(D\times[t_1,t_2])\), and the middle operator \(\mathcal{\tilde{T}}: L^2(D\times[t_1,t_2]) \rightarrow L^2(D\times[t_1,t_2])\) is again a multiplication operator defined by \[\begin{align} (\mathcal{\tilde{T}}\psi)(y,t):=\frac{1}{4\pi|x-y|}\int_{t_1}^{t_2} F(t-s)K(y,t)K(y,s)e^{{\rm i} k_0 (t-s)} e^{{\rm i} \tau_c (t+|x-y|)}\, {\rm d}s\, \psi(y,t). \nonumber \end{align}\] Note the adjoint of \(\tilde{L}:L^2(0,\tau_w)\rightarrow L^2(D\times[t_1,t_2])\) is defined by \[\begin{align} (\mathcal{\tilde{L}^*}\phi)(y,t)=\int_{t_1}^{t_2}\int_D e^{{\rm i} \gamma (|x-y|+t)} \phi (\gamma)\,{\rm d}\gamma,\qquad \phi\in L^2(0, \tau_w ). \end{align}\] Choose the test function: \[\tilde{\phi}_{y,\epsilon}^{(x)}(\tau)=\frac{1}{(t_2-t_1)|B_\epsilon(y)|} \int _{t_1}^{t_2} \int _{B_\epsilon(y)} \frac{e^{\rm i\tau(|x-z|+t)}}{4\pi|x-z|} \,{\rm d}z{\rm d}t\,\quad \tau \in (0, \tau_w).\] and introduce the indicator function \[\tilde{\mathcal{I}}^{(x)}(y):=\left[\sum_{n=1}^{n=\infty}\frac{|\langle \tilde{\phi}^{(x)}_{y,\epsilon}, \tilde{\psi}^{(x)}_n\rangle|^2_{L^2(0,\tau_w)}}{|\tilde{\lambda}_n^{(x)}|}\right]^{-1}, \quad y \in \mathbb{R}^3,\] where \((\tilde{\lambda}_n^{(x)}, \tilde{\psi}^{(x)}_n)\) is an eigensystem of the near-field operator \(\mathcal{N}^{(x)}_\#\). As a counterpart to \(\mathcal{I}\), we have the following result for the near-field case.

References↩︎