Loading [MathJax]/jax/output/SVG/jax.js

On a model of target detection in molecular communication networks

  • This paper is concerned with a target-detection model using bio-nanomachines in the human body that is actively being discussed in the field of molecular communication networks. Although the model was originally proposed as spatially one-dimensional, here we extend it to two dimensions and analyze it. After the mathematical formulation, we first verify the solvability of the stationary problem, and then the existence of a strong global-in-time solution of the non-stationary problem in Sobolev–Slobodetskiĭ space. We also show the non-negativeness of the non-stationary solution.

    Citation: Hirotada Honda. On a model of target detection in molecular communication networks[J]. Networks and Heterogeneous Media, 2019, 14(4): 633-657. doi: 10.3934/nhm.2019025

    Related Papers:

    [1] Hirotada Honda . On a model of target detection in molecular communication networks. Networks and Heterogeneous Media, 2019, 14(4): 633-657. doi: 10.3934/nhm.2019025
    [2] Eman S. Attia, Ashraf A. M. Khalaf, Fathi E. Abd El-Samie, Saied M. Abd El-atty, Konstantinos A. Lizos, Osama Alfarraj . Performance analysis of nanosystem based on cooperative relay for nanonetworks. Networks and Heterogeneous Media, 2023, 18(4): 1657-1677. doi: 10.3934/nhm.2023072
    [3] Matthieu Canaud, Lyudmila Mihaylova, Jacques Sau, Nour-Eddin El Faouzi . Probability hypothesis density filtering for real-time traffic state estimation and prediction. Networks and Heterogeneous Media, 2013, 8(3): 825-842. doi: 10.3934/nhm.2013.8.825
    [4] Manel Hmimida, Rushed Kanawati . Community detection in multiplex networks: A seed-centric approach. Networks and Heterogeneous Media, 2015, 10(1): 71-85. doi: 10.3934/nhm.2015.10.71
    [5] Yuntian Zhang, Xiaoliang Chen, Zexia Huang, Xianyong Li, Yajun Du . Managing consensus based on community classification in opinion dynamics. Networks and Heterogeneous Media, 2023, 18(2): 813-841. doi: 10.3934/nhm.2023035
    [6] Mary Luz Mouronte, Rosa María Benito . Structural properties of urban bus and subway networks of Madrid. Networks and Heterogeneous Media, 2012, 7(3): 415-428. doi: 10.3934/nhm.2012.7.415
    [7] Young-Pil Choi, Seung-Yeal Ha, Jeongho Kim . Propagation of regularity and finite-time collisions for the thermomechanical Cucker-Smale model with a singular communication. Networks and Heterogeneous Media, 2018, 13(3): 379-407. doi: 10.3934/nhm.2018017
    [8] Pau Erola, Albert Díaz-Guilera, Sergio Gómez, Alex Arenas . Modeling international crisis synchronization in the world trade web. Networks and Heterogeneous Media, 2012, 7(3): 385-397. doi: 10.3934/nhm.2012.7.385
    [9] Yacine Chitour, Frédéric Grognard, Georges Bastin . Equilibria and stability analysis of a branched metabolic network with feedback inhibition. Networks and Heterogeneous Media, 2006, 1(1): 219-239. doi: 10.3934/nhm.2006.1.219
    [10] Seung-Yeal Ha, Gyuyoung Hwang, Dohyun Kim . On the complete aggregation of the Wigner-Lohe model for identical potentials. Networks and Heterogeneous Media, 2022, 17(5): 665-686. doi: 10.3934/nhm.2022022
  • This paper is concerned with a target-detection model using bio-nanomachines in the human body that is actively being discussed in the field of molecular communication networks. Although the model was originally proposed as spatially one-dimensional, here we extend it to two dimensions and analyze it. After the mathematical formulation, we first verify the solvability of the stationary problem, and then the existence of a strong global-in-time solution of the non-stationary problem in Sobolev–Slobodetskiĭ space. We also show the non-negativeness of the non-stationary solution.



    Inverse problems for wave equations arise widely in non-destructive testing [1,2], underwater acoustics [3,4], antenna synthesis [5,6] and other application areas. Being two prototypical models in the field of inverse problems, the time-harmonic inverse obstacle scattering and inverse source problem have received significant investigations in the literature. For the nonlinear inverse obstacle scattering problem, the emanating source (incident wave) is given and the inhomogeneous scatterer serves as the target. Reciprocally, in the linear inverse source problem, the source is unknown while the passive scatterering inhomogeneity is vacant. Classical algorithms for inverse obstacle scattering problems involve the optimization methods [7,8,9], the singular source method [10] and the recursive linearization algorithm [11]. In addition, there are some sampling-type methods with qualitative features, for instance, the factorization method [12], the linear and direct sampling methods [13,14], the single-shot schemes [15,16] and the reverse time migration approach [17,18]. Recently, a variety of reconstruction methods have also been developed to deal with the inverse source problems. Bao et al. [19] applied a recursive algorithm to recover unknown sources of acoustic field with multi-frequency measurement data. Zhang et al. [20] used the direct sampling method to locate multipolar acoustic sources for the Helmholtz equation from near-field Cauchy data. Bousba et al. [21] utilized the direct sampling method to identify point sources for the Helmholtz equation from far-field data. The algebraic method has also been applied to the Helmholtz equation for determining monopolar sources [22].

    In this paper, we are concerned with an inverse problem to simultaneously reconstruct the acoustic obstacle and its excitation source points from the near-field measurements. This co-inversion problem can be mathematically formulated as an intrinsically coupled problem comprising the aforementioned inverse obstacle scattering and inverse source problem. In comparison with the vast studies on the conventional inverse obstacle/source problems, focuses on identification of both unknown source and scatterer are relatively rare. Over recent years, some interesting attempts have been made in dealing with such co-inversion problems. In [23], Liu et al. established sufficient conditions to recover both the sound speed of the medium and the source in thermo- and photo-acoustic tomography. The simultaneous recovery of an embedded obstacle and its surrounding medium by formally-determined scattering data was studied in [24]. In [25], the authors proved the uniqueness in simultaneously recovering the source functions and the obstacle for the wave equation. An adjoint state method for recovering both the source and attenuation in the attenuated X-ray transform was proposed in [26]. In addition, we also refer to [27,28,29] for some recent studies on the determination of both potential and source in the Schrödinger equations.

    In contrast to the standard problem of reconstructing either the source or obstacle solely, it is substantially more challenging to identify these unknown targets concurrently. It is worthwhile noting that only the total field is available and the contributions from inaccessible source and scatterer are superimposed in the total field measurements in the scenario of co-reconstruction. Hence, decoupling these incident and scattered components would inevitably suffer from inherent difficulties. As a result, the usual nonlinearity, ill-posedness as well as the additional lack of information constitute the mathematical and numerical challenges of our current study.

    In this work, we propose an optimization method to recover an obstacle and its illuminating point sources simultaneously from the measured total field. The salient features of our method can be summarized as follows. First, the proposed method is capable of simultaneously and quantitatively determining both the obstacle and source points since the reconstruction does not rely on any alternating iteration between the source and scatterer. Second, by the idea of the decomposition method and reformulating the inverse problem as a nonlinear optimization problem, the solution process for the direct problem is unnecessary. Third, by designing proper indicator functions, we provide an easy-to-implement and effective strategy for choosing the initial guesses for the optimization method, which avoids the a prior information of the unknown targets. Finally, the proposed method is insensitive to the noise-contaminated data and works well even with limited-aperture measurement data despite the severe ill-posedness.

    We now introduce the precise formulation of the co-inversion problem under consideration. Let $ D\subset \mathbb{R}^2 $ be an open and simply connected domain with $ C^2 $ boundary $ \partial D. $ Given the source point located at $ z\in\mathbb{R}^2\backslash\overline{D}, $ the incident field $ u^i $ due to $ z $ is given by

    $ ui(x;z)=Φ(x,z)=i4H(1)0(k|xz|),xR2(¯D{z}), $ (1.1)

    where $ k > 0 $ is the wave number, $ H_0^{(1)} $ is the Hankel function of the first kind of order zero. Then the forward scattering problem for the sound-soft obstacle $ D $ is to find the scattered field $ u^s\in H_{\rm loc}^1(\mathbb{R}^2\backslash\overline{D}) $ such that

    $ Δus+k2us=0in R2¯D, $ (1.2)
    $ u=0on D, $ (1.3)
    $ limr:=|x|r(usrikus)=0, $ (1.4)

    where $ u = u^i+u^s $ is the total field. The Sommerfeld radiation condition (1.4) holds uniformly with respect to all directions $ \hat{x} = x/|x|\in\mathbb{S}: = \{x\in\mathbb{R}^2:|x| = 1\} $, and (1.1) is also known as the fundamental solution to the Helmholtz equation (1.2). It is well known (see, e.g., [30]) that the direct problem (1.2)–(1.4) admits a unique solution $ u^s\in H^1_{\rm loc}\left(\mathbb{R}^2\backslash\overline{D}\right) $.

    In this paper, we assume that $ D\subset B_R: = \{x\in\mathbb{R}^2:|x| < R\} $ and $ B_R\backslash\overline{D} $ is connected. Let $ S: = \cup_{j = 1}^N\{z_j\}\subset B_R\backslash\overline{D} $ be a set of distinct source points with $ N\in\mathbb{N}_+ $ the number of source points. Denote by $ u^s(x; z_j) $ the scattered field corresponding to the incident field $ u^i(x; z_j) $ of the form (1.1), and we collect the total field $ u(x; z_j) = u^i(x; z_j)+u^s(x; z_j) $ on the measurement curve $ \Gamma_R: = \partial B_R = \{x\in\mathbb{R}^2: |x| = R\} $. Define the measured total field data by

    $ \mathbb{U}: = \{u(x; z):x\in \Gamma_R, z\in S\}, $

    then the inverse problem under consideration is to determine the obstacle-source pair $ (\partial D, S) $ from the measurements $ \mathbb{U} $, namely,

    $ U(D,S). $ (1.5)

    We refer to Figure 1 for an illustration of the geometry setting of the inverse obstacle and source problem (1.5). In addition, it deserves noting that, both the single source case ($ N = 1 $) and the multiple source case ($ N\ge 2 $) are admissible scenarios. The proposed algorithms in the subsequent sections can be applicable to both scenarios. We also want to point out that, the discrete total field measurement $ \bf{U} $ is a complex multistatic response matrix whose $ (i, j) $-th entry is the total field data due to the $ j $-th source point and the $ i $-th receiver. In other words, the total fields are collected separately for each source point one by one in our model and thus not directly related to the incident wave produced by the superposition of source points. Under the above configuration, the number $ N $ of the point sources is known as a priori, which is exactly the number of columns of the measurement matrix $ \bf{U} $.

    Figure 1.  Illustration of the co-inversion for imaging the obstacle and source points.

    The rest of this paper is organized as follows: In Section 2, an optimization method is proposed to reconstruct the source-obstacle pair and the corresponding convergence theory is established. In Section 3, an effective strategy based on the sampling-type method is provided to select the initial guess for the optimization method. Numerical experiments are given in Section 4 to verify the performance of the method. Finally, some conclusions are given in Section 5.

    Based on the decomposition method for inverse obstacle scattering problem [9], we propose a novel optimization method in this section by incorporating the source locations as new variables in the cost functional in order to reconstruct the obstacle-source pair $ (\partial D, S) $. The basic idea of the optimization method is to seek for the positions of source points and the boundary of obstacle as the locations where the Dirichlet boundary condition holds, namely, the total field $ u^i+u^s $ vanishes. We would like to emphasize that, since the admissible source locations and obstacles are radically distinct, the new optimization problem will thereby encounter nontrivial difficulties in the analysis and implementation.

    In this subsection, we present the optimization method for the inverse obstacle-source problem. According to the initial guess by the sampling results (see the next section for more details) about the obstacle-source pair $ (\partial D, S), $ we can place a closed $ C^2 $ curve $ \Lambda $ inside $ D $ such that $ k^2 $ is not a Dirichlet eigenvalue for the negative Laplacian in the interior of $ \Lambda $. Choosing some proper admissible set $ U\times V $ of the obstacle-source pair $ (\Gamma, z') $ such that $ \Gamma\in U, $ $ z' = (z'_1, z'_2, \cdots, z'_N)\in V, $ the scattered field $ u^s(x; z'_j) $ corresponding to $ z'_j $ can be represented by a single-layer potential

    $ us(x;zj)=ΛΦ(x,y)φj(y)ds(y), $ (2.1)

    with the unknown density $ \varphi_j\in L^2(\Lambda). $ Given the measured total field $ u(x; z_j), $ the density $ \varphi_j $ are supposed to satisfy the integral equation of the first kind

    $ (Sφj)(x)=u(x;zj)ui(x;zj),  xΓR, $ (2.2)

    where $ u^i(x; z'_j) $ is the incident wave due to source point $ z'_j $ and the compact integral operator $ \mathcal{S}:L^2(\Lambda)\to L^2(\Gamma_R) $ is defined by

    $ (\mathcal{S}\varphi)(x): = \int_\Lambda\Phi(x,y)\varphi(y)\mathrm{d}s(y),\quad x\in\Gamma_R. $

    The fact that the kernel of $ \mathcal{S} $ is analytic leads to the severe ill-posedness of (2.2), which motivates us to seek for a stable numerical solution by adopting the Tikhonov regularization strategy to replace (2.2) by

    $ αφαj+SSφαj=S(u(;zj)ui(;zj)),on  ΓR, $ (2.3)

    where $ \alpha > 0 $ is a regularization parameter and the adjoint operator $ \mathcal{S}^*:L^2(\Gamma_R)\to L^2(\Lambda) $ of $ \mathcal{S} $ is given by

    $ (Sg)(y):=ΓR¯Φ(x,y)g(x)ds(x),  yΛ, $

    where the overline signifies the complex conjugate.

    After obtaining the regularized density $ \varphi_j^\alpha $ by solving (2.3), we can calculate the approximation $ u^s_\alpha(x; z'_j) $ for the scattered field $ u^s(x; z'_j) $ by inserting the solution $ \varphi_j^\alpha $ of (2.3) into the single-layer potential representation (2.1). Then, the obstacle-source pair $ (\partial D, z) $ can be sought for by minimizing the defect

    $ Nj=1ui(x;zj)+usα(x;zj)L2(Γ) $ (2.4)

    over some suitable class $ U\times V $ of admissible obstacle-source pairs $ (\Gamma, z'). $ In (2.4), $ u^i(x; z'_j) $ is the incident field due to the source point $ z'_j. $

    Remark 2.1. We would like to point out that, as $ \alpha $ tends to zero, the solution $ \varphi_j^\alpha $ of (2.3) converges only if (2.2) is solvable. Following [7,Section 5.5] for the classic decomposition method, the solvability of (2.2) is related to regularity properties of the scattered field which, in general, cannot be known in advance for an unknown obstacle $ D $. In other words, $ u^s_\alpha(x; z'_j) $ may not converge to the $ u(x; z_j)-u^i(x; z'_j) $ since (2.2) may not be solvable.

    For the minimization problem, the admissible class $ U $ of $ \Gamma $ is chosen to be a compact set (with respect to the $ C^{1, \beta} $ norm, $ 0 < \beta < 1 $) of all star-like closed $ C^2 $ curves, described by

    $ Γ={r(ˆx)ˆx:ˆxS}, $ (2.5)

    where $ r\in C^2(\mathbb{S}) $ is the radial function satisfying the prior information

    $ 0 < a\le r(\hat{x})\le b,\ \ \hat{x}\in\mathbb{S}, $

    with given constants $ a $ and $ b. $ The admissible class $ V $ of $ z' = (z'_1, z'_2, \cdots, z'_N) $ is chosen to be a compact set (with respect to the Euclidean norm in $ \mathbb{R}^2 $) of all $ N $ ordered points such that $ z'_j\in B_R\backslash \overline{D_\Gamma}, j = 1, \cdots, N $, where $ D_\Gamma $ is the domain enclosed by $ \Gamma $ such that $ \Gamma = \partial D_\Gamma $.

    Remark 2.2. It deserves noting that, the purpose for choosing the admissible curve to be star-like is to provide convenience for the subsequent analysis and calculation, while the exact boundary curve $ \partial D $ does not necessarily require to be star-like. In fact, even if $ \partial D $ is not star-like, the obstacle-source pair may also be well reconstructed, which will be illustrated by the numerical experiments in Section 4.

    In this subsection, we will investigate the convergence properties of the optimization method with an analogous analysis to that of the decomposition methods [7,8].

    For later use, we define the acoustic single-layer operator $ \mathcal{S}_1:L^2(\Lambda)\to L^2(\Gamma) $ by

    $ (S1φ)(x):=ΛΦ(x,y)φ(y)ds(y),xΓ, $ (2.6)

    where $ \Gamma $ is a closed $ C^2 $ curve containing $ \Lambda $ in its interior. To start with, we state some properties of the operators $ \mathcal{S} $ and $ \mathcal{S}_1 $ (cf.[7]).

    Lemma 2.1. The single-layer operators $ \mathcal{S} $ and $ \mathcal{S}_1 $, defined by (2.1) and (2.6), respectively, are injective and have dense range provided that $ k^2 $ is not a Dirichlet eigenvalue for the negative Laplacian in the interior of $ \Lambda. $

    Combining the minimization of the Tikhonov functional for (2.4) and the defect minimization (2.2), we define the cost functional $ \mu:(L^2(\Lambda))^N\times U\times V\to \mathbb{R} $ by

    $ μ(φ,Γ,z;α)=Nj=1{u(x;zj)ui(x;zj)(Sφj)(x)2L2(ΓR)+αφj(x)2L2(Λ)+ui(x;zj)+(S1φj)(x)2L2(Γ)}, $ (2.7)

    where $ \varphi = (\varphi_1, \cdots, \varphi_N)\in(L^2(\Lambda))^N, $ $ z' = (z'_1, z'_2, \cdots, z'_N)\in V, $ $ \alpha > 0 $ is the regularization parameter in (2.3).

    Now, we turn to analyze the convergence properties of the optimization method. To this end, we first introduce the following definition of optimal obstacle-source pair.

    Definition 2.1. Given the measured total field $ u(x; z_j), j = 1, \cdots, N $ and a regularization parameter $ \alpha > 0, $ an obstacle-source pair $ (\Gamma^*, z^*)\in U\times V $ is called optimal if there exists $ \varphi^*\in (L^2(\Lambda))^N $ such that $ \varphi^* $ and $ (\Gamma^*, z^*) $ minimize the cost functional (2.7) simultaneously over all $ \varphi\in(L^2(\Lambda))^N $ and $ (\Gamma, z')\in U\times V, $ i.e.,

    $ μ(φ,Γ,z;α)=m(α), $ (2.8)

    where

    $ m(\alpha) = \inf\limits_{\substack{\varphi\in(L^2(\Lambda))^N,\\\Gamma\in U,\, z'\in V}}\mu(\varphi,\Gamma,z';\alpha). $

    In terms of Definition 2.1, the following theorem establishes the existence of the optimal solution.

    Theorem 2.1. For each $ \alpha > 0, $ there exists an optimal obstacle-source pair $ (\Gamma^*, z^*)\in U\times V. $

    Proof. Let $ (\varphi^n, \Gamma^n, z^n) $ be a minimizing sequence in $ \Big(L^2(\Lambda)\Big)^N\times U\times V $, i.e.,

    $ \lim\limits_{n\to\infty}\mu(\varphi^n,\Gamma^n, z^n;\alpha) = m(\alpha). $

    Since $ U, V $ are compact, we can assume that $ \Gamma^n\to\Gamma^*, $ $ z^n\to z^*, $ $ n\to\infty. $ In this paper, the convergence $ \Gamma^n\to\Gamma^* $ is understood in the sense that $ \|r^n-r^*\|_{C^{1, \beta}(\mathbb{S})}\to0, $ $ n\to\infty, $ with the functions $ r^n $ and $ r^* $ representing $ \Gamma^n $ and $ \Gamma^* $ via (2.5), respectively. From

    $ \alpha\|\varphi^n\|_{(L^2(\Lambda))^N}\le\mu(\varphi^n,\Gamma^n,z^n;\alpha)\to m(\alpha),\ \ n\to\infty, $

    with $ \alpha > 0, $ we conclude that the sequence $ \{\varphi^n\} $ is bounded, i.e., there exists some constant $ C $ such that $ \|\varphi^n\|_{(L^2(\Lambda))^N}\le C, \forall n. $ Further, we can assume that $ \varphi^n\rightharpoonup \varphi^*\in(L^2(\Lambda))^N, n\to\infty. $ Since $ \mathcal{S} $ is compact, it follows that $ \mathcal{S}\varphi_j^n\to\mathcal{S}\varphi_j^*, n\to\infty, j = 1, \cdots, N. $ The fact that the Hankel function $ H_0^{(1)} $ is continuous leads to the following convergence

    $ u^i(x; z_j^n) = \frac{\mathrm{i}}{4}H_0^{(1)}(k|x-z_j^n|)\to\frac{\mathrm{i}}{4}H_0^{(1)}(k|x-z_j^*|) = u^i(x; z_j^*),\ \ j = 1,\cdots, N. $

    Hence,

    $ \left\|u(x; z_j)-u^i(x; z_j^n)-(\mathcal{S}\varphi_j^n)(x)\right\|_{L^2(\Gamma_R)}\to \left\|u(x; z_j)-u^i(x; z_j^*)-(\mathcal{S}\varphi_j^*)(x)\right\|_{L^2(\Gamma_R)},\ \ n\to\infty. $

    By Taylor's formula, one can estimate that for all $ \hat{x}\in\mathbb{S}, $

    $ \left|\Phi(r^n(\hat{x})\hat{x},y)-\Phi(r^*(\hat{x})\hat{x},y)\right|\le K|r^n(\hat{x})-r^*(\hat{x})|, $

    with

    $ K = \sup\limits_{U\times\Lambda}\left|\nabla_x\Phi(x, y)\right|. $

    From the Schwarz inequality, it holds that

    $ \left|\int_\Lambda\left[\Phi(r^n(\hat{x})\hat{x},y)-\Phi(r^*(\hat{x})\hat{x},y)\right]\varphi_j^n(y)\mathrm{d}s(y)\right| \le CK|\Lambda||r_n(\hat{x})-r^*(\hat{x})|, $

    for all $ \hat{x}\in\mathbb{S}. $ Denote the dependence of $ \mathcal{S}_1:L^2(\Lambda)\to L^2(\Gamma^n) $ on $ n $ by writing $ \mathcal{S}_1^n, $ and the uniform convergence $ \mathcal{S}_1^n\varphi_j^n\to\mathcal{S}_1\varphi^*_j(n\to\infty $) holds. Therefore, we can deduce that

    $ \left\|u^i(x; z_j^n)+(\mathcal{S}_1^n\varphi_j^n)(x)\right\|^2_{L^2(\Gamma^n)}\to \left\|u^i(x; z_j^*)+(\mathcal{S}_1\varphi_j^*)(x)\right\|^2_{L^2(\Gamma^*)},\quad n\to\infty. $

    Further,

    $ αNj=1φnj2L2(Λ)m(α)Nj=1u(x;zj)ui(x;zj)(Sφj)(x)2L2(ΓR)Nj=1ui(x;zj)+(S1φj)(x)2L2(Γ)αNj=1φj2L2(Λ),  n. $

    Since $ \varphi_j^n\rightharpoonup \varphi_j^*, n\to\infty, j = 1, \cdots, N, $ it follows that

    $ \lim\limits_{n\to\infty}\sum\limits_{j = 1}^N\left\|\varphi_j^n-\varphi_j^*\right\|^2_{L^2(\Lambda)} = \lim\limits_{n\to\infty}\sum\limits_{j = 1}^N\left\|\varphi_j^n\right\|^2_{L^2(\Lambda)} -\sum\limits_{j = 1}^N\left\|\varphi_j^*\right\|^2_{L^2(\Lambda)}\le0. $

    Thus, we have the norm convergence $ \left\|\varphi_j^n-\varphi_j^*\right\|\to0, n\to\infty $ for $ j = 1, \cdots, N. $ Finally, it holds that

    $ \mu\Big(\varphi^*,\Gamma^*,z^*;\alpha\Big) = \lim\limits_{n\to\infty}\mu\Big(\varphi^n,\Gamma^n,z^n;\alpha\Big) = m(\alpha), $

    and the proof is complete.

    Before heading for the convergence analysis of the optimization method, we shall state the convergence property of the cost functional first.

    Theorem 2.2. Let $ z_j, $ $ j = 1, \cdots, N $ be the exact source points, $ u(x; z_j) $ be the exact total fields corresponding to the sound-soft obstacle $ D $ such that $ \partial D $ belongs to $ U. $ Then the cost functional $ m(\alpha) $ is convergent in the sense that

    $ limα0m(α)=0. $ (2.9)

    Proof. By Lemma 2.1, given $ \varepsilon_j > 0, $ there exists $ \varphi_j\in L^2(\Lambda) $ such that

    $ \left\|u^i(x;z_j)+(\mathcal{S}_1\varphi_j)(x)\right\|_{L^2(\partial D)}\le\varepsilon_j,\ \ j = 1,\cdots,N. $

    Since the scattered field of a radiating solution to the Helmholtz equation depends continuously on the boundary (cf. [7]), there exists some constant $ C_1 $ such that

    $ \left\|(\mathcal{S}\varphi_j)(x)-u^s(x;z_j)\right\|_{L^2(\Gamma_R)}\le C_1\left\|(\mathcal{S}_1\varphi_j)(x)-u^s(x;z_j)\right\|_{L^2(\partial D)}, $

    i.e.,

    $ u(x;zj)ui(x;zj)(Sφj)(x)L2(ΓR)C1u(x;zj)ui(x;zj)(S1φj)(x)L2(D)=C1ui(x;zj)+(S1φj)(x)L2(D)C1εj, $

    where we have used the fact that $ u(x; z_j) = 0, x\in\partial D $. The continuity of the Hankel function $ H_0^{(1)} $ implies that there exists $ \delta_j > 0 $ such that for $ z'_j $ satisfying $ |z_j-z'_j|\le\delta_j $, it holds that

    $ \left\|u^i(x; z_j)-u^i(x; z'_j)\right\|_{L^2(\Gamma_R)}\le \varepsilon_j. $

    Based on the above arguments, it holds that

    $ u(x;zj)ui(x;zj)(Sφj)(x)L2(ΓR)u(x;zj)ui(x;zj)(Sφj)(x)L2(ΓR)+ui(x;zj)ui(x;zj)L2(ΓR)(C1+1)εj:=Cεj. $

    Further, we deduce that

    $ \mu(\varphi,\partial D,z;\alpha)\le(1+C^2)\sum\limits_{j = 1}^N\varepsilon_j^2+\alpha\sum\limits_{j = 1}^N\left\|\varphi_j\right\|_{L^2(\Lambda)}^2\to (1+C^2)\sum\limits_{j = 1}^N\varepsilon_j^2,\ \ \alpha\to0. $

    The theorem is accomplished due to the arbitrariness of each $ \varepsilon_j, j = 1, 2, \cdots, N. $

    We are now in a position to state the main convergence result.

    Theorem 2.3. Let $ \{\alpha_n\} $ be a null sequence, $ (\Gamma^n, z^n) $ be the corresponding sequence of optimal obstacle-source pairs for the regularization parameter $ \alpha_n. $ Then there exists a convergent subsequence of $ \left\{(\Gamma^n, z^n)\right\}. $ Assume that $ u $ is the exact total field of a domain $ D $ such that $ \partial D $ is contained in $ U. $ Then every limit point $ (\Gamma^*, z^*) $ of $ (\Gamma^n, z^n) $ represents an optimal obstacle-source pair. In other words, if the incident wave due to the source point $ z^* $ is scattered by the sound-soft obstacle $ D^* $ bounded by $ \Gamma^* $, then the total field vanishes on $ \Gamma^* $.

    Proof. The existence of convergent subsequences of $ (\Gamma^n, z^n) $ follows from the compactness of $ U\times V $. Let $ (\Gamma^*, z^*) $ be the limit point. Without loss of generality, we assume that $ (\Gamma^n, z^n)\to(\Gamma^*, z^*), n\to\infty, $ i.e., $ \Gamma^n\to\Gamma^*, $ $ z^n\to z^*, $ $ n\to\infty. $ Let $ u_j^* $ be the solution to the direct scattering problem with the incident wave emitted by the source point $ z_j^* $ for the obstacle with boundary $ \Gamma^*, $ i.e., the boundary condition reads

    $ uj()+ui(;zj)=0,on Γ. $ (2.10)

    Since $ (\Gamma^n, z^n) $ is optimal for the parameter $ \alpha_n $, there exists $ \varphi^n\in (L^2(\Lambda))^N $ such that

    $ \mu(\varphi^n,\Gamma^n, z^n;\alpha_n) = m(\alpha_n), $

    for $ n = 1, 2, \cdots $. From Theorem 2.2, we conclude that the boundary data satisfies

    $ ui(x;znj)+(Sn1φnj)(x)L2(Γn)0,n. $ (2.11)

    From the condition that $ \Gamma^n\to\Gamma^*, z_j^n\to z_j^*, n\to\infty, $ we obtain that $ u_j^s = u_j^* $ with $ u_j^s = \lim\limits_{n\to\infty}\mathcal{S}_1^n\varphi_j^n, $ and the total field $ u^s_j(\cdot)+u^i(\cdot; z^*_j) $ must vanish on $ \Gamma^* $ as a result of (2.10).

    Remark 2.3. As far as we know, even for the standard decomposition method for the inverse obstacle problem, a uniqueness result for the optimal solution is not available. It is commented in [7] that

    "Since we do not have uniqueness either for the inverse scattering problem or for the optimization problem, in the above convergence analysis we cannot expect more than convergent subsequences."

    For the more complicated co-inversion problem considered here, the uniqueness for the optimal solution in Theorem 2.1 is still open. Therefore, the above comments from [7] are also applicable to our convergence analysis.

    Based on the previous discussions, the rudiments of our optimization procedure is summarized in Algorithm 1.

    Algorithm 1: Optimization method for the co-inversion problem.
    Step 1 Given the total field data $ u(x; z_j)|_{\Gamma_R}, j = 1, \cdots, N $, choose some suitable class $ U\times V $ of admissible obstacle-source pairs $ (\partial D, S) $;
    Step 2 Select an initial curve for the boundary $ \partial D $ and an initial guess for the $ N $ ordered points $ z' = (z'_1, z'_2, \cdots, z'_N) $;
    Step 3 Place a suitable auxiliary curve $ \Lambda $ inside $ D $;
    Step 4 Determine the obstacle-source pair $ (\partial D, S) $ by minimizing (2.4), which is formulated by the following implementations for $ j = 1, \cdots, N $:
    (a) Subtract the incident field $ u^i(x; z'_j) $ due to the source point $ z'_j $ from the measured total field $ u(x; z_j) $, where the incident field $ u^i(x; z'_j) $ is given by (1.1);
    (b) Solve the regularized Eq (2.3) with an appropriate regularization parameter $ \alpha $ to find the density $ \varphi_j^\alpha $;
    (c) Represent the approximate scattered field $ u_\alpha^s(x; z'_j) $ in terms of the single-layer representation (2.1) with density $ \varphi_j^\alpha $.
    Step 5 When the error tolerance is fulfilled, the minimizer $ (\Gamma^*, z^*) $ of (2.4) can be viewed as the desired reconstruction.

     | Show Table
    DownLoad: CSV

    Remark 2.4. For simplicity of implementation, the regularization parameter $ \alpha $ is chosen by the trial-and-error strategy and the error tolerance is independent of $ \alpha $ in this paper. To improve the performance of the algorithm, $ \alpha $ could be chosen by, for instance, the Morozov's discrepancy principle according to the noise level. Meanwhile, the error tolerance should dependent on the value of $ \alpha $, and may not be reached if $ \alpha $ is too large.

    At the end of this section, we would like to remark that, the nonlinear optimization problems generally suffer from the notorious difficulty of local minimums. Specifically, performance of any iterative solver for the optimization problem (2.8) would be unavoidably affected by the initial guess $ (\Gamma^0, z^0) $. Therefore, choosing an appropriate initial value is highly crucial for the success of our optimization method. Nevertheless, the a priori information of the rough $ (\Gamma^0, z^0) $ is not straightforwardly available. And it is thus vital to develop a fast and cheap strategy for resolving this crux, which will be the main concern of the next section.

    In the previous section, we develop an optimization method to reconstruct the obstacle-source pair $ (\partial D, S) $ from the total field $ \mathbb{U} $ and it has been pointed out that the initial guess $ (\Gamma^0, z^0) $ should play an important role in the quality of reconstruction. Noticing that the direct sampling methods (DSM) are independent of any a priori information of the unknown objects and easy to implement, we wish to utilize the sampling-type methods to seek the initial guess for the source $ S $ and the obstacle $ D $, respectively.

    Motivated by the direct sampling schemes for inverse obstacle/source scattering problems, we first propose the total field version of indicator functions to reconstruct the source points by knowledge of $ \mathbb{U} $. Let $ \Omega_1\subset B_R $ be a bounded sampling domain such that $ D\cup S\subset\Omega_1. $ For any sampling point $ y\in\Omega_1 $, we define $ N $ indicator functions as follows:

    $ Ij(y)=|ΓR¯u(x;zj)Φ(x,y)ds(x)|,  j=1,,N. $ (3.1)

    Note that the indicators rely on the inner product of the fundamental solution with respect to the total field data rather than the usual incident or scattered data. We expect to locate the peak of $ I_j(y) $ for each $ j $ and take the corresponding maximum point $ \tilde{z}_j $ as an approximation of the exact source $ z_j $.

    Before analyzing the characteristics of $ I_j(y), $ we first claim two lemmas which will play a key role in the subsequent analysis.

    Lemma 3.1. [17,Lemma 3.1] Let $ g \in H^{1 / 2}\left(\partial{D}\right). $ Then the scattering problem

    $ Δw+k2w=0in  R2¯D, $ (3.2)
    $ w=gon  D, $ (3.3)
    $ r(wrikw)0as  r, $ (3.4)

    admits a unique solution $ w \in H_{\mathrm{loc}}^{1}\left(\mathbb{R}^2 \backslash \overline{D}\right) $. Moreover, there exists a constant $ C > 0 $ such that

    $ \left\|\frac{\partial w }{\partial \nu}\right\|_{H^{-1 / 2}\left(\partial{D}\right)} \leq C\|g\|_{H^{1 / 2}\left(\partial{D}\right)}, $

    with $ \nu $ the unit outer normal to the boundary $ \partial D $.

    For any $ g \in H^{1 / 2}\left(\partial{D}\right) $, the Dirichlet-to-Neumann mapping $ \mathbb{T}: H^{1 / 2}\left(\partial{D}\right) \rightarrow H^{-1 / 2}\left(\partial{D}\right) $ for the scattering problems (3.2)–(3.4) is defined as

    $ \mathbb{T}(g) = \frac{\partial w}{\partial \nu}. $

    By Lemma 3.1, $ \mathbb{T} $ is a bounded linear operator and we will denote $ \|\mathbb{T}\| $ its operator norm in this paper.

    Analogous to [17,Lemma 3.2], one can immediately obtain the following estimate concerning the scattered field.

    Lemma 3.2. Let $ L: = \min_{1\le j\le N}{\rm dist}\left(z_j, \overline{D}\right) $. Then for each $ j = 1, \cdots, N $, the following estimate holds

    $ |u^s(x; z_j)|\le C_{js}(1+\|\mathbb{T}\|)k^{-1}(RL)^{-1/2}, $

    where $ x\in\Gamma_R $ is the receiver and the constant $ C_{js} $ is dependent of $ k, R $ and $ L $.

    The following theorem sheds light on the indicating properties of $ I_j(y), j = 1, \cdots, N $.

    Theorem 3.1. Let $ I_j(y), j = 1, \cdots, N, $ be defined by (3.1) and assume that

    $ Lj:=dist(zj,ΓR)>supyΩ1|yzj|. $ (3.5)

    Then there exists a constant $ C_j > 0 $ such that for each $ y\in \Omega_1, \ \ y\ne z_j, $

    $ Ij(y)Cjk1|yzj|1/2,j=1,,N. $ (3.6)

    Proof. We first recall the following estimate for the Hankel function [17,(3.8)]:

    $ |H(1)0(t)|2πt,  t>0. $ (3.7)

    Then from the representation of incident field $ u^i(x; z_j), $ it can be easily derived that

    $ ΓR|¯ui(x;zj)Φ(x,y)|ds(x)=ΓR|¯i4H(1)0(k|xzj|)i4H(1)0(k|xy|)|ds(x)116ΓR2πk|xzj|2πk|xy|ds(x)=18kπΓR1|xzj||xy|ds(x). $ (3.8)

    Further, by Lemma 3.2, we deduce that there is a positive constant $ C_{js} $ such that

    $ ΓR|¯us(x;zj)Φ(x,y)|ds(x)=ΓR|us(x;zj)||i4H(1)0(k|xy|)|ds(x)Cjs24k3/2(πRL)1/2ΓR1+T|xy|ds(x). $ (3.9)

    Combining (3.8) and (3.9), together with the aid of $ |x-y|\ge\left||x-z_j|-|y-z_j|\right| $, we have

    $ Ij(y)=|ΓR¯ui(x;zj)Φ(x,y)ds(x)+ΓR¯us(x;zj)Φ(x,y)ds(x)|ΓR|¯ui(x;zj)Φ(x,y)|ds(x)+ΓR|¯us(x;zj)Φ(x,y)|ds(x)Cj1(ΓRk1(|xzj||xy|)1/2ds(x)+ΓRk3/2(1+T)|xy|1/2ds(x))Cj1(ΓRk1(|xzj|||xzj||yzj||)1/2ds(x)+ΓRk3/2(1+T)||xzj||yzj||1/2ds(x)), $

    where

    $ C_{j1} = \max\left\{\frac{1}{8\pi}, \frac{C_{js}\sqrt{2}}{4(\pi R L)^{1/2}}\right\}. $

    Assumption (3.5) implies that there exists a positive constant $ C_{j2} $ such that $ L_j\ge (1+C_{j2})\sup_{y\in\Omega_1}|y-z_j| $. Then

    $ Ij(y)2πCj1RkLjCj2|yzj|+2Cj1π(1+T)Rk3/2Cj2|yzj|Cjk1|yzj|1/2, $

    where

    $ C_j = \frac{2\pi R C_{j1}}{\sqrt{C_{j2}}}\left(\frac{1}{\sqrt{L_j}}+\frac{1+\|\mathbb{T}\|}{\sqrt{k}}\right), $

    and this completes the proof.

    In virtue of the above theorem, each function $ I_j(y) $ should decay as the sampling point $ y $ recedes from the corresponding source point $ z_j $. Hence, the initial detection can be conveniently chosen by locating the point $ \tilde{z}_j $ where $ I_j(y) $ reaches its apex. Then in the optimization process, the exact source location $ z_j $ could be effortlessly and accurately found in the vicinity of $ \tilde{z}_j $. This indicating behavior will be numerically verified by the experiments in the next section.

    The indicator function proposed in the previous subsection has the capability of locating the approximate source points, but it usually fails to image the shape of obstacle. Thus a further sampling step for identifying the rough obstacle should be taken into account. By integrating the approximate source points into the reverse time migration (RTM) approach [18], we propose an approximate reverse time migration method to reconstruct the obstacle in this subsection.

    To reconstruct the obstacle, information about the scattered field is indispensable in the classical RTM. However, the scattered field can not be measured directly due to the existence of unknown source points. So we are going to approximate the scattered field by subtracting the incident field $ u^i(x; \tilde{z}_j) $ due to the numerical source point $ \tilde{z}_j $ from the measured total field $ u(x; z_j), $ i.e.,

    $ usj(x)=u(x;zj)ui(x;˜zj),xΓR, $ (3.10)

    and develop an imaging functional $ I_D(y) $ using the approximate scattered field $ u_j^s(x) $. To this end, confining the imaging domain to $ \Omega_2\subset\Omega_1 $ such that $ \overline{D}\subset \Omega_2, $ $ S\subset\Omega_1\backslash\overline{\Omega}_2, $ and we define the imaging function for $ y\in\Omega_2 $ by

    $ ID(y):=k2Im(Nj=1Φ(y,zj)ΓR¯usj(x)Φ(x,y)ds(x)). $ (3.11)

    Note that, in (3.11), instead of the exact scattered field $ u^s(x; z_j) $ corresponding to the $ j $-th source point $ z_j, $ we use an approximate scattered field $ u^s_j(x) $ for $ z_j $. That is why this method is called the approximate reverse time migration method. Once the source points are retrieved, the indicating behavior of the approximate RTM immediately follows the original RTM [17,18]. The imaging function $ I_D(y) $ provides a qualitative information of the obstacle $ D $, which will be used as a useful clue to discover the obstacle in the optimization method. We emphasize that, if this sampling procedure for $ D $ is not implemented, then it would be difficult to find a good initial guess for $ D $ or a suitable internal curve $ \Lambda $. In this case, such an apriori information about the obstacle is usually necessary to facilitate the realization of Algorithm 1.

    We end this section with a brief description of Algorithm 2, which can be used to determine the initial guess for the optimization method proposed in Section 2.

    Algorithm 2: Determine the initial guess by the direct sampling methods.
    Step 1 For each $ j = 1, \cdots, N $, collect the total field data $ u(x; z_j)|_{\Gamma_R} $;
    Step 2 Choose a proper sampling domain $ \Omega_1\subset B_R $ such that $ D\cup S\subset \Omega_1 $ and generate the sampling grid $ \mathcal{T}_1 $ over $ \Omega_1 $;
    Step 3 Locate the maximum of $ I_j(y) $ for each $ j $ and take the corresponding $ \tilde{z}_j $ as an approximation of $ z_j $. The locations $ \{\tilde{z}_j\}_{j = 1}^N $ could be found one by one in this way and will be reused as the initial guess for the source points in the optimization method;
    Step 4 For $ j = 1, \ldots, N $, compute the approximate scattered field $ u^s_j $ by (3.10);
    Step 5 Choose a proper imaging domain $ \Omega_2\subset\Omega_1 $ such that $ \overline{D}\subset\Omega_2 $, $ S\subset \Omega_1\backslash\overline{\Omega}_2 $ and generate the sampling grid $ \mathcal{T}_2 $ over $ \Omega_2 $;
    Step 6 Evaluate $ I_D(y) $ for the imaging point $ y\in\mathcal{T}_2 $ by (3.11) and accordingly select an approximate initial guess for the boundary curve in the optimization method.

     | Show Table
    DownLoad: CSV

    In this section, we will present several numerical experiments to demonstrate the performance of the proposed method. In our experiments, synthetic total field data are generated by reformulating the direct problems (1.2)–(1.4) as a boundary integral equation and the integral equation is solved by the Nyström method [7]. Here we use a numerical quadrature rule with $ 64 $ equidistant grid points on $ [0, 2\pi] $. The receivers are set to be $ x_r = 4(\cos\theta_r, \sin\theta_r), $ $ \theta_r = {r\theta}/{N_R}, $ $ r = 1, \cdots, N_R $ with $ \theta\in (0, 2\pi] $ the observation aperture. For the full-aperture case, $ \theta = 2\pi, $ $ N_R = 120 $. For the limited-aperture case, we consider two cases with $ \theta = 3\pi/2, $ $ N_R = 90 $ and $ \theta = \pi, $ $ N_R = 60 $, respectively. The forward solver provides us with the total field data

    $ u(x_r; z_j),\quad r = 1,\cdots,N_R,\quad j = 1,\cdots,N. $

    To test the stability of the numerical method, random noise of different levels is added to the total field data $ u(x; z_j) $ to produce the noisy total field data

    $ u^\varepsilon: = u+\varepsilon r_1|u|\mathrm{e}^{\mathrm{i}\pi r_2}, $

    where $ r_1, $ $ r_2 $ are two uniformly distributed random numbers ranging from $ -1 $ to $ 1 $ and $ \varepsilon > 0 $ is the noise level. We use $ \varepsilon = 10\% $ in Examples 1–3 and $ \varepsilon = 5\% $ in Example 4.

    In all the subsequent examples, the sampling domains $ \Omega_1 $ are chosen as squares centered at the origin with $ \mathcal{T}_1 $ being a $ 200\times200 $ uniform grid. The imaging domains $ \Omega_2 $ are individualized case by case, according to the sampling results for the point sources. The imaging points $ \mathcal{T}_2 $ are selected as $ 100\times100 $ uniform grids distributed over $ \Omega_2 $. For the optimization method, we choose the auxiliary curve $ \Lambda $ to be a circle centered at the origin with radius $ r_\Lambda = 0.6. $ The integral over the unit circle $ \mathbb{S} $ is numerically approximated by the trapezoidal rule with $ 64 $ grid points. To solve the nonlinear least-squares problem, the Levenberg-Marquardt algorithm is adopted with both functional value stopping criterion and successive iterate stopping criterion chosen to be $ 10^{-6} $ for the full-aperture case and $ 10^{-4} $ for the limited-aperture problems. The starting curve $ \Gamma_0 $ for the Levenberg-Marquardt algorithm is a circle centered at the origin with radius $ r_{\Gamma_0} = \left|\arg\max_{y\in\Omega_2}I_D(y)\right| $. For $ j = 1, \cdots, N $, the starting guesses for the source points $ \{z_j^0\}_{j = 1}^N $ are chosen to be $ \{\tilde{z}_j\}_{j = 1}^N $. The regularization parameter $ \alpha $ is selected by trial and error. The ansatz function $ r(t) $ is chosen in the $ (2M+1) $-dimensional subspace $ U_M\subset U $ and can be expanded by trigonometric polynomials of degree less than or equal to $ M, $ i.e.,

    $ rM(t)=rΛ+a0+Mm=1(amcosmt+bmsinmt), $ (4.1)

    with some $ M\in\mathbb{N}_+ $. Unless otherwise specified, $ M = 8 $ is used in what follows.

    To evaluate the reconstruction quantitatively, we introduce an equidistant set of knots on $ [0, 2\pi] $ by $ t_i: = 2\pi i/N_1, i = 0, 1, \cdots, N_1-1, $ with $ N_1 $ the number of knots chosen to be 256 in our numerical experiments. The discrete relative $ L^2 $ error of the boundary curve is defined by

    $ ED:=(N1i=1|rM(τi)r(τi)|2)1/2(N1i=1|r(τi)|2)1/2, $ (4.2)

    with $ r^* $ standing for the exact boundary curve and $ \tau_i = \tau(t_i) $ being $ t_i $ or some transformation of $ t_i $ depending on whether the exact boundary is star-like or not, respectively.

    In the following figures regarding geometrical settings of the problem, the red and blue points denote respectively the exact source locations and receivers, the green solid line denotes the boundary of the obstacle, the black and purple dashed lines denote the boundary of the sampling domains $ \Omega_1 $ and $ \Omega_2, $ respectively. In the figures about the results of the numerical experiments, the small black '+' markers are the exact source points, the small red circles stand for the reconstructed source points, the green solid line denotes the exact boundary of the obstacle, the blue dashed line denotes the reconstructed boundary. In addition, for ease of comparison, the indicator functions $ I_j(y), $ $ j = 1, \cdots, N, $ and the imaging function $ I_D(y) $ are normalized so that the maximum value of each function is $ 1. $

    Example 1. In the first example, we consider the reconstruction of a starfish-shaped obstacle and the influence of source quantities. The exact boundary of the starfish-shaped obstacle is parameterized by

    $ x_S(t) = (1+0.2\cos 5t)(\cos t,\sin t),\quad0\le t\le2\pi. $

    We first consider the case where the five source points are respectively located at

    $ S_1 = \{(2,2),(-0.5,1.8),(-1.6,-1),(0,-2.2),(2,-1.2)\}. $

    The sampling domains are chosen to be $ \Omega_1 = [-2.5, 2.5]\times[-2.5, 2.5] $, $ \Omega_2 = [-1.5, 1.5]\times[-1.5, 1.5] $. The regularization parameter and the wave number are set to be $ \alpha = 10^{-8} $ and $ k = 5 $, respectively.

    Figure 2(a) exhibits the geometry setting. Figure 2(b) shows the reconstructed locations of $ S_1 $ by the DSM and we can find that the reconstructed locations of the source points are very close to the exact locations. Figure 2(c) shows the image of the imaging function $ I_D(y) $ over the sampling domain $ \Omega_2. $ We find that in Figure 2(c), the shape of the obstacle can be roughly captured.In Figure 2(d), reconstruction of the optimization method is depicted, where the initial guess for the source points are chosen to be the numerical source points obtained by the DSM and the initial curve is represented by the black dashed circle in Figure 2(c). Further, we compute the error of boundary curve by (4.2) with $ \tau_i = t_i, i = 1, \cdots, N $ and in this case $ E_D = 5.73\% $. In Table 1, we list the exact and reconstructed locations of source points $ S_1 $ by the DSM and the optimization method, respectively. According to Table 1 and Figure 2, we want to point out that both the DSM and optimization method could produce satisfactory reconstructions of the source points. A promising advantage of the optimization method is that the obstacle can also be reconstructed quantitatively, i.e., a significantly more sharper profile can be identified by the optimization method.

    Figure 2.  Geometry setting and the reconstruction of the starfish and $ S_1 $. (a) Problem geometry; (b) exact and reconstructed locations of the source points by the DSM; (c) image of $ I_D(y) $; (d) reconstruction by the optimization method.
    Table 1.  Reconstruction of the source points $ S_1 $.
    Exact Location Reconstructed Location
    DSM Optimization
    Point 1 $ (2, 2) $ $ (1.98, 1.98) $ $ (2.01, 1.98) $
    Point 2 $ (-0.5, 1.8) $ $ (-0.52, 1.83) $ $ (-0.50, 1.83) $
    Point 3 $ (-1.6, -1) $ $ (-1.58, -0.98) $ $ (-1.56, -0.98) $
    Point 4 $ (0, -2.2) $ $ (0.02, -2.19) $ $ (-0.01, -2.19) $
    Point 5 $ (2, -1.2) $ $ (2.03, -1.23) $ $ (1.99, -1.23) $

     | Show Table
    DownLoad: CSV

    Next we consider the reconstruction from a different set of source points, which is given by $ S_2 = \{(1.3, 0), (1.35, 0.1), (0, 2), (0, -2), (-2, 0)\} $. In this case, the sampling domains are adapted to be $ \Omega_1 = [-2.2, 2.2]\times[-2.2, 2.2], $ $ \Omega_2 = [-1.2, 1.2]\times[-1.2, 1.2] $. The other parameters remain unchanged. We refer to Figure 3(a) for the geometry setting.

    Figure 3.  Geometry setting and the reconstruction of the starfish and $ S_2 $. (a) Problem geometry; (b) image of $ I_1(y) $ (white "$ + $": exact location of the first source point); (c) exact and reconstructed locations of the source points with $ I_j(y), j = 1, \cdots5 $; (d) image of $ I_D(y) $; (e) reconstruction by the optimization method; (f) optimization result in a local domain.

    Figure 3(b) shows the indicator $ I_1(y) $ over the sampling domain $ \Omega_1 $. One can observe that the maximizer of $ I_1(y) $ coincides well with the exact source location. Figure 3(c) shows that all the reconstructed source locations match well with the exact positions. In Figure 3(d), the imaging function $ I_D(y) $ is plotted. From Figure 3(d), we find that the maximizer of $ I_D(y) $ incurs on the boundary of the obstacle.

    Compared with Figure 2(c), it seems that the current result deteriorates slightly. A reason accounting for this fact is the specific distribution of the source points where two points are almost adjacent to each other.In Figure 3(e), we show the reconstruction of the starfish by the optimization method. Figure 3(e) demonstrates that the quality of reconstruction for the left side of the obstacle is better than that of the right side. Physically speaking, this is due to the fact that the right side of the obstacle is affected by the nearby source points 1 and 2. Moreover, the adjacent source points 1 and 2 could be individually identified as long as the sampling grid is sufficiently fine. To see the local reconstruction more clearly, we intercept a part $ [1, 1.4]\times[-0.1, 0.2] $ of Figure 3(e) and it is zoomed in as Figure 3(f).

    In Table 2, we list the exact and reconstructed source locations. To evaluate the reconstruction of the obstacle quantitatively, we compute the errors between the exact boundary curve and the reconstructed boundary curve as (4.2) with $ \tau_i = t_i $. In this case, the error of the boundary curve is $ E_D = 7.90\% $.

    Table 2.  Reconstruction of the source points $ S_2 $.
    Exact source Reconstruction
    Point 1 $ (1.3, 0) $ $ (1.25, 0.02) $
    Point 2 $ (1.35, 0.1) $ $ (1.30, 0.07) $
    Point 3 $ (0, 2) $ $ (-0.00, 1.98) $
    Point 4 $ (0, -2) $ $ (-0.00, -1.98) $
    Point 5 $ (-2, 0) $ $ (-1.99, 0.02) $

     | Show Table
    DownLoad: CSV

    Comparing the above two cases, it can be found that our method could simultaneously produce a rather accurate reconstruction of the obstacle and sources, no matter whether the source points are relatively sparse or dense, or even quite close to the obstacle.

    In Figure 4, we consider the reconstruction with $ N( = 2, 4, 6, 8) $ source points. The wave number and the regularization parameter is chosen to be $ k = 8 $ and $ \alpha = 10^{-6} $, respectively.The sampling domains are chosen to be $ \Omega_1 = [-2.5, 2.5]\times[-2.5, 2.5], $ $ \Omega_2 = [-1.5, 1.5]\times[-1.5, 1.5]. $ Further, we take $ k = 5, 8, $ respectively, and list the relative $ L^2 $ errors for the reconstruction of the starfish-shaped obstacle in Table 3. Table 3 shows that the error of reconstruction decreases as the number of source points increases.

    Figure 4.  Reconstruction of the starfish and different numbers of source points. The 1st column: problem geometry; the 2nd column: exact and reconstructed locations of the source points with DSM; the 3rd column: images of $ I_D(y) $; the 4th column: reconstruction by the optimization method.
    Table 3.  Relative $ L^2 $ errors for reconstructions of the starfish.
    $ N $ 2 4 6 8
    $ k=5 $ $ 8.60\% $ $ 7.45\% $ $ 8.01\% $ $ 5.84\% $
    $ k=8 $ $ 15.46\% $ $ 6.63\% $ $ 5.93\% $ $ 4.68\% $

     | Show Table
    DownLoad: CSV

    Example 2. In this example, we aim to investigate the accuracy of reconstruction by considering a circular obstacle as well as different number of source points. The parametric boundary of the obstacle is given by $ x_C(t) = (\cos t, \sin t), \ 0\leq t < 2\pi, $ and $ \alpha = 10^{-8}, $ $ \Omega_1 = [-2.5, 2.5]\times[-2.5, 2.5], $ $ \Omega_2 = [-1.2, 1.2]\times[-1.2, 1.2]. $

    To analyze the reconstruction quantitatively, the relative errors of the boundary curves are computed as (4.2) with $ \tau_i = t_i $. Table 4 shows the reconstructions in different cases by varying the wave numbers and the number of source points. In particular, a suitable choice of the parameters (for instance, $ k = 8 $ and $ N = 6 $) would produce a stunningly accurate reconstruction. Figure 5 exhibits the reconstruction of the circle and $ N( = 1, 2, 4, 6, 8) $ source points with wave number $ k = 5 $. In addition, it deserves noting that, even if we take $ N = 1, $ both the source and the obstacle can be reconstructed, though the error is slightly large (see Figure 5(a) and the second column in Table 4).

    Table 4.  Relative $ L^2 $ errors for reconstructions of the circle with different numbers of source points.
    $ N $ 1 2 4 6 8
    $ k=5 $ $ 14.28\% $ $ 7.32\% $ $ 2.30\% $ $ 2.47\% $ $ 1.24\% $
    $ k=8 $ $ 11.65\% $ $ 5.44\% $ $ 1.38\% $ $ 0.72\% $ $ 0.76\% $

     | Show Table
    DownLoad: CSV
    Figure 5.  Reconstructions of the circle and $ N $ source points. (a) $ N = 1 $; (b) $ N = 2 $; (c) $ N = 4 $; (d) $ N = 6 $ (e) $ N = 8 $.

    Further, we test the effectiveness of the DSM proposed in Section 3 for choosing the initial guess. The parameters are chosen to be $ \alpha = 10^{-9}, k = 8, N = 8 $. The exact source points are given by

    $ S_3 = \{(2,1),(1.8,1.4),(0,1.6),(-1.6,1.8),(-2,0),(-2,-2),(0,-2.1),(1.8,-1.9)\}. $

    For comparison purpose, we consider the following two sets of artificially chosen initial guess

    $ S_4\! = \!\{(2.1,0.9), (1.85,1.3), (-0.1,1.75), (-1.7,1.7), (-1.9, 0.05),(-1.94,-1.9),(0,-2.13),(1.8,-1.8)\}, $

    and

    $ S_5 = \{(2.3,1.3),(2.1,1.7),(0.3,1.3),(-1.3,1.5),(-2.3,-0.3),(-2.3,-2.3),(-0.3,-1.8),(1.5,-1.6)\}, $

    respectively. Taking $ S_4 $ and $ S_5 $ as the initial guesses for the source points respectively, we then compare the numerical results with the case where the initial guess is determined by the DSM. The relative $ L^2 $ errors for $ k = 5 $ are listed in Table 5. Furthermore, we present the reconstructions of the circle as well as the source points in Figure 6. From Table 5 and Figure 6, one can observe that the DSM plays a significant role in the reconstruction by providing a good initial guess for the optimization process, which is a key step for achieving a satisfactory reconstruction. If the initial source points somewhat deviate from the true ones, then the reconstruction would drastically deteriorate, see Figure 6(b). Therefore, results of the DSM indeed exert an immense impact on the effectiveness of the optimization. And a good initial knowledge generated by a high-resolution sampling grid $ \mathcal{T}_1 $ is the essential prerequisite for achieving the high-quality reconstruction.

    Table 5.  Relative $ L^2 $ errors for reconstructions of the circle subject to different initial guesses.
    $ S_4 $ $ S_5 $ DSM
    $ k=5 $ $ 5.26\% $ $ 48.87\% $ $ 1.78\% $
    $ k=8 $ $ 3.35\% $ $ 31.41\% $ $ 1.22\% $

     | Show Table
    DownLoad: CSV
    Figure 6.  Reconstructions of the circle and $ 8 $ source points with different initial guesses. (a) $ S_4; $ (b) $ S_5; $ (c) DSM.

    Example 3. This example is designed to check the validity of our method when the boundary of obstacle is not star-like. Now we consider the reconstruction of two kite-shaped obstacles with boundaries $ \partial D $ described by the parametric representations

    $ kiteI: x(1)K:=(x(1)1,x(1)2)=(cost+0.15sint+0.35cos2t0.35,1.2sint),kiteII: x(2)K:=(x(2)1,x(2)2)=(cost+0.65cos2t0.65,1.5sint), $

    respectively. In this example, the sampling domains are taken to be $ \Omega_1 = [-2.5, 2.5]\times[-2.5, 2.5], $ $ \Omega_2 = [-1.5, 1.5]\times[-1.5, 1.5] $. To compute the accuracy of reconstruction, the relative errors of the boundary curves are computed by (4.2) with

    $ \tau_i^{(\ell)} = \mathit{\rm{arc}}\tan\frac{x_2^{(\ell)}(t_i)}{x_1^{(\ell)}(t_i)},\quad \ell = 1,2. $

    In Figure 7, we take the ansatz function space of different dimensions and present the reconstruction of the kite-I with $ k = 5 $ and $ \alpha = 10^{-6}. $ By taking different $ M( = 2, 4, 6, 8, 12, 20) $ in (4.1), we compute the reconstruction errors defined by (4.2). From the errors in Table 6, we observe that the accuracy of inversion is insensitive to moderately large $ M $. Except for the cases $ M = 2 $ and $ M = 4 $, the reconstruction is almost negligibly influenced by the dimension of $ U_M $.

    Figure 7.  Reconstructions of the kite-I and $ 12 $ source points with subspace $ U_M $ of different dimensions: (a) $ M = 2; $ (b) $ M = 4; $ (c) $ M = 12; $ (d) $ M = 20 $.
    Table 6.  Relative $ L^2 $ errors for reconstructions of the kite-I with 12 source points.
    $ M $ 2 4 6 8 12 20
    $ E_D $ $ 14.39\% $ $ 6.68\% $ $ 5.97\% $ $ 5.77\% $ $ 5.75\% $ $ 5.75\% $

     | Show Table
    DownLoad: CSV

    Next, we fix $ M = 8 $ and reconstruct the two kites. In Table 7, we provide the errors for reconstructing the kites with different number of source points. These results illustrate that the reconstructions of kite-I are superior to those of kite-II, which is probably due to the fact that the boundary fluctuation of kite-I is relatively mild in comparison to that of kite-II.

    Table 7.  Relative $ L^2 $ errors for reconstructions of the kite-shaped obstacles.
    Error of kite-I Error of kite-II
    $ N $ $ k=3, \alpha=10^{-10} $ $ k=5, \alpha=10^{-9} $ $ k=3, \alpha=10^{-9} $ $ k=5, \alpha=1\times10^{-8} $
    $ 2 $ $ 10.29\% $ $ 7.00\% $ $ 10.57\% $ $ 8.27\% $
    $ 4 $ $ 7.72\% $ $ 4.81\% $ $ 11.30\% $ $ 9.96\% $
    $ 6 $ $ 7.33\% $ $ 5.07\% $ $ 11.46\% $ $ 12.98\% $
    $ 8 $ $ 7.33\% $ $ 4.63\% $ $ 11.87\% $ $ 10.91\% $
    $ 10 $ $ 6.93\% $ $ 4.84\% $ $ 12.04\% $ $ 10.56\% $

     | Show Table
    DownLoad: CSV

    In Figure 8, we plot the reconstruction of kite-II with $ N(N = 2, 6, 8, 10) $ source points. It can be found that both the source points and the obstacle are reconstructed well. So the proposed method also works for the recovery of non-starlike scatterers.

    Figure 8.  Reconstructions of kite-II and source points of different number: (a) $ N = 2; $ (b) $ N = 6; $ (c) $ N = 8; $ (d) $ N = 10 $.

    Example 4. In the last example, we plan to reconstruct the obstacle-source pair with limited-aperture measurements. Motivated by realistic applications where only partial data is available, we test the proposed method when the total field is only accessed on a portion of the measurement curve. In this example, we collect the total field recorded by the receivers located on a part of $ \Gamma_R $ with the aperture chosen to be $ \theta = 3\pi/2 $ and $ \pi, $ respectively. The regularization parameter is chosen to be $ \alpha = 10^{-6} $. The sampling domains are $ \Omega_1 = [-2.5, 2.5]\times[-2.5, 2.5] $ and $ \Omega_2 = [-1.5, 1.5]\times[-1.5, 1.5] $.

    In Table 8, we list the relative $ L^2 $ errors of the boundary curves with respect to different parameters $ k $, $ \theta $ and $ N $. Figure 9 shows the reconstruction for the starfish together with $ N( = 4, 8) $ source points in the limited-aperture case with $ k = 8 $. Figure 9 illustrates that the part of boundary curve encompassed by the receivers can always be well reconstructed. Typically, due to the lack of information, recovery of the opposite side oriented towards the open aperture is less accurate. Anyhow, all the numerical results show that, even with limited-aperture data, the reconstruction method still reasonably works.

    Table 8.  Relative $ L^2 $ errors for the reconstructions of the starfish with limited angles.
    $ k=8 $ $ k=5 $
    $ \theta $ $ N=8 $ $ N=4 $ $ N=8 $ $ N=4 $
    $ \pi $ $ 14.30\% $ $ 16.17\% $ $ 8.39\% $ $ 12.04\% $
    $ 3\pi/2 $ $ 7.06\% $ $ 7.33\% $ $ 7.99\% $ $ 7.92\% $

     | Show Table
    DownLoad: CSV
    Figure 9.  Reconstructions of the starfish and different numbers of source points from limited aperture data. (The 1st column: problem geometry; the 2nd column: exact and reconstructed locations of the source points with DSM; the 3rd column: image of $ I_D(y) $; the 4th column: reconstruction by the optimization method. The 1st row: $ N = 4, \theta = 3\pi/2 $; the 2nd row: $ N = 8, \theta = 3\pi/2 $; the 3rd row: $ N = 4, \theta = \pi $; the 4th row: $ N = 8, \theta = \pi $.).

    In this work, we develop a novel optimization method to simultaneously reconstruct a sound-soft obstacle and multiple source points from near field scattering data. To obtain a good initial guess for the optimization method, a two-step sampling method is also developed to respectively image the rough source points and obstacle. Theoretically, we analyze the convergence of the optimization method and the indicating behaviors of the sampling scheme. Numerical results are provided to demonstrate the feasibility and promising features of the proposed method. Our future work consists of the theoretical justification of the uniqueness issue, which may be established by exploiting the nonlinearity of the DtN map and the Complex Geometrical Optics (CGO) solutions as in [29]. We will also consider the extension to three-dimensional model and the other boundary conditions, as well as the applicability in the scenarios of electromagnetic or elastic waves.

    This work was supported by the NSFC grant under NO. 11971133 and the Fundamental Research Funds for the Central Universities. The authors would like to thank Professor Xianchao Wang for valuable discussions. We would also like to thank the anonymous referees for their insightful and constructive suggestions on improving the quality of this paper.

    The authors declare there is no conflicts of interest.



    [1] On a Keller-Segel system with logarithmic sensitivity and non-diffusive chemical. Discrete Contin. Dyn. Syst. (2014) 34: 5165-5179.
    [2] Large-time regularity of viscous surface waves. Arch. Ration. Mech. Anal. (1983/84) 84: 307-352.
    [3] Global solutions of some chemotaxis and angiogenesis systems in high space dimensions. Milan J. Math. (2004) 72: 1-28.
    [4]

    A. Einolghozati, M. Sardari, A. Beirami and F. Fekri, Capacity of discrete molecular diffusion channels, Proc. IEEE International Symposium on Information Theory, (2011).

    [5]

    A. Einolghozati, M. Sardari and F. Fekri, Capacity of diffusion-based molecular communication with ligand receptors, Proc. IEEE Information Theory Workshop, (2011).

    [6] Maximum principles for the primitive equations of the atmosphere. Discrete Contin. Dynam. Systems (2001) 7: 343-362.
    [7] Mathematical analysis of a model for the initiation of angiogenesis. SIAM J. Math. Anal. (2002) 33: 1330-1355.
    [8] Stability of solutions of chemotaxis equations in reinforced random walks. J. Math. Anal. Appl. (2002) 272: 138-162.
    [9] Global existence and uniqueness of solutions for multidimensional weakly parabolic systems arising in chemistry and biology. Comm. Pure and Appl. Anal. (2007) 6: 287-309.
    [10] Nonlinear transmission problems for quasilinear diffusion systems. Networks and Heterogeneous Media (2007) 2: 359-381.
    [11] Some boundedness of solutions for the primitive equations of the atmosphere and the ocean. ZAMM Journal of Applied Mathematics and Mechanics (2015) 95: 38-48.
    [12] Local-in-time solvability of target detection model in molecular communication network. International Journal of Applied Mathematics (2018) 31: 427-455.
    [13] From 1970 until present: The Keller-Segel model in chemotaxis and its consequences. I. Jahresber Dtsch. Math.-Verein. (2003) 105: 103-165.
    [14] Convergence of solutions to simplified self-organizing target-detection model. Sci. Math. Japnonicae (2016) 81: 115-129.
    [15] A mathematical model of mon-diffusion-based mobile molecular communication networks. IEEE Comm. Lettr. (2017) 21: 1967-1972.
    [16] The stability and dynamics of a spike in the 1D Keller-Segel Model. IMA J. Appl. Math. (2007) 72: 140-162.
    [17] Model for chemotaxis. J. Theor. Biol. (1971) 30: 225-234.
    [18]

    O. A. Ladyženskaja, V. A. Solonnikov and N. N. Ural'ceva, Linear and Quasilinear Equations of Parabolic Type, Translations of Mathematical Monographs, Vol. 23, American Mathematical Society, Society, Providence, R.I., 1968.

    [19] (1968) Linear and Quasilinear Elliptic Equations. New York-London: Academic Press.
    [20]

    J.-L. Lions and E. Magenes, Non-Homogeneous Boundary Value Problems and Applications. Vol. I, Die Grundlehren der Mathematischen Wissenschaften, Band 181. Springer-Verlag, New York-Heidelberg, 1972.

    [21] (2013) Molecular Communication. Cambridge: Cambridge University Press.
    [22]

    T. Nakano and et al., Performance evaluation of leader-follower-based mobile molecular communication networks for target detection applications, IEEE Trans. Comm., 65 (2017), 663–676.

    [23] Remarks on strongly elliptic partial differential equations. Comm. Pure Appl. Math. (1955) 8: 649-675.
    [24]

    Y. Okaie and et al., Modeling and performance evaluation of mobile bionanocensor networks for target tracking, Proc. IEEE ICC, (2014), 3969–3974.

    [25]

    Y. Okaie and et al., Cooperative target tracking by a mobile bionanosensor network, IEEE Trans. Nanobioscience, 13 (2014), 267–277.

    [26] Atsushi Finite dimensional attractor for one-dimensional Keller-Segel equations. Funkcialaj Ekvacioj (2001) 44: 441-469.
    [27] Stationary solutions of chemotaxis systems. Trans. Amer. Math. Soc. (1985) 292: 531-556.
    [28] Some structures of the solution set for a stationary system of chemotaxis. Adv. Math. Sci. Appl. (2000) 10: 191-224.
    [29] On multivortex solutions in Chern-Simons gauge theory. Boll. Unione Mat. Ital. Sez. B Artic. Ric. Mat. (1998) 1: 109-121.
    [30] Global solutions to a chemotaxis system with non-diffusive memory. J. Math. Anal. Appl. (2014) 410: 908-917.
    [31] Instability of Turing patterns in reaction-diffusion-ODE systems. J. Math. Biol. (2017) 74: 583-618.
    [32] Surface waves for a compressible viscous fluid. J. Math. Fluid Mech. (2003) 5: 303-363.
    [33] Steady state solutions of a rReaction-diffusion systems modeling chemotaxis. Math. Nachr. (2002) 233/234: 221-236.
    [34]

    J. Wloka, Partielle Differentialgleichungen, B. G. Teubner, Stuttgart, 1982,500 pp.

  • This article has been cited by:

    1. Shoutao Tan, Zhanfeng Fang, Yanyi Liu, Zhe Wu, Hang Du, Renjie Xu, Yunfei Liu, An SSD-MobileNet Acceleration Strategy for FPGAs Based on Network Compression and Subgraph Fusion, 2022, 14, 1999-4907, 53, 10.3390/f14010053
    2. O. V. Maksymuk, V. V. Sobchuk, I. P. Salanda, Yu. V. Sachuk, A system of indicators and criteria for evaluation of the level of functional stability of information heterogenic networks, 2020, 7, 23129794, 285, 10.23939/mmc2020.02.285
    3. Hirotada Honda, 2023, 10.5772/intechopen.112599
  • Reader Comments
  • © 2019 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(2685) PDF downloads(385) Cited by(3)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog