Quantitative inverse scattering via reduced order modeling Liliana Borcea Mathematics, University of Michigan Ann Arbor Joint work with: • Vladimir Druskin Worcester Polytechnic Institute • Alexander Mamonov University of Houston • Mikhail Zaslavsky Schlumberger • J¨ orn Zimmerling University of Michigan Support: U.S. Office of Naval Research N00014-17-1-2057. 1
Inverse scattering for generic hyperbolic system • Primary wave P ( s ) ( t, x ) and dual wave � P ( s ) ( t, x ) satisfy � � � � � � P ( s ) ( t, x ) P ( s ) ( t, x ) 0 − L q t > 0 , x ∈ Ω ⊂ R d ∂ t = , L T P ( s ) ( t, x ) � P ( s ) ( t, x ) � 0 q with homogeneous boundary conditions and initial conditions P ( s ) (0 , x ) = b ( s ) ( x ) , P ( s ) (0 , x ) = 0 . � • Ω is half space ( x d > 0) and measurements on ∂ Ω + at x d = 0 + . • Array of sensors on ∂ Ω + . Source excitation b ( s ) ( x ) supported near ∂ Ω + , where ( s ) counts sensor index and polarization. • Finite duration of measurements � can truncate Ω to compact cube Ω c ⊂ Ω with accessible boundary ∂ Ω ac ⊂ ∂ Ω and inaccessi- c ble boundary ∂ Ω inac ⊂ Ω. c 2
Nonlinear reflectivity to data mapping • Unknown medium modeled by reflectivity q . First order partial differential operator L q is affine in q . Kinematics is assumed known! • Inverse scattering problem: Find q in Ω c from array mea- surements of reflected primary wave � � � � � � � D ( r,s ) b ( r ) , P ( s ) ( kτ, · ) b ( r ) , cos L q L T b ( s ) := = kτ q k for r, s = 1 , . . . , m and time instants kτ, k = 0 , . . . , 2 n − 1 . � � • L q is affine in q but map q �→ 0 ≤ k ≤ 2 n − 1 to invert is nonlinear D k 3
Outline of talk • Most imaging assumes reflectivity to data map is linear (Born approximation). • Nonlinear methods: Qualitative (linear sampling, factoriza- tion,...) mostly at single frequency. Optimization is difficult. Goal 1: Use reduced order model (ROM) to approximate the Data to Born (DtB) map ( D k ) 0 ≤ k ≤ 2 n − 1 → ( D Born ) 0 ≤ k ≤ 2 n − 1 k D Born defined using Fr´ echet derivative of map q �→ D k at q = 0. k Goal 2: Use the ROM to obtain quantitative estimate of q . 4
ROM for Data to Born (DtB) transformation Data are m × m matrices � � � � � � b (1) , . . . b ( m ) � L q L T D k = b , cos kτ b , b = , 0 ≤ k ≤ 2 n − 1 q � � � L q L T ROM of propagator ∗ P q = cos τ q • Wave at time kτ is Chebyshev polynomial T k of 1 st kind of P q � � � P ( s ) ( kτ, x ) = cos L q L T b ( s ) ( x ) = T k ( P q ) b ( s ) ( x ) kτ q • ROM propagator P P gives exact data D k , 0 ≤ k ≤ 2 n − 1. P ROM q It is constructed from the data and inherits properties of P q to allow DtB transformation. 5 ∗ P ( s ) (( k + 1) τ, x ) = 2 P q P ( s ) ( kτ, x ) − P ( s ) (( k − 1) τ, x )
Definition of reduced order model (ROM) Algebraic setting: from continuum to fine grid discretization • Operator L q � lower block bidiagonal matrix L q ∈ R N × N � � � L q L T P • Propagator P P q = cos τ is N × N matrix. q � � P ( s ) ( kτ, · ) P q ) b ∈ R N × m • Snapshots P P P k = 1 ≤ s ≤ m = T k ( P P P P P ROM obtained by projection on range of P P := ( P P 0 , . . . ,P P n − 1 ) ROM = V T b ∈ R nm × n = V T P P q V ∈ R nm × nm P P P P ROM b q Here V ∈ R N × nm satisfies V V T = orthogonal projector on range( P V T V = I nm , P P ) 6
Data interpolation Theorem: Projection ROM satisfies D k = b T T k ( P ROM ) T T k ( P P P P q ) b = ( b P q ) b ROM , 0 ≤ k ≤ 2 n − 1 . ROM Proof: ROM for k = 0 , . . . , n − 1. Step 1: Prove P P P k = V T k ( P P q ) b P ROM Step 2: This gives ROM = ( b D k = b T P P k = b T V T k ( P ROM ) T T k ( P P P P q ) b P P q ) b ROM , 0 ≤ k ≤ n − 1 ROM ROM Step 3: For n ≤ k ≤ 2 n − 1 use the above and the recursion T k ( x ) = 2 T n − 1 ( x ) T k − n +1 ( x ) − T | 2 n − 2 − k | ( x ) � 7
Proof that P P P k = V T k ( P P ) b for k = 0 , . . . , n − 1 P ROM ROM q ROM = V T b = V T P Using definition P P P P P q V and b ROM q ROM = V T 0 ( P P 0 = b = V V T b = V b • For k = 0 we have P P P q ) b P ROM ROM • Hypothesis: true for k < n − 1. • For k + 1 use T k +1 ( x ) = 2 x T k ( x ) − T k − 1 ( x ) ROM = 2 V P ROM − V T k − 1 ( P V T k +1 ( P P P q ) b P P q T k ( P P P q ) b P P q ) b ROM ROM ROM ROM ROM ROM − P = 2 V V T P P P P P q V T k ( P P q ) b ROM P k − 1 = V V T 2 P P P P P q P P k − P P k − 1 = V V T ( P P P k +1 + P P P k − 1 ) − P P P k − 1 = P P P k +1 8
Our choice of V Note: Any V satisfies the data interpolation. Which V is best? • Define V by Gram-Schmidt (QR factorization) P P P = V R P • Causality and finite speed of propagation make P P ≈ block upper-tridiagonal with coordination of temporal and spatial mesh. This requires knowing kinematics! Basis that transforms P P P to block upper-tridiagonal R is almost the canonical one � V is approximate identity. = V T P Theorem: Matrix V from QR factorization makes P P P P P q V ROM q block tridiagonal . This result proved using recursion relations of polynomials be- comes important in inversion. 9
Illustration for sound waves in 1-D 10
Illustration for sound waves in 2-D σ c P P P o P P P q V o V Array with m = 50 sensors × Snapshots plotted for a single source ◦ 11
From data to ROM • Start with P = ( P P P P n − 1 ) = V R and use P P j = T j ( P P P q ) b P P 0 , . . . ,P 2 b T � � P q ) b = 1 ( P T P ) jk = b T T j ( P P P P P P q ) T k ( P T j + k ( P P q ) + T | j − k | ( P P q ) b � � = 1 = ( R T R ) jk , D j + k + D | j − k | 0 ≤ j, k ≤ n − 1 . 2 Block Cholesky decomposition to get R is ill-conditioned part of P T P computation � spectral truncation of Gramian P P P P . P q V = R − T � � = V T P P T P R − 1 P P P P P • ROM propagator: P P P P q P P ROM q � � � � j,k = 1 P T P P P P P P q P P D j + k +1 + D | k − j +1 | + D | k − j − 1 | + D | k + j − 1 | 4 ROM = V T b = V T P P 0 = V T V R E 1 = R E 1 • ROM sensor function: b P 12
Properties of ROM propagator • Propagator factorization � � �� � � � 2 = 2 L q L T = L q L T P I − cos I − P P q τ q , L q ≈ L q q τ 2 τ 2 L q = block lower bidiagonal (discretized 1 st order operator L q ). • By construction ROM propagator is symmetric, block tridiag- onal with factorization � � 2 q ) T = V T L q L T P = L q ( L I nm − P P ROM ROM ROM q V . q τ 2 = V T L q � • Cholesky factor L V is block lower bidiagonal. ROM q � Galerkin approximation on spaces of primary and dual snap- shots with orthogonal bases in V and � V . 13
Data to Born transformation • Approximate Fr´ echet derivative of ROM ) T T k ( P P q �→ D k = ( b P q ) b ROM ROM using = V T L q � - L V is approximately affine in q . ROM q ROM = V T b is independent of q . - b • For a scaled down reflectivity ǫq , with ε ≪ 1, εq := I mn − τ 2 � � εq ) T L εq := L 0 + ε L − L , P P P 2 L εq ( L ROM ROM ROM ROM ROM ROM ROM q 0 • The transformed (to Born) data: � ROM ) T d � D Born � := D 0 ,k + ( b dε T k ( P P P εq ) ROM , 0 ≤ k ≤ 2 n − 1 ROM b � ε =0 k 14
DtB transformation: Sound waves 1-D 15
DtB transformation: Sound waves 2-D σ ( x ) c ( x ) Scattered field True Born DtB transform Axes in km. Colorbars show σ , c normalized by values at array. 16
Robustness of transformation to background velocity True wave speed Incorrect wave speed Wrong velocity model induces artifacts due to domain boundary 17
DtB transformation: Sound waves - 2D Model c ( x ) Scattered data DtB • Here we considered constant ρ and variable velocity. Only the constant background c o is assumed known. • Note how the echo from small reflector, masked by a multiple, is revealed by the DtB transformation. Results for 2-D isotropic elasticity are in our paper. 18
Quantitative inversion: 2 possibilities • Use DtB output in linear least-squares Born data fit: 2 n − 1 � � D Born − F Born ( q s ) � 2 q = arg min q s F k =0 ≈ affine in q ( x ) ≈ � • Match ROM instead. Since L j q j φ j ( x ) ROM q � � � q � 2 q s q = arg min q s �L q s − L F , L q s = L 0 + L φ j − L ROM ROM ROM ROM ROM ROM j 0 j 0 ) T (tridiagonal matrix discretization of ∆) Grid from L 0 ( L ROM ROM 19
Quantitative inversion Linear LS data fit without the DtB transformation. 20
Quantitative inversion Linear LS data fit with the DtB transformation. 21
Quantitative inversion: ROM match Iteration 1 and 6 (top and middle) and true medium bottom. 22
Quantitative inversion: ROM match (iteration 3) 23
Conclusions • We introduced a linear algebraic algorithm for transforming the scattered wave measured by an active array of sensors to the single scattering (Born) approximation which is linear in the unknown reflectivity. • We showed that ROM can be used for quantitative inversion. Lots left to do: • Synthetic aperture setup; transmission setup; time harmonic waves , anisotropic and attenuating media . • Approach can be extended to select multiple scattering effects. 24
Recommend
More recommend
