Scattering by impenetrable scatterer in a stratified ocean waveguide

In this work, the direct and inverse scattering problems of wave impenetrable scatterers in the three-layered ocean waveguide are under investigation. We have established the well-posedness of forward problem and proposed a novel direct sampling method for the inverse problem. The direct recovery approach only applies the matrix-vector operations to approximate the wave impenetrable obstacle from the received partial data. The method is capable of reconstructing the objects of different shapes and locations, computationally quite cheap and easy to carry out. The theoretical analysis and the novel direct recovery algorithm are expected to have wide applications in the direct and inverse scattering problems of submerged acoustics.

1. Introduction. The forward and inverse scattering problems of submerged acoustics have wide applications in identification of submerged wreckages, navigational scatterers, submarines, mineral deposits and so on, and have received great attention in recent years; see [10,14,18,19] and references therein. One of the popular models employed for acoustics in a finite depth ocean is the one-layered waveguide bounded by two horizontal planes (see Figure 1(a)), and the corresponding theoretical and numerical achievements are presented in [2,12,17,20,21]. As we know that the ocean environment is consisted of multiple layers, thus a two-layered ocean waveguide (see Figure 1(b)) is under investigation recently. In the two-layered model, the well-posedness of direct problem, the uniqueness of inverse problem and several effective numerical algorithms are exhibited in [6,9,11,17].
However, it is well known that a typical sound velocity profile consists of the surface channel, the thermocline and isothermal layers [20]. Therefore, the two-layered model is not appropriate to represent the true ocean waveguide. For the purpose of maintaining the physical features of ocean structure, the horizontal stratified waveguide in [1,5] is brought in as a simple but reasonably realistic model, see Figure 2 for the geometrical illustration. In this model, the sound speed c is assumed to be depending fundamentally on the depth of the waveguide. A distinctive feature of the submerged acoustic model is the transmission loss (propagation loss), which is because of the absorption of acoustic energy by the geometric spreading (divergence effect) and the propagation medium [5]. Moreover, due to the reflection and refraction of the surface, bottom and the interfaces of each two adjacent layers, sound wave propagates in a special way. Especially, owing to the presence of two parallel infinite boundaries of the waveguide, only a finite number of wave modes are able to propagate at a long distance, while the other modes decay exponentially [22]. All the above phenomenons increase the ill-posedness of the inverse problem evidently, therefore the detection of buried scatterers D in the stratified system is more complicated compared with that in the homogeneous one. Nevertheless, a few efficient numerical methods are established for the reconstructions of wave penetrable scatterers in the three-layered ocean waveguide, for example, the multilevel sampling method (MSM) [14,18], the direct sampling method (DSM) [15,19], the simple method (SM) [13] and so on. However, as far as the author has noticed, the theoretical analysis and the numerical approaches for the detection of wave impenetrable obstacles in the three-layered ocean waveguide is seldomly addressed and considered in the literature. In addition, the total acoustic pressure field is entirely different for the wave penetrable obstacles and the wave impenetrable objects. Generally, in the stratified ocean waveguide model, a large mount of theoretical analysis and numerical methods of the wave penetrable scatterers are no longer applicable for the wave impenetrable inclusions. Consequently, we shall investigate the well-posedness of forward problem for the wave impenetrable objects embedded in the three-layered ocean waveguide in this paper. Moreover, a novel direct sampling method would be constructed for the corresponding inverse scattering problem in this work. Among all the existing techniques, one popular and effective type of algorithms [14,18,19] is relied on some indicator functions that exhibit extremely different behaviors for the interior and exterior of scatterers D. Based on the property of indicators, we can easily verify if a sampling point locates inside or outside the inhomogeneous inclusions, and thus recovering their positions and shapes. Motivated by the characteristic of the index function, we shall design a novel indicator to reconstruct the buried impenetrable scatterers in the stratified ocean waveguide.
The initial key step of the inverse scattering problem is to effectively estimate some domains that contain all the unknown obstacles. Otherwise, a fairly large computational region would bring in a huge additional computational burden for the entire reconstruction process when the initial key step is neglected. Our newly proposed direct sampling method (NDSM) is devoted to provide an appropriate initial computational areas for the recovery of wave impenetrable scatterers in the stratified ocean waveguide. Whereafter, any existing effective but computationally more demanding algorithms can be applied in the initial region for the accurate recoveries of physical features (i.e. density, refractive index, etc.) of different impenetrable objects. The direct approach only involves matrix-vector multiplications, and does not require any solutions of large-scale ill-posed linear systems, matrix inversions or optimization procedures. Accordingly, it is easy to carry out and have advantage in computation for both 2D and 3D.
We would like to mention that the simple method [13], the multilevel sampling method [14] and the direct sampling method [19] are proposed for imaging the embedded wave penetrable scatterers in our early work. Furthermore, these methods are fundamentally different from the novel direct sampling method we introduce in this work, as the method in [13] depends mainly on the indicator functions of [14,19], and the algorithm in [14] needs to update the values of an index function by solving integral equations recursively and is computationally quite expensive. In addition, the method in [19] applies the inner product of scattered field and the Green's function as the index function, and the behaviors of indicator in [19] is different from the one we propose in this paper. This paper is organized as follows. In Section 2, we shall describe the forward problem of three-layered ocean waveguide in the three dimensions. The wellposedness of direct problem would be stated as well. Moreover, we will suggest a numerical technique for solving the forward problem at the end of this section. Section 3 shows the construction of novel direct sampling method for the total field with different boundary conditions. Finally, some concluding remarks would be presented in Section 4.

2.
The description of direct scattering in a stratified ocean waveguide. We shall present the forward problem of our interest in the three dimension in this section. A three-layered ocean waveguide of the following form is under consideration where d is the height of waveguide. The bottom plane, surface plane and the three interior layers of stratified ocean waveguide are separately represented as follows with the third coordinate axis is singled out as the one orthogonal to the horizontal surface, and the heights satisfy 0 < d 1 < d 2 < d. Moreover, the interface between M 1 and M 2 is expressed as Γ 1 = x ∈ R 3 d x 3 = d 1 , and the other one between M 2 and M 3 is denoted as We suppose an impenetrable scatterer D of any geometric shape is simply connected and bounded C 2 domain in the ocean waveguide R 3 d , and the other region of waveguide is connected. We exhibit the entire physical waveguide R 3 d in Figure 2 for the geometrical demonstration, where the impenetrable object D is situated in the thermocline M 2 . The Green's function G(x; x s ) of the stratified ocean waveguide R 3 d is written as [1,18], , ϕ 1 and ϕ 2 are the normal modes of the propagating waves, W p (ϕ 1 , ϕ 2 ) represents the corresponding Wronskian, ζ p stands for the zeros of the Wronskian W and H (1) 0 is the zero-order Hankel function of a first kind.
The total acoustic pressure field u without convection in the three-layered ocean waveguide is of the following form By the separation of variables, we can obtain the mode expansion of propagating wave u as follows [1,18] (3) and the following decay property is satisfied by the mode expansion coefficient v p ( x), where r = | x|. In addition, the total acoustic pressure field u satisfies the Helmholtz system as follows, where x s = ( x s , x s3 ) expresses a source point located in the ocean waveguide, k i , n i and i (i = 1, 2, 3) are respectively the wavenumber, the refractive index and the density in the layer M i , ν means the normal vector, and Bu is depending on the physical feature of the impenetrable scatterer D which represents the following different boundary conditions, where ρ > 0 and ρ ∈ C(∂D).
We bring in the following piecewise functions for the ease of depiction, With the aid of (9), the equations (5) can be simplified as (10) ∆u And we have the following equation for the scattered field u s , (11) ∆u s + nu s = 0 in R 3 d \D . We consider an acoustic source point located at ξ and G(x; ξ) is the Green's function of the stratified ocean waveguide. Multiplying both sides of (10) by G(x; ξ) and integrating over Ω L , we are readily derive where Ω L is a cubic region which consists of vertical plane boundary Γ L (L is the length of cube) and two horizontal plane boundaries Γ + (part of ocean surface) and Γ − (part of ocean bottom). Accordingly, we are able to present the scattered acoustic pressure field u s as With the aid of the boundary conditions on Γ + and Γ − , we acquire (1) and (3), the scattered field possesses the following representation when L approaches to infinity [1,6] As we know that the incident field u i satisfies the Helmholtz equation in the ocean waveguide, thus we obtain Therefore, the total acoustic field u := u s + u i can be expressed as In addition, from (12) and (13), we can represent the scattered fields for different boundary conditions as follows (D) : sound-soft boundary condition (N) : sound-hard boundary condition (R) : impedance boundary condition In order to establish the well-posedness of the direct problem, we bring in some operators from C 0,α (∂D) into C 0,α (∂D) (0 < α < 1) [3] as follows and define the operator The total acoustic pressure field u can be expressed in the form of the above operators [7] (22) where u satisfies the Dirichlet boundary condition, say, u = f with f ∈ C(∂D), λ means a nonzero real coupling parameter and φ ∈ C(∂D) is a density function which satisfies Moreover, the total field u with Nuemann boundary condition is able to be written as where ∂u ∂ν = g with g ∈ C(∂D), S 0 denotes the potential theoretic limit case k = 0 and the continuous density function φ satisfies In addition, we can represent the total field u of Robin boundary condition as in (24), and the corresponding density function φ is the solution of following equation where ∂u ∂ν +iρ(x)u = h with h ∈ C(∂D). Now we would establish the well-posedness of the direct scattering problem by the following theorem Theorem 2.1. The density function φ can be uniquely determined by (23), (25) and (26) for different boundary conditions, and thus the total acoustic pressure field u is uniquely represented by (22) and (24) separately with the corresponding boundary condition. Moreover, the direct scattering problem of wave impenetrable scatterers in the stratified ocean waveguide is well-posed.
The proof of this theorem is similar to the Theorem 3.11-3.13 in [7]. Based on the above well-posedness theorem, we suggest the Nyström method [7] and the finite element method to simulate the total field and the corresponding scattered field in the practical applications.

3.
A novel direct sampling method for the inverse scattering problem. In this section, we shall propose a novel direct sampling method (NDSM) to provide a stable and reliable approximation for the physical features (i.e. locations, shapes) of impenetrable scatterers Ω in the stratified ocean waveguide R 3 d . The essential part of NDSM is to design an indicator function which has fairly different behaviors between the boundary and the other part of the impenetrable obstacles, namely, the value of index function is relatively large on the boundary of objects while relatively small in the other part of it. As far as the author concerns about this problem, the existing methods are only applicable for the wave penetrable inhomogeneities in the stratified waveguide. Hence, a novel index function to detect the wave impenetrable inclusions is considerably significant and indispensable.
We know that the Green's function G(x; x m ) of stratified ocean waveguide satisfies, where Ω R is a cylindrical region: for R and L are respectively the radius and the height of Ω R , we refer to Figure 3 for the geometrical illustration. 3.1. A priori estimate. In this section, a priori estimate would be attained which is the critical part for the construction of indicator function. We multiply (27) by the conjugate G(x, x n ) of Green's function G(x, x n ) and integrating over the domain Ω R , we readily derive Now we replace x m in (27) by point x n ∈ Ω R , and take its conjugate. Multiplying both sides by G(x, x m ) and take the integral over Ω R , we obtain With the help of (28) and (29), we attain the following equality by the Green's formula, where {·} denotes the imaginary part of the complex function, Γ R means the surface of cylinder Ω R and ν represents the normal vector. Recall the expression of the Green's function in (1), we are able to rewrite G(x, x m ) as follows, Wp(ϕ1,ϕ2) . We can now calculate the normal derivative of G as follows, From the orthogonal property of Ψ p (·), the following approximation can be acquired for p ≥ 1 [19], where L means the height of the cylinder Ω R . From (30), (31) and the approximation (34), a priori estimate is acquired as follows Assume that x ∈ Γ R , the total field u(x) has the following representation We shall investigate the inner product of the total field u and the Green's function G over the surface Γ R in the following two subsections.

3.2.
Indicator function for the total field u with the (D) boundary condition. In this subsection, the construction of indicator function for the total field u with the Dirichlet boundary condition is under consideration. We first study the following integration over the cylindrical surface Γ R : From (32), (35) and the orthogonal property of Ψ p (·), the integral (36) is able to be simplified as follows, With help of (31), (35) and the orthogonal property of Ψ p (·), we can predigest the integral of (37) in the following form, Accordingly, we can possess the following approximation where We would like to mention that φ is a continuous function at the boundary, Ψ p (·) are the combinations of sine and cosine functions, W p is the corresponding Wronskians and H (1) 0 (·) is the zero-order Hankel function of the first kind. Therefore, it can be concluded that the integration in (40) possesses a relative large value when x m approaches to the boundary of scatterer D and decays quickly as x m moves away from it. This typical behavior motivates us with the following indicator function, where Ω denotes the sampling region and the inner product is of the following form In the practical applications, the integral (40) is estimated as follows where the boundary ∂D is discretized into a set of tiny segments, the point ξ j is situated in the jth segment, and s j means the area of jth segment. We would like to mention that the identity (43) is the Helmholtz-Kirchhoff type identity, and it plays a key role in understanding the resolution limit in imaging with waves [4]. From (41) and (43), we can observe that the sharper the behavior of the imaginary part of function (41) around the boundary of scatterer is, the higher is the resolution.

3.3.
Indicator function for the total field u with the (N) or (R) boundary condition. In this subsection, we shall design an indicator function for the total field u with the Nuemann or Robin boundary condition. Likewise, the following integration over the cylindrical surface Γ R is under investigation: From (38) and (39), the integrations (44) and (45) can be expressed respectively in the following forms Consequently, the following estimation is able to be obtained where As we know that c 1 ≤ |ζ p | ≤ c 2 for p ∈ N, where c 1 and c 2 are positive constants [6]. We then obtain where It is known that Ψ p (·) are the combinations of sine and cosine functions, W p is the corresponding Wronskians and H (1) 0 (·) is the zero-order Hankel function of the first kind. Hence, we can deduce that the integral in (46) has a relative large value when x m moves to the boundary of scatterer D and decays quickly as x m tends away from it. Similarly, we can define the index function as follows, where Ω denotes the sampling region and the inner product is of the following form The same technique in (43) can be applied for the discretization of (46) in practice.
We can observe that the computation of indicator function I(x m ) only involves matrix-vector multiplications, and does not require any solutions of large-scale illposed linear systems, matrix inversions or optimization procedures. Accordingly, it is explicit, simple and straightforward. Based on the indicator functions that exhibit considerably different behaviors for the boundary and the other parts of scatterers D, the novel direct sampling method would be effective and robust for the recovery of wave impenetrable objects of different physical features in the practical applications.
In practice, we can directly calculate the value of index function at every sampling point, and the computation complexity of this procedure is O(N 2 ) in 2D and O(N 3 ) in 3D, where N is the number of sampling points along one direction. However, the computation complexity is quite demanding when N is large. In order to get rid of this drawback, we suggest implementing the simple method [13] with the aforementioned indicator function to recover the buried impenetrable scatterers since the computation complexity of the method is O(log N ) for both 2D and 3D, and thus the computation complexity is sharply reduced.
Finally, we would like to remark that many effective direct imaging functionals are introduced in the book [4]. And we apply the Green function (31) into those direct imaging functionals would lead to direct imaging functionals in stratified waveguide in the multiscale configuration. Moreover, the estimations of signal to noise ratio for these algorithms can be derived as well. 4. Numerical experiments. In this section, we shall demonstrate some numerical experiments to verify the robustness and effectiveness of the proposed novel direct sampling method (NDSM). All the programs in our examples are written in MATLAB and run on a 2.83 GHz PC with 16GB memory.
Initially, we present the parameters that are employed in our numerical simulations. The ocean depth is set to be 100, and the heights of the first interface Γ 1 and the second interface Γ 2 of the ocean waveguide are selected as d 1 = 100/3 and d 2 = 200/3 separately. In each layer M i (i = 1, 2, 3), the wavenumber is chosen as k i = 2πf /c i , where f=75, c 1 = 1450, c 2 = 1500 and c 3 = 1550, and the refractive indies are set to be n 1 = 1, n 2 = 2 and n 3 = 1/2, while the densities are selected as 1 = 1400, 2 = 1500, 3 = 1600 and 4 = 1550 separately. For the purpose of better observation, we only keep the sampling point with the value of index function larger than 0.95 in the following numerical tests. And the finite element method is applied to generate the total field u. In addition, we add some random noises into the exact data as follows, whereρ = 10% denotes the noise level, Θ 1 (x r ) and Θ 2 (x r ) are two random numbers varying from -1 to 1. We would like to mention that when computing the total field of the forward problem, we have made the results as accurate as our machine can do, and the remaining error serves as noise in the inverse problem. The source points (the red star points) and the receivers (cyan circle points) are located on Γ s and Γ r respectively, and the obstacle D is embedded in the sampling domain Ω, see Figure 4 for the geometrical exhibition.   In order to optimize the computation complexity of NDSM, we apply the parallel bisection technique [13] and the NDSM to recover the obstacles, and the approximations under 10% and 15% random noises are shown in Figure 5 (c) and (d) respectively. The black cubes are the exact scatterers while the blue star points are the reconstructions. One can observe that the impenetrable inclusions can be separated clearly. Considering the partial detected data is only available in the upper part of layers M 2 and M 3 , thus the recovered region of D 1 is not as good as the one of D 2 . Experiment 2. In this example, we consider two impenetrable inhomogeneities of different shapes and scales, which are situated at We can find that the estimations are also quite impressive and reliable for the scatterers of different shapes, scales and boundary conditions. Furthermore, the two inclusions can be clearly separated when the distance between them is quite small. Accordingly, we can conclude that the proposed approach is able to deal with the impenetrable components of different physical features.  We could see that the reconstruction is quite satisfactory and good, and the estimate is acceptable under even 15% random noise. Furthermore, we can conclude that the proposed algorithm is able to deal with the impenetrable component of special shapes. Consequently, the approach can be viewed as a feasible and effective alternative to provide impressive recoveries.

5.
Conclusions. This work states the well-posedness of direct scattering problem for the wave penetrable scatterers embedded in the stratified ocean waveguide, and designs a novel direct sampling approach which aims to provide an initial computational domain for the numerical recovery of impenetrable obstacles in the threelayered ocean waveguide. The proposed method only applies the matrix-vector operations and does not implement any matrix inversions or optimization procedures, and thus it can be considered as an efficient and simple approach. The theoretical analysis and the novel direct recovery algorithm are expected to have wide applications in the direct and inverse scattering problems of submerged acoustics.