Reduced Order Model Approach to Inverse Scattering orn Zimmerling J - PowerPoint PPT Presentation

Reduced Order Model Approach to Inverse Scattering orn Zimmerling † J¨ Liliana Borcea † , Vladimir Druskin ∗ , Mikhail Zaslavsky & ,Alexander Mamonov ⊗ † University of Michigan, ∗ WPI, & Schlumberger, ⊗ UH 22 April 2020, ICERM & J¨ orn’s Couch

Introduction Consider inverse scattering for a hyperbolic equation ac ( ∂ 2 t + L q L T q ) u ( t , x ) = 0 u (0 , x ) = b ( x ) ∂ t u (0 , x ) = 0 with first order waveop. L q affine in reflectivity q as q u = √ c ∇ [ √ c u ] + c L T 2 ∇ q u Ω inac UMICH (J¨ orn Zimmerling) 2 / 32

Introduction Consider inverse scattering for a hyperbolic equation ( ∂ 2 t + L q L T q ) u ( t , x ) = 0 ac u (0 , x ) = b ( x ) ∂ t u (0 , x ) = 0 with first order waveop. L q affine in reflectivity q as q u = √ c ∇ [ √ c u ] + c L T 2 ∇ q u Ω inac Coinciding sources and receivers yield data at times t = j τ for j = 0 , 1 , . . . , 2 n − 1. � � � � � L q L T D j = � b ( x ) , u ( j τ, x ) � = b ( x ) , cos j τ b ( x ) q In a MIMO setting b ( x ) = ( b (1) ( x ) , ..., b ( m ) ( x )) with m sources/receivers Problem 1. Find reflectivity q ( x ) from backscattering data D , i.e. invert q �→ D UMICH (J¨ orn Zimmerling) 2 / 32

Outline of approach 1. Construct a ROM L ROM from the measured data that has the q structure of a wave operator ◮ Data interpolation � � � � � D = b ( x ) , cos L q L T b ( x ) j τ q � � � � � T b ROM , cos b ROM = j τ L ROM L ROM q q ◮ Sparsity pattern of L ROM should resemble discretization q ◮ Regularization on level of the ROM 2. Interpret ROM L ROM as a wave operator (affine in q ) q 3. Invert for the reflectivity by minimizing ROM mismatch O LS ( q s ) = ||L ROM − L ROM || + regularization q s q UMICH (J¨ orn Zimmerling) 3 / 32

Comparison to standard imaging (FWI) ROM Approach FWI 1. Objective function data 1. Objective function ROM mismatch missmatch 2. Highly nonlinear (many 2. Close to linear due to ROM iterations) transform 3. Difficult to regularize 3. Intrinsic regularization UMICH (J¨ orn Zimmerling) 4 / 32

The Propagator � � � ◮ We define the Propagator P q = cos L q L T τ q ◮ Field can be expressed in Chebyshev polynomials of first kind of P q � � � � � � � ��� L q L T L q L T u ( j τ, x ) = cos j τ b = cos j arccos τ b cos q q = cos ( j arccos [ P q ]) b = T j ( P q ) b ◮ We move from continuous time ∂ 2 t u ( t , x ) = − L q L T q u ( t , x ) u (0 , x ) = b ( x ) ∂ t u (0 , x ) = 0 (1) to a discrete time equation u j = u ( j τ, x ) (Chebyshev recursion) u j +1 = 2 P q u j − u j − 1 u 0 = b ( x ) u 1 = u − 1 (2) ◮ Note the similarity of (2) to discretizing ∂ 2 t in (1) u j +1 − 2 u j + u j − 1 = 2 τ 2 ( P q − I ) u j = −L q L T q u j (3) τ 2 UMICH (J¨ orn Zimmerling) 5 / 32

Galerkin projection - ROM ◮ Project (2) onto the solutions at the first n times U ( x ) = ( u 0 ( x ) , u 1 ( x ) , . . . , u n − 1 ( x )) ◮ Expand u j ≈ U ( x )g j and use the Galerkin condition to find � U , U � g j +1 = 2 � U , P q U � g j − � U , U � g j − 1 D ROM = g T g 0 = e 1 g 1 = g − 1 0 � U , U � g j j ◮ This defines our ROM with Mg j +1 = 2Sg j − Mg j − 1 Propositions 1. The data obtained from the ROM interpolates the true data D ROM = D j j ≤ 2 n − 1 j 2. We can find M and S from the data D j UMICH (J¨ orn Zimmerling) 6 / 32

Reduced order model propagator operator I Mg j +1 = 2Sg j − Mg j − 1 ◮ Cholseky factorize M = R T R to orthogonalize the basis V = U R − 1 ◮ Introduce ROM field u ROM U R − 1 , u j � � = Rg j = j = 2R − T SR − 1 u ROM u ROM − u ROM u ROM = Re 1 = b ROM j +1 j − 1 0 j = R − T SR − 1 ◮ This defines the ROM propagator operator P ROM q ◮ P ROM is a projection of the propagator P q onto the orthogonalized q snapshots V ( x ) = U ( x )R − 1 UMICH (J¨ orn Zimmerling) 7 / 32

Reduced order model propagator operator I Mg j +1 = 2Sg j − Mg j − 1 ◮ Cholseky factorize M = R T R to orthogonalize the basis V = U R − 1 ◮ Introduce ROM field u ROM U R − 1 , u j � � = Rg j = j = 2R − T SR − 1 u ROM u ROM − u ROM u ROM = Re 1 = b ROM j +1 j − 1 0 j = R − T SR − 1 ◮ This defines the ROM propagator operator P ROM q ◮ P ROM is a projection of the propagator P q onto the orthogonalized q snapshots V ( x ) = U ( x )R − 1 ◮ Elements of U are strongly dependent on q ◮ Elements of V are localized and weakly dependent on q ◮ ⇒ Known background ( q = 0) snapshots V 0 ( x ) ≈ V ( x ) provide embedding of ROM into physical space UMICH (J¨ orn Zimmerling) 7 / 32

Reduced order model propagator operator II u ROM = 2 P ROM u ROM − u ROM u ROM = Re 1 = b ROM j +1 q j j − 1 0 ◮ Defines data D ROM = (b ROM ) T T j ( P ROM )b ROM j q ◮ Interpretation as discrete-time, discrete-space WEQ u ROM − 2u ROM − u ROM = − 2 j +1 j − 1 j τ 2 (I − P ROM )u ROM = −L ROM ( L ROM ) T u ROM q j q q j τ 2 UMICH (J¨ orn Zimmerling) 8 / 32

Reduced order model propagator operator II u ROM = 2 P ROM u ROM − u ROM u ROM = Re 1 = b ROM j +1 q j j − 1 0 ◮ Defines data D ROM = (b ROM ) T T j ( P ROM )b ROM j q ◮ Interpretation as discrete-time, discrete-space WEQ u ROM − 2u ROM − u ROM = − 2 j +1 j − 1 j τ 2 (I − P ROM )u ROM = −L ROM ( L ROM ) T u ROM q j q q j τ 2 ◮ Compare to to continuous-time, continuous-space WEQ ∂ 2 t u ( t , x ) = − L q L T q u ( t , x ) ◮ L ROM can be interpreted as a discretization of L q absorbing O ( τ 2 ) q error from time discretion UMICH (J¨ orn Zimmerling) 8 / 32

Interpolation Proof, D ROM = D j j ≤ 2 n − 1 j j ≤ n − 1 ◮ For j ≤ n − 1 the true field is in the span of U ( x ) ◮ The Galerkin condition is unique and g j = e j +1 (block identity matrix) j ≤ 2 n − 2 ◮ The field is a (Chebyshev) polynomial in P q and P q = P T q u j = T j ( P q ) b ◮ For l = 1 , . . . , n − 1 we use the Chebyshev identity T n − 1+ l = 2 T n − 1 T l − T | n − 1 − l | D n − 1+ l = � b , T n − 1+ l ( P q ) b � � � = 2 �T n − 1 ( P q ) b , T l ( P q ) b � − b , T | n − 1 − l | ( P q ) b � � = 2 � u n − 1 , u l � − b , u | n − 1 − l | = 2(u ROM n − 1 ) T u ROM − b ROM u ROM l | n − 1 − l | � � b , T n − 1+ l ( P ROM = D ROM = ) b n − 1+ l q ◮ Similar prove for j = 2 n − 1 UMICH (J¨ orn Zimmerling) 9 / 32

Finding M and S from the Data ◮ The Mass matrix is defined as M = � U , U � (M) i +1 , j +1 = � u i , u j � = �T i ( P q ) b , T j ( P q ) b � = 1 2 � b , 2 T i ( P q ) T j ( P q ) b � = 1 � � b , [ T i + j ( P q ) + T | i − j | ( P q )] b 2 = 1 2( D i + j + D | i − j | ) ◮ Similar result can be obtained for S ◮ With noise M and S need to be regularized such that 1. M is positive definite 2. The pencil ( M , S ) has a maximum eigenvalue ≤ 1 (Guarantees that P ROM is contractive operator) ◮ Basically L¨ owner inner products in the time domain ◮ Next: V ( x ) = U ( x )R − 1 independence of q UMICH (J¨ orn Zimmerling) 10 / 32

Orthogonalized Snapshot matrix - 1D ◮ Orthogonalized Snapshots approximately independent of medium UMICH (J¨ orn Zimmerling) 11 / 32

Orthogonalized Snapshot matrix - 2D q c U | q =0 V | q =0 U V Array with m = 50 sensors × Snapshots plotted for a single source ◦ UMICH (J¨ orn Zimmerling) 12 / 32

Summary ROM ◮ From the data we can obtain a ROM that explains the data u ROM − 2u ROM − u ROM j +1 j j − 1 = −L ROM ( L ROM ) T u ROM q q j τ 2 ◮ Size of the ROM is dictated by the data ◮ The ROM lives in a basis of orthogonalized snapshots U R − 1 (localization) ) T block tridiagonal ◮ L ROM ( L ROM q q ◮ L ROM is nearly affine in the reflectivity q q ◮ L ROM − L ROM close to linear in q [1] q 0 √ 2 τ 2 LL T + O ( τ 4 ≈ − 1 � � LL T I − P q = I − cos τ 4! ) [1] Liliana Borcea, Vladimir Druskin, Alexander Mamonov and Mikhail Zaslavsky , Untangling the nonlinearity in inverse scattering with data-driven reduced order models , Inverse Problems, Volume 34, Number 6 UMICH (J¨ orn Zimmerling) 13 / 32

Inversion Method ◮ Minimize ROM mismatch rather than data mismatch O LS ( q s ) = ||L ROM − L ROM || F + regularization q q s ◮ We assume a kinematic model c 0 and an initial guess q 0 = 0 to obtain L ROM 0 ◮ Parametrize the reflectivity as N s q s ( x ) = � q s j φ j ( x ) j =1 ◮ Use (approximate) affine relationship (approximate as V 0 ≈ V ) N s � L ROM ≈ L ROM q s j ( L ROM − L ROM + ) q s 0 φ j 0 j =1 ◮ The coefficients q j follow from a least-squares problem ◮ Approach can be iterated ◮ How to choose φ j ( x )? ⇒ Based on Resolution UMICH (J¨ orn Zimmerling) 14 / 32

Choosing φ j ( x ) ◮ The resolution of imaging method depends on location in Ω ◮ A point like perturbation δ j ( x ) at the point x j causes a local perturbation of the diff.op. ∆ L ( δ j ) such that ∆ L ( δ j )∆ L ( δ j ) T ϕ ( x ) = c 2 ( x ) |∇ δ j ( x ) | 2 ϕ ( x ) 4 ◮ Idea: Lift the perturbation ∆ L ROM that δ j causes in the ROM into δ j the physical space using the orthogonalized basis function of the background V 0 ( x ) � � V 0 ( x )∆ L ROM (∆ L ROM ) T V T Ψ j ( x ) = 0 ( x ) δ j δ j ◮ We do a partition of unity using these point spread function ◮ If we need to regularize our ROM due to noisy data, it directly impacts this resolution UMICH (J¨ orn Zimmerling) 15 / 32