FACTORIZATION METHOD IN INVERSE INTERACTION PROBLEMS WITH BI-PERIODIC INTERFACES BETWEEN ACOUSTIC AND ELASTIC WAVES

Consider a time-harmonic acoustic wave incident onto a doubly periodic (biperiodic) interface from a homogeneous compressible inviscid fluid. The region below the interface is supposed to be an isotropic linearly elastic solid. This paper is concerned with the inverse fluid-solid interaction (FSI) problem of recovering the unbounded periodic interface separating the fluid and solid. We provide a theoretical justification of the factorization method for precisely characterizing the region occupied by the elastic solid by utilizing the scattered acoustic waves measured in the fluid. A computational criterion and a uniqueness result are presented with infinitely many incident acoustic waves having common quasiperiodicity parameters. Numerical examples in 2D are demonstrated to show the validity and accuracy of the inversion algorithm.


Introduction
Consider a time-harmonic acoustic wave incident onto an unbounded doubly periodic (or bi-periodic) surface from above.The medium above the surface is supposed to be filled with homogeneous compressible inviscid fluid with a constant mass density, whereas the region below is occupied by an isotropic linearly elastic solid characterized by the Lamé constants.Due to the external incident acoustic field, an elastic wave propagating downward is incited inside the solid, while the incident acoustic wave is scattered back into the fluid (see Figure 1).This leads to the fluid-solid interaction (FSI) problem with unbounded bi-periodic interfaces between acoustic and elastic waves.Our paper concerns the inverse FSI problem of recovering the unbounded periodic interface separating the fluid and solid from knowledge of the scattered acoustic waves measured in the fluid.The periodic interfaces exist for many reasons, e.g., grain structure, lamination and fiber reinforcement as well as in the manufacturing of material surfaces.In particular, the fluid-solid (for example, water-brass, water-perspex) interface is very common in the real world, and the investigation of Rayleigh waves generated by periodic structures can be important in developing new surface acoustic wave devices and planar actuators ( [8]).We refer to [8,5,21,30] and references therein for the discussions of the forward scattering model with practical applications in underwater acoustics, sonic and photonic crystals as well as in the field of ultrasonic non-destructive evaluation.Fluid-solid interfaces are also important in seismology and exploration geophysics for the acoustical characterization of the first layer of the sea floor.This is because the propagation of Stoneley-Scholte waves, another type of surface waves that exit on a fluid-solid interface (see [21] for a detailed discussion), is strongly related to the shear wave velocity of the sedimentary bottom [14].Inversion schemes for identifying the shape of a (non-periodic) bounded elastic body immersed in fluid can be found in [11,12] where an optimization-based technique is applied and in [31,32] using the reciprocity gap (RG) and linear sampling methods (LSM).The factorization method established in [25] provides not only a necessary and sufficient condition for characterizing the scatterer but also a rigorous justification of the modified linear sampling method.Compared to the optimization method, the previously mentioned sampling-type methods require neither computation of direct solutions nor initial guesses.Uniqueness with far-field patterns for all incident directions at a single frequency is a side product of the factorization method [25] and has also been verified in [31] based on a technical boundary regularity result.
To investigate the inverse problem with an unbounded interface, we combine the novel functional framework of [25] with the original ideas from [1,2] for inverse grating diffraction problems modeled by the Helmholtz equation.Our studies are carried out analogously to [25] but modified to be applicable to the case of unbounded periodic interfaces.See also [28,18] for applications of the factorization method to inverse electromagnetic scattering from a penetrable periodic layer and to inverse elastic scattering from a rigid periodic interface.Due to the exponentially decaying of the surface Rayleigh waves and the quasi-periodicity of the wave fields, we use near-field data measured in a single periodic cell other than far-field data to recover the scattering interface.In particular, the Rayleigh coefficients are utilized to establish a proper factorization of the near-field operator.Such an idea has recently led to the inversion scheme [19] of imaging non-periodic bounded obstacles by using the Fourier coefficients of radiating waves over a sphere.Our arguments could be readily extended to inverse fluid-solid interactions problems in a waveguide, which might have import applications in underwater acoustics.
We set up the mathematical formulation of the forward scattering problem following the grating diffraction problems for pure acoustic and elastic waves (see e.g., [23,3,13,9,10]).Coupling conditions on the interface (see e.g., [17,29]) are employed to model the interacting of acoustic and elastic waves.The scattered acoustic fields and transmitted elastic waves are required to satisfy the upward and downward quasi-periodic Rayleigh expansions with a finite number of propagating plane waves, respectively.Existence and uniqueness of quasiperiodic solutions were established in [20] via the variational approach for all frequencies excluding possibly a discrete set with the only accumulating point at infinity.The variational argument applies directly to the second auxiliary quasi-periodic boundary value problems (QBVP) introduced in Section 4 for all but a discrete set of frequencies.The role of the third QBVP (III) (see (19)) is the same as the interior transmission eigenvalue problem discussed as in [31,25].In the Appendix, we shall verify the unique solvability of the problem (III) in a periodic strip for small frequencies or for all frequencies except for a possible discrete set.
In contrast to the bounded obstacle case, the incident angles have to be restricted to the upper half-space in order to identify the unbounded interface from above.However, it seems not suitable to employ incident waves with all possible angles, since the quasiperiodicity parameters of the scattered fields vary with the direction of incidence.Inspired by earlier factorization methods for diffraction gratings [1] and for bounded obstacles in a half-space [24, Chapter 2.6], we utilize a set of incident acoustic waves with common quasiperiodicity parameters to factorize the near-field operator arising from our inverse FSI problem.In comparison with the far-field operator (see [25]), this near-field operator fails to be normal.Thanks to properties of the middle operator (see Lemma 6.1), we can still apply a proper range identity to the resulting factorization of the near-field operator.The main technical difficulty in our analysis is the denseness proof of the solution operator which is carried out in a non-trivial way and can be straightforwardly extended to the nonperiodic case; see Lemma 5.3.In this paper we define the adjoint of the middle operator in two different ways (cf.(34) and (38)).We employ the first definition given in (34) to prove the injectivity of the middle operator in an easy way and to provide a simple compactness proof in Lemma 6.1 (iii) (cf.[25,Theorem 2.4]) as well.The second definition (see Remark 3 (i)) turns out to be equivalent to that for bounded elastic bodies.In addition, a small mistake in the proof of [25,Lemma 2.5] is corrected in Lemma 5.1 of the current paper.
The paper is organized as follows.In Section 2 we rigorously formulate the direct and inverse interaction problems with bi-periodic Lipschitz surfaces between acoustic and elastic waves.Sections 3 and 4 are devoted to introducing an admissible set of incident acoustic waves and several auxiliary quasiperiodic boundary values problems for facilitating a factorization of the near-field operator to be defined later.The solution operator is defined in Section 4 and its properties are verified in Section 5. We introduce the near-field operator and derive its factorization in Section 6.In Section 7, a theoretical justification of the factorization method is presented for characterizing the elastic body in terms of the spectrum of the nearfield operator.Numerical experiments are reported in Section 8 where the sensitivity of the inversion algorithm to several parameters is tested.
We end up this section by introducing some notation that will be used throughout the paper.Denote by (•) the transpose of a vector or a matrix, and by (•) * the adjoint of an operator.The symbols e j , j = 1, 2, 3, denote the Cartesian unit vectors in R 3 .For a ∈ C, let |a| denote its modulus, and for a ∈ R 3 , let |a| denote its Euclidean norm.The notation a • b stands for the inner product we write x = (x 1 , x 2 ) so that x = (x, x 3 ).We shall use the symbol ∂ j to denote ∂/∂x j .

Mathematical formulations
Assume that a bi-periodic Lipschitz surface Γ divides the space R 3 into two unbounded connected components Ω + and Ω − which are above and below Γ, respectively; see Figure 2. Without loss of generality, we suppose Γ to be 2π-periodic in x 1 and x 2 , i.e., Note that in our studies Γ is not necessarily the graph of some bi-periodic function.It is supposed that Ω + is filled with a homogeneous compressible inviscid fluid with the constant mass density ρ f > 0, and that Ω − is occupied by an isotropic linearly elastic solid characterized by the real-valued constant mass density ρ > 0 and the Lamé constants λ, µ ∈ R satisfying µ > 0, 3λ + 2µ > 0.
Assume an acoustic wave is incident onto Γ ⊂ R 3 from Ω + .The incident wave is supposed to be a time-harmonic plane wave of the form v in (x) exp(−iωt) with the angular frequency ω > 0 and speed of sound c 0 > 0, where the spatially dependent function v in takes the form In (1), θ ∈ S 2 := {x ∈ R 3 : |x| = 1} denotes the incident direction with the incident angles θ 1 ∈ [0, π/2), θ 2 ∈ [0, 2π), and k := ω/c 0 is the wave number in the fluid.
Under the hypothesis of small amplitude oscillations both in the solid and fluid, the direct or forward scattering problem can be formulated as the following boundary value problem: Find the scalar total acoustic field v = v in + v sc and the transmitted elastic field u = (u 1 , u 2 , u 3 ) generated from a known (prescribed) incident wave v in such that (see e.g., [29,17]) Here, the notation ν = (ν 1 , ν 2 , ν 3 ) ∈ S 2 denotes the unit normal vector on Γ pointing into Ω − and ∂ ν u = ν •∇ u.In (2), T u stands for the stress vector or traction having the form Let us Betti's formula (see e.g., [26]) in a bounded Lipschitz domain D ⊂ R 3 : 3 , with It can be observed that the role of the stress operator (3) in the Lamé equation is the same as that of the normal derivative in the scalar Helmholtz equation.
Throughout the paper, we define the quasi-periodic parameter Obviously, the incident field v in is α-quasiperiodic in the sense that v in (x) exp(−iα• x) is 2π-periodic with respect to x.The periodicity of the interface together with the form of the incident wave implies that the solution (v, u) must also be α-quasiperiodic.Hence, for w = v in Ω + and w = u in Ω − it holds that Since the domain Ω ± is unbounded in the ±x 3 -direction, a radiation condition must be imposed at infinity to ensure well-posedness of the boundary value problem (2).Set We require the scattered acoustic field v sc to satisfy a weighted upward Rayleigh expansion (see [23,13,3] for the non-weighted version where w n = 1) with the parameters α n = (α n ) ∈ R 2 , η n ∈ C and the weights w n ∈ R given respectively by for all n = (n 1 , n 2 ) ∈ Z 2 and some fixed b > h + .In this paper, the constants v n ∈ C in ( 5) are referred to as the Rayleigh coefficients of the scattered acoustic wave and the weights w n are introduced only for the convenience of efficiently computing the Rayleigh coefficients when |α n | > k.To see the Rayleigh expansion of the elastic displacement, we decompose u into the compressional and shear parts, where k p and k s are the compressional and shear wave numbers defined as k p := ω ρ/(2µ + λ) , k s := ω ρ/µ, respectively.Analogously, the longitudinal and shear wave speeds are defined as follows: The scalar functions ϕ and the vector function ψ in (7) satisfy the homogeneous Helmholtz equations Applying the downward (non-weighted) Rayleigh expansion for the scalar Helmholtz equation to ϕ and ψ, we obtain an expansion of u into downward propagating elastic waves (see [10] for the upward Rayleigh expansion) The parameters β n = β n (ω) and γ n = γ n (ω) involved in (9) are defined analogously to η n with k replaced by k p and k s , respectively.Since η n , β n and γ n are real for at most finitely many indices n ∈ Z 2 , we observe that only a finite number of plane waves, namely, the modes corresponding to |η n | ≤ k in (5) and those corresponding to (9), propagate into the far field.The remaining part consists of evanescent (or surface) waves decaying exponentially as |x 3 | → +∞.Thus, the Rayleigh expansion (5) converges uniformly with all derivatives in the upper half-space x 3 > b for all b > h + , while (9) converges uniformly in the lower half-space x 3 < a for all a < h − .Now, we can formulate our direct scattering problem (DP) as the following boundary value problem.

Inverse Problems and Imaging
Volume 9, No. 1 (2015), X-XX (DP) Given a bi-periodic Lipschitz surface Γ ⊂ R 3 (which is 2π-periodic in x 1 and x 2 ) and an incident field v in , find a scalar function v = v in + v sc ∈ H 1 loc (Ω + ) and a vector function u ∈ H 1 loc (Ω − ) 3 that satisfy the equations and transmission conditions in (2), the quasiperiodic condition (4) and the radiation conditions ( 5) and ( 9).Below we recall from [20] the solvability results to (DP) established via variational arguments.Assume that The above condition is only a technical condition in the variational approach [20], under which the DtN mapping for the Lamé operator is well-defined.In particular, (11) can be guaranteed if ω is sufficiently small or if the relation λ + 2µ ≤ ρc 2 0 ( equivalently k ≤ k p ) holds.
Lemma 2.1.Let Γ be a bi-periodic Lipschitz surface.(i): For the incident plane wave of the form there always exists a weak solution to (DP).(ii): Assume λ + 2µ ≤ ρc 2 0 and let v in be a solution to the Helmholtz equation in Ω + .Then there exists a small frequency ω 0 > 0 such that (DP) is uniquely solvable for all ω ∈ (0, ω 0 ].Moreover, there admits a uniqueness solution for all frequencies excluding possibly a discrete set D with the only accumulation point at infinity. Note that the so-called Jones frequency arising from the FSI problem for bounded elastic bodies may exist in periodic structures (see [20]).Throughout the paper, the incident frequency is supposed to be such that (DP) is always solvable for any v in satisfying the Helmholtz equation (∆ + k 2 )v in = 0 in Ω + .Our numerical solutions show that the finite element method based on the variational formulation in a bounded periodic cell involving truncated DtN maps may be still stable even if the condition λ + 2µ ≤ ρc 2 0 is not true.This condition means that the acoustic wave speed of the fluid is bigger than that in the solid, i.e., c 0 ≥ c p > c s , which however does not apply to mostly encountered fluid-solid-interaction cases (where c < c s < c p ).
Since the evanescent (or surface) waves in ( 5) can be hardly measured in the fluid far away from the interface, we shall use near-field rather than far-field data to recover the interface.Our concern on the inverse problem is to detect Γ from knowledge of the scattered acoustic field for some b > h + .In our inversion algorithms, we shall send several incident waves from the admissible set for some M ∈ N, and then record the corresponding near-field data for each incident wave.In (12), the coefficients w j are given in (6) and the functions v in,j share the same quasiperiodicity parameter α for all j ∈ Z 2 .If |α j | ≤ k, v in,j is a downward (homogeneous) plane wave with the direction (α j , −η j )/k ∈ S 2 .In ultrasonics, such kind of waves can be generated by an ultrasound source of large dimensions (i.e., a large-diameter transducer as compared with the wavelength).When |α j | > k, v in,j is a downward evanescent wave mode taking the form of an inhomogeneous plane wave.It may be physically produced at the prism face by total internal reflection (see e.g.[4,7,15] and references therein for the practical use in near-field optics).Denote by v sc,j (x, x 3 ) the scattered acoustic field in the fluid generated by the incident wave v in,j ∈ A, and by (v The inverse problem (IP) under consideration can be formulated as follows (IP): Recover the scattering interface Γ from the set {v We remark that the Rayleigh coefficients v (j) n can be computed straightforwardly from the near-field data v sc,j (x, b) via an integral over Γ b , i.e., For the inverse problem (IP), we suppose there is a priori information that h + > 0 and that the unknown interface Γ lies between the flat surfaces Γ b and Γ Below we introduce notation to be used in subsequent sections.Since most of our discussions will be restricted to a single periodic cell, we set (see Figure 2) Clearly, Ω + b and Ω − 0 are both bounded Lipschitz domains in R3 .For simplicity we still use Γ to denote one period of the scattering surface.Let H s α (•) (s ∈ R) denote the Sobolev space of scalar functions which are α-quasiperiodic with respect to x 1 and x 2 .For s = 0, we write

A new admissible set of incident acoustic waves
Inspired by earlier factorization methods for diffraction gratings [1] as well as for bounded obstacle scattering in a half-space [24, Chapter 2.6], we define a new admissible set of incident acoustic waves for facilitating a factorization of the nearfield operator to be defined in Section 6.
We first recall the free space Green function Φ (k) (x, y) for the Helmholtz equation and the free space α-quasiperiodic Green function G (k) (x, y) for the Helmholtz equation: Inverse Problems and Imaging Volume 9, No. 1 (2015), X-XX for x−y = (2nπ, 0),n ∈ Z 2 with α n , η n defined as in (6).In the subsequent sections, we suppose η n = η n (ω) = 0 for all n ∈ Z 2 , so that ( 13) is well-defined.It is wellknown that the difference of Φ (k) and G (k) is an analytic function in R3 .Moreover, the function D (x, y), with y := (ỹ, −y 3 ), is the α-quasiperiodic Green function of the Helmholtz equation in the half space x 3 > 0 satisfying the homogeneous Dirichlet boundary condition on x 3 = 0.For x 3 > y 3 , we can rewrite with w n defined in (6) and The expression (14) implies that, for x 3 > y 3 , the function x → G (k) D (x, y) is actually a weighted upward-travelling solution with the Rayleigh coefficients given by g n (y)− g n (y ).In contrast with the incident plane wave (1), we introduce a 'new' admissible set of α-quasiperiodic incident acoustic waves: By the definitions of η j (see ( 6)) and w j , one can derive from (15) that = ṽin,j (y) for all j ∈ Z 2 , that is, the new incident wave ṽin,j coincides with the conjugate of the j-th Rayleigh coefficient of the function x → G (k) D (x, y) for x 3 > y 3 .
Remark 1.The incident wave ṽin,j ∈ Ã has been used in [1] for the Helmholtz equation.It consists of two parts: ṽin,j = ṽin,j 1 + ṽin,j 2 , ṽin,j One can observe that the first part is a downward wave mode, whereas the second part is an upward mode that is not physically meaningful as an incoming wave.By uniqueness, the solution (v sc , u) to the forward problem (DP) with the 'artificial' incoming wave v in = ṽin,j 2 is nothing than (−ṽ in,j 2 , 0).Thus the scattered acoustic field ṽsc,j in the fluid incited by the incident wave ṽin,j can be obtained by linear superposition.More precisely, the n-th Rayleigh coefficient of ṽsc,j is given by ṽ where v (j) n stands for the n-th Rayleigh coefficient of v sc,j and δ nj is the Kronecker delta function.

Auxiliary boundary value problems and DtN maps
The aim of this subsection is to introduce several auxiliary boundary value problems that will be used later for establishing the factorization method.For h ∈ H 1/2 α (Γ), consider the boundary value problem of finding an α-quasiperiodic Recall that the region Ω − 0 denotes the periodic layer between Γ 0 and the scattering interface in one periodic cell.The above problem (I) is uniquely solvable for each h ∈ H , where D1 consists of α-quasiperiodic Dirichlet eigenvalues of the negative Laplace in Ω − 0 .It was shown in [16] that D1 has only a finite number of eigenvalues in the finite interval (−N, N ) for any N > 0. Now we assume k 2 = ω 2 /c 2 0 / ∈ D1 .Then, the normal derivative of v on Γ defines the Dirichlet-to-Neumann (DtN) operator , where v is the unique solution to problem (I).In particular, for each incident wave ṽin,j ∈ Ã, it holds that Applying Green's identity, one can readily prove that Lemma 4.1.
−α (Γ), there holds the relation where the operator T −α : is defined analogously to T α with α replaced by −α.
From Lemma 4.1, we see T α is a self-adjoint operator with respect to the duality between H and that v sc and u satisfy the upward and downward Rayleigh expansions (5) resp.(9), respectively.By Lemma 2.1, problem (II) admits a unique solution for every ϕ ∈ H , where D is the discrete set appearing in Lemma 2.1 (ii) and D 1 := {ω : ω 2 /c 2 0 ∈ D1 }.Analogously to the data-to-pattern operator for bounded obstacle scattering problems, we define the solution operator G : where v n denotes the n-th Rayleigh coefficient of v sc solving problem (II).Note that we have G(ϕ) ∈ l 2 due to the fact that v sc (x, b) ∈ L 2 α (Γ b ) for any b > h + .Obviously, our scattering problem (DP) with v in = ṽin,j ∈ Ã can be reformulated as the boundary value problem (II) with ϕ = ṽin,j .This implies that G(ṽ in,j | Γ ) = (ṽ (j) n ) n∈Z 2 for all ṽin,j ∈ Ã, j ∈ Z 2 , Inverse Problems and Imaging Volume 9, No. 1 (2015), X-XX with ṽ(j) n given by ( 17).Below we introduce another DtN map defined on Γ 0 which is equivalent to the downward Rayleigh expansion of u for the Navier equation.α (Γ 0 ) 3 , the DtN operator T α is defined as , where T is the stress operator and u is the unique α-quasiperiodic solution to the homogeneous Navier equation subject to the Dirichlet boundary value u = w on Γ 0 and the downward Rayleigh expansion (9) in x 3 < 0.

T −
α is also called the Dirichlet-to-Neumann map, because it maps the Dirichlet data u| Γ0 of a radiation condition to the traction of u on Γ 0 .From the numerical point of view, one may extend the solution from Ω − 0 to the region x 3 < 0 by using the Rayleigh coefficients of u| Γ0 .The DtN map given by Definition 4.2 ensures the continuity of (T u)| Γ0 across the interface Γ 0 .The explicit expression of T − α (see [20]) shows that T − α is a bounded operator from H s α (Γ 0 ) 3 to H s−1 α (Γ 0 ) 3 for any s ∈ R. To justify the factorization method, we still need to consider the problem of with 3 .Problem (III) plays the same role of the interior transmission eigenvalue problem for the interaction problem with a nonperiodic bounded elastic body (cf.[25]).Since the variational argument for proving well-posedness of (III) differs from that for (DP), we shall verify the existence and uniqueness of solutions to (III) for all ω ∈ R + \D 2 with some discrete set D 2 ; see Lemma 9.1 in the Appendix.Further, one can observe that, if v is a solution to (I) In the subsequent sections we assume that the above problems (I), (II), (III) and (DP) are always uniquely solvable in both α-quasiperiodic and −α-quasiperiodic Sobolev spaces.In particular, the mapping (f, g)

Properties of solution operator
This section concerns properties of the solution operator G.We first show that the range Range(G) of G can be utilized to characterize the domain Ω − beneath Γ, and then prove the denseness of Range(G) in l 2 and the compactness of G.
Lemma 5.1.Let g n (y) be given by (15) with y ∈ R3 .The sequence (g n (y)) n∈Z 2 lies in the range of G if and only if y ∈ Ω− .Proof.Assume y ∈ Ω− .Let (w, u) be the unique solution to problem (III) with where G (k) (•, y) is the free-space α-quasiperiodic Green function to the Helmholtz equation.By the definition of T α , we see T α (w| Γ ) = (∂ ν w)| Γ .Hence, the solution (v sc , u) := (G (k) (•, y), u) solves problem (II) with ϕ = w| Γ .Together with the definition of G, this implies that where g n (y) is the n-th Rayleigh coefficient of the outgoing solution x → G (k) (x, y) in Ω+ .
Our choice of w in (20) corrects a mistake in the proof of [25,Lemma 2.5].Before proving the denseness of G, we derive an explicit expression of its adjoint G * .

Lemma 5.2. The explicit expression of G
with the over bar denoting the complex conjugate.Here, (ṽ sc , ũ) is the unique -αquasiperiodic solution pair to (II) with ϕ(y) = φ(y) = q(y)| Γ , q(y) := where Remark 2. We remark that w is α-quasiperiodic and converges absolutely in any compact set of x 3 > b.Moreover, since q is −α-quasiperiodic in Ω− , it follows from the transmission conditions in problem (II) that (ν Therefore, the operator G * is well-defined. Proof of Lemma 5.2.For ϕ ∈ H 1/2 α (Γ), denote by (v sc , u) the unique solution to problem (II) and by (v n ) n∈Z 2 the Rayleigh coefficients of v sc .By the definitions of G (see ( 18)), we see where •, • l 2 denotes the inner product in l 2 .Using the α-quasiperiodic Green's function in a half-space, we may represent v sc as Inserting the above expression into (23) and changing the integration order yield with q given by (22).Let (ṽ sc , ũ) be given as in Lemma 5.2, with the Rayleigh expansion ṽsc = Using the coupling conditions it follows from (24) that Applying Green's formula and using the upward Rayleigh expansions of v sc and ṽsc , one can straightforwardly verify that Combining ( 26) and ( 27) and making use of the coupling boundary conditions where in deriving the last equality we have used Lemma 4.1 and the relation which can be proved analogously to (27).Therefore, we obtain the expression of the adjoint operator G * shown as in (21).
Lemma 5.3.The operator G : is compact and has dense range.
Proof.Recalling that ω ∈ R is supposed not to be the exceptional frequency, the problem (II) is well-posed for any ϕ ∈ H 1/2 α (Γ).Let (v sc , u) be the unique solution with the stability estimate for some c > 0. We further observe that

Inverse Problems and Imaging
Volume 9, No. 1 (2015), X-XX To prove the denseness of G, it suffices to verify the injectivity of G * : l 2 → H −1/2 α (Γ).Supposing that G * (d) = 0, we claim that d = 0 in l 2 .Indeed, by the expression of G * (see ( 21)) we get ũ • ν = −T −α (ν • T ũ) on Γ.Hence, by (25), where ṽsc and ũ are given in Lemma 5.2.Let w be the unique −α-quasiperiodic solution to problem (I) with h = −(q + ν • T ũ)| Γ , and define the -α-quasiperiodic function The relations in (28) combined with the definition of T −α imply that where the superscripts + and − denote respectively the limits from above and below.Thus, w is a −α-quasiperiodic solution to the Helmholtz equation in the half-space x 3 > 0 satisfying the homogeneous Dirichlet boundary condition w = 0 on Γ 0 and the upward Rayleigh expansion.By uniqueness, we obtain w = 0 in x 3 > 0 and in particular, ṽsc = w = 0 in Ω + .Consequently, it follows from (28) and the transmission conditions between ũ and ṽsc that Hence, the solution pair (q, ũ) is the unique −α-quasiperiodic solution to problem (III) with f = 0, g = 0.By uniqueness it holds that q = 0 in Ω − 0 , and by unique continuation, q = 0 in Ω + b .Hence, we get q = 0 on Γ b and q = 0 in x 3 > b due to the uniqueness of α-quasiperiodic solutions to the exterior boundary value problem for a flat surface.Recalling the jump relation 0 = (∂ ν q) − − (∂ ν q) + = w/2 on Γ b , we obtain w = 0 on Γ b .This implies that d n = 0 for all n ∈ Z 2 by the definition of w.Therefore, d = 0 in l 2 .The proof is thus complete.

Near-field operator and its factorization
Analogously to the Herglotz wave function in bounded obstacle scattering problems, we define the operator H : The operator H is the restriction of a linear superposition of α-quasiperiodic waves ṽin,n with the weight a n to Γ.In view of the relation ( 16), one may derive the adjoint operator with g n given in (15).The near-field operator N : l 2 → l 2 is then defined as By the definition of G (see (18)), the near-field operator maps a superposition of incident waves to the set of Fourier coefficients of the corresponding scattered field on x 3 = b.
In the rest of the paper we shall use the quasiperiodic single layer potential whose kernel is the Dirichlet Green's function in the half-space x 3 > 0. Since (Sψ)| Γ0 = 0, we have T α ((Sψ)| Γ ) = ∂ ν (Sψ)| Γ by the definition of T α .The jump relations for the normal derivatives of Sψ on Γ lead to where the normal ν to Γ has been assumed to point into Now, we introduce the operator J : where w solves (33).Since Sψ satisfies the upward Rayleigh expansion (5) for x 3 > h + , rearranging the terms in (33) yields Together with (14) and the definition of G, this implies that for all ψ ∈ H −1/2 α (Γ), from which the relation H = J * G * /2 follows.Hence, we obtain from (30) a factorization of the near-field operator in which J * is always referred to as the middle operator.
that is, (p, u) solves problem (III) with f = ψ, g = 0. Hence, the operator J has an alternative representation as follows: This expression of J is equivalent to the one defined in [25,Section 2] for the inverse FSI problem with a bounded elastic body.

Inverse Problems and Imaging
Volume 9, No. 1 (2015), X-XX In the following we present properties of the operator J.

(Γ). (ii):
The operator J and its adjoint J * are both one-to-one.(iii): There exists a coercive operator J 0 : where c > 0 is a positive constant independent of ψ.
Proof.(i) For ψ ∈ H −1/2 α (Γ), let (p, u) be the unique solution to (37), so that the relation (38) holds.Using the jump relation (32) and the coupling conditions between u and p in (37), we arrive at Recalling the boundary condition w = 0 on Γ 0 and applying Green's formula to p and Betti's formula to u in Ω − 0 yield Similarly, from the vanishing data Sψ = 0 on Γ 0 and Green's formulae applied to Sψ in Ω + b and Ω − 0 , it follows that where the integral over Ω b,0 := {x : 0 < x 3 < b, 0 < x 1 , x 2 < 2π} is understood as the sum of the integrals over Ω − 0 and Ω + b .Inserting (40)-( 42) into (39) and taking the imaginary part of the resulting expression, we get Inverse Problems and Imaging Volume 9, No. 1 (2015), X-XX Employing the upward and downward Rayleigh expansions for Sψ and u, one can verify where ψ n is the n-th Rayleigh coefficient of Sψ for x 3 > h + , and with A p,n and A s,n being the corresponding Rayleigh coefficients of u.We refer to [10] for the proof of (43) where u satisfies a corresponding upward Rayleigh expansion.Consequently, we get Im ψ, Jψ ≥ 0 for all ψ ∈ H 1/2 α (Γ).(ii) To prove the injectivity of J, we employ the definition (34).Let w, u and Sψ be solutions to (33) and (35).Assume w| Γ = Jψ = 0 on Γ.By uniqueness to problem (I), we get w ≡ 0 in Ω − 0 , and in particular, ∂ ν w = 0 on Γ.Therefore, the solution pair (Sψ, u) solves problem (II) with ϕ = 0, implying that Sψ = 0 in Ω + b and u = 0 in Ω − 0 .In particular it holds Sψ = 0 on Γ. Recalling that Sψ also vanishes on Γ 0 and ω / ∈ D 1 , we obtain Sψ = 0 in Ω − 0 by uniqueness.Finally, applying the jump relation (32) yields ψ = 0, i.e., J is injective.In view of Remark 3 (ii), we can verify the injectivity of J * in the same manner.
(iii) Define the Green function G D and the single layer potential (Sψ) (i) in the same way as G (k) D and (Sψ) (k) with k = i (see (31)).We observe that Setting p 0 = Sψ (i) + w 0 , we have The coercivity of J 0 can be treated in the same manner as in the proof of [25,Theorem 2.4] for the fluid-solid interaction problem with a bounded elastic body.Note that the non-periodic bounded elastic body and its exterior in [25] correspond respectively to our periodic layer Ω − 0 and the region Ω+ in one periodic cell.In the present situation, the scalar function w 0 (or p 0 ) is additionally required to satisfy the homogeneous Dirichlet boundary condition on Γ 0 .The upward Rayleigh expansion for w 0 has been used in place of the Sommerfeld radiation condition arising from bounded elastic body scattering.For brevity we omit the proof here.

Inverse Problems and Imaging
Volume 9, No. 1 (2015), X-XX It remains to prove the compactness of Now, the compactness of J 1 : and the compactness of

Inversion algorithm
In this subsection we report the inversion algorithm for finding the bi-periodic interface separating the fluid and solid in one-periodic cell.By Lemma 5.1, the sequence (g n (y)) n∈Z 2 can be used to characterize the domain Ω − beneath Γ.However, we still need to bridge Range(G) with Range(N ), since the near-field operator N other than G can be straightforwardly computed from knowledge of the Rayleigh coefficients due to the admissible incident waves v in,j ∈ A. For this purpose, we shall apply the following abstract range identity (see [28]) to the factorization of the near-field operator established in (36).Lemma 7.1 (Range Identity).Let X ⊂ U ⊂ X * be a Gelfand triple with Hilbert space U and reflexive Banach space X such that the embedding is dense.Furthermore, let Y be a second Hilbert space and F : Y → Y , G : X → Y and T : X * → X be linear and bounded operators with F = GT G * .Suppose further that (a): G is compact and has a dense range.(b): There exists t ∈ (0, 2π) with cos t = 0 such that Re [exp(it)T ] has the form Re [exp(it)T ] = T 0 + T 1 with some compact operator T 1 and some coercive operator T 0 : X * → X, that is there exists c > 0 with ϕ, T 0 ϕ ≥ c ϕ 2 for all X * .(44) (c): Im (T ) is non-negative on X, that is, Im (T )ϕ, ϕ ≥ 0 for all ϕ ∈ X.
Moreover, we assume that one of the following conditions is fulfilled.The above range identity generalizes the one contained in [27, Chapter 1] and its proof was essentially based on the approach of Kirsch and Grinberg [24, Theorem 2.15] (cf.[28]).It plays a crucial role in various versions of the factorization method in wave scattering from impenetrable and penetrable scatterers.To apply Lemma 7.1, we set In our settings, all the conditions in Lemma 7.1 are satisfied.In fact, conditions (a) and (b) follow from Lemma 5.3 and Lemma 6.1 (iii), respectively.Conditions (c) and (d) are guaranteed by Lemma 6.1 (i) and (ii).Combining Lemmas 5.1 and 7.1 yields Theorem 7.2.Set N := |Re N | + |Im N |, and let g n (y) be given by (15).Then (i): The sequence (g n (y)) n∈Z 2 belongs to Range(N 1/2 ) if and only if y ∈ Ω− .
Note that the uniqueness described in Theorem 7.2 (ii) is only a corollary of the first assertion, and that Theorem 7.2 holds if k is not in any of the exceptional sets mentioned in Section 4. By Picard's theorem (see e.g, [6,Theorem 4.8]), the region Ω− below Γ can be characterized through the eigensystem of the near-field operator as follows.
Corollary 1.Let (λ j , ψ j ) be an eigensystem of the (positive) operator N .We have the following characterization of Ω − 0 : or equivalently, Thus, the interface Γ can be identified by first selecting sampling points from the set {(ỹ, y 3 ) ∈ R3 : 0 < y 3 < b} and then computing the values of the indicator function I(y).The values I(y) for y lying below Γ will be relatively larger than those above Γ which are actually zero.

Numerical experiments in two dimensions
In this section we report numerical experiments to test the validity and accuracy of the factorization method.For simplicity we supposed the bi-periodic interface to be invariant in the x 2 -direction, and its cross-section in the (x 1 , x 3 )-plane to be represented by the function x 3 = f (x 1 ) which is periodic in x 1 .Moreover, all elastic waves are assumed to be propagating perpendicular to the x 2 -axis (this implies that the incident angle θ 2 = 0), so that the interaction problem can be treated as a two-dimensional problem in the ox 1 x 3 -plane.The theoretical analysis of the factorization method in 3D carries over to this case straightforwardly.

Inverse Problems and Imaging
Volume 9, No. 1 (2015), X-XX To generate the synthetic scattering data for downward incoming waves v in,j ∈ A, we use the finite element method established in [20].The near-field operator N can be discretized by the (2M + 1) × (2M + 1) matrix , where the positive integer M is given as in (IP) and the Rayleigh coefficients ṽ(j) n are defined by (17) for −M ≤ j, n ≤ M .In our numerical experiments, we consider two smooth interfaces parameterized by (see Figures 3 (i) and (ii)) (i): f (x) = 0.8 + 0.4 sin(2x), x ∈ (0, 2π).
(ii): f (x) = 0.8 + 0.3 sin(x) + 0.2 sin(2x), x ∈ (0, 2π) and the piecewise linear profile given by Figures 3 (iii).Note that the periodicity of the surfaces (i) and (ii) are π and 2π, respectively.Our numerical tests show that the inversion algorithm can be applied to recover more oscillating profiles (e.g., surface (i)), if the periodicity of the interface is unknown in advance.
In the following numerical experiments, unless otherwise stated we always set This implies the quasiperiodicity parameters Hence, there are totally 14 upward acoustic propagating wave modes corresponding to the indexes {n ∈ Z : |n+α (1) | ≤ k} scattered back into the fluid in the x 1 x 3 -plane, given by exp(iα (1)  n We shall examine the sensitivity of the factorization method to the number M ∈ N, the acoustic wave number k > 0, the height of the measurement position b > h + as well as to the noise level δ of the near-field data.
Experiment 1.We apply the factorization method to reconstructing surfaces (i)-(iii) with different number of scattered wave modes.For fixed M ∈ N, there are exactly 2M + 1 wave modes involved in inversion algorithms.By (47), all upward propagating plane waves (i.e., the far-field data) are taken into account in the case M = 10, 15, 20 and k = 7. Figures 4 and 5 show that better numerical results can be achieved in recovering smooth profiles if we increase the number of surface (evanescent) wave modes.The same can be seen by comparing the pictures (c),(g), (h) and their contour lines (f),(i),(j) in Figure 6.However, the depth and corner points of the piecewise linear profile are not well-reconstructed.
Experiment 2. We test the imaging scheme with different acoustic wave numbers.Fixing the total number of scattered wave modes, we find that increasing the number of propagating wave modes gives rise to more stable and accurate images; see Figure 7 for recovering surface (i) and the pictures (a)-(f) in Figure 6 for surface (iii).The increased stability can be observed from the contour lines of the pictures.Note that a higher wave number implies more propagating plane wave modes, for instance, only six propagating wave modes exist when k = 3, less than the case of k = 5 or k = 7.This can be explained by the fact that the near-field measurement of propagating wave modes contain more information of the scattering surface than that of surface modes.The latter propagate only along the grating profile and decay exponentially in the x 3 -direction, whereas the propagating wave modes may travel along the x 3 -direction.
Experiment 3. We fix the number M = 20 and take the measurement of acoustic field at x 3 = b with b = 1.3, 1.5, 1.7.Since surface waves are exponentially decaying as b increases, lowering the height of the measurement position contributes to better imaging quality (see Figure 8).
Experiment 4. Since unpolluted data are used in the previous experiments.in the final experiment we investigate the sensitivity of the factorization method to the noisy data.The near-field acoustic data are perturbed by the multiplication of (1 + δ ξ) with the noise level δ, where ξ is an independent and uniformly distributed random variable generated between −1 and 1. Figure 9 illustrates the reconstructions from different noise levels at δ = 1%, 3%, 5%, respectively.The factorization method for solving inverse FSI problem can tolerate noisy data at a relatively low level, perhaps due to the large errors produced by our FEM code for solving the direct scattering problem.

Appendix
In this section we verify via the variational approach the unique solvability of the boundary value problem (III) for small frequencies or for all frequencies except for a possible discrete set.Arguing analogously to [10, Theorem 3], we can show the existence of a small frequency ω 1 such that the quasi-periodic boundary value problem 3 .Further, there holds the stability estimate for some constant C > 0 uniformly in all ω ∈ (0, ω 1 ].In [10], similar results were proved for upward propagating elastic waves in the upper half-space with the first, second, third or fourth kind boundary value condition on Γ, which carry over to the downward propagating waves as considered here. 3.Then there exists a small frequency ω 2 < ω 1 and a discrete set D 2 ⊂ R + such that the quasi-periodic coupling problem

Figure 1 .
Figure 1.Scattering of plane waves from an egg-crate shaped biperiodic interface separating the regions of fluid (above) and solid (below) in R 3 .

Figure 2 .
Figure 2. Geometry settings of our fluid-solid-interaction problem.

1 / 2 α
(Γ), and let (v sc , u) be the unique solution to (II) with the same ϕ.Then it holds thatG (k) (•, y) = v sc on Γ b .Furthermore, we have G (k) (•, y) = v sc in Ω+ \{y},due to the uniqueness to the Dirichlet boundary value problem in the half space x 3 > b and the unique continuation of solutions to the Helmholtz equation.If y ∈ Ω+ , the boundedness of lim x→y v sc (x) contradicts the singularity of G (k) (x, y) at x = y.If y ∈ Γ, we can always find an open domain D ⊂ Ω+ such that y ∈ ∂D.Clearly,

Figure 3 .
Figure3.The one-dimensional periodic interfaces to be reconstructed.

Figure 4 .
Figure 4. Reconstructions of the profile (i) with sensitivity to the number of wave modes (2M + 1) involved in inversion algorithms.We set k = 7.

Figure 5 .
Figure 5. Reconstructions of the profile (ii) with sensitivity to the number of wave modes (2M + 1) involved in inversion algorithms.We set k = 7.

Figure 6 .
Figure 6.Reconstructions of the piecewise linear profile.We set b = 1.3.

Figure 7 .
Figure 7. Reconstructions of the profiles (i) with sensitivity to the wave number k.We set M = 20, b = 1.3.