On an inverse source problem for the Biot's equations in electro-seismic imaging
aa r X i v : . [ m a t h . A P ] F e b ON AN INVERSE SOURCE PROBLEM FOR THE BIOTEQUATIONS IN ELECTRO-SEISMIC IMAGING
YIXIAN GAO, PEIJUN LI, AND YANG YANG
Abstract.
Electro-seismic imaging is a novel hybrid imaging modality ingeophysical exploration. This paper concerns an inverse source problem forBiot’s equations that arise in electro-seismic imaging. Using the time reversalmethod, we derive an explicit reconstruction formula, which immediately givesthe uniqueness and stability of the reconstructed solution. Introduction
Electro-Seismic (ES) imaging and Seismo-Electric (SE) imaging are emergingmodalities in geophysical exploration. Compared with traditional modalities, theyprovide images of high accuracy with low cost, which has led to wide applicationsin locating groundwater aquifers and petroleum hydrocarbon reservoir. The un-derlying physical phenomena of ES and SE imaging are known as electro-seismicconversion and seismo-electric conversion, respectively. These conversions usuallyoccur in a porous medium and have been used as bore-hole logging and cross-holelogging tools [1, 11, 15, 27].Both ES and SE imaging involve conversions of electromagnetic energy and seis-mic energy. The governing equations were derived by Pride [21]: ∇ × E = − µ∂ t H, (1) ∇ × H = ( ǫ∂ t + σ ) E + L ( −∇ p − ρ ∂ t u s ) , (2) ρ ∂ t u s + ρ ∂ t u f = div τ, (3) ρ ∂ t u s + ρ ∂ t u f + ηκ ∂ t u f + ∇ p = ηκ LE, (4) τ = ( λ div u s + q div u f ) I + µ ( ∇ u s + ( ∇ u s ) T ) , (5) − p = q div u s + r div u f , (6)where E is the electric field, H is the magnetic field, u s is the solid displacement, u f is the fluid displacement, µ is the magnetic permeability, ǫ is the dielectricconstant, σ is the conductivity, L is the electro-kinetic mobility parameter, ρ isa linearly combined density of the solid and the fluid, ρ is the density of the porefluid, ρ is the mass coupling coefficient, κ is the fluid flow permeability, η is theviscosity of pore fluid, λ and µ are the Lam´e parameters, q and r are the Biot moduliparameters, τ is the bulk stress tensor, p is the pore pressure, I is the 3 × The research of YG was supported in part by NSFC grant 11671071, JLSTDP 20160520094JHand FRFCU2412017FZ005. The research of YY was partly supported by NSF Grant DMS-1715178, an AMS Simons travel grant, and a start-up fund from Michigan State University. wave propagation. Equations (3)–(6) are Biot’s equations, describing the seismicwave propagation in the porous medium [4, 5].In this work, we focus on the ES imaging. The inverse problem in ES imagingconcerns recovery of the physical parameters in (1)–(6) from boundary measurementof the displacements u s and u f . This problem can be studied in two mutuallyrelevant steps. The first step concerns an inverse source problem for the Biotequations (3)–(6) to recover the coupling term ηκ LE . The second step utilizesthis term as internal measurement to retrieve physical parameters in the Maxwellequations (1)–(2). The second step has been considered in [7, 8].The first step is not well understood mathematically. In an unpublished workof Chen and de Hoop [6], they suggested using Gassmann’s approximation [13] toreduce Biot’s equations to the elastic equation and then applying the result in [26].This approach has the limitation that Gassmann’s approximation is valid only whenthe fluid permeability is small and the wave frequency is sufficiently low. In [3], aninverse source problem was investigated in a different context by assuming accessto internal measurement rather than boundary measurement.Our goal is to demonstrate a general approach to the first step and completethe two-step approach towards the coupled physics inverse problem in ES imaging.We study the inverse source problem from boundary measurement of u s and u f directly for Biot’s equations, without resorting to Gassmann’s approximation orany internal measurement. We derive an explicit reconstruction formula in termsof a Neumann series. Uniqueness and stability of the reconstructed solution areimmediate consequences of the explicit formula.2. Problem formulation and main result
We make two simplifications to the Biot equations (3)–(6). First we ignore theattenuation term ηκ ∂ t u f in (4). This simplification is not essential since the firstorder time derivative can be handled by adding a perturbation argument to ourproof. Second, we assume that the source term ηκ LE is instantaneous in time andhas the form η ( x ) κ ( x ) L ( x ) E ( t, x ) = δ ( t ) f ( x ) , where f is an unknown source function. This assumption is necessary to gainuniqueness in the inverse problem since our boundary measurement has merelythree degrees of freedom and cannot be expected to recover four degrees uniquely.Let us rewrite Biot’s equations (without attenuation) in different forms for sub-sequent analysis. It follows from (5)–(6) that we havediv τ = ∆ µ,λ u s + ∇ ( q div u f ) , ∇ p = −∇ ( q div u s ) − ∇ ( r div u f ) , where the elastic operator ∆ µ,λ is defined by(7) ∆ µ,λ u s := div ( µ ( ∇ u s + ( ∇ u s ) T )) + ∇ ( λ div u s ) , = µ ∆ u s + ( µ + λ ) ∇ (div u s ) + (div u s ) ∇ λ + ( ∇ u s + ( ∇ u s ) T ) ∇ µ. Then the Biot equations (without attenuation) can be written as(8) ( ρ ∂ t u s + ρ ∂ t u f − ∆ µ,λ u s − ∇ ( q div u f ) = 0 ρ ∂ t u s + ρ ∂ t u f − ∇ ( q div u s ) − ∇ ( r div u f ) = δ ( t ) f ( x ) . N INVERSE PROBLEM IN ELECTRO-SEISMIC IMAGING 3
One can further write the system (8) into a matrix equation. Let(9) M = (cid:18) ρ I ρ I ρ I ρ I (cid:19) , P ( D ) = (cid:18) − ∆ µ,λ · −∇ ( q div · ) −∇ ( q div · ) −∇ ( r div · ) (cid:19) . The matrix form of (8) is
M ∂ t u + P ( D ) u = δ ( t ) f . which, by Duhamel’s principle, is equivalent to the initial value problem:(10) M ∂ t u + P ( D ) u = in R × R , u | t =0 = f ( x ) ,∂ t u | t =0 = . Hereafter, a vector with 6 components is written in bold face such as u . Itis often split into two 3-vectors such as u = ( u s , u f ) ∈ R × R . The namesof the two 3-vectors are chosen to be consistent with the splitting for solutionsof the Biot equations. Following this notation, the right hand side is defined as f = ( f s , f f ) := (0 , f ). Hence we will not refer to the specific form any longer butconsider a general f . Well-posedness of the initial value problem (10) can be provedin a similar way as [3, Theorem 1.1].The following hypothesis are crucial and assumed throughout the paper.( H1 ): The functions ρ , ρ , ρ , µ , λ , q , r are bounded from below by aconstant, say c >
0, and have bounded derivatives.( H2 ): ρ ρ − ρ > R . ( H3 ): λr − q > R . Here ( H1 ) states that all the physical parameters are positive and sufficientlysmooth. ( H2 ) ensures that the matrix M is positive definite. ( H3 ) guaranteesthat the energy functional defined in Section 5 is non-negative.In ES imaging, the measurement is the solid displacement u s and the fluid dis-placement u f on the boundary of a domain-of-interest Ω which is a smooth boundedopen subset in R . Introduce the source-to-measurement operator Λ as follows: Λf := u | [0 ,T ] × ∂ Ω where u is the solution of the initial value problem (10) and T > f is supported in the interior of Ω. We areinterested in recovering information on f from the boundary measurement Λf .A quick look at the equations (8) suggests that there is a certain gauge transformfor the recovery of u f . In fact, adding to u f by any vector field g that is divergencefree (i.e, div g = 0) does not affect (8). This is a reflection of the simple fact that u f appears in the equations only in the form of div u f . Taking such gauge intoaccount, we raise the following question which is the central topic of this paper. Inverse source problem:
Suppose that M and P ( D ) are known and satisfythe hypotheses ( H1 )-( H3 ), can one reconstruct the initial source f = ( f s , f f ), com-pactly supported in Ω, from the boundary measurement Λf up to a pair of vectorfields (0 , g ) with g divergence free?We give an affirmative answer to the question, including a reconstruction for-mula. Our main result can be summarized as follows. A rigorous restatement ofthis theorem is given in Theorem 10. Y. GAO, P. LI, AND Y. YANG
Theorem 1.
Under appropriate assumptions, f can be uniquely and stably recon-structed from Λf by a convergent Neumann series, up to a pair of vector fields (0 , g ) with g divergence free. Our proof is based on the modified time reversal method proposed by Stefanovand Uhlmann [23] for the thermo-acoustic tomography (TAT) [18]. TAT is a hybridmodality in medical imaging where optical or electromagnetic waves are exerted totrigger ultrasound wave in tissue through thermo-elastic conversion. Conventionaltime reversal method is known to provide an approximate reconstruction of thesource [16]. Stefanov and Uhlmann improved it and obtained an accurate recon-struction of the source by a Neumann series [23] which was numerically implementedin [9, 20]. This improved time reversal method has since been adapted and gener-alized to many other models [10, 14, 17, 19, 24, 26]. We refer to [2] for a survey onmore hybrid modalities in the context of medical imaging.A slight variation of our approach can be used to reconstruct f ( x ) when thesource term takes the form L ηκ E = δ ′ ( t ) f ( x ). Indeed, it follows from Duhamel’sprinciple that we have M ∂ t u + P ( D ) u = in R × R , u | t =0 = ,∂ t u | t =0 = f ( x ) . If u solves this problem, the ∂ t u solves (10) by [25]. Hence one can first reconstruct ∂ t u using our approach and then integrate to obtain u .The rest of the paper is organized as follows. In Section 3 we convert Biot’sequations (without attenuation) to two hyperbolic systems. These systems areused in Section 4 to prove finite speed of propagation and unique continuationresults. Section 5 is devoted to discussion of function spaces. The main theorem isstated and proved in Section 6.3. Biot’s equations
In this section we transform Biot’s equations into two hyperbolic systems: aprincipally scalar system and a symmetric hyperbolic system.3.1.
The principally scalar system.
Let us start with the principally scalarsystem. Recall the definition of a principally scalar system [12]. For a function a = a ( x ), define the scalar wave operator ✷ a := a∂ t − ∆. Definition 2.
A principally scalar system refers to (11) ✷ a j u j + b j ( t, x, ∇ u ) + c j ( t, x, u ) = f j , j = 1 , . . . , m. Here u = ( u , . . . , u m ) , a j = a j ( x ) ∈ C are real-valued functions, b j and c j arelinear functions with L ∞ loc -coefficients of ∇ u and u . The system is called principally scalar as the principal part of each equation isa scalar wave operator. We shall see in Section 4 that a principally scalar systemscan be uniquely continued towards inside from proper boundary data. This uniquecontinuation property is crucial to obtain Proposition 7.We can write the equation in (10) as a principally scalar system by following theprocedures in [3]. Let ρ ( x ) = ρ ( x ) ρ ( x ) − ρ ( x ). We have ρ > H2 ) and M − = 1 ρ (cid:18) ρ I − ρ I − ρ I ρ I (cid:19) . N INVERSE PROBLEM IN ELECTRO-SEISMIC IMAGING 5
Multiplying the equation in (10) by M − to get ∂ t (cid:18) u s u f (cid:19) + M − P ( D ) (cid:18) u s u f (cid:19) = (cid:18) (cid:19) , where M − P ( D ) = 1 ρ (cid:18) − ρ ∆ µ,λ · + ρ ∇ ( q div · ) − ρ ∇ ( q div · ) + ρ ∇ ( r div · ) ρ ∆ µ,λ · − ρ ∇ ( q div · ) ρ ∇ ( q div · ) − ρ ∇ ( r div · ) . (cid:19) . Now, we write the matrix equation as a system of equations and move all thefirst and zeroth order derivatives to the right hand side to get(12) (cid:26) ∂ t u s − µ ∆ u s − ( µ + λ ) ∇ div u s − q ∇ div u f = P ( u s , u f ) ,∂ t u f + µ ∆ u s − r ∇ (div u f ) − q ∇ (div u s ) = Q ( u s , u f ) , where the new coefficients are µ = ρ − ρ µ, λ = ρ − ( ρ λ − ρ q ) , q = ρ − ( ρ q − ρ r ) ,µ = ρ − ρ µ, r = ρ − ( ρ r − ρ q ) , q = ρ − ( ρ q − ρ ( µ + λ )) . Here and below, the script letters P , Q , R , . . . denote various first-order lineardifferential operators.A straightforward calculation shows that first order derivatives in P ( u s , u f ) and Q ( u s , u f ) appear in terms of ∇ u s , ( ∇ u s ) T , div u s , and div u f . The right hand sidecan be simplified by the substitutions as indicated below.Introduce the substitutions:(13) v s = div u s , v f = div u f , v s = curl u s . Applying div to (12) and utilizing the relationsdiv ∇ u s = ∆ u s = ∇ (div u s ) − curl curl u s , div ( ∇ u s ) T = ∇ (div u s ) , we get (cid:26) ∂ t v s − a ∆ v s − a ∆ v f = R ( v s , v f , u s , v s ) ∂ t v f − a ∆ v s − a ∆ v f = S ( v s , v f , u s , v s );or equivalently ∂ t (cid:18) v s v f (cid:19) − A ( x )∆ (cid:18) v s v f (cid:19) = L ( v s , v f , u s , v s ) , where the coefficient A ( x ) is A ( x ) := (cid:18) a a a a (cid:19) = 1 ρ (cid:18) ρ − ρ − ρ ρ (cid:19) . (cid:18) µ + λ qq r (cid:19) . The hypotheses ( H1 )–( H3 ) ensures that A is a symmetric positive definite matrix.Let a , a be its eigenvalues, then there exists a non-singular matrix Q ( x ) such that( Q − AQ )( x ) = Diag ( a , a )( x ) . Making the change of variable(˜ v s , ˜ v f ) = Q − ( v s , v f )yields that (˜ v s , ˜ v f ) solves (cid:26) ∂ t ˜ v s − a ∆˜ v s = M (˜ v s , ˜ v f , u s , v s ) ,∂ t ˜ v f − a ∆˜ v f = N (˜ v s , ˜ v f , u s , v s ) . Y. GAO, P. LI, AND Y. YANG
Applying curl to the first equation of (12) and using curl ∇ u = 0 gives ∂ t v s − µ ∆ v s = T (˜ v s , ˜ v f , u s , v s ) . The first equation of (12) can be written as ∂ t u s − µ ∆ u s = U (˜ v s , ˜ v f , u s , v s ) . Thus we obtain the following principally scalar system in the variables (˜ v s , ˜ v f , u s , v s ):(14) ✷ a ˜ v s − a M (˜ v s , ˜ v f , u s , v s ) = 0 , ✷ a ˜ v f − a N (˜ v s , ˜ v f , u s , v s ) = 0 , ✷ µ u s − µ T (˜ v s , ˜ v f , u s , v s ) = 0 , ✷ µ u f − µ U (˜ v s , ˜ v f , u s , v s ) = 0 . This is the desired principally scalar system. Note that a , a , µ are smooth andstrictly positive by ( H1 ), hence their reciprocals exist and are smooth as well.Another observation is that (˜ v s , ˜ v f ) = 0 if and only if ( v s , v f ) = 0.3.2. The symmetric hyperbolic system.
We proceed to write the principallyscalar system (14) as a first order symmetric hyperbolic system [22], which will beexploited to show that the Biot equations have finite speed of propagation. Werestrict to the case where the coefficients are time-independent matrices.For x ∈ R , let A ( x ) , . . . , A ( x ) be matrix-valued functions. Denote by L thepartial differential operator(15) L ( x, ∂ t , ∂ x ) = A ( x ) ∂ t + X j =1 A j ( x ) ∂ x j + A ( x ) , where the coefficient matrices A , . . . , A are assumed to have uniformly boundedderivatives, i.e.,sup x ∈ R k ∂ αx ( A ( x ) , . . . , A ( x )) k < ∞ for any multi-index α. Definition 3. L is called symmetric hyperbolic if the following conditions hold:(i) A ( x ) , . . . , A ( x ) are symmetric;(ii) A is strictly positive, i.e., there is a constant C > such that for all x ∈ R , A ( x ) ≥ CI, where I is the identity matrix. Our goal is to convert the principally scalar system (14) to a symmetric hyper-bolic system. We begin with a scalar wave equation to demonstrate the procedures.The principally scalar system is treated afterwards.Let v ( x ) be a scalar function defined in R and satisfy the equation(16) ∂ t v − c ( x )∆ v = L ( ∂ x v, ∂ x v, ∂ x v, v ) , where L ( ∂ x v, ∂ x v, ∂ x v, v ) is linear in each derivative. We define a vector V ( v ): V ( v ) = ( V , V , V , V , V ) := ( ∂ t v, ∂ x v, ∂ x v, ∂ x v, v ) . N INVERSE PROBLEM IN ELECTRO-SEISMIC IMAGING 7
The scalar equation (16) can be written in terms of V , . . . , V as follows: ∂ t V − a ∂ x V − a ∂ x V − a ∂ x V − L ( V , . . . , V ) = 0 ,a ∂ t V − a ∂ x V = 0 ,a ∂ t V − a ∂ x V = 0 ,a ∂ t V − a ∂ x V = 0 ,a ∂ t V − a V = 0 . In the matrix form, this system reads B ( x ) ∂ t V + B ( x ) ∂ x V + B ( x ) ∂ x V + B ( x ) ∂ x V + B ( x ) V = 0 , where B , . . . , B are 5 × B = Diag(1 , a , a , a , a ), B is determined by the concrete form of L , and the other matrices are B i ( x ) := (cid:26) − a for (1 , i + 1) and ( i + 1 ,
1) entry0 for other entries , i = 1 , , . Note this is a symmetric hyperbolic system since B , . . . , B are symmetric matricesand B is strictly positive.Now we turn to the principally scalar system (14). As each equation in thesystem takes the form (16), we define a vector U , which has 40 components and isobtained by juxtaposing V ( v ) with v replaced successively by the components of(˜ v s , ˜ v f , u s , v s ), i.e., U = ( U , . . . , U ) := ( V (˜ v s ) , V (˜ v f ) , V ( u s1 ) , V ( u s2 ) , V ( u s3 ) , V ( v s1 ) , V ( v s2 ) , V ( v s3 )) . Since the principal part of each equation in (14) is uncoupled, one can write (14)as a first order system in a similar manner:(17) A ( x ) ∂ t U + A ( x ) ∂ x U + A ( x ) ∂ x U + A ( x ) ∂ x U + A ( x ) U = 0 . Here each A i is a 40 ×
40 matrix, A represents the diagonal matrixDiag ((1 , a , a , a , a ) , (1 , a , a , a , a ) , (1 , µ , µ , µ , µ ) , . . . , (1 , µ , µ , µ , µ ))which is strictly positive, and A , A , A are symmetric. This is the desired sym-metric hyperbolic system that is equivalent to the principally scalar system (14)and the Biot equations in (10).4. Finite speed of propagation and unique continuation
We derive some results for the Biot equations regarding finite speed of propaga-tion and unique continuation. It is crucial to have the symmetric hyperbolic system(17) and the principally scalar system (14).4.1.
Finite speed of propagation.
Let L be the symmetric hyperbolic operatordefined as in (15). For any ξ ∈ R \{ } , letΛ( ξ ) = inf { ℓ : A ( x ) − X j =1 A j ( x ) ξ j A ( x ) − ≤ ℓI } , which is the smallest upper bound of the eigenvalues for A − (cid:16)P j =1 A j ξ j (cid:17) A − .We state the following result, which shows that the solution of a symmetric hyper-bolic system has finite speed of propagation. Y. GAO, P. LI, AND Y. YANG
Proposition 4. [22, Theorem 2.3.2] Suppose s ∈ R and u ∈ C ([0 , ∞ ); H s ( R )) satisfies a symmetric hyperbolic system L u = 0 . Ifsupp u (0 , x ) ⊂ { x ∈ R : x · ξ ≤ } , then for t ≥ , supp u ( t, x ) ⊂ { x ∈ R : x · ξ ≤ Λ( ξ ) t } . In particular, set c max := max | ξ | =1 Λ( ξ ) , ifsupp u (0 , x ) ⊂ { x ∈ R : | x | ≤ R } , then supp u ( t, x ) ⊂ { x ∈ R : | x | ≤ R + c max | t |} . Here | ξ | denotes the Euclidean norm of ξ . If u is a vector, supp u stands for theunion of the supports of its components.Using this proposition, one can deduce finite speed of propagation for the solu-tions of Biot’s equation. Let c max be the number in Proposition 4 with L definedon the left hand side of (17). Corollary 5.
Let u = ( u s , u f ) solve M ∂ t u + P ( D ) u = . Ifsupp ( u s (0 , · ) , div u f (0 , · )) ⊂ { x ∈ R : | x | ≤ r } , then supp ( u s ( t, · ) , div u f ( t, · )) ⊂ { x ∈ R : | x | ≤ r + c max | t |} . Proof.
Given ( u s , u f ), one can construct the functions (˜ v s , ˜ v f , u s , v s ) as in (13),and further the solution U of (17). The proof is completed by applying Proposition4 to the symmetric hyperbolic system (17). (cid:3) Unique continuation.
A principally scalar system satisfies certain uniquecontinuation property [12], which will be used as an intermediate step towards themain theorem. Let T ′ > D ⊂ R be a C -domaincontaining the origin. Denote by B (0; R ) = { x ∈ R : | x | < R } the ball of radius R centered at the origin. We state the following unique continuation result. Proposition 6. [12, Corollary 3.5] Let u = ( u , . . . , u m ) be a solution of a generalprincipally scalar system (11) . Suppose that there exists θ > such that D ⊂ B (0; θT ′ ) and the coefficients a j in (11) satisfy the constraints: θ a j ( a j + a − j | t ∇ a j | ) < a j + x · ∇ a j in ( − T ′ , T ′ ) × Dθ a j ≤ in D. Then u = ∂ ν u = 0 on ( − T ′ , T ′ ) × ∂D implies u = 0 on { ( t, x ) ∈ ( − T ′ , T ′ ) × D : | x | > θt } . The inequality constraints in the theorem justify the pseudo-convexity of a cer-tain phase function with respect to the wave operators ✷ a j [12]. Observe that if a j are positive constants, the above constraints reduce to θ a j < D . In this case,any θ > { √ a j : j = 1 , · · · , m } fulfills the inequalities.Next proposition is the main result of this section. Choose R > ⊂ B (0; R ). Recall that c max is defined in Corollary 4. N INVERSE PROBLEM IN ELECTRO-SEISMIC IMAGING 9
Proposition 7.
Let u = ( u s , u f ) be the solution of the forward problem (10) with f compactly supported in Ω . Suppose there exists θ > such that B (0; R + c max T ) ⊂ B (0; θT ) and the following inequality constraints hold when a is replaced by a ( x ) , a ( x ) , µ ( x ) respectively: θ a ( a + a − | t ∇ a | ) < a + x · ∇ a, in ( − T , T ) × B (0; R + c max T ) θ a ≤ in B (0; R + c max T ) . If u s ( T, x ) = div u f ( T, x ) = 0 for x ∈ R \ Ω , then u s (0 , x ) = div u f (0 , x ) = 0 for x ∈ Ω , x = 0 , Proof.
In view of the initial conditions in (10), we can extend the solution u as aneven function of t to ( − T, T ) × Ω. This extension is denoted by u again.Given u s ( T, x ) = div u f ( T, x ) = 0 for x ∈ R \ Ω, it follows from Corollary 5 that u s ( t, x ) = div u f ( t, x ) = 0 in { ( t, x ) : | x | ≥ R + c max | t − T |} . As u is an even function of t , we also have u s ( t, x ) = div u f ( t, x ) = 0 in { ( t, x ) : | x | ≥ R + c max | t + T |} . On the other hand, the fact that f has compact support in Ω implies u s (0 , x ) =div u f (0 , x ) = 0 for x ∈ R \ Ω. Hence by Corollary 5, we get u s ( t, x ) = div u f ( t, x ) = 0 in { ( t, x ) : | x | ≥ R + c max | t |} . Combing the above equations gives u s ( t, x ) = div u f ( t, x ) = 0 on { ( t, x ) : | x | = R + c max T , − T < t < T } . Using Proposition 6 with the choice T ′ := T , D := B (0 , R + c max T ) yields u s (0 , x ) = div u f (0 , x ) = 0 for x ∈ Ω , x = 0 , which completes the proof. (cid:3) Energy Conservation
Define a Sobolev space H (div ; Ω) := { v f ∈ ( L (Ω)) : div v f ∈ L (Ω) } , which is equipped with the norm k v f k H (div ;Ω) := k v f k L (Ω) + k div v f k L (Ω) . Consider the space V := { v = ( v s , v f ) ∈ ( H (Ω)) × H (div ; Ω) } , which has the norm k ( v s , v f ) k V := k v s k H (Ω) + k v f k H (div ;Ω) . Introduce a symmetric bilinear form on V : B ( v , w ) = Z Ω (cid:2) λ div v s · div w s + 2 µ ( ǫ ( v s ) : ǫ ( w s )) + r div v f · div w f (cid:3) dx + Z Ω q (cid:2) div v f · div w s + div v s · div w f (cid:3) dx, where ǫ ( v s ) = ( ∇ v s + ( ∇ v s ) T ), A : B = tr( AB ⊤ ) is the Frobenius inner productof matrices A and B .It is clear to note that B ( v , v ) ≥ v ∈ V in view of ( H3 ); moreover, B ( v , v ) = 0 if and only if v s = 0 and div v f = 0. Thus B induces a semi-norm k v k B := ( v , v ) B . We relate the bilinear form B to the differential operator P ( D ). Lemma 8.
Suppose v , w ∈ V with v | ∂ Ω = or w | ∂ Ω = , then B ( v , w ) = ( P ( D ) v , w ) L (Ω) = ( v , P ( D ) w ) L (Ω) . Proof.
By symmetry, we only need to prove the first equality with the assumptionthat w | ∂ Ω = . Recalling w = ( w s , w f ) T and P ( D ) v = ( − ∆ µ,λ v s − ∇ ( q div v f ) , −∇ ( q div v s ) − ∇ ( r div v f )) T , we have ( P v , w ) L (Ω) = Z Ω (cid:2) − ∆ µ,λ v s · w s − ∇ ( q div v f ) · w s − ∇ ( q div v s ) · w f − ∇ ( r div v f ) · w f (cid:3) dx. (18)To deal with the first integrand, we expand ∆ µ,λ using (7) to have Z Ω [ − ∆ µ,λ v s · w s ] dx = − Z Ω div ( µǫ ( v s )) · w s dx − Z Ω ∇ ( λ div v s ) · w s dx. We claim that div ( µǫ ( v s )) · w s = div ( µǫ ( v s )) w s − µǫ ( v s ) : ǫ ( w s ) . To justify this, we write ǫ ( v s ) := ( ǫ , ǫ , ǫ ) T , where ǫ , ǫ , ǫ are the three rows;write w s = ( w s1 , w s2 , w s3 ) T . Thendiv ( µǫ ( v s )) · w s = X j =1 div ( µǫ j ) w sj = X j =1 (cid:2) div (cid:0) µ w sj ǫ j (cid:1) − µǫ j · ∇ w sj (cid:3) = div ( µǫ ( v s )) w s − µǫ ( v s ) : ǫ ( w s ) . Using the integration by parts yields Z Ω [ − ∆ µ,λ v s · w s ] dx = Z Ω [2 µǫ ( v s ) : ǫ ( w s ) + λ div v s · div w s ] dx + Z ∂ Ω [ − µ ( ǫ ( v s ) w s ) − ( λ div v s ) w s ] · ν dx. (19)The boundary term vanishes owing to the compact support of w .The remaining three integrands in (18) can be treated using the standard inte-gration by parts: Z Ω (cid:2) −∇ ( q div v f ) · w s − ∇ ( q div v s ) · w f − ∇ ( r div v f ) · w f (cid:3) dx = Z Ω (cid:2) q div v f · div w s + q div v s · div w f + r div v f · div w f (cid:3) dx + Z ∂ Ω (cid:2) − ( q div v f ) w s − ( q div v s ) w f − ( r div v f ) w f (cid:3) · ν dx. (20)The boundary term again vanishes. This completes the proof. (cid:3) N INVERSE PROBLEM IN ELECTRO-SEISMIC IMAGING 11
Let L (Ω; M ) be a weighted L -space with measure M ( x ) dx where M ( x ) isthe positive definite matrix defined in (9); in other words, for any v ∈ ( L (Ω)) , k v k L (Ω; M ) := ( v , M v ) L . The space L ( R ; M ) and the norm k · k L ( R ; M ) isdefined similarly with the domain of integration replaced by R .Given a time-dependent function u ( t, x ), define its total energy over the domainΩ at time t : E Ω ( t, u ) = k ∂ t u ( t, · ) k L (Ω; M ) + k u ( t, · ) k B = Z Ω (cid:0) M ( x ) ∂ t u · ∂ t u + λ | div u s | + 2 µ | ǫ ( u s ) | + r | div u f | + 2 q (div u f )(div u s ) (cid:1) dx. This quantity is conservative on a bounded domain Ω for the solution of Biot’sequations when imposed with appropriate boundary conditions.
Lemma 9.
Let u satisfy Biot’s equations with zero Dirichlet boundary condition: (cid:26) M ∂ t u + P ( D ) u = in (0 , T ) × Ω , u | [0 ,T ] × ∂ Ω = , then E Ω ( t, u ) = E Ω (0 , u ) 0 ≤ t ≤ T. Proof.
We briefly sketch the proof since it is similar to that of Lemma 8. Takingthe inner product of the equation with ∂ t u = ( ∂ t u s , ∂ t u f ) T , we obtain0 = ( ∂ t u , ∂ t u ) L (Ω; M ) + ( P ( D ) u , ∂ t u ) L (Ω) = ( ∂ t u , ∂ t u ) L (Ω; M ) + B ( u , ∂ t u ) + b.t. = 2 ddt E Ω ( t, u ) + b.t.. Here the second equality is justified by Lemma 8; “ b.t. ” represents the arisingboundary terms, which are the boundary integrals in (19) and (20) with v and w replaced by u and ∂ t u respectively. The zero Dirichlet boundary conditionannihilates b.t. , which completes the proof of conservation of energy. (cid:3) The proof verifies the well known fact that zero Dirichlet boundary conditionpreserves energy. In fact each of the following boundary conditions(i) u s ( t, x ) = 0 and u f · ν ( t, x ) = 0 on (0 , T ) × ∂ Ω;(ii) u s · ν ( t, x ) = 0 and u f ( t, x ) = 0 on (0 , T ) × ∂ Ω;(iii) u s · ν ( t, x ) = 0 and u f · ν ( t, x ) = 0 on (0 , T ) × ∂ Ω.is energy preserving as well, since each of them is sufficient to annihilate the bound-ary term “ b.t. ” in the proof.Let u be the solution of the direct problem (10). We can also consider the energyover the entire space E R ( t, u ). This global energy is conservative as well, i.e., E R ( t, u ) = E R (0 , u ) = k f k B ≤ t ≤ T. To show this, one can take a large ball B (0; R ) so that the solution u ( t, · ) issupported inside B (0; R ) for any t ∈ [0 , T ]. Then Lemma 9 is applicable since u | [0 ,T ] × ∂B (0; R ) = and E R ( t, u ) = E Ω ( t, u ) for any t ∈ [0 , T ]. Main Theorem
Let v be the solution of(21) M ∂ t v + P ( D ) v = 0 in (0 , T ) × Ω , v | [0 ,T ] × ∂ Ω = h, v ( T, · ) = φ ,∂ t v ( T, · ) = , where φ is the function satisfying P ( D ) φ = 0 , φ | ∂ Ω = h ( T, · ) . The solution φ exists since the analysis in Section 3 manifests that P ( D ) φ = 0 canbe transformed into an elliptic system, which is the time-independent counterpartof (14). Define the time reversal operator Ah := v (0 , · ) . Now we take h = Λf as the boundary measurement and expect A Λf to be a rea-sonable approximation of f . The rationale, from a microlocal viewpoint, is that thehyperbolic operator in the forward problem propagates microlocal singularities of f to ∂ Ω, while the time-reversal process tends to send back these singularities. Thissuggests a possible reconstruction of f , as least on the level of principal symbols.The microlocal viewpoint also suggests the necessity of an additional assumptionto make sure that all the microlocal singularities of u , the solution of (10), are nottrapped, i.e., all of them are able to reach ∂ Ω in a finite time. This leads to thenon-trapping condition: there exists a maximal escaping time T ( M, P ( D ) , Ω) > M , P ( D ) and the geometry of Ω, such that all the (microlocal)singularities of u are out of Ω whenever t > T ( M, P ( D ) , Ω); in other words, u ( t, x )is smooth for x ∈ Ω whenever t > T ( M, P ( D ) , Ω). Violation of the non-trappingcondition would cause loss of singularities in the measured data Λf , leading tounstable reconstruction of f .Now we are in the position to state and prove the main theorem. We show that A is the inverse of Λ up to a compact operator and the compact operator becomesa contraction on a suitable function space.Recall that k · k B is merely a semi-norm on V : it is not positive definite since k w k B = 0 implies only w s = 0 and div w f = 0. We can take V modulo the closedsubspace { w ∈ V : w s = 0 , div w f = 0 } to make it a genuine norm. Denote by B the quotient space, then ( B, k · k B ) is a Banach space and k w k B = 0 implies w = 0in B . We also project Ah = v (0 , · ) to this quotient space and abuse the notationto call A composed with the canonical projection as A . So A maps into ( B, k · k B ). Theorem 10.
Let Ω be non-trapping and T > T ( P ( D ) , Ω) . Suppose that thehypotheses ( H1 )-( H3 ) and the assumption of Proposition 7 are satisfied. Then K := I − A Λ is compact and contractive on B in the sense that k K k B → B < . Asa consequence, I − K is invertible on B and f = ∞ X j =0 K j A Λf = A Λf + KA Λf + K A Λf + . . . . in B. Proof.
The proof is divided into two claims. We first show the inequality k K k B → B ≤
1, and then prove by a contra-positive argument that the inequality is strict. Givena time-dependent function g ( t, x ), we abbreviate g ( t ) for the spatial function g ( t, · ). N INVERSE PROBLEM IN ELECTRO-SEISMIC IMAGING 13
Claim 1: k K f k B < k f k B unless f = in B . Let us give another representation of K f . Let u be the solution of (10); let v bethe solution of (21) with h replaced by Λf . Denote w := u − v , then w satisfies(22) M ∂ t w + P ( D ) w = in (0 , T ) × Ω , w | [0 ,T ] × ∂ Ω = , w ( T ) = u ( T ) − φ ,∂ t w ( T ) = ∂ t u ( T ) . Moreover, we have(23) K f = f − A Λf = u (0) − v (0) = w (0) . On the other hand, it is clear to note that ( u ( T ) − φ ) | ∂ Ω = by the constructionof φ . It follows from Lemma 8 that( u ( T ) − φ , φ ) B = ( u ( T ) − φ , P ( D ) φ ) L (Ω) = 0 , which gives k u ( T ) − φ k B = k u ( T ) k B − k φ k B . It is easy to verify that E Ω ( T, w ) = k ∂ t w ( T ) k L (Ω; M ) + k w ( T ) k B = k ∂ t u ( T ) k L (Ω; M ) + k u ( T ) − φ k B = k ∂ t u ( T ) k L (Ω; M ) + k u ( T ) k B − k φ k B = E Ω ( T, u ) − k φ k B ≤ E Ω ( T, u ) . Combining Lemma 9 and conservation of energy in R yields E Ω (0 , w ) = E Ω ( T, w ) ≤ E Ω ( T, u ) ≤ E R ( T, u ) = E R (0 , u ) = k f k B . Thus we have from (23) that k K f k B = k w (0) k B ≤ E Ω (0 , w ) ≤ k f k B . Suppose the equality holds for some f ∈ B , then all the above inequalities becomeequalities. In particular E Ω ( T, u ) = E R ( T, u ), which implies u s ( T, x ) = div u f ( T, x ) = 0 for x ∈ R \ Ω . By Proposition 7, we obtain f s ( x ) = 0 and div f f ( x ) = 0 for x = 0. Changing thevalue at the single point x = 0 does not affect a function in B . This completes theproof of Claim 1. Claim 2: K : B → B is compact and k K k B → B < . Claim 1 alone implies k K k B → B ≤
1. To prove the strict inequality, we show K is a compact operator on B . The spectrum of a compact operator consistsof countably many eigenvalues which may accumulate only at 0. Since Claim 1excludes eigenvalues of modulus 1, the spectral radius of K must be strictly lessthan 1, proving that k K k B → B < K . Using the representation (23), we decom-pose K into composition of bounded operators: f π ∗ ( f , ) U ( u | t = T − φ , ∂ t u | t = T )= ( w ( T ) , ∂ t w ( T )) U ( w (0) , ∂ t w (0)) π w (0) . Here π : ( f , g ) f is the natural projection onto the first component; π ∗ is itsadjoint; U is the solution operator of the forward problem (10), mapping the state t = 0 to the state t = T ; and U is the solution operator of (22) sending t = T to t = 0. These are all bounded operators.Consider U : ( f , ) ( u | t = T − φ , ∂ t u | t = T ) . In view of the assumption that
T > T ( P ( D ) , Ω), all the microlocal singularitiesof f have escaped from Ω at the moment T , hence ( u | t = T , ∂ t u | t = T ) is a pair ofsmooth functions. On the other hand, the function φ , as a solution of the ellipticequations P ( D ) φ = 0, is smooth by elliptic regularity. We conclude that U is asmoothing operator, hence compact. This means K is compact as well, since it isthe composition of U with other bounded operators.We know that K is a contraction on B , ( I − K ) − exists as a bounded operator.Applying ( I − K ) − to the identity ( I − K ) f = A Λf and expanding it in termsof Neumann series, we obtain the reconstruction formula in the statement of thetheorem. (cid:3) The following stability estimate shows that the faster the energy escapes fromΩ, the faster the convergence of the Neumann series is.
Corollary 11.
Under the assumption of Theorem 10, the following stability esti-mate holds k K f k B ≤ (cid:18) E Ω ( T, u ) E Ω (0 , u ) (cid:19) k f k B for f = 0 in B. Proof.
A simple calculation yields that k K f k B k f k B = k w (0) k B E Ω (0 , u ) ≤ E Ω (0 , w ) E Ω (0 , u ) ≤ E Ω ( T, u ) E Ω (0 , u ) . (cid:3) Acknowledgement . The third author would like to thank Prof. Plamen Ste-fanov for many helpful discussions and Prof. Gunther Uhlmann for bringing thereference [17]. He is also grateful to the Institute of Computational and Exper-imenta Research in Mathematics (ICERM), where part of this research was con-ducted during his participation in the program “Mathematical Challenges in Radarand Seismic Imaging” in 2017.
References [1] A. H. Araji, A. Revil, A. Jardani, B. J. Minsley, and M. Karaouslis,
Imaging with cross-holeseismoelectric tomography , Geophys. J. Int., , (2012), 1285-1302.[2] G. Bal,
Hybrid inverse problems and internal functionals , Inside Out, Cambridge UniversityPress, Cambridge, UK, G. Uhlmann, Editor, 2012.[3] M. Bellassoued, M. Yamamoto,
Carleman estimates and inverse source problem of Biotsequations describing wave propagation in porous media , Inverse Problems, , (2013), 115002.[4] M. A. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid. I. Lowfrequency range , J. Acoust. Soc. Am. , (1956), 168-178.[5] M. A. Biot, Theory of propagation of elastic waves in a fluid-saturated porous solid. II. Highfrequency range , J. Acoust. Soc. Am. , (1956), 179-191.[6] J. Chen, The inverse problem of electroseismic conversion (a talk at BIRS) , (2015), [7] J. Chen and Y. Yang,
Inverse problem of electro-seismic conversion , Inverse Problems, (2013), 115006.[8] J. Chen and M. V. de Hoop, Inverse problem of electroseismic conversion I: Inversion ofMaxwells equations with internal data , preprint, (2014), arXiv:1406.0367.
N INVERSE PROBLEM IN ELECTRO-SEISMIC IMAGING 15 [9] E. Chung, C. Y. Lam, and J. Qian,
A Neumann series based method for photoacoustictomography on irregular domains , Contemporary Mathematics, (2014), 89-104.[10] M. V. de Hoop and J. Tittelfitz,
An inverse source problem for a variable speed wave equationwith discrete-in-time sources , Inverse Problems, (7) (2015), 075007.[11] J. C. Dupuis, K. E. Butler, A. W. Kepic, and B. D. Harris, Anatomy of a seismoelectricconversion: measurements and conceptual modeling in boreholes penetrating a sandy aquifer ,J. Geophys. Res., , B10306, (2009) DOI: 10.1029/2008JB005939.[12] M. Eller, V. Isakov, G. Nakamura, and D. Tataru,
Uniqueness and stability in the Cauchyproblem for Maxwell’s and elasticity systems , College de France Seminar, 14, Studies in Math.Appl., Vol.31, (2002), North-Holland, Elsevier Science, 329-349.[13] F. Gassmann, ¨Uber die elastizit¨at por¨oser Medien , Vierteljahrschrift der NaturforschendenGessellschaft in Zurich 96, (1951), 1856-1999.[14] A. Homan,
Multi-wave imaging in attenuating media , Inverse Problems Imaging, (2013),1235-1250.[15] H. Hu, W. Guan, J. Harris, Theoretical simulation of electroacoustic borehole logging inuidsaturated prous formation , J. Acoust. Soc. Am., , (2007), 135-145.[16] Y. Hristova, P. Kuchment, and L. Nguyen,
Reconstruction and time reversal in thermoacous-tic tomography in acoustically homogeneous and inhomogeneous media , Inverse Problems, (2008), 055006.[17] V. Katsnelson and L. Nguyen, Time reversal method for thermoacoustic acous-tic tomography in elastic media: convergence with sharp observation time ,https://arxiv.org/pdf/1707.05822.pdf, (2017).[18] P. Kuchment and L. Kunyansky,
Mathematics of thermoacoustic tomography , European Jour-nal of Applied Mathematics, (2) (2008), 91224.[19] B. Palacios, Reconstruction for multi-wave imaging in attenuating media with large dampingcoefficient , Inverse Problems, (12) (2016), 125008.[20] J. Qian, P. Stefanov, G. Uhlmann, and H. Zhao, An efficient Neumann series-based algorithmfor thermoacoustic and photoacoustic tomography with variable sound speed , SIAM J. ImagingSci., (2011), 850–883.[21] S. R. Pride, Governing equations for the coupled electro-magnetics and acoustics of porousmedia , Phys. Rev. B, (1994), 15678-15696.[22] J. Rauch, Hyperbolic partial differential equations and geometric optics , Graduate Studies inMathematics, vol. 133, American Mathematical Society, Providence, RI, (2012), MR 2918544.[23] P. Stefanov and G. Uhlmann,
Thermoacoustic tomography with variable sound speed , InverseProblems, (2009), 075011.[24] P. Stefanov and G. Uhlmann, Thermoacoustic tomography arising in brain imaging , InverseProblems, (2011), 045004.[25] P. Stefanov and G. Uhlmann, Instability of the linearized problem in multiwave tomography ofrecovery both the source and the speed , Inverse Problems and Imaging, (4) (2013), 1367-1377.[26] J. Tittelfitz, Thermoacoustic tomography in elastic media , Inverse Problems, (2012),055004.[27] Z. Zhu and M. N. Toks¨oz, Seismo-electric and seismo-magnetic measurements in fracturedbore-hole models , Geophysics, vol. 70 no. 4, (2005) F45-F51.
School of Mathematics and Statistics, Center for Mathematics and InterdisciplinarySciences, Northeast Normal University, Changchun, Jilin 130024, P. R. China
E-mail address : [email protected] Department of Mathematics, Purdue University, West Lafayette, IN 47907, USA.
E-mail address : [email protected] Department of Computational Mathematics, Science and Engineering, Michigan StateUniversity, East Lansing, MI 48824, USA.
E-mail address ::