A hybrid adaptive phase space method for reflection traveltime tomography
AA hybrid adaptive phase space method for reflectiontraveltime tomography
Hongkai Zhao Yimin ZhongMarch 8, 2018
Abstract
We present a hybrid imaging method for a challenging travel time tomographyproblem which includes both unknown medium and unknown scatterers in a boundeddomain. The goal is to recover both the medium and the boundary of the scatterersfrom the scattering relation data on the domain boundary. Our method is composedof three steps: 1) preprocess the data to classify them into three different categories ofmeasurements corresponding to non-broken rays, broken-once rays, and others, respec-tively, 2) use the the non-broken ray data and an effective data-driven layer strippingstrategy–an optimization based iterative imaging method–to recover the medium ve-locity outside the convex hull of the scatterers, and 3) use selected broken-once raydata to recover the boundary of the scatterers–a direct imaging method. By numericaltests, we show that our hybrid method can recover both the unknown medium andthe not-too-concave scatterers efficiently and robustly.
Traveltime tomography is an important class of inverse problems which appear in variousapplications such as global seismology [1–5], ocean acoustic tomography [6–9], ultrasoundtomography [10–12] in biomedical imaging and so on. It determines the internal velocity ofthe medium by measuring the wave traveltime between points on the boundary.Theoretically, the traveltime tomography is very closely related to boundary rigidity andlens rigidity problems in differential geometry [13–22]. The boundary rigidity problem is todetermine the metric of compact Riemannian manifold up to a diffeomorphism from firstarrival time information, and the traveltime is the length of geodesic, which is also called rayin geometric optics context, connecting two points on boundary. The lens rigidity problemutilizes multiple arrival times information to determine the Riemann metric. The multiplearrival times are encoded in scattering relation which consists of incoming and outgoingpoints and directions as well as the traveltime.For boundary rigidity, the uniqueness of reconstruction (up to an action of a diffeomor-phism) is known for simple metrics (see [13, 16, 23–25] and references therein) and manyother cases [26–30]. A compact Riemannian manifold (
M, ∂M, g ) is simple if the boundary ∂M is strictly convex with respect to its metric g and there are no conjugate points along any1 a r X i v : . [ m a t h . NA ] M a r eodesic. Moreover, for simple manifolds, the knowledge of scattering relation does not pro-vide more information than boundary distance function. See [17, 18, 31–33] and referencestherein for recent progress on lens rigidity for non-simple manifolds. Numerically, there aremany numerical algorithms motivated by the theoretical progress in boundary rigidity andlens rigidity problems, see [34–41] for algorithmic developments.When there are strong scattering effects or impenetrable obstacles inside the medium,then the geodesics could be broken. In [19], Kurylev, Lassas and Ulhmann established auniqueness result for reconstructing Riemannian metric from the broken scattering relation.For reflective obstacles, we consider an incident ray jointed with its corresponding reflectedray as a broken geodesic by imposing reflection condition at the joint point. When there isonly one strictly convex obstacle inside the manifold, then under certain conditions such assimple manifolds of dimension ≥ ≥ O (1) mismatch in the data.In [34–36], the authors have developed a phase-space approach for transmission and re-flection traveltime tomography for acoustic and elastic media by using the Stefanov-Ulhmannidentity formulated in [45]. The method is advantageous over traditional methods in inversekinematic problems [3, 46–49], because it uses multiple arrival times systematically and hasthe potential to handle anisotropic metrics as well, while these traditional methods can onlyrecover isotropic metrics by utilizing first arrival times.However, for the challenging case where both the medium and the scatterers are unknown,adaptive phase space method developed in [34] only uses non-broken ray to recover themedium outside the convex hull of the scatterers. In this work, we combine the adaptivephase space method, which is an optimization based iterative method, with a direct imagingmethod using selected broken-once ray data. This will give us the possibility to recover non-convex part of the boundary of the unknown scatterers. We also make several improvementsthat include preprocessing of the scattering relation data to classify data correspondingto non-broken rays, broken-once rays, and the others respectively and improvements inefficiency and robustness for the adaptive phase method.Although Stefanov-Ulhmann identity can be used as feedback from computed metric tothe exact metric, however, the identity itself is nonlinear. The linearization of Stefanov-Ulhmann identity will require the two metrics to be close enough, therefore the initial guessis critical for stable reconstruction. In our method, we first consider those geodesics withshort traveltimes, by taking Taylor expansion for such geodesics, we can obtain Dirichlet andNeumann data of the metric on boundary, then we can extrapolate the initial guess of metricfrom these boundary data and the initial guess should be quite close to the exact solutionnear boundary. For the construction, our method also follows the layer stripping idea, butquite different from [34], which selects the rays according to smallness in mismatch and could2nd up with some long rays which may deviate from layer stripping process. In our method,we introduce an auxiliary fidelity function to guide the layer stripping process. It can seenfrom our numerical experiments in Section 4 that the iteration number can be reduced andthe reconstruction process is very stable. For the reflection traveltime tomography, themethod in [34] will take more iterations and more time due to its trial and error strategyin distinguishing broken and nonbroken rays, while our method first preprocess the dataand directly detect non-broken rays from the scattering relation by scanning discontinuitiesin derivatives. The non-broken rays can immediately be used to reconstruct the metricoutside the convex hull of obstacles by Helgason support theorem [44]. Furthermore, whenthe obstacles are not large and not too concave or the metric does not vary too much nearobstacle, then our method can be used to capture non-convex shape of the obstacles bytracking those rays which hit the obstacle in normal direction. Such rays will reverse theirtrace back to their initial location after reflection and provide a direct and stable way oflocating points on boundaries of obstacles by tracing the ray to half of the traveltime (seethe numerical experiments in Section 4).The paper is organized as follows: we introduce the mathematical formulation for re-flection traveltime tomography and broken geodesics in Section 2. Then we describe ournumerical algorithm and the hybrid method in Section 3. Test results of our method fordifferent setups are presented in Section 4. Let (
M, g ) be a compact Riemann manifold with boundary dimension of d , and denote S ( M )its unit tangent bundle. The scattering relation or lens relation [19] is L = { (( x , ξ ) , ( y , ζ ) , t ) ∈ S ( M ) × S ( M ) × R + ∪ { } : x , y ∈ ∂M, ( γ x , ξ ( t ) , ˙ γ x , ξ ( t )) = ( y , ζ ) for some t ≥ } , (1)where γ x , ξ is the geodesic of ( M, g ) starts from x with direction ξ at t = 0.As defined in [19], a broken-once geodesic is a path α = α x , ξ , z , η ( t ) where z = γ x , ξ ( s ) ∈ M for some s ≥ η ∈ S z ( M ), and α x , ξ , z , η ( t ) = (cid:40) γ x , ξ ( t ) , for t < s,γ z , η ( t − s ) , for t ≥ s. (2)The entering and exiting points of broken geodesics define the broken scattering relation [19] R = { (( x , ξ ) , ( y , ζ ) , t ) ∈ S ( M ) × S ( M ) × R + ∪ { } : ( x , ξ ) ∈ Γ + , ( y , ζ ) ∈ Γ − ,t = l ( α x , ξ , z , η ) , and ( α x , ξ , z , η ( t ) , ˙ α x , ξ , z , η ( t )) = ( y , ζ ) for some ( z , η ) ∈ S ( M ) } , (3)where l ( α x , ξ , y , η ) ∈ R + ∪ {∞} denotes the smallest l > α x , ξ , y , η ∈ ∂M . Let ν bethe interior unit normal vector of ∂M and we define the following incoming and outgoing3ubbundles: Γ + = { ( x , ξ ) ∈ S ( M ) : x ∈ ∂M, (cid:104) ξ , ν (cid:105) g > } , Γ − = { ( x , ξ ) ∈ S ( M ) : x ∈ ∂M, (cid:104) ξ , ν (cid:105) g < } . (4)Note that the scattering relation does not contain any information about the point z or itscorresponding direction η where the broken ray α x , ξ , z , η changes its direction [19]. Let Ω ⊆ R d be a compact domain and let ( g ij ) be a Riemann metric on it. We define theHamiltonian H g by H g ( x , ξ ) = 12 (cid:32) (cid:88) ≤ i,j ≤ d g ij ( x ) ξ i ξ j − (cid:33) , (5)where ( g ij ) = ( g ij ) − . Let X (0) = ( x (0) , ξ (0) ) be the initial condition belonging to the inflowset: S − = { ( x , ξ ) | x ∈ ∂ Ω , H g ( x , ξ ) = 0 , (cid:88) ≤ i,j ≤ d g ij ξ i ν j < } , (6)where ν ( x ) is the unit outward normal vector at x ∈ ∂ Ω and ν j ( x ) is the j -th componentof ν ( x ). The geodesic X g ( s, X (0) ) satisfies following Hamiltonian system: d x ds = ∂ ξ H g , d ξ ds = − ∂ x H g , (7)with initial condition ( x (0) , ξ (0)) = X (0) . Then the solution X g ( s, X (0) ) = ( x ( s ) , ξ ( s ))defines a geodesic (or a ray) in phase space, where x ( s ) is the projection onto physical spaceΩ, ξ ( s ) is the cotangent vector at x ( s ), and s denotes the traveltime.In the following sections, we consider the case in which there are obstacles inside thedomain Ω. Rays are broken and reflected at the boundary when they hit the obstacles.Then the Hamiltonian system (7) needs to impose a jump condition at the reflection point.The jump condition for the case when there is only one obstacle strictly lying in Ω has beenderived in [34]. Let Γ be the interface where rays are reflected, for each reflected ray, there isa unique time s ∗ > x ( s ∗ ) ∈ Γ hits the interface at an incoming direction ξ in := ξ ( s ∗ ).The ray will be reflected to an outgoing direction ξ out := ( I − nn T ) ξ in , where n is unitoutward normal vector at x ( s ∗ ) of the obstacle. The Hamiltonian system for a broken-onceray will be d x ds = ∂ ξ H g , d ξ ds = − ∂ x H g , < s ≤ s ∗ , ( x (0) , ξ (0)) = X (0) ,d x ds = ∂ ξ H g , d ξ ds = − ∂ x H g , s > s ∗ , ( x ( s ∗ ) , ξ ( s ∗ )) = ( x ( s ∗ ) , ξ out ) . (8)To derive Stefanov-Ulhmann identity, we need the following Jacobian matrix with respectto the initial condition J g ( s, X (0) ) := ∂X g ∂X (0) ( s, X (0) ) = (cid:32) ∂ x ∂ x (0) ∂ x ∂ ξ (0) ∂ ξ ∂ x (0) ∂ ξ ∂ ξ (0) (cid:33) . (9)4et M = (cid:18) H ξ , x H ξ , ξ − H x , x − H x , ξ (cid:19) , (10)then J g ( s, X (0) ) satisfies system dJds = M J, J (0) = I for 0 < s < s ∗ ,dJds = M J, J ( s ∗ ) = B for s > s ∗ . (11)where B = (cid:18) J ( s ∗ ) J ( s ∗ ) ∂ ξ ξ out J ( s ∗ ) + ∂ x ξ out J ( s ∗ ) ∂ ξ ξ out J ( s ∗ ) + ∂ x ξ out J ( s ∗ ) (cid:19) . (12)Similar to [34], we consider function F ( s ) = X g ( t − s, X g ( s, X (0) )) , (13)where t = t g and (cid:90) t F (cid:48) ( s ) ds = X g ( t, X (0) ) − X g ( t, X (0) ) . (14)The left-hand side is (cid:90) t F (cid:48) ( s ) ds = (cid:90) t J g ( t − s, X g ( s, X (0) )) × ( V g − V g )( X g ( s, X (0) )) ds, (15)where V g = ( ∂ ξ H g , − ∂ x H g ). Linearize above integral’s right-hand side at metric g , then weapproximately have X g ( t, X (0) ) − X g ( t, X (0) ) (cid:39) (cid:90) t J g ( t − s, X g ( s, X (0) )) × ∂ g V ( g − g , g , X g ( s, X (0) )) ds. (16)For simplicity, we only consider isotropic metric in following sections. Let X = ( x , ξ ), then g ij = c − δ ij , ∂ g V ( λ, g, X ) = (2 cλ ξ , − ( λ ∇ c + c ∇ λ ) | ξ | ) . (17) In [34], the authors have introduced the adaptive phase space method. The numericalmethod is an iterative algorithm based on linearized Stefanov-Uhlmann identity (16), andthe algorithm automatically follows layer-stripping process by choosing those rays with smallmismatches on exiting phase measurements. However, using only mismatch informationcould deviate from layer-stripping since small mismatches do not guarantee small errors onpaths. Hence we propose a stabilized iterative method to overcome the “false picking”.5or simplicity, we first introduce the method for the medium without interior obstacle.The metric g is discretized over an underlying Eulerian grid in the physical domain Ω.The linearized Stefanov-Ulhmann identity (16) is discretized along the ray for each initialcondition X (0) in phase space, the Jacobian matrix along the ray is computed by (11), thevalue of metric g on non-grid points are linearly interpolated from neighborhood grid values.Therefore each integral equation along the ray X g ( s, X (0) ) represents a linear equation forneighboring grid values.Let X (0) k , k = 1 , , . . . , m , be initial coordinates in phase space of those measurements X g ( t k , X (0) k ), where t k is the traveltime of corresponding ray starts from X (0) k . We theniteratively construct a sequence of metric g n as follows.First, we need to construct a good initial guess g for the linearized problem, which isimportant for convergence of the metric.For short geodesic X g ( s, X (0) ) = ( x ( s ) , ξ ( s )) , ≤ s (cid:28)
1, with assumption on differentia-bility and local analyticity of g , we can easily deduce that x ( s ) = x (0) + sH g, ξ + s H (0) g, ξ , x H (0) g, ξ − H (0) g, ξ , ξ H (0) g, x ) + O ( s ) , ξ ( s ) = ξ (0) − sH (0) g, x − s H (0) g, x , x H (0) g, ξ − H (0) g, x , ξ H (0) g, x ) + O ( s ) , (18)where H (0) g = H g ( X (0) ). In case of isotropic metric, g ij = c ( x ) δ ij , choose X (0) = ( x (0) , c − ( x (0) ) v )for a unit direction v that is closest to tangential direction, and we have the following Talyorexpansions x ( s ) = x (0) + sc v + s c (2 v ( ∇ c ) T v − ∇ c ) + O ( s ) , ξ ( s ) = ξ (0) − s ∇ cc − s ∇ c ∇ c T + c ∇ c ) c v − ∇ c v T ∇ cc ) + O ( s ) . (19)By ignoring the O ( s ) term of the first approximation above at a boundary point x ∈ ∂ Ω,we can solve ∇ c ( x ) from the following linear equation x ( s ) = x (0) + sc v + s c (2 v ( ∇ c ) T v − ∇ c ) . (20)After solving ∇ c for all boundary points, we can construct a smooth initial guess g ij = c ( x ) δ ij by solving following minimization problem c = arg min ˜ c (cid:90) ∂ Ω (cid:0) |∇ c − ∇ ˜ c | + | c − ˜ c | (cid:1) ds + γ (cid:90) Ω |∇ ˜ c | dx, (21)where γ is a small regularization parameter. Define measurement mismatch d nk of k -th rayas d nk = X g ( t k , X (0) k ) − X g n ( t k , X (0) k ) . (22)By linearized Stefanov-Ulhmann identity, we define linear operator F nk along the k -th ray F nk ˜ g = J g n ( t, X (0) k ) (cid:90) t J − g n ( s, X (0) k ) ∂ g V g n (˜ g, X g n ( s, X (0) k )) ds. (23)6he linear operators F nk are matrices of size 4 × N , where N is the number of unknowns inphysical domain D , d nk are vectors of size 4 ×
1. For each g n , we define block matrix A n ofsize m × N , with each block of size 4 × A n = F n F n ... F nm , and A n g (cid:39) b n := d n d n ... d nm . (24)From a geometrical viewpoint, if the matrix F nk ’s l -th column is nonzero, then it meansthe k -th ray passes nearby the l -th grid, the total number of nonzero columns in F nk roughlyrepresents the length of this ray. During each iteration, suppose the k -th ray’s mismatch isnegligible, then it is highly possible that the matrix F nk is close to the correct one. In otherwords, the velocity field at those grid points used in the computation of the linear operator F nk in (23) along the k -th ray is likely to be accurate. In order to characterize this property,we define a fidelity function 0 ≤ p n ( x ) ≤ x ∈ Ω at n -th iteration, whichapproximately represents the confidence of the current value of g at x and initially p ( x ) ≡ F nk of matrix A n , res ( F nk ) = nnz ( F nk ) − (cid:88) i : F nk (: ,i ) (cid:54) =0 p ( x i ) , (25)where nnz ( F nk ) is the number of nonzero columns in F nk and x i , i = 1 , , . . . , N are the gridpoints. Geometrically, res ( F nk ) approximately reveals effective length of the part of the k -thray along which the velocity field is unknown. Or simply an indicator of the accuracy ofthe linearized Stefanov-Ulhmann identity (16) along the k -th ray. At n -th iteration, we use S n to represent a subset of row blocks, which corresponds to all indices k that res ( F nk ) areunder some effective length threshold r min , S n = { k ∈ { , , . . . , m }| res ( F nk ) ≤ r min } . (26)We notice that the rays belonging to S n will provide more stable reconstruction than otherrays do.At n -th iteration, we construct a perturbation ˜ g by minimizing functional H (˜ g ) = 12 (cid:107) A n ( S n , :)˜ g − b ( S n ) (cid:107) + β (cid:107)∇ ˜ g (cid:107) . (27)where β is a regularization parameter. We then update g n +1 by g n +1 = g n + ˜ g, (28)until the relative error of measurement is below certain tolerance level. The fidelity function p n +1 is updated by the following steps.1. For each k ∈ S n , we calculate the residual error along the k -th ray e k = (cid:107) A n ( k, :)˜ g − b ( k ) (cid:107) .
7. If e k is under some threshold τ , we update the fidelity for grid points involved in thelinear operator F nk by p n +1 ( x i ) = max( p n ( x i ) , − αe k ) , when A n ( k, i ) (cid:54) = × . (29)where the threshold τ is chosen to be small and α is a parameter to control the fidelitydecay. Here we have made assumption that if an effectively “short” ray has smallermismatch with measurement, then metric g ’s value along the ray has higher fidelity.This adaptive method follows a layer-stripping process by using shorter rays near-boundaryfirst and stripping them by updating fidelity function, which corresponds to the foliationprocess introduced in [43]. The layers are implicitly charaterized by p n at each iteration,and only depend on the measurements/data. We can see from the numerical experimentsin Section 4 that this method not only improves stability and efficiency over the previousadaptive phase space method introduced in [34], but also achieves a better accuracy onreconstructed metric. When there is a reflective obstacle lying inside the medium, and neither the obstacle nor themetric is known, then the reconstruction of both can be very challenging. In [34], the authorsused the adaptive phase space method to distinguish most of the unbroken rays. These rayscan help to recover the convex hull (under the same metric) of the unknown obstacle andthe medium outside the convex hull. However, for concave part of the boundary of theobstacle, one also needs to use those broken rays. In the following we propose a directimaging method to find points on the boundary of the obstacle using broken-once raysand the metric g recovered from the adaptive phase method based on non-broken rays asdescribed in the previous section.Let’s consider a simple case in this scenario. Suppose the physical domain Ω is convex and interior reflector D ⊆ Ω have C boundaries, which are homotopic equivalent to S d . Wecall a point y ∈ ∂D a radiative boundary point if a ray starts from y with direction along theoutward normal n ( y ) of ∂D only intersects with ∂ Ω. Denote the set of all radiative boundarypoints by R , assume R is nonempty, then consider a continuous mapping L g : R → ∂ Ω that L g ( y ) = P x X g ( t, ( y , c − n ( y ))) , (30)where t is the traveltime such that X g ( t, ( y , c − n ( y ))) ∈ Γ − and P x is projection mappingfrom phase space to physical space. Let the range of L g be T g , then for each x (0) ∈ T g , thereexists a direction ξ (0) and travel time s ∗ ( x (0) , ξ (0) ) such that x ( s ∗ ) ∈ ∂D, ξ out = − ξ ( s ∗ ),which means the ray hits the obstacle in the inward normal direction. Such reflected rayswill go back to its initial physical location and have exact opposite directions. Thoughnumerically such rays are rare, but they are very stable and can be used to parametrize theobstacle implicitly, because the reflection occurs exactly at halfway of a broken-once ray, seeFigure 1.We remark that once such rays exist, then an obstacle will be detected because otherwisethese non-broken rays (or geodesics) are not unique along the their initial phases.8igure 1: The broken rays hit the obstacle in (nearly) normal direction.Numerically, let X (0) k , k = 1 , , . . . , m be initial coordinates in phase space as previoussection, t k is the traveltime of corresponding ray ( x k ( s ) , ξ k ( s )) starts from X (0) k = ( x (0) k , ξ (0) k ).We use T to represent a subset of index { , , . . . , m } of rays satisfying T = (cid:110) k (cid:12)(cid:12)(cid:12) (cid:107) x (0) k − x k ( t k ) (cid:107) L + (cid:107) ξ (0) k + ξ k ( t k ) (cid:107) L < (cid:15) (cid:111) . (31)from the given data, i.e., the scattering relations for certain numerical tolerance (cid:15) > s ∗ k = t k /
2. Comparing to other broken rays, these aremore predictable and easier to use to locate the boundary due to the knowledge of reflectiontime and angle.
The adaptive phase space method [34] has shown its potential to distinguish most of thenon-broken rays during the layer-stripping process, if a broken ray or a non-broken ray isfalsely predicted, then it is likely to produce an O (1) mismatch in scattering relation dueto the jump condition (12). But this method would have difficulties in critical cases, e.g.distinguishing near-tangent broken and non-broken rays.For each physical location x (0) on the boundary of the domain, ∂ Ω, we define the set ofdirections for broken-rays (including tangential rays) by B ( x (0) ) = { ξ | ( x (0) , ξ ) ∈ S − , X g ( t, x (0) , ξ ) ∩ ∂D (cid:54) = ∅ for some t > } . (32)Let C ( x (0) ) be the smallest simply connected set containing B ( x (0) ), then we only have tofind out the directions ξ ∈ ∂ C ( x (0) ), since the rays with directions outside C ( x (0) ) are all non-broken. In order to be able to detect the set ∂ C (0) with both medium and object unknown,we need to make a few further assumptions on the metric and obstacle. For simplicity, wefocus our study on the case where the medium is isotropic, i.e., H ( x , ξ ) = ( c ( x ) (cid:107) ξ (cid:107) − c ( x ) is a C function and non-trapping.2. Both boundaries ∂ Ω and ∂D are C . Without loss of generality, we assume the bound-aries ∂ Ω and ∂D are represented by G ( x ) = 0 and F ( x ) = 0 respectively, where G and F are C functions.3. If a broken ray X g ( s, x (0) , ξ (0) ) = ( x ( s ) , ξ ( s )) is tangential to ∂D and its exiting phaseis ( y , ζ ) = ( x ( t ) , ξ ( t )), where t is traveltime, then there is a constant δ > | ζ · n ( y ) | > δ , where n ( y ) is the outward unit normal vector of ∂ Ω at y . Physicallyspeaking, this means for those the rays that are tangential to the obstacle’s boundary,they exit in non-tangential directions.4. If a broken ray X g ( s, x (0) , ξ (0) ) = ( x ( s ) , ξ ( s )) is tangential to ∂D at p = x ( t p ) ∈ ∂D , n p · ∂ ξ ∂s (cid:12)(cid:12)(cid:12) t p + (cid:18) ∂ x ∂s ∇ ∇ F (cid:107)∇ F (cid:107) ξ (cid:19) (cid:12)(cid:12)(cid:12) t p (cid:54) = 0 , (33)where n p is outward unit normal at p . This term stands the interaction between thegeometry of the obstacle and the medium. When this term is nonzero, reflection at theobstacle boundary would impose a significant change on an impinging ray directioncompared to the change of direction due to the medium variation. In the isotropiccase, this requirement is equivalent to − n p · ∇ c ( p ) c ( p ) + k · Hess( F ) (cid:107)∇ F (cid:107) (cid:12)(cid:12)(cid:12) p k (cid:54) = 0 , (34)where k = ξ (cid:107) ξ (cid:107) . The second term is bounded below by the minimal principal curvature λ min , therefore, if we assume that (cid:107)∇ c ( p ) (cid:107) c ( p ) < λ min , (35)then (33) is satisfied. Simply speaking, if the obstacle boundary’s minimal principalcurvature is not small, or the speed c varies little, then reflection’s is strong enoughto be observed. In addition, we require n p not parallel to ξ (0) and n p · ∂ x ∂ ξ (0) (cid:12)(cid:12)(cid:12) t p (cid:54) = .This condition guarantees one to find a ray starting at x (0) with a direction in a smallneighborhood of ξ (0) whose reflection by the obstacle can be observed in the scatteringrelation (see Lemma 3.4 and 3.6).By the differentiability theorem of initial value problems, we can easily show the followinglemma. Lemma 3.1.
If metric g ij = c − δ ij satisfies that c ( x ) is C k , k ≥ in Ω , then the Hamiltoniansystem’s solution X g is C k − before and after hitting the obstacle. Lemma 3.2.
For any ξ ∈ ∂ C ( x (0) ) , X g ( · , x (0) , ξ ) is tangential to ∂D . roof. We prove by contradiction. If X g ( s, x (0) , ξ ) impinges on the obstacle D with non-tangential direction, then there is an open neighborhood U at ξ such that ∀ ξ b ∈ U , X g ( s, x (0) , ξ b ) still intersects with D , which contradicts with the assumption of C ( x (0) ) con-tains B ( x (0) ). On the other hand, if X g ( s, x (0) , ξ ) is non-broken, then there is also an openneighborhood V at ξ such that ∀ ξ n ∈ V , X g ( s, x (0) , ξ n ) are non-broken, which contradictsthe smallness assumption of C ( x (0) ). Lemma 3.3.
If a broken ray X g ( s, x (0) , ξ (0) ) = ( x ( s ) , ξ ( s )) is tangential to ∂D and ξ (0) ∈ ∂ C ( x (0) ) . U is an open neighborhood of ( x (0) , ξ (0) ) , if the phase ( x (0) n , ξ (0) n ) ∈ U and the ray X g ( s, x (0) n , ξ (0) n ) = ( x n ( s ) , ξ n ( s )) is non-broken with (cid:107) ( x (0) , ξ (0) ) − ( x (0) n , ξ (0) n ) (cid:107) (cid:28) , denotethe traveltimes of X g ( s, x (0) , ξ (0) ) and X g ( s, x (0) n , ξ (0) n ) are t and t n respectively, then t n − t = − u · ( x (0) n − x (0) ) − v · ( ξ (0) n − ξ (0) ) + o ( (cid:107) x (0) n − x (0) (cid:107) ) + o ( (cid:107) ξ (0) n − ξ (0) (cid:107) ) . (36) where u and v are defined as u = (cid:18) ∇ G · ∂ x ∂s (cid:12)(cid:12)(cid:12) t (cid:19) − (cid:18) ∇ G · ∂ x ∂ x (0) (cid:12)(cid:12)(cid:12) t (cid:19) , v = (cid:18) ∇ G · ∂ x ∂s (cid:12)(cid:12)(cid:12) t (cid:19) − (cid:18) ∇ G · ∂ x ∂ ξ (0) (cid:12)(cid:12)(cid:12) t (cid:19) . (37) Proof.
The existence of such a non-broken ray in U is directly from previous lemma. At theexiting locations, G ( x ( t )) = 0 and G ( x n ( t n )) = 0. Since G is twice differentiable, we havethe following expansion, G ( x n ( t n )) = G ( x ( t )) + ∇ G · ∂ x ∂s (cid:12)(cid:12)(cid:12) t ( t n − t )+ ∇ G · ∂ x ∂ x (0) (cid:12)(cid:12)(cid:12) t ( x (0) n − x (0) ) + ∇ G · ∂ x ∂ ξ (0) (cid:12)(cid:12)(cid:12) t ( ξ (0) n − ξ (0) )+ O (( t n − t ) ) + o ( (cid:107) x (0) n − x (0) (cid:107) ) + o ( (cid:107) ξ (0) n − ξ (0) (cid:107) ) . According to the assumption 3, there exists a constant η > (cid:12)(cid:12)(cid:12) ∇ G · ∂ x ∂s (cid:12)(cid:12)(cid:12) t (cid:12)(cid:12)(cid:12) > η .Therefore we can find vectors u and v that t n − t = − u · ( x (0) n − x (0) ) − v · ( ξ (0) n − ξ (0) ) + o ( x (0) n − x (0) ) + o ( ξ (0) n − ξ (0) ) , (38)where u = (cid:18) ∇ G · ∂ x ∂s (cid:12)(cid:12)(cid:12) t (cid:19) − (cid:18) ∇ G · ∂ x ∂ x (0) (cid:12)(cid:12)(cid:12) t (cid:19) , v = (cid:18) ∇ G · ∂ x ∂s (cid:12)(cid:12)(cid:12) t (cid:19) − (cid:18) ∇ G · ∂ x ∂ ξ (0) (cid:12)(cid:12)(cid:12) t (cid:19) . (39) Lemma 3.4.
Suppose ξ ∈ S d − is a fixed vector and v (cid:54) = is not parallel to ξ . For anyopen set B ⊂ S d − sup ζ ∈ B (cid:107) ξ − ζ (cid:107) < ε (cid:28) , there exists ˆ ζ ∈ B such that | v · ( ξ − ˆ ζ ) | = O ( (cid:107) ξ − ˆ ζ (cid:107) ) . (40)11 roof. Without loss of generality, we assume (cid:107) v (cid:107) = 1. Let r = ξ · v , | r | <
1. Since B is anopen set in S d − , there exists ˆ ζ ∈ B , ˆ ζ (cid:54) = ξ such that | ( ξ − ˆ ζ ) · ( r ξ − v ) | ≥ ˆ γt √ − r for someˆ γ >
0, where t = (cid:107) ξ − ˆ ζ (cid:107) < (cid:15) . Since | ( ξ − ˆ ζ ) · ξ | = O ( t ), we have | v · ( ξ − ˆ ζ ) | = O ( (cid:107) ξ − ˆ ζ (cid:107) ). Lemma 3.5.
For any initial phase X (0) = ( x (0) , ξ (0) ) , suppose the ray X g ( s, x (0) , ξ (0) ) = ( x ( s ) , ξ ( s )) , x (0) = x (0) , ξ (0) = ξ (0) is non-broken on s ∈ (0 , T ) . If vector n satisfies ∂ x ( t ) ∂ ξ (0) n = , ∂ ξ ( t ) ∂ ξ (0) n = , ∀ t ∈ (0 , T ) (41) then n = .Proof. Since the determinant of the Jacobian matrix J = (cid:32) ∂ x ∂ x (0) ∂ x ∂ ξ (0) ∂ ξ ∂ x (0) ∂ ξ ∂ ξ (0) (cid:33) is 1, if the vector n (cid:54) = , then J (cid:18) (cid:19) = (42)which contradicts to the non-degeneracy of J . Lemma 3.6.
From a fixed point x (0) ∈ ∂ Ω , we denote the ray with initial direction ζ as X g ( s, x (0) , ζ ) = ( x ( s, ζ ) , ξ ( s, ζ )) and let T ( ζ ) be the traveltime. Suppose ξ (0) ∈ ∂ C ( x (0) ) ,then either ∂ ζ T or ∂ ζ x or ∂ ζ ξ is discontinuous at ξ (0) .Proof. Consider a small open neighborhood V ⊂ S − of ξ (0) , then the following sets N = { ζ ∈ V | X g ( · , x (0) , ζ ) is nonbroken } ,B = { ζ ∈ V | X g ( · , x (0) , ζ ) is broken } . (43)are both nonempty sets. Then by Lemma 3.3, for any (cid:15) >
0, we can select ξ (0) n ∈ N that (cid:107) ξ (0) n − ξ (0) (cid:107) < (cid:15) and | T ( ξ (0) ) − T ( ξ (0) n ) | = O ( (cid:107) ξ (0) − ξ (0) n (cid:107) ) . (44)and the differences between the exiting locations and phases are also of order O ( (cid:107) ξ (0) − ξ (0) n (cid:107) ).On the other hand, suppose the broken ray X g ( s, x (0) , ξ (0) ) is tangential to ∂D at point p = x ( t p , ξ (0) ), let z = ∇ F ( p ) · ∂ x ∂ ζ (cid:12)(cid:12)(cid:12) t p , ξ (0) , then by the assumption 4, z (cid:54) = and not parallel to ξ (0) , using the Lemma 3.4, For the same (cid:15) , we can select ξ (0) b (cid:54) = ξ (0) with ξ (0) b ∈ B , (cid:107) ξ (0) b − ξ (0) (cid:107) < (cid:15) and z · ( ξ (0) b − ξ (0) ) = O ( (cid:107) ξ (0) b − ξ (0) (cid:107) ) . (45)12igure 2: The illustration of the tangential ray and broken ray, both rays start from the samephysical locations but different directions. The tangential ray intersects with the obstacleat p and the other broken ray intersects with the obstacle at q .Then the ray with initial phase X (0) b = ( x (0) , ξ (0) b ) is broken. Note that the tangent ray X g ( · , x (0) , ξ (0) ) satisfies ∇ F ( p ) · ∂ x ∂s (cid:12)(cid:12)(cid:12) t p , ξ (0) = 0 ,F ( p ) = 0 . (46)Assume the broken ray X g ( s, x (0) , ξ (0) b ) intersects ∂D at q = x ( t q , ξ (0) b ) as illustrated inFigure 2, then take Taylor expansion at ( t p , ξ (0) ),0 = F ( x ( t q , ξ (0) b )) = F ( x ( t p , ξ (0) ))+ ∇ F ( p ) · ∂ x ∂s (cid:12)(cid:12)(cid:12) ( t p , ξ (0) ) ( t q − t p ) + ∇ F ( p ) · ∂ x ∂ ζ (cid:12)(cid:12)(cid:12) ( t p , ξ (0) ) ( ξ (0) b − ξ (0) )+ o ( (cid:107) ξ (0) b − ξ (0) (cid:107) ) + O (( t q − t p ) ) . (47)The first two terms on the right hand side are zero according (46). By using (45) ∇ F ( p ) · ∂ x ∂ ζ (cid:12)(cid:12)(cid:12) ( t p , ξ (0) ) ( ξ (0) b − ξ (0) ) = z · ( ξ (0) b − ξ (0) ) = O ( (cid:107) ξ (0) b − ξ (0) (cid:107) ) , we then conclude | t p − t q | = O (cid:18)(cid:113) (cid:107) ( ξ (0) b − ξ (0) ) (cid:107) (cid:19) . (48)After reflection at time t q , the direction of the broken ray X g ( s, x (0) , ξ (0) b ) has been reflectedto ξ b, out = ( I − n q ⊗ n q ) ξ ( t − q , ξ (0) b ) , n q = ∇ F ( q ) (cid:107)∇ F ( q ) (cid:107) is the outward unit normal vector at q . Then by taking Taylor expan-sion at ( t p , ξ (0) ), the difference between the directions of the two rays X g ( t + q , x (0) , ξ (0) ) and X g ( t + q , x (0) , ξ (0) b ) is ξ ( t q , ξ (0) ) − ξ b, out = 2 n p (cid:18) n p · ∂ ξ ∂s (cid:12)(cid:12)(cid:12) t p , ξ (0) + (cid:18) ∂ x ∂s ∇ ∇ F (cid:107)∇ F (cid:107) ξ (cid:19) (cid:12)(cid:12)(cid:12) t p , ξ (0) (cid:19) ( t q − t p )+ O (( t q − t p ) ) . (49)where n p = ∇ F ( p ) (cid:107)∇ F ( p ) (cid:107) is the outward normal vector at p . According to assumption 4, the firstterm on right hand side does not vanish, therefore ξ ( t q , ξ (0) ) − ξ b, out = n p · O ( t q − t p ) = n p · O (cid:18)(cid:113) (cid:107) ( ξ (0) b − ξ (0) ) (cid:107) (cid:19) . (50)Regarding the phases at time t + q as initial phases, then we can define the rest of the tangentialray as X g ( s, x ( t q , ξ (0) ) , ξ ( t q , ξ (0) )) = ( y ( s ) , θ ( s )) , y (0) = x ( t q , ξ (0) ) , θ (0) = ξ ( t q , ξ (0) ) , (51)and denote the traveltime of this partial ray as t y . For the other ray, we define the rest ofthe broken ray as X g ( s, q , ξ b, out ) = ( w ( s ) , η ( s )) , w (0) = q , η (0) = ξ b, out , (52)and similarly, the traveltime is denoted as t w . Then use the Lemma 3.3 and (50) T ( ξ (0) ) − T ( ξ (0) b ) = ( t q + t y ) − ( t q + t w )= − u · ( y (0) − w (0)) − v · ( θ (0) − η (0))+ o ( (cid:107) y (0) − w (0) (cid:107) ) + o ( (cid:107) θ (0) − η (0) (cid:107) )= − ( v · n p ) · O (cid:18)(cid:113) (cid:107) ( ξ (0) b − ξ (0) ) (cid:107) (cid:19) + O ( (cid:107) ( ξ (0) b − ξ (0) ) (cid:107) ) , (53)where u and v are defined as following, u = (cid:18) ∇ G · ∂ y ∂s (cid:12)(cid:12)(cid:12) t y (cid:19) − (cid:18) ∇ G · ∂ y ∂ y (0) (cid:12)(cid:12)(cid:12) t y (cid:19) , v = (cid:18) ∇ G · ∂ y ∂s (cid:12)(cid:12)(cid:12) t y (cid:19) − (cid:18) ∇ G · ∂ y ∂ θ (0) (cid:12)(cid:12)(cid:12) t y (cid:19) . (54)We consider following three cases,1. If v · n p (cid:54) = 0, then ∂ ζ T is discontinuous at ξ (0) . Otherwise, T ( ξ (0) ) − T ( ξ (0) b ) = t y − t w = O ( (cid:107) ( ξ (0) b − ξ (0) ) (cid:107) ) and we consider the next case.14. Use (50) and (cid:107) x ( t q , ξ (0) ) − q (cid:107) = O ( (cid:107) ( ξ (0) b − ξ (0) ) (cid:107) ), the difference between exitingphysical locations is y ( t y ) − w ( t w ) = ∂ y ∂s (cid:12)(cid:12)(cid:12) t y ( t y − t w )+ ∂ y ∂ y (0) (cid:12)(cid:12)(cid:12) t y ( x ( t q , ξ (0) ) − q ) + ∂ y ∂ θ (0) (cid:12)(cid:12)(cid:12) t y ( ξ ( t q , ξ (0) ) − ξ b, out )+ O (( t y − t w ) ) + O ( (cid:107) x ( t q , ξ (0) ) − q (cid:107) ) + O ( (cid:107) ( ξ ( t q , ξ (0) ) − ξ b, out (cid:107) )= (cid:18) ∂ y ∂ θ (0) (cid:12)(cid:12)(cid:12) t y · n p (cid:19) · O (cid:18)(cid:113) (cid:107) ( ξ (0) b − ξ (0) ) (cid:107) (cid:19) + O ( (cid:107) ( ξ (0) b − ξ (0) ) (cid:107) ) . (55)If ∂ y ∂ θ (0) (cid:12)(cid:12)(cid:12) t y · n p is nonzero vector, then the derivative ∂ ζ x will suffer from a discontinuityat y ( t y ), otherwise if ∂ y ∂ θ (0) (cid:12)(cid:12)(cid:12) t y · n p = , we consider the next case.3. Similarly, the difference between exiting directions is θ ( t y ) − η ( t w ) = ∂ θ ∂s (cid:12)(cid:12)(cid:12) t y ( t y − t w )+ ∂ θ ∂ y (0) (cid:12)(cid:12)(cid:12) t y ( x ( t q , ξ (0) ) − q ) + ∂ θ ∂ θ (0) (cid:12)(cid:12)(cid:12) t y ( ξ ( t q , ξ (0) ) − ξ b, out )+ O (( t y − t w ) ) + O ( (cid:107) x ( t q , ξ (0) ) − q (cid:107) ) + O ( (cid:107) ( ξ ( t q , ξ (0) ) − ξ b, out (cid:107) )= (cid:18) ∂ θ ∂ θ (0) (cid:12)(cid:12)(cid:12) t y · n p (cid:19) · O (cid:18)(cid:113) (cid:107) ( ξ (0) b − ξ (0) ) (cid:107) (cid:19) + O ( (cid:107) ( ξ (0) b − ξ (0) ) (cid:107) ) . (56)If ∂ y ∂ θ (0) (cid:12)(cid:12)(cid:12) t y · n p is nonzero vector, then the derivative ∂ ζ ξ will suffer from a discontinuityat θ ( t y ). Otherwise we will have following equations ∂ y ∂ θ (0) (cid:12)(cid:12)(cid:12) t y · n p = 0 ,∂ θ ∂ θ (0) (cid:12)(cid:12)(cid:12) t y · n p = 0 , (57)by Lemma 3.5, we must have n p = , which is a contradiction. Therefore either ∂ ζ T or ∂ ζ x or ∂ ζ ξ must have a discontinuity at ξ (0) .With these assumptions, we can directly detect non-broken rays from the measurementsby scanning the traveltimes, exiting directions and exiting locations for jumps in the deriva-tives with respect to initial directions, see Figure 3. And using these non-broken rays enablesus to recover the metric outside of the convex hull of the obstacle [42, 50].15igure 3: The plots of traveltimes, exiting directions, exiting locations from numerical ex-ample 4.3.2. Left: From top to bottom, the plots are the traveltimes, exiting directions,exiting locations corresponding to the initial phases with varying directions at the 15thboundary point, the shadowed spots are placed at the detected jumps near 52nd and 195thrays respectively. Right: From top to bottom, the plots are the traveltimes, exiting direc-tions, exiting locations corresponding to the initial phases with varying directions at the30th boundary point, the shadowed spots are placed at the detected jumps near 98th and175th rays respectively. In this section, we present the hybrid method for reconstructing both the metric and includedobstacles. First we find out if there exists an obstacle inside the medium by checking if theset T defined in (31) is empty.If no obstacle is detected, then we can use the improved adaptive phase space methodin Section 3.1 to recover the metric efficiently and stably with the layer-stripping strategy.16nce an obstacle is detected, then we distinguish the broken rays and non-broken onesby scanning the exiting traveltimes, locations and directions in scattering relation explainedin Section 3.3 and use the improved adaptive phase space method in Section 3.1 on non-broken rays to recover the metric. In this case, our layer-stripping reconstruction strategywill be able to recover the metric starting from the boundary and continuing inward all theway to the convex hull of the obstacle. Since our adaptive phase space method is based onoptimization formulation with a regularization (27) for the metric at all grid points in thedomain, the numerically reconstructed metric in the whole domain can be viewed as a goodapproximation of the true metric outside the convex hull of the obstacle plus a harmonicextension to the interior of the convex hull.If the obstacle is convex (under the metric), one can reconstruct both the metric and theobstacle using non-broken rays as we can see from the numerical experiments in Section 4.Since non-broken rays only contain information outside the convex hull of the obstacle, it isimpossible to reconstruct a concave obstacle with only non-broken rays.On the other hand, since our adaptive phase method using non-broken rays reconstructsthe metric on the whole domain, then by tracing back the rays in T (defined in (31)) tohalf of the traveltime, we will get the approximated reflection points of such rays, if the truemetric varies slowly inside the convex hull. This gives us a direct imaging method for theboundary of the obstacle, see the numerical experiments in Section 4.For the cases that the metric has large variations inside the convex hull, our methodthen can not recover the obstacle and metric inside the convex hull without using otherbroken rays. One possible way is to introduce an representation of the obstacle’s boundaryand iteratively reconstruct the metric inside the convex hull as well as morph the boundarysimultaneously to minimize the mismatch, e.g., using the result from our hybrid method asan initial guess. However, this will be a daunting task due to the highly non-convex andcoupled optimization problem.Finally, we briefly discuss the computational cost of our hybrid method in 2D. Suppose wehave N s sources and each source probes N a directions, then there are N s N a scattering relationmeasurements. For obstacle detection, it will take O ( N s N a ) complexity at the worst case.For non-broken rays detection, it will take O ( N s N a ) complexity due to linear scan. At eachiteration of the stabilized adaptive phase space method, we have to solve the Hamiltoniansystem for X g ( s, X (0) ) in (7) and Jacobian matrix (11), which has the worst complexity as O ( T max N s N a ), where T max is length of the longest geodesic. These solutions are then usedto calculate mismatch and formulate the linearized Stefanov-Ulhmann identity (16) over theEulerian grid in (24) with the worst complexity O ( T max N s N a ). Then we use the standardmultifrontal solver umfpack to solve the perturbation for the minimization problem (27),the complexity is O ( n ) in general, where n is the number of unknowns. Therefore the totaltime complexity is O ( K ( T max N s N a + n )), where K is the number of iterations. All numerical experiments are implemented in
Julia and performed on a dual-core laptop of2 . GHz
CPU and 16
GByte memory. Source code is hosted on https://github.com/lowrank/ray.We take the physical domain Ω as unit disk for all examples. The discretization of metric17 over a uniform grid is parametrized by Q4 element. If a ray passes through a grid, thenit will involve 12 surrounding grid values, under such situation we can set rank threshold r min = 12. For other parameters, our selections are conservative, we take τ = 5% and α = 10for fidelity function updating, and regularization parameter β = 0 . (cid:15) = 0 .
5% for obstacle detection (see Section 3.2). We keep them fixedfor all of the examples.
In this scenario, we experiment our improved adaptive phase space method described inSection 3.1 on simple cases without interior obstacle.
The exact solution is c ( x, y ) = 1 + 0 . πx ) sin( πy ) . The grid’s resolution is h = 1 /
15. We put 50 equispaced sources and each source probes 100uniformly distributed directions. The method converges to a solution with relative L error2 . × − at 11th iteration. We plot the numerical and the exact solutions in Figure 4.Figure 4: Left: the exact solution. Middle: the numerical solution at 11th iteration. Right:the error between numerical solution and exact solution. The exact solution is c ( x, y ) = 1 + 0 . . πx ) sin(1 . πy ) . The grid’s resolution is h = 1 /
25. We put 100 equispaced sources and each source probes 100uniformly distributed directions. The method converges to a solution with relative L error3 . × − at 22th iteration. We plot the numerical and the exact solutions in Figure 5.We also plot the auxiliary fidelity function at three different iterations to illustrate the layerstripping process in Figure 6. 18igure 5: Left: the exact solution. Middle: the numerical solution at 22th iteration. Right:the error between numerical solution and exact solution.Figure 6: The fidelity function p n in different iterations. From left to right: 6th, 11th, 17thiteration. In this scenario, we experiment our hybrid method for imaging an unknown convex obstacleinside an unknown metric.
We consider the obstacle as a circle at center x + y = 116 , and the exact solution is given by c ( x, y ) = 1 + 0 . (cid:16) π (cid:112) ( x − . + ( y − . (cid:17) + 0 . (cid:16) π (cid:112) ( x + 0 . + ( y + 0 . (cid:17) . The grid’s resolution is 1 /
15. We put 50 equispaced sources and each source probes 300uniformly distributed directions. We first detect the obstacle by checking the scatteringrelation as in Section 3.2, and then distinguish the non-broken rays from the broken raysas in Section 3.3, see Figure 7. And then we use all the non-broken rays to reconstruct themetric by the stabilized adaptive phase space method, the method converges to a solutionwith relative L error 4 . × − at 9th iteration. We plot the numerical and exact solutionsin Figure 8. From the experiment, we can see that the metric outside of the obstacle has19igure 7: Left: The broken rays in T , which hit the obstacle in (nearly) normal direction inExample 3. Right: Detected tangent rays in Example 3.Figure 8: Left: exact solution. Middle: numerical solution at 9th iteration. Right: errorbetween numerical solution and exact solution.been recovered well. Then the obstacle’s convex hull can be approximated by the envelopeof all the tangent rays computed through the recovered metric, see Figure 9. After thereconstruction of the convex hull of the obstacle and the metric outside the convex hull, wetrace the rays in collection T (defined in (31)) to half of the traveltime to get the reflectionpoints on the boundary, see also in Figure 9. We can see that the computed reflection pointsare quite close to the boundary. In this scenario, we will use our hybrid method to recover a non-convex unknown interiorobstacle and the underlying metric.
In this example, we consider an easier case. The obstacle’s boundary is parameterized inpolar coordinate ( r, θ ) as r ( θ ) = 0 .
25 + 0 .
05 sin(3 θ ) , T to half traveltime, the dashed blue circle at center is theexact boundary.which is a slightly concave shape, the exact solution is given by c ( x, y ) = 1 + 0 . (cid:16) π (cid:112) ( x − . + ( y − . (cid:17) + 0 . (cid:16) π (cid:112) ( x + 0 . + ( y + 0 . (cid:17) . The grid’s resolution is 1 /
15. We put 50 equispaced sources and each source probes 300uniformly distributed directions. From the scattering relation, we can directly extract therays that hit the obstacle in almost normal direction, and also distinguish the non-brokenrays by detecting the jumps in scattering relation. We plot those rays in Figure 10. Then weFigure 10: Left: The rays that hit the obstacle in (nearly) normal direction in Example 4.Right: The tangent rays detected from scattering relation.follow the method in Section 3.3 to distinguish the non-broken rays and broken rays. Thenwe use all the non-broken rays to recover the metric outside the obstacle by the stabilizedadaptive phase space method. The method converges to a solution with relative L error5 . × − at 9th iteration. We plot the numerical and exact solutions in Figure 11.Since the error of metric is small outside the obstacle, then the obstacle’s convex hullcan be approximated well by all the tangent rays computed through the recovered metric,see Figure 12. By tracing back the rays in T , we approximately obtain the reflection points21igure 11: Left: exact solution. Middle: numerical solution at 9th iteration. Right: errorbetween the numerical solution and the exact solution.Figure 12: Left: The tangent rays computed from recovered metric in Example 4. Right:The traced rays in collection T to half of traveltime. The blue dashed line is the exactboundary of obstacle.on the obstacle, also see Figure 12. However, since no information is available inside theconvex hull, the error of reflection points can be large in general. In this example, we take a more challenging obstacle. The obstacle’s boundary is parame-terized in polar coordinate ( r, θ ) as r ( θ ) = 0 . . θ ) , which is a more concave shape than previous example, and the exact solution is again givenby c ( x, y ) = 1 + 0 . (cid:16) π (cid:112) ( x − . + ( y − . (cid:17) + 0 . (cid:16) π (cid:112) ( x + 0 . + ( y + 0 . (cid:17) . The grid’s resolution is 1 /
15. We put 50 equispaced sources and each source probes 300uniformly distributed directions. From the scattering relation, we can directly extract therays that hit the obstacle in almost normal direction, and also distinguish the non-brokenrays by detecting the jumps in scattering relation. We plot such rays in Figure 13. Then22igure 13: Left: The rays that hit the obstacle in (nearly) normal direction in Example 5.Right: The tangent rays detected from scattering relation.we use the stabilized adaptive phase space method in Section 3.1 on all the non-brokenrays to recover the metric as much as possible. The method converges to a solution withrelative L error of 2 . × − at 9th iteration. We plot the numerical and exact solutions inFigure 14. After having recovered the metric from the non-broken rays’ scattering relation,Figure 14: Left: exact solution. Middle: numerical solution at 9th iteration. Right: errorbetween the numerical solution and the exact solution.we can approximate the convex hull of the obstacle by the recovered non-broken rays, seeFigure 15. And by tracing the rays in collection T , we can approximately obtain the reflectionpoints on the boundary of the obstacle, also see Figure 15. In this work, we proposed a hybrid phase space method for traveltime tomography whichincludes both an unknown medium and unknown scatterer. The underlying medium outsidethe convex hull of the scatterer is reconstructed by a optimization based iterative method.The newly developed method is more stable than the previous adaptive phase method pro-posed in [34] due to the introduction of an auxiliary fidelity function to guide the layerstripping process and a direct detection of all non-broken rays. To image the boundary ofthe scatterer, we use a direct imaging method that can locate points on the boundary of23igure 15: Left: The tangent rays computed from recovered metric in Example 5. Right:The traced rays in collection T to half of traveltime. The blue dashed line is the exactboundary of obstacle.the scatterer by selecting those broken-once rays that hit the scatterer almost normally andtracing back those rays to half traveltime in the reconstructed medium. Acknowledgement
H. Zhao is partially supported by NSF grant DMS-1418422. Both authors would like tothank ICERM 2017 Fall program on Mathematical and Computational Challenges in Radarand Seismic Reconstruction, where this project was started. The authors also would like tothank Kui Ren for valuable discussions.
References [1] BLN Kennett and ER Engdahl. Traveltimes for global earthquake location and phaseidentification.
Geophysical Journal International , 105(2):429–465, 1991.[2] Hiroshi Inoue, Yoshio Fukao, Kunio Tanabe, and Yosihiko Ogata. Whole mantle p-wavetravel time tomography.
Physics of the Earth and Planetary Interiors , 59(4):294–328,1990.[3] TN Bishop, KP Bube, RT Cutler, RT Langan, PL Love, JR Resnick, RT Shuey,DA Spindler, and HW Wyld. Tomographic determination of velocity and depth inlaterally varying media.
Geophysics , 50(6):903–923, 1985.[4] Richard A Clarke, Bertrand Alazand, Laure Pelle, Delphine Sinoquet, Patrick Lailly,Florence Delprat-Jannaud, and Lionel Jannaud. 3d traveltime reflection tomographywith multi-valued arrivals. In
SEG Technical Program Expanded Abstracts 2001 , pages1875–1878. Society of Exploration Geophysicists, 2001.[5] Susan E Minkoff. A computationally feasible approximate resolution matrix for seismicinverse problems.
Geophysical Journal International , 126(2):345–359, 1996.246] Walter Munk and Carl Wunsch. Ocean acoustic tomography: A scheme for large scalemonitoring.
Deep Sea Research Part A. Oceanographic Research Papers , 26(2):123–161,1979.[7] Walter Munk, Peter Worcester, and Carl Wunsch.
Ocean acoustic tomography . Cam-bridge University Press, 2009.[8] MD Collins and WA Kuperman. Inverse problems in ocean acoustics.
Inverse Problems ,10(5):1023, 1994.[9] Finn B Jensen, William A Kuperman, Michael B Porter, and Henrik Schmidt.
Com-putational ocean acoustics . Springer Science & Business Media, 2011.[10] Ali Hormati, Ivana Jovanovic, Olivier Roy, and Martin Vetterli. Robust ultrasoundtravel-time tomography using the bent ray model. In
Proceedings of the SPIE MedicalImaging , number LCAV-CONF-2010-001. Spie-Int Soc Optical Engineering, Po Box 10,Bellingham, WA 98227-0010 USA, 2010.[11] Hermann Schomberg. An improved approach to reconstructive ultrasound tomography.
Journal of Physics D: Applied Physics , 11(15):L181, 1978.[12] Xing Jin and Lihong V Wang. Thermoacoustic tomography with correction for acousticspeed variations.
Physics in Medicine and Biology , 51(24):6437, 2006.[13] Christopher B Croke et al. Rigidity and the distance between boundary points.
Journalof Differential Geometry , 33(2):445–464, 1991.[14] Bela Frigyik, Plamen Stefanov, and Gunther Uhlmann. The x-ray transform for ageneric family of curves and weights.
Journal of Geometric Analysis , 18(1):89–108,2008.[15] Plamen Stefanov, Gunther Uhlmann, et al. Stability estimates for the x-ray transformof tensor fields and boundary rigidity.
Duke Mathematical Journal , 123(3):445–467,2004.[16] Plamen Stefanov and Gunther Uhlmann. Boundary rigidity and stability for genericsimple metrics.
Journal of the American Mathematical Society , 18(4):975–1003, 2005.[17] Plamen Stefanov and Gunther Uhlmann. Boundary and lens rigidity, tensor tomographyand analytic microlocal analysis.
Algebraic Analysis of Differential Equations, Fetschriftin Honor of Takahiro Kawai, edited by T. Aoki, H. Majima, Y. Katei and N. Tose , pages275–293, 2008.[18] Colin Guillarmou. Lens rigidity for manifolds with hyperbolic trapped sets.
Journal ofthe American Mathematical Society , 30(2):561–599, 2017.[19] Yaroslav Kurylev, Matti Lassas, and Gunther Uhlmann. Rigidity of broken geodesicflow and inverse problems.
American journal of mathematics , 132(2):529–562, 2010.2520] Leonid Pestov and Gunther Uhlmann. On characterization of the range and inversionformulas for the geodesic x-ray transform.
International mathematics research notices ,2004(80):4331–4347, 2004.[21] Plamen Stefanov and Gunther Uhlmann. The geodesic x-ray transform with fold caus-tics.
Analysis & PDE , 5(2):219–260, 2012.[22] L Pestov and G Uhlmann. Two dimensional simple compact manifolds with boundaryare boundary rigid.
Annals of Math , 161:1089–1106, 2005.[23] Ravil G Mukhometov. A problem of reconstructing a riemannian metric.
SiberianMathematical Journal , 22(3):420–433, 1981.[24] Ren´e Michel. Sur la rigidit´e impos´ee par la longueur des g´eod´esiques.
Inventionesmathematicae , 65(1):71–83, 1981.[25] Plamen Stefanov and Gunther Uhlmann. Recent progress on the boundary rigidityproblem.
Electronic research announcements of the American Mathematical Society , 11(8):64–70, 2005.[26] Mikhael Gromov et al. Filling riemannian manifolds.
Journal of Differential Geometry ,18(1):1–147, 1983.[27] G´erard Besson, Gilles Courtois, and Sylvestre Gallot. Entropies et rigidit´es des espaceslocalement sym´etriques de courbure strictement n´egative.
Geometric and functionalanalysis , 5(5):731–799, 1995.[28] Christopher B Croke. Rigidity for surfaces of non-positive curvature.
CommentariiMathematici Helvetici , 65(1):150–169, 1990.[29] Matti Lassas, Vladimir Sharafutdinov, and Gunther Uhlmann. Semiglobal boundaryrigidity for riemannian metrics.
Mathematische Annalen , 325(4):767–793, 2003.[30] Vladimir Altafovich Sharafutdinov.
Integral geometry of tensor fields , volume 1. Walterde Gruyter, 1994.[31] Plamen Stefanov and Gunther Uhlmann. Local lens rigidity with incomplete data fora class of non-simple riemannian manifolds. arXiv preprint math/0701595 , 2007.[32] Plamen Stefanov. Microlocal approach to tensor tomography and boundary and lensrigidity.
Serdica Mathematical Journal , 34(1):67p–112p, 2008.[33] Plamen Stefanov and Gunther Uhlmann. Integral geometry of tensor fields on a class ofnon-simple riemannian manifolds.
American journal of mathematics , 130(1):239–268,2008.[34] Eric Chung, Jianliang Qian, Gunther Uhlmann, and Hongkai Zhao. An adaptive phasespace method with application to reflection traveltime tomography.
Inverse problems ,27(11):115002, 2011. 2635] Eric Chung, Jianliang Qian, Gunther Uhlmann, and Hong-Kai Zhao. A phase-spaceformulation for elastic-wave traveltime tomography. In
Journal of Physics: ConferenceSeries , volume 124, page 012018. IOP Publishing, 2008.[36] Eric Chung, Jianliang Qian, Gunther Uhlmann, and Hongkai Zhao. A new phase spacemethod for recovering index of refraction from travel times.
Inverse Problems , 23(1):309, 2007.[37] Shingyu Leung, Jianliang Qian, et al. An adjoint state method for three-dimensionaltransmission traveltime tomography using first-arrivals.
Communications in Mathemat-ical Sciences , 4(1):249–266, 2006.[38] Shingyu Leung and Jianliang Qian. Transmission traveltime tomography based onparaxial liouville equations and level set formulations.
Inverse Problems , 23(2):799,2007.[39] Wenbin Li and Shingyu Leung. A fast local level set adjoint state method for firstarrival transmission traveltime tomography with discontinuous slowness.
GeophysicalJournal International , 195(1):582–596, 2013.[40] Wenbin Li, Shingyu Leung, and Jianliang Qian. A level-set adjoint-state method forcrosswell transmission-reflection traveltime tomography.
Geophysical Journal Interna-tional , 199(1):348–367, 2014.[41] Roland Glowinski, Shingyu Leung, and Jianliang Qian. A penalization-regularization-operator splitting method for eikonal based traveltime tomography.
SIAM Journal onImaging Sciences , 8(2):1263–1292, 2015.[42] Venkateswaran P Krishnan. A support theorem for the geodesic ray transform onfunctions.
Journal of Fourier Analysis and Applications , 15(4):515–520, 2009.[43] Gunther Uhlmann and Andr´as Vasy. The inverse problem for the local geodesic raytransform.
Inventiones mathematicae , 205(1):83–120, 2016.[44] Sigurdur Helgason. The radon transform on r n. In
Integral Geometry and RadonTransforms , pages 1–62. Springer, 2011.[45] Plamen Stefanov and Gunther Uhlmann. Rigidity for metrics with the same lengths ofgeodesics.
Mathematical Research Letters , 5:83–96, 1998.[46] Vladimir Gavrilovich Romanov.
Inverse problems of mathematical physics . Brill, 1987.[47] Alain Sei and William W Symes. Gradient calculation of the traveltime cost functionwithout ray tracing. In
SEG Technical Program Expanded Abstracts 1994 , pages 1351–1354. Society of Exploration Geophysicists, 1994.[48] Alain Sei and William W Symes. Convergent finite-difference traveltime gradient fortomography. In
SEG Technical Program Expanded Abstracts 1995 , pages 1258–1261.Society of Exploration Geophysicists, 1995.2749] John K Washbourne, James W Rector, and Kenneth P Bube. Crosswell traveltimetomography in three dimensions.
Geophysics , 67(3):853–871, 2002.[50] Joonas Ilmavirta and Mikko Salo. Broken ray transform on a riemann surface with aconvex obstacle. arXiv preprint arXiv:1403.5131arXiv preprint arXiv:1403.5131