Convexification method for a coefficient inverse problem and its performance for experimental backscatter data for buried targets
CConvexification method for a coefficient inverseproblem and its performance for experimentalbackscatter data for buried targets ∗ Michael V. Klibanov † Aleksandr E. Kolesov †‡ Dinh-Liem Nguyen § Abstract
We present in this paper a novel numerical reconstruction method forsolving a 3D coefficient inverse problem with scattering data generatedby a single direction of the incident plane wave. This inverse problemis well-known to be a highly nonlinear and ill-posed problem. Therefore,optimization-based reconstruction methods for solving this problem wouldtypically suffer from the local-minima trapping and require strong a prioriinformation of the solution. To avoid these problems, in our numericalmethod, we aim to construct a cost functional with a globally strictlyconvex property, whose minimizer can provide a good approximation forthe exact solution of the inverse problem. The key ingredients for theconstruction of such functional are an integro-differential formulation ofthe inverse problem and a Carleman weight function. Under a (partial)finite difference approximation, the global strict convexity is proven usingthe tool of Carleman estimates. The global convergence of the gradientprojection method to the exact solution is proven as well. We demon-strate the efficiency of our reconstruction method via a numerical studyof experimental backscatter data for buried objects.
Keywords.
Carleman weight function, Carleman estimates, reconstructionmethod, convexification, global convergence, coefficient inverse problem, exper-imental data
AMS subject classification.
We develop in this paper a novel numerical method for solving a coefficientinverse problem (CIP) for the 3D Helmholtz equation with scattering data gen- ∗ Supported by US Army Research Laboratory and US Army Research Office grantW911NF-15-1-0233 and by the Office of Naval Research grant N00014-15-1-2330. In addition,the work of Kolesov A.E. was partially supported by Mega-grant of the Russian FederationGovernment (N14.Y26.31.0013) and RFBR (N17-01-00689A) † Department of Mathematics & Statistics, University of North Carolina at Charlotte, Char-lotte, NC 28223, USA ([email protected]) ‡ Institute of Mathematics and Information Science, North-Eastern Federal University,Yakutsk, Russia ([email protected]) § Department of Mathematics, Kansas State University, Manhattan, KS 66506 ([email protected]) a r X i v : . [ m a t h . NA ] M a y rated by a single direction of the incident plane wave at multiple frequencies.More precisely, the goal of this CIP is to recover a coefficient in the Helmholtzequation from boundary measurements of its solutions for a single direction ofthe incident plane wave at multiple frequencies.This CIP arises in a wide range of applications including non-destructivetesting, detection of explosives, medical imaging, geophysics, etc. It is alsowell-known that any CIP is a highly nonlinear and ill-posed problem caus-ing substantial challenges in the design of numerical algorithms for solving it.Optimization-based reconstruction methods can be considered as the most stud-ied approach for solving CIPs in general. However, these methods suffer fromthe fact that they might converge to a local minimum, which is not the truesolution of the CIP. Moreover, these methods typically require strong a prioriinformation of the solution, which is not always available in practice.The goal of the so-called globally convergent method (GCM), recently devel-oped by the first author and coauthors (see e.g. [6]) is to overcome the drawbacksmentioned above when solving CIPs. This method aims to provide a point ina sufficiently small neighborhood of the true solution of the CIP without anyadvanced knowledge of this neighborhood. The size of this neighborhood shoulddepend only on approximation errors and the level of noise in the data.The numerical method we develop in this paper can be considered as thesecond type of GCMs, which has certain advantages compared with the firsttype of GCMs in [6,35,36]. More precisely, we do not impose in the convergenceanalysis here the assumption on a small interval of wavenumbers. Neither we donot iterate here with respect to the so-called tail functions. The combination ofthe latter two features with the globally strictly convex cost functional (below)are the main improvements of the convexification over the first type of globallyconvergent methods.This second type of GCMs is also called convexification methods, whichwas studied for the 1D case in [27]. The present work can be considered as ageneralization to the 3D case of the cited 1D version. Convexification methodsare based on the minimization of the weighted cost functional with a Carlemanweight function (CWF) in it. The CWF is the function which is involved inthe Carleman estimate for the corresponding PDE operator. The CWF can bechosen in such a way that the cost functional becomes strictly convex. Notethat the majority of known numerical methods of solutions of nonlinear ill-posed problems minimize conventional least squares cost functionals (see, e.g.[11, 13, 14]), which are usually non convex and have multiple local minima andravines, see, e.g. [39] for a good numerical example of multiple local minima.We work in this paper with a semidiscrete version of the convexification,which is more realistic for computations than continuous versions used in previ-ous works on the convexification of the first author with coauthors [3, 20, 24, 25,27, 29]. “Semidiscrete” means that we develop the theory for the case when thedifferential operator we work with is written in finite differences with respect totwo out of three spatial variables and in the continuous form with respect to thethird variable. We impose a computationally reasonable assumption that thegrid step size does not tend to zero (unlike the case of some forward problems).The fully discrete case, i.e. when derivatives with all three variables are writtenin finite differences, is not investigated yet. Indeed, it is well known that thiscase is quite a complicated one for ill-posed problems for PDEs, especially innonlinear cases, such as we work with. There are known only a few results for2he fully discrete cases of linear ill-posed problems, see, e.g. [8, 19]. We alsorefer to the recent publication [29] of the first two authors about a 3D versionof the convexification method. In [29] convexification was numerically tested onsome computationally simulated data. This is unlike the current paper in whichtesting is done for a significantly more challenging case of experimental data.The theory in [29] is developed for the continuous case. Although the idea ofthe semidiscrete version is briefly outlined in [29], corresponding theorems areneither formulated nor proved there, unlike the current paper.We point out that the CIP considered in this paper is also called a inversescattering problem in some contexts. There is a vast literature on both the-oretical and numerical studies on this inverse problem and its variations, see,e.g. [1, 2, 10, 12, 15–18, 30–34]. These cited papers have considered the cases ofmultiple measurements and/or shape reconstructions. We recall that we con-sider in this paper the CIP with a single measurement which is both differentand more challenging than the configurations considered in those cited papers.In the next section, we provide a statement of the forward and inverse prob-lems. In Section 3 we present an integro-differential equation formulation ofthe CIP. Section 4 involves the approximation of the tail function which is animportant component in the integro-differential equation. We introduce in Sec-tion 5 the partial finite difference approximation and related function spaces forthe integro-differential formulation. We describe in Section 6 the weighted costfunctional with the Carleman weight function in it. Section 7 is dedicated to thetheoretical analysis, including a Carleman estimate and proofs of global strictconvexity of that functional as well as convergence results for the optimizationproblem. Finally, our numerical study is presented in Section 8. Let x = ( x, y, z ) ∈ R and consider positive numbers b > d >
0. For theconvenience for our numerical study (Section 8), we define from the beginningthe domain of interest Ω and the backscatter part Γ of its boundary asΩ = { ( x, y, z ) : | x | , | y | < b, z ∈ ( − ξ, d ) } , Γ = { ( x, y, z ) : | x | , | y | < b, z = − ξ } . (2.1)Let the function c ( x ) be the spatially distributed dielectric constant and k bethe wavenumber. We consider the forward scattering problem for the Helmholtzequation: ∆ u + k c ( x ) u = 0 , x ∈ R , (2.2) u ( x , k ) = u s ( x , k ) + u i ( x , k ) , (2.3)lim r →∞ r ( ∂u s /∂r − iku s ) = 0 , r = | x | , (2.4)where u ( x , k ) is the total wave, u i ( x , k ) is the incident wave and u s ( x , k ) is thescattered wave satisfying the Sommerfeld radiation condition. This conditionmeans that the scattered field behaves like a outgoing spherical wave far awayfrom the scattering medium.Here we consider u i ( x , k ) as the incident plane wave propagating along thepositive direction of the z − axis: u i ( x , k ) = e ikz . (2.5)3lso, the function c ( x ) satisfies with the following conditions: c ( x ) = 1 + β ( x ) , β ( x ) ≥ , x ∈ R , and c ( x ) = 1 , x / ∈ Ω . (2.6)The assumption of (2.6) c ( x ) = 1 in R \ Ω means that we have vacuum outsideof the domain Ω . Finally, we assume that c ( x ) ∈ C ( R ). This smoothnesscondition was imposed to derive the asymptotic behavior of the solution ofthe Helmholtz equation (2.2) (see [26]). We also note that extra smoothnessconditions are usually not of a significant concern when a CIP is considered,see, e.g. Theorem 4.1 in [37]. Also, it follows from Lemma 3.3 of [28] thatthe derivative ∂ k u ( x , k ) exists for all x ∈ R , k > u ( x , k ) . Coefficient Inverse Problem (CIP).
Let Ω and Γ ⊂ ∂ Ω be as in (2.1).Let the wavenumber k ∈ [ k, k ] , where [ k, k ] ⊂ (0 , ∞ ) is an interval of wavenum-bers. Determine the function c ( x ) , x ∈ Ω , given the boundary data g ( x , k ) as u ( x , k ) = g ( x , k ) , x ∈ Γ , k ∈ [ k, k ] . (2.7)In addition to the data (2.7) we can obtain the boundary conditions for thederivative of the function u ( x , k ) in the z − direction using the data propagationprocedure (see [35]), u z ( x , k ) = g ( x , k ) , x ∈ Γ , k ∈ [ k, k ] . (2.8)Even though we use the data propagation procedure in our computations below,we do not describe it here for brevity. Instead, we refer to detailed descriptionsin [35, 36]. In fact, this procedure is widely used in Optics under the name the angular spectrum representation .In addition, we complement Dirichlet (2.7) and Neumann (2.8) boundaryconditions on Γ with the heuristic Dirichlet boundary condition at the rest ofthe boundary ∂ Ω as: u ( x , k ) = e ikz , x ∈ ∂ Ω \ Γ , k ∈ [ k, k ] . (2.9)The boundary condition (2.9) coincides with the one for the uniform mediumwith c ( x ) ≡ . To justify (2.2), we recall that, using the tail functions method, itwas demonstrated in sections 7.6 and 7.7 of [36] that (2.2) does not affect muchthe reconstruction accuracy as compared with the correct Dirichlet boundarycondition. Besides, (2.2) has always been used in works [36] with experimen-tal data, where accurate results were obtained by the tail functions globallyconvergent method.The uniqueness of the solution of this CIP is an open and long standingproblem. In fact, uniqueness of a similar coefficient inverse problem can becurrently proven only in the case if the right hand side of equation (2.2) isa function which is not vanishing in Ω . This can be done by the Bukhgeim-Klibanov method [9], also see, e.g. [7,21,22] and references cited therein for thismethod. Hence, for the computational purpose, we assume below the uniquenessof our CISP.In this last part of this section we want to briefly describe the travel time τ ( x ) which is important in our analysis. The Riemannian metric generated bythe function c ( x ) is: dτ ( x ) = (cid:112) c ( x ) | d x | , | d x | = (cid:112) ( dx ) + ( dy ) + ( dz ) . a >
0, consider the plane P a = { ( x, y, − a ) : x, y ∈ R } . We assume that Ω ⊂ { z > − a } and impose everywhere below the followingcondition on the function c ( x ): Regularity Assumption . For any point x ∈ R there exists a uniquegeodesic line Γ( x , a ) , with respect to the metric dτ , connecting x with the plane P a and perpendicular to P a near the intersection point. A sufficient condition of the regularity of geodesic lines is [38]: (cid:88) i,j =1 ∂ ln( c ( x )) ∂x i ∂x j ξ i ξ j ≥ , for all x ∈ Ω , ξ ∈ R . We introduce the travel time τ ( x ) from the plane P a to the point x as [26] τ ( x ) = (cid:90) Γ( x ,a ) (cid:112) c ( ξ ) dσ. In this section we reformulate our coefficient inverse problems as an integro-differential equation, which is one of the main ingredients in our reconstructionmethod. To this end, we first need a result on (high frequency) asymptoticbehavior of the total field u ( x , k ) in [26]. It was shown in this cited paper that u ( x , k ) = A ( x ) e ik ( τ ( x ) − a ) [1 + s ( x , k )] , x ∈ Ω , k → ∞ , (3.1)where A ( x ) > s ( x , k ) satisfies s ( x ,k ) = O (cid:18) k (cid:19) , ∂ k s ( x ,k ) = O (cid:18) k (cid:19) , x ∈ Ω , k → ∞ . (3.2)Here τ ( x ) is the length of the geodesic line generated by the function c ( x ) inthe Riemannian metric. Define w ( x , k ) = u ( x , k ) u i ( x , k ) . (3.3)From (2.5), (3.1) and (3.3), we have w ( x , k ) = A ( x ) e ik ( τ ( x ) − z − a ) [1 + s ( x , k )] , x ∈ Ω , k → ∞ . (3.4)From (3.1) and (3.4), and for x ∈ Ω, k ∈ [ k, k ], we can uniquely define thefunction log w ( x , k ) for sufficiently large values of k aslog w ( x , k ) = ln A ( x ) + ik ( τ ( x ) − z − a ) + ∞ (cid:88) n =1 ( − n − n ( s ( x , k )) n . (3.5)It is clear that, with log w ( x , k ) defined as above, exp[log w ( x , k )] equals to theright hand side of (3.4). Thus, we assume below that the number k is sufficientlylarge. 5ow we are ready to derive the integro-differential equation. For x ∈ Ω , k ∈ [ k, k ] we define the function v ( x , k ) ,v ( x , k ) = log w ( x , k ) k . (3.6)Then ∆ v + k ∇ v · ∇ v = − c ( x ) . (3.7)Setting q ( x , k ) as q ( x , k ) = ∂ k v ( x , k ) , (3.8)we obtain v ( x , k ) = − k (cid:90) k q ( x , κ ) dκ + V ( x ) . (3.9)Here we call V ( x ) the tail function, V ( x ) = v ( x , k ) . (3.10)Combining (2.2), (2.5), (2.6) and (3.3), we obtain∆ w + k βw + 2 ik ∂w∂z = 0 . (3.11)Taking into account (3.6), equation (3.11) becomes∆ v + k ∇ v · ∇ v + 2 ik ∂v∂z + β ( x ) = 0 . (3.12)To eliminate the function β ( x ) we differentiate (3.12) with respect to k, ∆ q + 2 k ∇ v · ( k ∇ q + ∇ v ) + 2 i (cid:18) k ∂q∂z + ∂v∂z (cid:19) = 0 . (3.13)Substituting (3.9) into (3.13) leads to the following integro-differential equation L ( q ) = ∆ q + 2 k (cid:32) ∇ V − (cid:90) kk ∇ q ( κ ) dκ (cid:33) · (cid:32) k ∇ ( q + V ) − (cid:90) kk ∇ q ( κ ) dκ (cid:33) +2 i (cid:32) kq z + V z − (cid:90) kk q z ( κ ) dκ (cid:33) = 0 . (3.14)This equation is complemented with the overdetermined boundary conditions: q ( x , k ) = φ ( x , k ) , q z ( x , k ) = φ ( x , k ) , x ∈ Γ , k ∈ [ k, k ] ,q ( x , k ) = 0 , x ∈ ∂ Ω \ Γ , k ∈ [ k, k ] , (3.15)where the functions φ and φ are computed from the functions g and g in(2.7), (2.8). The third boundary condition (3.15) follows from (2.5), (3.3), (3.6)and (3.8). 6ote that in (3.14) we have two unknowns q ( x , k ) and V ( x ). Hence, we willsolve the problem (3.14), (3.15) using a predictor-corrector method. Here wefind some approximation of V ( x ) first and use it as a predictor, and then solvefor q ( x , k ). One can see that if certain approximations of q ( x , k ) and V ( x ) arefound, then an approximation for the unknown coefficient c ( x ) can be foundvia (3.9) and (3.7) for a certain value of k ∈ [ k, k ]. In our computations we use k = k for that value. Therefore, we focus below on approximating functions q ( x , k ), V ( x ). In this section we present a method for finding an approximation of the tailfunction V ( x ). We note that this method is different the one studied in [27].It follows from (3.5) and (3.10) that there exists a function p ( x ) such that v ( x , k ) = p ( x ) k + O (cid:18) k (cid:19) , q ( x , k ) = − p ( x ) k + O (cid:18) k (cid:19) , k → ∞ , x ∈ Ω . (4.1)For sufficiently large k , we drop O (1 /k ) and O (1 /k ) in (4.1) and set v ( x , k ) = p ( x ) k , q ( x , k ) = − p ( x ) k , k ≥ k, x ∈ Ω . (4.2)Next, substituting (4.2) in (3.14) and setting k = k , we obtain∆ V ( x ) = 0 , x ∈ Ω . (4.3)This equation is supplemented by the following boundary conditions: V ( x ) = ψ ( x ) , V z ( x ) = ψ ( x ) , x ∈ Γ , V ( x ) = 0 , x ∈ ∂ Ω \ Γ , (4.4)where functions ψ and ψ are computed using (2.7) and (2.8). Boundary con-ditions (4.4) are over-determined ones. Due to the approximate nature of (4.2),we have observed that the obvious approach of finding the function V ( x ) bydropping the second boundary condition (4.4) and solving the resulting Dirich-let boundary value problem for Laplace equation (4.3) with the boundary data(4.4) does not provide satisfactory results. The same observation was madein [27] for the 1D case.We show in section 5 how do we approximately solve the problem (4.3)–(4.4). We now write differential operators in (3.14) and (4.3) in finite differences withrespect to x, y . Let the domain Ω ⊂ R be the orthogonal projection of thedomain Ω ⊂ R in (2.1) on the plane { z = 0 } , Ω = { ( x, y, z ) : | x | < b, | y | < b, z = 0 } . with the uniform grid step size h. Thisgrid consists of points { ( x j , y s ) } N h j,s =1 ⊂ Ω . DenoteΩ h = { ( x j , y s , z ) : ( x j , y s ) N h j,s =1 ⊂ Ω , z ∈ ( − ξ, d ) } . (5.1)For every interior point ( x j , y s , z ) ∈ Ω \ ∂ Ω four neighboring points are:( x j +1 , y s , z ) = ( x j + h, y s , z ) , ( x j − , y s , z ) = ( x j − h, y s , z ) , ( x j , y s +1 , z ) = ( x j , y s + h, z ) , ( x j , y s − , z ) = ( x j , y s − h, z ) . The corresponding Laplace operator written in partial finite differences is∆ h u = u zz + u hxx + u hyy , (5.2)where u hxx and u hyy are finite difference analogs of continuous derivatives u xx and u yy ,u hxx ( x j , y s , z ) = u ( x j − h, y s , z ) − u ( x j , y s , z ) + u ( x j + h, y s , z ) h (5.3)and similarly for u hyy . Next, ∇ h u = ( ∂ hx u, ∂ hy u, ∂ z u ) , (5.4)where ∂ hx u ( x j , y s , z ) = u ( x j + h, y s , z ) − u ( x j − h, y s , z )2 h and similarly for ∂ hy u ( x j , y s , z ) . We now rewrite problem (3.14)–(3.15) in partial finite differences. To this end,we keep in mind that only interior grid points are involved in differential oper-ators below. Using (5.1)–(5.4), we obtain for x ∈ Ω h L h ( q ) = ∆ h q + 2 k (cid:32) ∇ h V − (cid:90) kk ∇ h q ( κ ) dκ (cid:33) · (cid:32) k ∇ h ( q + V ) − (cid:90) kk ∇ h q ( κ ) dκ (cid:33) +2 i (cid:32) kq z + V z − (cid:90) kk q z ( κ ) dκ (cid:33) = 0 , (5.5) q ( x , k ) = φ ( x , k ) , q z ( x , k ) = φ ( x , k ) , x ∈ Γ , k ∈ [ k, k ] ,q ( x , k ) = 0 , x ∈ ∂ Ω \ Γ , k ∈ [ k, k ] . (5.6)Similarly, problem (4.3)–(4.4) becomes∆ h V ( x ) = 0 , x ∈ Ω . (5.7) V ( x ) = ψ ( x ) , V z ( x ) = ψ ( x ) , x ∈ Γ , V ( x ) = 0 , x ∈ ∂ Ω \ Γ . (5.8)8 emark 5.1.
1. From now on functions q ( x , k ) and V ( x ), and other functions we con-sider are semidiscrete, i.e. they are defined on Ω h . This means that, e.g. q ( x , k ) = { q ( x j , y s , z ) } N h j,s =1 , V ( x ) = { V ( x j , y s , z ) } N h j,s =1 , etc. Boundaryconditions at ∂ Ω for the functions q and V are also defined only on gridpoints which belong to the boundary ∂ Ω.2. Since the grid step size h is not changing in our arrangement, we willnot indicate below for brevity the dependence of some parameters on h ,although they do depend on h . Thus, for example below C = C ( ξ, d ) > ξ , d and h . Denote by z the complex conjugate of z ∈ C . It is convenient for us to considerany complex valued function U = Re U + i Im U = U + iU as the 2D vector func-tion U = ( U , U ) . Furthermore, each component U j of this vector function is,in turn, another vector function defined on the above grid, U j = U j ( x j , y s , z, k ) . Hence, below any Banach space of complex valued functions is actually the spaceof these real valued vector functions with the well known definitions of normsand scalar products (if in Hilbert spaces). For brevity we do not differentiatebelow between complex valued functions and corresponding vector functions.These things are always clear from the context.We introduce the Hilbert spaces H ,h (Ω h ) , L h (Ω h ) and H hn of semidiscretecomplex valued functions as H ,h (Ω h ) = { f ( x j , y s , z ) : (cid:107) f (cid:107) H n,h (Ω h ) = N h (cid:88) j,s =1 2 (cid:88) r =0 h d (cid:90) − ξ | ∂ rz f ( x j , y s , z ) | dz < ∞} ,L h (Ω h ) = { f ( x j , y s , z ) : (cid:107) f (cid:107) L h (Ω h ) = N h (cid:88) j,s =1 h d (cid:90) − ξ | f ( x j , y s , z ) | dz < ∞} ,H hn = { f ( x j , y s , z, k ) : (cid:107) f (cid:107) H hn = k (cid:90) k (cid:107) f ( x , k ) (cid:107) H n,h (Ω h ) dk < ∞} , n = 2 , . Denote [ , ] the scalar product in the space H ,h (Ω h ). We also define subspaces H ,h (Ω h ) ⊂ H ,h (Ω h ) and H h , ⊂ H h as H ,h (Ω h ) = { f ( x j , y s , z ) ∈ H ,h (Ω h ) : f ( x ) | ∂ Ω = 0 , f z ( x ) | Γ = 0 } ,H h , = { f ( x j , y s , z, k ) ∈ H h : f ( x , k ) | ∂ Ω = 0 , f z ( x , k ) | Γ = 0 , ∀ k ∈ [ k, k ] } . Note that since, for all f ∈ H h , , f ( x j , y s , z, k ) = z (cid:90) − ξ f z ( x j , y s , ρ, k ) dρ, f z ( x j , y s , z, k ) = z (cid:90) − ξ f zz ( x j , y s , ρ, k ) dρ, H h , is equivalent with (cid:107) f ( x ) (cid:107) H h , (Ω h ) = N h (cid:88) j,s =1 h d (cid:90) − ξ (cid:12)(cid:12) ∆ h f ( x j , y s , z ) (cid:12)(cid:12) dz. (5.9)In addition, for l = 0 , C l (Ω h ) = { f ( x j , y s , z ) : (cid:107) f (cid:107) C l (Ω h ) = max j,s (cid:107) f ( x j , y s , z ) (cid:107) C l [ − ξ,d ] < ∞} ,C hl = { f ( x j , y s , z, k ) : (cid:107) f (cid:107) C hl = max k ∈ [ k,k ] (cid:107) f ( x j , y s , z, k ) (cid:107) C l (Ω h ) < ∞} . By embedding theorem H ,h (Ω h ) ⊂ C (Ω h ) , H hn ⊂ C hn − and (cid:107) f (cid:107) C (Ω h ) ≤ C (cid:107) f (cid:107) H ,h (Ω h ) , for all f ∈ H ,h (Ω h ) , (5.10) (cid:107) f (cid:107) C hn − ≤ C (cid:107) f (cid:107) H hn , for all f ∈ H hn . (5.11) It is our computational experience for the 1D case [3, 25, 27] that one shoulduse for computations such a CWF which would be a simple one. A similarconclusion can be found on page 1581 of [5]. Thus, the CWF we use in thispaper is: ϕ λ ( z ) = e − λz . (6.1) First, we present the cost functional for the solution of problem (5.7)–(5.8)which is about the tail function. Non-zero boundary conditions in (5.8) areinconvenient for us. Hence, we assume that there exists a function Q ( x ) ∈ H ,h (Ω h ) such that Q ( x ) = ψ ( x ) , ∂ z Q ( x ) = ψ ( x ) , x ∈ Γ; Q ( x ) = 0 , x ∈ ∂ Ω \ Γ , (6.2)Define W ( x ) = V ( x ) − Q ( x ) ∈ H ,h (Ω h ) . (6.3)Hence, we consider the following minimization problem: Minimization Problem 1 . For W ∈ H ,h (Ω h ) , minimize the functional I µ ( W ) ,I µ ( W ) = e µd N h (cid:88) j,s =1 h d (cid:90) − ξ (cid:12)(cid:12) (∆ h W + ∆ h Q )( x j , y s , z ) (cid:12)(cid:12) ϕ µ ( z ) dz. (6.4)The multiplier e µd is introduced here to ensure that e µd min [ − ξ,d ] ϕ µ ( z ) =1 . Remark 6.1.
Since the operator ∆ h is linear, then, in principle at least, onecan apply straightforwardly the quasi-reversibility method to find an approximatesolution of the problem ∆ W + ∆ Q = 0 for W ∈ H ,h (Ω h ) [23]. This means thatone can use λ = 0 in (6.4). However, it was observed in [3] that the involvementof the CWF like in (6.4) leads to a better solution accuracy.
10e now follow the classical Tikhonov regularization concept [4, 40]. By thisconcept, we should assume that there exists an exact solution V ∗ ( x ) of theproblem (5.7)–(5.8) with the noiseless data ψ ∗ ( x ) , ψ ∗ ( x ) . Below the subscript“ ∗ ” is related only to the exact solution. In fact, however, the data ψ ( x ) and ψ ( x ) contain noise. Let δ ∈ (0 ,
1) be the level of noise in the data ψ ( x )and ψ ( x ). Again, following the same concept, we should assume that thenumber δ ∈ (0 ,
1) is sufficiently small. Assume that there exists the function Q ∗ ( x ) ∈ H ,h (Ω h ) such that Q ∗ ( x ) = ψ ∗ ( x ) , ∂ z Q ∗ ( x ) = ψ ∗ ( x ) , x ∈ Γ; Q ∗ ( x ) = 0 , x ∈ ∂ Ω \ Γ , (6.5) (cid:107) Q − Q ∗ (cid:107) H ,h (Ω h ) < δ, (6.6)where Q is defined in (6.2). We will choose in Theorem 7.2 of section 6 acertain dependence µ = µ ( δ ) of the parameters µ on the noise level δ . Denote W µ ( δ ) ( x ) = W min ( x ) the unique minimizer of the functional I µ ( δ ) ( W ) (Theorem7.2) and by (6.3) let V µ ( δ ) ( x ) = W µ ( δ ) ( x ) + Q ( x ) = W min ( x ) + Q ( x ) . (6.7) Suppose that there exists a function F ( x , k ) ∈ H h such that (see (3.15)): F ( x , k ) = φ ( x , k ) , F z ( x , k ) = φ ( x , k ) , x ∈ Γ , F ( x , k ) = 0 , x ∈ ∂ Ω \ Γ . (6.8)Also, assume that there exists an exact solution c ∗ ( x ) of our CIP satisfyingthe above conditions imposed on the coefficient c ( x ) and generating the noise-less boundary data φ , ∗ and φ , ∗ in (3.15). Also, assume that there exists thefunction F ∗ ( x , k ) ∈ H h satisfying the following analog of boundary conditions(6.8): F ∗ ( x , k ) = φ , ∗ ( x , k ) , ∂ z F ∗ ( x , k ) = φ , ∗ ( x , k ) , x ∈ Γ , F ∗ ( x , k ) = 0 , x ∈ ∂ Ω \ Γ . (6.9)We assume that (cid:107) F − F ∗ (cid:107) H h < δ. (6.10)Let q ∗ ∈ H h be the function q generated by the exact coefficient c ∗ ( x ). Wedefine functions p and p ∗ as p ( x , k ) = q ( x , k ) − F ( x , k ) , p ∗ ( x , k ) = q ∗ ( x , k ) − F ∗ ( x , k ) . (6.11)Hence, the functions p , p ∗ ∈ H h , . Let
R > B ( R ) ⊂ H h , of the radius R , B ( R ) = { r ∈ H h , : (cid:107) r (cid:107) H h < R } . (6.12)Using the integro-differential equation (3.14), boundary conditions (3.15) forit, (6.8), (6.9) and (6.11), we construct our cost functional J λ ( p ) with the CWF116.1) in it as: J λ ( p ) = e λd N h (cid:88) j,s =1 h k (cid:90) k d (cid:90) − ξ | L h ( p + F )( x j , y s , z, κ ) | ϕ λ ( z ) dzdκ, p ∈ B ( R ) , (6.13)where the tail function in L h is defined in (6.7). Similarly with (6.4), themultiplier e λd is introduced to balance two terms in the right hand side of(6.13). We consider the following minimization problem: Minimization Problem 2 . Minimize the functional J λ ( p ) on the set p ∈ B ( R ) . In this section we formulate theorems about the minimization problems 1 and2 of section 6. First, we are concerned with the Carleman estimate with theCWF (6.1).
Theorem 7.1 (Carleman estimate).
For λ > let B h ( u, λ ) = M h (cid:88) j,s =1 h d (cid:90) − ξ | ∆ h u ( x j , y s , z ) | ϕ λ ( z ) dz. Then there exists a sufficiently large number λ = λ ( ξ, d ) > such that for all λ ≥ λ the following estimate is valid for all functions u ∈ H ,h (Ω h ) B h ( u, λ ) ≥ C M h (cid:88) j,s =1 h d (cid:90) − ξ | u zz ( x j , y s , z ) | ϕ λ ( z ) dz + Cλ M h (cid:88) j,s =1 h d (cid:90) − ξ | u z ( x j , y s , z ) | ϕ λ ( z ) dz + Cλ M h (cid:88) j,s =1 h d (cid:90) − ξ | u ( x j , y s , z ) | ϕ λ ( z ) dz. (7.1) Proof.
Recall that we do not indicate the dependence of neither constants C nor other constants on h (second item in Remarks 5.1). Since | v | = (Re v ) +(Im v ) for v ∈ C , then it is sufficient to prove estimate (7.1) for real valuedfunctions u ∈ H ,h (Ω h ) . The following Carleman estimate was proven in lemma3.1 of [27] for all real valued functions w ( z ) ∈ H ( − ξ, d ) such that w ( − ξ ) = w (cid:48) ( − ξ ) = 0: d (cid:90) − ξ ( w (cid:48)(cid:48) ) ϕ λ ( z ) dz ≥ (cid:101) C d (cid:90) − ξ ( w (cid:48)(cid:48) ) ϕ λ ( z ) dz + (cid:101) Cλ d (cid:90) − ξ ( w (cid:48) ) ϕ λ ( z ) dz + (cid:101) Cλ d (cid:90) − ξ w ϕ λ ( z ) dz, (7.2)for all λ ≥ λ ( ξ, d ) , where the number (cid:101) C = (cid:101) C ( ξ, d ) > ξ and12 . Next, it follows from (5.2) and (5.3) that B h ( u, λ ) = M h (cid:88) j,s =1 h d (cid:90) − ξ [( u zz + u hxx + u hyy )( x j , y s , z )] ϕ λ ( z ) dz ≥ M h (cid:88) j,s =1 h d (cid:90) − ξ [ u zz ( x j , y s , z )] ϕ λ ( z ) dz − C M h (cid:88) j,s =1 d (cid:90) − ξ [ u ( x j , y s , z )] ϕ λ ( z ) dz. Hence, using (7.2), we obtain B h ( u, λ ) ≥ C M h (cid:88) j,s =1 h d (cid:90) − ξ [ u zz ( x j , y s , z )] ϕ λ ( z ) dz + Cλ M h (cid:88) j,s =1 h d (cid:90) − ξ [ u z ( x j , y s , z )] ϕ λ ( z ) dz + Cλ M h (cid:88) j,s =1 h d (cid:90) − ξ [ u ( x j , y s , z )] ϕ λ ( z ) dz − C M h (cid:88) j,s =1 d (cid:90) − ξ [ u ( x j , y s , z )] ϕ λ ( z ) dz. (7.3)Now choosing λ so large that λ h > C/
2, we obtain from (7.3) the targetestimate (7.1) for all λ > λ . (cid:3) The next theorem is about the functional I µ ( W ) in (6.4). Theorem 7.2.
Assume that there exists a function Q ∈ H ,h (Ω h ) satisfyingconditions (6.5). Introduce the function W ∈ H ,h (Ω h ) via (6.3). Then for each µ > there exists unique minimizer W µ ∈ H ,h (Ω h ) of the functional (6.4).Suppose now that there exists an exact solution V ∗ ∈ H ,h (Ω h ) of equation(4.3) with the boundary data ψ ∗ ( x ) and ψ ∗ ( x ) in (4.4). Also, assume thatthere exists a function Q ∗ ∈ H ,h (Ω h ) satisfying conditions (6.5) and such thatinequality (6.6) holds, where δ ∈ (0 , is the noise level in the data. Let λ > be the number of Theorem 6.1. Choose a number δ ∈ (0 , e − d + ξ ) λ ). For any δ ∈ (0 , δ ) let µ = µ ( δ ) = ln( δ − / (2( d + ξ )) ) . (7.4) Let the function V µ ( δ ) ( x ) be defined via (6.7) . Then the following convergenceestimate of V µ ( δ ) ( x ) to the exact solution V ∗ ( x ) holds as δ → (cid:13)(cid:13) V µ ( δ ) − V ∗ (cid:13)(cid:13) H ,h (Ω h ) ≤ C √ δ. (7.5) In addition, V µ ( δ ) ∈ C (Ω h ) and C (cid:13)(cid:13) ∇ V µ ( δ ) (cid:13)(cid:13) C (Ω h ) ≤ (cid:13)(cid:13) V µ ( δ ) (cid:13)(cid:13) H ,h (Ω h ) ≤ C (cid:104) (cid:107) V ∗ (cid:107) H ,h (Ω h ) (cid:105) . (7.6) Proof.
It follows from (6.4) and the variational principle that the vectorfunction W min = ( W , min , W , min ) ∈ H ,h (Ω h ) is a minimizer of the functional13 µ,α ( W ) if and only if e µd M h (cid:88) j,s =1 h d (cid:90) − ξ (∆ h W , min ∆ h r + ∆ h W , min ∆ h r )( x j , y s , z ) ϕ µ ( z ) dz = − e µd M h (cid:88) j,s =1 h d (cid:90) − ξ (∆ h Q ∆ h r + ∆ h Q ∆ h r )( x j , y s , z ) ϕ µ ( z ) dz, (7.7)for all r = ( r , r ) ∈ H ,h (Ω h ). For any vector function P = ( P , P ) ∈ H ,h (Ω h ) consider the expression in the left hand side of (7.7) in which thevector function ( W , min , W , min ) is replaced with ( P , P ). Then (5.9) impliesthat this expression defines a new scalar product { P, r } in the space H ,h (Ω h ) , and the corresponding norm { P, P } / is equivalent to the norm in the space H ,h (Ω h ). Next, for all r = ( r , r ) ∈ H ,h (Ω h ), we have (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − e µd M h (cid:88) j,s =1 h d (cid:90) − ξ (∆ h Q ∆ h r + ∆ h Q ∆ h r )( x j , y s , z ) ϕ µ ( z ) dz (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ≤ D (cid:107) Q (cid:107) H ,h (Ω h ) (cid:107) r (cid:107) H ,h (Ω h ) ≤ D (cid:112) { Q, Q } (cid:112) { r, r } , where the constants D, D do not depend on Q and r . Hence, by Riesz theoremthere exists unique vector function (cid:98) Q = (cid:16) (cid:98) Q , (cid:98) Q (cid:17) = (cid:98) Q ( Q ) ∈ H ,h (Ω h ) suchthat { (cid:98) Q, r } = − e µd M h (cid:88) j,s =1 h d (cid:90) − ξ (∆ h Q ∆ h r + ∆ h Q ∆ h r )( x j , y s , z ) ϕ µ ( z ) dz, for all r = ( r , r ) ∈ H ,h (Ω h ). Hence, by (7.7) { W min , r } = { (cid:98) Q, r } , ∀ r ∈ H ,h (Ω h ). This implies that W min = (cid:98) Q . Thus, both existence and uniquenessof the minimizer of the functional I µ ( W ) are established.We now prove convergence estimate (7.5). Let W ∗ = ( W ∗ , , W ∗ , ) = V ∗ − Q ∗ . Then W ∗ ∈ H ,h (Ω h ) . Denote (cid:102) W = W min − W ∗ and (cid:101) Q = Q − Q ∗ . Since e µd M h (cid:88) j,s =1 h d (cid:90) − ξ (∆ h W ∗ , ∆ h r + ∆ h W ∗ , ∆ h r )( x j , y s , z ) ϕ µ ( z ) dz = − e µd M h (cid:88) j,s =1 h d (cid:90) − ξ (∆ h Q ∗ ∆ h r + ∆ h Q ∗ ∆ h r )( x j , y s , z ) ϕ µ ( z ) dz, (7.8)14hen subtracting (7.8) from (7.7) and setting r = (cid:102) W , we obtain e µd M h (cid:88) j,s =1 h d (cid:90) − ξ (∆ h (cid:102) W ) ϕ µ ( z ) dz = − e µd M h (cid:88) j,s =1 h d (cid:90) − ξ (∆ h (cid:101) Q ∆ h (cid:102) W + ∆ h (cid:101) Q ∆ h (cid:102) W ) ϕ µ ( z ) dz. Using the Cauchy-Schwarz inequality, taking into account (6.6) and (7.4), weobtain e µd M h (cid:88) j,s =1 h d (cid:90) − ξ (∆ h (cid:102) W ( x j , y s , z )) ϕ µ ( z ) dz ≤ Ce µ ( d + ξ ) δ . (7.9)By (7.4) µ = ln( δ − / (2( d + ξ )) ), which implies e µ ( d + ξ ) δ = δ. Hence, (7.9) impliesthat e µd M h (cid:88) j,s =1 h d (cid:90) − ξ (∆ h (cid:102) W ( x j , y s , z )) ϕ µ ( z ) dz ≤ Cδ. (7.10)We now apply Theorem 7.1 to the left hand side of (7.10). We obtain for all µ ≥ µ δ ≥ Ce µd M h (cid:88) j,s =1 h d (cid:90) − ξ ( (cid:102) W zz ( x j , y s , z )) ϕ µ ( z ) dz + Ce µd µ M h (cid:88) j,s =1 h d (cid:90) − ξ ( (cid:102) W z ( x j , y s , z )) ϕ µ ( z ) dz + µ M h (cid:88) j,s =1 h d (cid:90) − ξ ( (cid:102) W ( x j , y s , z )) ϕ µ ( z ) dz . (7.11)Since e µd ϕ µ ( z ) ≥ e µd ϕ µ ( d ) = 1 , then (7.11) implies that (cid:107) (cid:102) W (cid:107) H ,h (Ω h ) ≤ C √ δ. Hence, (6.6), (6.7) and triangle inequality imply that (cid:13)(cid:13) V µ ( δ ) − V ∗ (cid:13)(cid:13) H ,h (Ω h ) = (cid:13)(cid:13)(cid:13)(cid:102) W + ( Q − Q ∗ ) (cid:13)(cid:13)(cid:13) H ,h (Ω h ) ≤ C √ δ, (7.12)which proves (7.5). Next, by (7.5) and triangle inequality imply the right esti-mate (7.6). The left estimate (7.6) follows from (5.10). (cid:3) The main analytical result of this paper is Theorem 7.3.
Theorem 7.3 (globally strict convexity).
Assume that conditions of The-orem 6.2 hold. Let λ ≥ λ be the number defined below in the formulation ofthis theorem. Assume that there exist functions F ( x , k ), F ∗ ( x , k ) ∈ H h satisfy-ing conditions (6.8)–(6.10), where δ ∈ (0 , δ ) and δ ∈ (0 , e − d + ξ ) λ ). Set in(6.13) V = V µ ( δ ) , where the function V µ ( δ ) is defined in Theorem 7.2. First, the unctional J λ ( p ) has the Frech´et derivative J (cid:48) λ ( p ) ∈ H h , at any point p ∈ H h , . Second, there exist numbers λ = λ (Ω h , R, (cid:107) F ∗ (cid:107) H h , (cid:107) V ∗ (cid:107) H ,h (Ω h ) , k, k ) ≥ λ ,C = C (Ω h , R, (cid:107) F ∗ (cid:107) H h , (cid:107) V ∗ (cid:107) H ,h (Ω h ) , k, k ) > , depending only on listed parameters, such that for any λ ≥ λ the functional J λ ( p ) is strictly convex on B ( R ) . In other words, the following estimate holds: J λ ( p + r ) − J λ ( p ) − J (cid:48) λ ( p )( r ) ≥ C (cid:107) r (cid:107) H h , for all p, p + r ∈ B ( R ) . (7.13) Proof.
In this proof C = C (Ω h , R, (cid:107) F ∗ (cid:107) H h , (cid:107) V ∗ (cid:107) H ,h (Ω h ) , k, k ) > V ( x ) = V µ ( δ ) ( x ) and also sometimes we do not indicate the dependencies on( x j , y s , z, k ). Note that (5.10), (5.11), (6.10)-(6.12) and (7.6) imply that (cid:107)∇ V (cid:107) C (Ω h ) ≤ C , (cid:107) F (cid:107) C h ≤ C , (7.14) (cid:107)∇ h r (cid:107) C h ≤ C . (7.15)Consider an arbitrary vector function p = ( p , p ) ∈ B ( R ) and an arbitraryfunction r = ( r , r ) ∈ H h , such that p + r ∈ B ( R ). By (6.13) we need toconsider A , where A = | L h ( p + r + F ) | − | L h ( p + F ) | . (7.16)First, we will single out such a part of A , which is linear with respect to r . Thiswill lead us to the Frech´et derivative J (cid:48) λ . Next, we will single out | ∆ h r | . Basedon this, we will apply the Carleman estimate of Theorem ?? . For all z , z ∈ C ,we have | z | − | z | = ( z − z ) z + ( z − z ) z . (7.17)Denote z = L h ( p + r + F ) , z = L h ( p + F ) . (7.18)Then A = ( z − z ) z , A = ( z − z ) z , A = A + A . (7.19)Using (5.5), (6.11) and (7.18), we obtain z − z = ∆ h r − k ∇ h r · (cid:34) ∇ h V − k (cid:90) k ∇ h ( p + F ) dκ (cid:35) + 2 i (cid:34) r z − k (cid:90) k r z dκ (cid:35) . +2 k k (cid:90) k ∇ h r dκ · (cid:34) ∇ h V − k (cid:90) k ∇ h ( p + F ) dκ + k ∇ h ( p + F ) (cid:35) (7.20)16ext, z = ∆ h ( r + p + F ) − k (cid:34) ∇ h V − k (cid:90) k ∇ h ( r + p + F ) dκ (cid:35) · (cid:34) k ∇ h ( r + p + F )+ ∇ h V − k (cid:90) k ∇ h ( r + p + F ) dκ (cid:35) − i (cid:34) k (cid:0) r z + p z + F z (cid:1) + V z − k (cid:90) k (cid:0) r z + p z + F z (cid:1) dκ (cid:35) . Hence, A = ( z − z ) z = | ∆ h r | + Z l, ( r, k ) + Z ( r, k ) , (7.21)where Z l, ( h, k ) is linear with respect to the vector function r = ( r , r ) ,Z l, ( r, k ) = ∆ h r · Y + ∇ h r ∇ h Y · Y + ∇ h r ∇ Y · Y + k (cid:90) k ∇ h rdκ ∇ Y · Y + k (cid:90) k ∇ h rdκ ∇ h Y · Y + r z − k (cid:90) k r z dκ Y + r z − k (cid:90) k r z dκ Y , (7.22)where explicit expressions for functions Y j ( x , k ) , j = 1 , ...,
11 can be written inan obvious way. Also, it follows from those formulae as well as from (7.14) that Y , Y , Y , Y ∈ C h and Y , Y , Y , Y , Y , Y ∈ C h . In addition, (cid:40) (cid:107) Y (cid:107) C h , (cid:107) Y (cid:107) C h , (cid:107) Y (cid:107) C h , (cid:107) Y (cid:107) C h ≤ C , (cid:107) Y (cid:107) C h , (cid:107) Y (cid:107) C h , (cid:107) Y (cid:107) C h , (cid:107) Y (cid:107) C h , (cid:107) Y (cid:107) C h , (cid:107) Y (cid:107) C h ≤ C . (7.23)The term Z ( r, k ) in (7.21) is nonlinear with respect to r . Applying the Cauchy-Schwarz inequality and (7.15), we obtain Z ( r, k ) ≥ | ∆ h r | − C |∇ h r | − C k (cid:90) k |∇ h r | dκ. (7.24)Similarly with (7.21)–(7.24) we obtain A = ( z − z ) z = Z l, ( r, k ) + Z ( r, k ) , (7.25)where the term Z l, ( r, k ) is linear with respect to r and has the form similarwith the one in (7.22), although with different functions Y j , which still satisfydirect analogs of estimates (7.23). As to the term Z ( r, k ) , it is nonlinear withrespect to r and, as in (7.24), Z ( r, k ) ≥ − C |∇ h r | − C k (cid:90) k |∇ h r | dκ. (7.26)17n addition, the following upper estimate is valid | Z ( r, k ) | + | Z ( r, k ) | ≤ C | ∆ h r | + |∇ h r | + k (cid:90) k |∇ h r | dκ . (7.27)Thus, it follows from (6.13) and (7.18)-(7.26) that J λ ( p + r ) − J λ ( p ) = e λd N h (cid:88) j,s =1 h k (cid:90) k d (cid:90) − ξ (cid:2) ( S ∆ h r + S · ∇ h r )( x j , y s , z ) (cid:3) ϕ λ ( z ) dzdκ + e λd N h (cid:88) j,s =1 h k (cid:90) k d (cid:90) − ξ Z ( r, k )( x j , y s , z ) ϕ λ ( z ) dzdκ, (7.28)where Z ( r, k )( x j , y s , z ) = Z ( r, k )( x j , y s , z ) + Z ( r, k )( x j , y s , z ) . (7.29)The second line of (7.28),Lin( r ) = e λd N h (cid:88) j,s =1 h k (cid:90) k d (cid:90) − ξ (cid:2) ( S ∆ h r + S · ∇ h r )( x j , y s , z ) (cid:3) ϕ λ ( z ) dzdκ (7.30)is linear with respect to r . Also, the vector functions S ( x , k ) and S ( x , k ) aresuch that | S ( x , k ) | , | S ( x , k ) | ≤ C in Ω h × [ k, k ] . (7.31)As to the third line of (7.28), it can be estimated from the below as e λd M h (cid:88) j,s =1 h k (cid:90) k d (cid:90) − ξ Z ( r, k )( x j , y s , z ) ϕ λ ( z ) dzdκ ≥ e λd M h (cid:88) j,s =1 h k (cid:90) k d (cid:90) − ξ | ∆ h r ( x j , y s , z ) | ϕ λ ( z ) dzdκ (7.32) − C e λd k (cid:90) k d (cid:90) − ξ |∇ h r ( x j , y s , z ) | ϕ λ ( z ) dzdκ .
18n addition, using (7.27) and (7.29), we obtain e λd N h (cid:88) j,s =1 h k (cid:90) k d (cid:90) − ξ | Z ( r, k )( x j , y s , z ) | ϕ λ ( z ) dzdκ ≤ C e λd N h (cid:88) j,s =1 h k (cid:90) k d (cid:90) − ξ (cid:32) | ∆ h r | + |∇ h r | + k (cid:90) k |∇ h r | dκ (cid:33) ( x j , y s , z ) ϕ λ ( z ) dzdκ. (7.33)The functional Lin( h ) in (7.30) is linear with respect to r . Also, by (7.30) and(7.31) | Lin( r ) | ≤ C e λ ( d + ξ ) (cid:107) r (cid:107) H h , for all r ∈ H h , . Hence, Lin( r ) : H h , → R is a bounded linear functional. Hence, by Riesztheorem for each pair λ > X λ ∈ H h , independenton r such that Lin( r ) = [ X λ , r ] , for all r ∈ H h , . (7.34)In addition, (7.28), (7.33) and (7.34) imply that | J λ ( p + r ) − J λ ( p ) − [ X λ , r ] | ≤ C e λ ( d + ξ ) (cid:107) r (cid:107) H h , for all r ∈ H h , . (7.35)Thus, (7.28)–(7.35) imply that X λ ∈ H h , is the Frech´et derivative of the func-tional J λ ( p ) at the point p, i.e. X λ = J (cid:48) λ ( p ).Next, using (7.28) and (7.32), we obtain J λ ( p + r ) − J λ ( p ) − J (cid:48) λ ( p )( r ) ≥ e λd N h (cid:88) j,s =1 h k (cid:90) k d (cid:90) − ξ | ∆ h r ( x j , y s , z ) | ϕ λ ( z ) dzdκ − C k (cid:90) k d (cid:90) − ξ |∇ h r ( x j , y s , z ) | ϕ λ ( z ) dzdκ . We now apply Carleman estimate of Theorem 7.1 for λ ≥ λ ,e λd N h (cid:88) j,s =1 h k (cid:90) k d (cid:90) − ξ (cid:12)(cid:12) ∆ h r ( x j , y s , z ) (cid:12)(cid:12) ϕ λ ( z ) dzdκ − C k (cid:90) k d (cid:90) − ξ (cid:12)(cid:12) ∇ h r ( x j , y s , z ) (cid:12)(cid:12) ϕ λ ( z ) dzdκ ≥ e λd N h (cid:88) j,s =1 h k (cid:90) k d (cid:90) − ξ | r zz ( x j , y s , z ) | ϕ λ ( z ) dzdκ + Cλ d (cid:90) − ξ [ r z ( x j , y s , z )] ϕ λ ( z ) dz + λ d (cid:90) − ξ [ r ( x j , y s , z )] ϕ λ ( z ) dz − C e λd k (cid:90) k d (cid:90) − ξ |∇ h r ( x j , y s , z ) | ϕ λ ( z ) dzdκ. Hence, from these two equations it follows that for sufficiently large λ λ = λ (Ω h , R, (cid:107) F ∗ (cid:107) H h , (cid:107) V ∗ (cid:107) H ,h (Ω h ) , k, k ) ≥ λ λ ≥ λ J λ ( p + r ) − J λ ( p ) − J (cid:48) λ ( p )( r ) ≥ e λd N h (cid:88) j,s =1 h k (cid:90) k d (cid:90) − ξ | r zz ( x j , y s , z ) | ϕ λ ( z ) dzdκ + C λ d (cid:90) − ξ [ r z ( x j , y s , z )] ϕ λ ( z ) dz + C λ d (cid:90) − ξ [ r ( x j , y s , z )] ϕ λ ( z ) dz ≥ C (cid:107) r (cid:107) H h , which establishes (7.13). (cid:3) Theorem 7.4 . Assume that the conditions of Theorems 7.2 and 7.3 re-garding the tail function V = V µ ( δ ) and the functions F and F ∗ are satis-fied. Then the Frech´et derivative J (cid:48) λ of the functional J λ satisfies the Lip-schitz continuity condition in any ball B ( R (cid:48) ) as in (6.12) with an arbitrary R (cid:48) > . More precisely, the following inequality holds with the constant M = M (Ω h , R (cid:48) , (cid:107) F ∗ (cid:107) H h , (cid:107) V ∗ (cid:107) H ,h (Ω h ) , λ, k, k ) > depending only on listed parame-ters: (cid:107) J (cid:48) λ ( p ) − J (cid:48) λ ( p ) (cid:107) H h ≤ M (cid:107) p − p (cid:107) H h , for all p , p ∈ B ( R (cid:48) ) . The proof of this theorem is completely similar with that of theorem 3.1of [3] and is, therefore, omitted.Denote P B : H h , → B ( R ) the projection operator of the Hilbert space H h , on B ( R ) ⊂ H h , . Let p ∈ B ( R ) be an arbitrary point of the ball B ( R ). Let thenumber γ ∈ (0 , p n = P B ( p n − − γJ (cid:48) λ ( p n − )) , n = 1 , , . . . (7.36)The following theorem follows immediately from the combination of Theorems7.3 and 7.4 with lemma 2.1 and Theorem 2.1 of [3]. Theorem 7.5.
Assume that conditions of Theorems 7.2 and 7.3 are satis-fied. Let λ ≥ λ , where λ is defined in Theorem 7.3. Then there exists uniqueminimizer p min ,λ ∈ B ( R ) of the functional J λ ( p ) on the set B ( R ) and J (cid:48) λ ( p min ,λ )( y − p min ,λ ) ≥ , for all y ∈ H h , . (7.37) Also, there exists a sufficiently small number γ = γ (Ω h , R, (cid:107) F ∗ (cid:107) H h , (cid:107) V ∗ (cid:107) H ,h (Ω h ) , k, k, λ ) ∈ (0 , depending only on listed parameters such that for any γ ∈ (0 , γ ) the se-quence (7.36) converges p min ,λ , (cid:107) p min ,λ − p n (cid:107) H h ≤ θ n (cid:107) p min ,λ − p (cid:107) H h , n = 1 , , . . . (7.38) where the number θ = θ (Ω h , R, (cid:107) F ∗ (cid:107) H h , (cid:107) V ∗ (cid:107) H ,h (Ω h ) , k, k, λ, γ ) ∈ (0 , de-pends only on listed parameters. Thus, (7.38) estimates the convergence rate of the sequence (7.36) to theminimizer p min ,λ . We now need to estimate the convergence rate of this sequenceto the exact solution. To do this, we follow the Tikhonov regularization concept[4, 40] in Theorem 7.6 via assuming that the exact solution p ∗ ∈ B ( R ) . Theorem 7.6.
Assume that conditions of Theorems 7.2 and 7.3 are satis-fied. Let λ be the number of Theorem 7.3, δ ∈ (0 , e − d + ξ ) λ ) and δ ∈ (0 , δ ) . et λ = λ ( δ ) = ln( δ − / (4( d + ξ )) ) > λ . Furthermore, assume that the function p ∗ ∈ B ( R ) . Then there exists a number C = C (Ω h , R, (cid:107) F ∗ (cid:107) H h , (cid:107) V ∗ (cid:107) H ,h (Ω h ) , k, k ) > depending only on listed parameters such that (cid:13)(cid:13) p ∗ − p min ,λ ( δ ) (cid:13)(cid:13) H h ≤ C δ / , (7.39) (cid:13)(cid:13) c ∗ − c min ,λ ( δ ) (cid:13)(cid:13) L h (Ω h ) ≤ C δ / , (7.40) In addition, the following convergence estimates hold (cid:107) p ∗ − p n (cid:107) H h ≤ C δ / + θ n (cid:13)(cid:13) p min ,λ ( δ ) − p (cid:13)(cid:13) H h , n = 1 , , . . . (7.41) (cid:107) c ∗ − c n (cid:107) L h (Ω h ) ≤ C δ / + θ n (cid:13)(cid:13) p min ,λ ( δ ) − p (cid:13)(cid:13) H h , n = 1 , , . . . (7.42) where θ ∈ (0 , is the number of Theorem 7.5 and functions c min ,λ ( δ ) ( x ) and c n ( x ) is reconstructed from functions p min ,λ ( δ ) and p n ( x , k ) respectively using(3.7)–(3.10) and (6.11). Remark 7.1.
Since
R > is an arbitrary number and p is an arbitrary pointof the ball B ( R ) , then Theorems 7.5 and 7.6 ensure the global convergence ofthe gradient projection method for our case, see section 1. We note that if afunctional is non convex, then the convergence of a gradient-like method of itsminimization can be guaranteed only if the starting point of iterations is locatedin a sufficiently small neighborhood of its minimizer.Proof. We temporarily denote J λ ( p ) := J λ ( p + F ), see (6.13). We have J λ ( p ∗ + F ∗ ) = e λd N h (cid:88) j,s =1 h k (cid:90) k d (cid:90) − ξ | L h ( p ∗ + F ∗ ) ( x j , y s , z, κ ) | ϕ λ ( z ) dzdκ = 0 . (7.43)It follows from (5.5), (5.6), (6.8)-(6.13), (7.5) (7.43) that J λ ( p ∗ + F ) ≤ C δe λ ( d + ξ ) . (7.44)Next, using (7.13) and , we obtain J λ ( p ∗ + F ) − J λ ( p min ,λ ( δ ) + F ) − J (cid:48) λ ( p min ,λ ( δ ) + F )( p ∗ − p min ,λ ) ≥ C (cid:107) p ∗ − p min ,λ (cid:107) H h . Hence, since − J λ ( p min ,λ ( δ ) + F ) ≤ − J (cid:48) λ ( p min ,λ + F )( p ∗ − p min ,λ ( δ ) ) ≤ , we obtain, using (7.44) and recalling that λ = ln( δ − / (4( d + ξ )) ): (cid:13)(cid:13) p ∗ − p min ,λ ( δ ) (cid:13)(cid:13) H h ≤ C √ δ, which implies (7.39). Estimate (7.40) follows immediately from (3.7)-(3.10),(6.11), (7.5) and (7.39).We now prove (7.41) and (7.42). Using (7.38), (7.39) and the triangle in-equality, we obtain for n = 1 , , . . . (cid:107) p ∗ − p n (cid:107) H h ≤ (cid:13)(cid:13) p ∗ − p min ,λ ( δ ) (cid:13)(cid:13) H h + (cid:13)(cid:13) p min ,λ ( δ ) − p n (cid:13)(cid:13) H h ≤ C δ / + (cid:13)(cid:13) p min ,λ ( δ ) − p n (cid:13)(cid:13) H h ≤ C δ / + θ n (cid:107) p min ,λ − p (cid:107) H h , (cid:107) c ∗ − c n (cid:107) L h (Ω h ) ≤ (cid:13)(cid:13) c ∗ − c min ,λ ( δ ) (cid:13)(cid:13) L h (Ω h ) + (cid:13)(cid:13) c min ,λ ( δ ) − c n (cid:13)(cid:13) L h (Ω h ) ≤ C δ / + C (cid:13)(cid:13) p min ,λ ( δ ) − p n (cid:13)(cid:13) H h ≤ C δ / + C (cid:13)(cid:13) p min ,λ ( δ ) − p n (cid:13)(cid:13) H h ≤ C δ / + C θ n (cid:13)(cid:13) p min ,λ ( δ ) − p (cid:13)(cid:13) H h . The latter establishes (7.42). (cid:3)
We present in this section a numerical study of the application of our convexifi-cation method to microwave experimental backscatter data for buried objects.One of possible applications is in the standoff detection of explosives. We notethat these data were treated in [35] by a different globally convergent method.We first describe very briefly the measured data and its preprocessing which isimportant for the application of our convexification method. We refer to [35]for all the details of data collection and preprocessing.
The experimental data were measured by a scattering facility at the Universityof North Carolina at Charlotte. We have measured the backscatter data forobjects buried in a sandbox. This sandbox was filled with dry sand and containsno moisture, see Figure 1. The data were measured on a rectangular surface ofdimensions 1 m × x − axis and the y − axis are respectively the horizontal and the vertical axis,while the z − axis is orthogonal to the measurement surface. The direction fromthe measurement surface to the target is the positive direction of the z − axis.Figure 1: A schematic diagram of the collection of our experimental data.The measurements consist of multi-frequency backscatter data associatedwith 300 frequency points uniformly distributed over the range from 1 GHzto 10 GHz. However, we work with the preprocessed data which are stable onnarrow intervals of frequencies centered at 2.6 GHz, 3.01 GHz or 3.1 GHz. Sincethe corresponding wavelength for 2.6 GHz is 11.5 cm, the distance between the22ource and the buried targets was about at least 6.17 wavelengths. This distanceis sufficiently large in terms of wavelengths, and therefore justifies our modelingof the source as a plane wave. The backscatter data were generated by a singledirection of the incident plane wave.Recall that these experimental data were preprocessed in [35] and we willstudy the performance of our inversion method on that preprocessed data in-stead of the raw ones. The preprocessing developed in the cited paper comprisestwo main goals: distill the signals reflected by our buried targets from signalsreflected by the sandbox and other unwanted objects, and reduce the noise inthe data as well as the computational domain.For the convenience of the readers we briefly summarize the main steps ofthe data preprocessing developed in [35].Step 1. Subtract the reference data from the measured data for buried objects.The reference data are the ones measured in the case when the sandboxcontains no buried objects. This subtraction helps us to sort of extractthe signals of the buried targets from the total signal and also to reducethe noise.Step 2. The data obtained after Step 1 were back propagated to the sandbox usingthe data propagation process. This process aims to “move” the data closerto the target. As a result, we obtain reasonable estimates for the locationof the buried targets, particularly in the ( x, y ) − plane, see Figure 1. Inaddition, this step helps us reduce the computational domain.Step 3. Determine an interval of frequencies on which the data obtained after Step2 are stable. (a) (b) Figure 2: a) Absolute value of measured experimental data for the buried target4 (sycamore, see table 1). b) Absolute value of the propagated data of a).
In this section we present the results of reconstructions from experimental datafor the objects buried in a sandbox in Table 1 using our convexification method.23able 1: Buried objectsNumber Description Size in x × y × z directions (in cm)1 Bamboo 3 . × . × .
82 Geode 8 . × . × .
83 Rock 10 . × . × .
04 Sycamore 3 . × . × .
85 Wet wood 9 . × . × .
86 Yellow pine 9 . × . × . k, k ] for our objects. We refer to [35] for the details of thedetermination of these intervals.The objects with their directly measured dielectric constant c meas and com-puted coefficient c comp along with corresponding measurement ε meas and com-putational errors ε comp = | c comp − c meas | /c meas ∗ c comp in Table 3 are the maximal values of the re-constructed functions c ( x ). In all our numerical tests we have used reasonablevalues of parameters µ = λ = 3 . . Considering the significant amount of noise in the measured data, the com-putational errors ε comp of reconstructed coefficients are sufficiently small. Thecomputed dielectric constant of object 3 (a piece of rock) has the biggest error ε comp = 9 . . d [35], estimated locationof objects and location of the reconstructed objects, i.e. the location of themaximum value of computed coefficient max( c comp ( x )). Errors of locations aresmall comparable with the size of the computational domain where we solve ourinverse problemFig. 3 and 4 illustrate the exact and computed images for the objects 2 and4, respectively. Images are obtained using the contour filter in Paraview. Table 3 and Figures 3, 4 demonstrate that our numerical method accuratelyreconstructs both dielectric constants and locations of targets in a quite chal-lenging case of backscatter experimental data collected for buried targets.
References [1]
H. Ammari, J. Garnier, W. Jing, H. Kang, M. Lim, K. Solna, andH. Wang , Mathematical and Statistical Methods for Multistatic Imaging ,vol. 2098 of Lecture Notes in Mathematics, Springer, Cham, 2013.24 a) (b)
Figure 3: Reconstruction result for target 2: (a) exact image, (b) computedimage (a) (b)
Figure 4: Reconstruction result for the target 4: (a) exact image, (b) computedimage 25able 2: Optimal frequencies and interval of wavenumbersNumber Optimal frequency, GHz Interval of wavenumbers [ k, k ]1 3.10 [6.322, 6.638]2 3.01 [6.133, 6.448]3 3.01 [6.070, 6.385]4 3.10 [6.322, 6.638]5 2.62 [5.313, 5.691]6 2.62 [5.313, 5.691]Table 3: Measured and reconstructed coefficients of objectsNumber c meas ε meas c comp ε comp .
99% 4.69 4 . .
13% 5.28 3 . .
3% 5.07 9 . .
89% 4.95 1 . .
69% 8.06 6 . .
54% 5.22 8 . x, y, z ) Computed location in ( x, y, z )1 (0.80, -0.11, 0.19) (0.83, 0.03, -0.05)2 (0.58, -0.14, 0.44) (0.63, 0.03, 0.16)3 (0.62, -0.14, 0.20) (0.63, 0.08, -0.20)4 (0.80, -0.04, 0.19) (1.04, 0.08, -0.30)5 (0.57, -0.42, 0.29) (0.53, -0.08, 0.16)6 (0.54, -0.33, 0.29) (0.53, -0.03, 0.21)[2] H. Ammari, Y.T. Chow, and J. Zou , The concept of heterogeneousscattering and its applications in inverse medium scattering , SIAM J. Math.Anal., 46 (2014), 2905-2935.[3]
A. B. Bakushinskii, M. V. Klibanov, and N. A. Koshev , Carlemanweight functions for a globally convergent numerical method for ill-posedCauchy problems for some quasilinear PDEs , Nonlinear Analysis: RealWorld Applications, 34 (2017), pp. 201–224.[4]
A. B. Bakushinskii, M.Yu. Kokurin and M.M. Kokurin,
Regular-ization Algorithms for Ill-Posed Problems, De Guyter, Berlin, 2018.[5]
L. Baudouin, M. de Buhan and S. Ervedoza , Convergent algorithmbased on Carleman estimates for the recovert of a potential in the waveequation , SIAM J. on Numerical Analysis, 55 (2017), 1578-1613.266]
L. Beilina and M. V. Klibanov , Approximate global convergence andadaptivity for coefficient inverse problems , Springer, 2012.[7]
M. Bellassoued and M. Yamamoto , Carleman Estimates and Appli-cations to Inverse Problems for Hyperbolic Systems , Springer Japan KK,2017.[8]
E. Burman, J. Ish-Horowicz and L. Oksanen , Fully discrete finiteelement data assimilation method for the heat equation, arxiv : 1707.06908,2017.[9]
A. Bukhgeim and M. Klibanov , Uniqueness in the large of a classof multidimensional inverse problems , Soviet Math. Doklady, 17 (1981),pp. 244–247.[10]
F. Cakoni and D. Colton , Qualitative Methods in Inverse ScatteringTheory. An Introduction , Springer, Berlin, 2006.[11]
G. Chavent , Nonlinear Least Squares for Inverse Problems - TheoreticalFoundations and Step-by-Step Guide for Applications , Springer, 2009.[12]
D. Colton and R. Kress , Inverse Acoustic and Electromagnetic Scat-tering Theory , Springer, New York, 3rd ed., 2013.[13]
A. Goncharsky and S. Romanov , Supercomputer technologies in in-verse problems of ultrasound tomography , Inverse Problems, 29 (2013),p. 075004.[14]
A. V. Goncharsky and S. Y. Romanov , Iterative methods for solvingcoefficient inverse problems of wave tomography in models with attenuation ,Inverse Problems, 33 (2017), p. 025003.[15]
K. Ito, B. Jin, and J. Zou , A direct sampling method to an inversemedium scattering problem , Inverse Problems 28 (2012), 025003.[16]
K. Ito, B. Jin, and J. Zou , A direct sampling method for inverse elec-tromagnetic medium scattering , Inverse Problems 29 (2013), 095018.[17]
S. I. Kabanikhin, A. D. Satybaev, M. Shishlenin , Direct Methods ofSolving Multidimensional Inverse Hyperbolic Problem , VSP, Utrecht, 2004.[18]
S. Kabanikhin, K. Sabelfeld, N. Novikov, M. Shishlenin , Numer-ical solution of the multidimensional Gelfand-Levitan equation , J. Inverseand Ill-Posed Problems 23 (2015) 439-450.[19]
M. V. Klibanov and F. Santosa , A computational quasi-reversibilitymethod for Cauchy problems for Laplace’s equation , SIAM J. Applied Math-ematics, 51 (1991), pp. 1653-1675.[20]
M. V. Klibanov , Global convexity in a three-dimensional inverse acousticproblem , SIAM Journal on Mathematical Analysis, 28 (1997), pp. 1371–1388.[21]
M. V. Klibanov and A. Timonov , Carleman Estimates for CoefficientInverse Problems and Numerical Applications , de Gruyter, Utrecht, 2004.2722]
M. V. Klibanov , Carleman estimates for global uniqueness, stability andnumerical methods for coefficient inverse problems , Journal of Inverse andIll-Posed Problems, 21 (2013), pp. 477–560.[23]
M. V. Klibanov , Carleman estimates for the regularization of ill-posedCauchy problems,
Applied Numerical Mathematics, 94 (2015), pp. 46–74.[24]
M. V. Klibanov , Carleman weight functions for solving ill-posed Cauchyproblems for quasilinear PDEs , Inverse Problems, 31 (2015), p. 125007.[25]
M. V. Klibanov and N. T. Th`anh , Recovering dielectric constants ofexplosives via a globally strictly convex cost functional , SIAM Journal onApplied Mathematics, 75 (2015), pp. 518–537.[26]
M. V. Klibanov and V. Romanov , Two reconstruction procedures fora 3-D phaseless inverse scattering problem for the generalized Helmholtzequation , Inverse Problems, 32 (2016), p. 0150058.[27]
M. V. Klibanov, A. E. Kolesov, L. Nguyen, and A. Sullivan , Globally strictly convex cost functional for a 1-D inverse medium scatteringproblem with experimental data , SIAM J. Appl. Math., 77 (2017), 1733-1755.[28]
M. V. Klibanov, D.-L. Nguyen, L. H. Nguyen, and H. Liu , A globallyconvergent numerical method for a 3D coefficient inverse problem with asingle measurement of multi-frequency data , Inverse Problems and Imaging ,
12 (2018), 493-523.[29]
M. V. Klibanov and A. E. Kolesov,
Convexification of a 3-D coef-ficient inverse scattering problem , Computers and Mathematics with Ap-plications, published online, https://doi.org/10.1016/j.camwa.2018.03.016,2018.[30]
A. Lakhal , A decoupling-based imaging method for inverse medium scat-tering for Maxwell’s equations , Inverse Problems, 26 (2010), 015007.[31]
J. Li, H. Liu, and Q. Wang , Enhanced multilevel linear sampling methodsfor inverse scattering problems , J. Comput. Phys., 257 (2014), pp. 554–571.[32]
J. Li, P. Li, H. Liu, and X. Liu , Recovering multiscale buried anomaliesin a two-layered medium , Inverse Problems, 31 (2015), 105006.[33]
L. A. Nazarova, L. A. Nazarov, A. L. Karchevsky, M. Vandamme , Determining kinetic parameters of a block coal bed gas by solving inverseproblem based on data of borehole gas measurements , Journal of MiningScience, 2015, Vol. 51, No. 4, pp. 666–672.[34]
A. A. Duchkov, A. L. Karchevskii,
Application of temperature mon-itoring to estimate the heat flux and thermophysical properties of bottomsediments,
Doklady Earth Sciences , 2014, Vol. 458, Part 2, p. 1285–1288.[35]
D.-L. Nguyen, M. V. Klibanov, L. H. Nguyen, and M. A. Fiddy , Imaging of buried objects from multi-frequency experimental data using aglobally convergent inversion method , J. Inverse and Ill-Posed Problems,accepted for publication (2017), available online of this journal, DOI:10.1515/jiip-2017-0047; also available at arxiv: 1705.01219, 2017.2836]
D.-L. Nguyen, M. V. Klibanov, L. H. Nguyen, A. E. Kolesov,M. A. Fiddy, and H. Liu , Numerical solution of a coefficient inverseproblem with multi-frequency experimental raw data by a globally convergentalgorithm , Journal of Computational Physics, 345 (2017), pp. 17–32.[37]
R. G. Romanov , Inverse Problems of Mathematical Physics , VSP,Utrecht, 1986.[38]
V. G. Romanov , Inverse problems for differential equations with memory ,Eurasian J. of Mathematical and Computer Applications, 2, issue 4, pp.51-80, 2014.[39]
J. A. Scales, M. L. Smith, and T. L. Fischer , Global optimiza-tion methods for multimodal inverse problems , Journal of ComputationalPhysics, 103 (1992), pp. 258–268.[40] A.N. Tikhonov, A.V. Goncharsky, V.V. Stepanov and A.G. Yagola,