On the Effectiveness of Joint Inversion Jodi Mead James Ford - PowerPoint PPT Presentation

On the Effectiveness of Joint Inversion Jodi Mead James Ford Department of Mathematics Clearwater Analytics Diego Domenzain John Bradford Department of Geosciences Department of Geophysics Colorado School of Mines Thanks to: NSF DMS-1418714

Outline • Joint inversion of Electrical Resistivity and Ground Penetrating Radar for near subsurface imaging – Additional data as prior information or regularization – Simultaneous joint inversion • Method to assess data effectiveness in joint inversion

Single Inversion

Joint Inversion - additional data for initial estimates

Near subsurface imaging Boise Hydrogeophysical Research Site (BHRS) • Field laboratory on a gravel bar adjacent to the Boise River, 15 km southeast of downtown Boise. • Consists of coarse cobble and sand. Braided stream fluvial deposits overlie a clay layer at about 20 m depth. Difference in retention properties in a lenticular sand feature yields signifi- cantly different geophysical properties.

Electrical Resistivity Tomography (ERT) • 2D grid of observations pro- vides a 2.5-D inverted model that emphasizes the sand lenticular feature. • BHRS survey consisted of 12 electrodes spaced 1 meter apart acquired with a dipole- dipole configuration. BHRS survey acquired at surface when subsurface achieved saturation.

Electrical Resistivity Model −∇ · σ ∇ ϕ = i ( δ ( x − s + ) − δ ( x − s − )) ϕ - electric potential, σ - conductiv- ity, i - current intensity, s ± - source- sink position. L dc ϕ = s dc , d dc = M dc ϕ, Forward model: finite volume method to solve for d dc . 1 Inverse model: adjoint method to invert for σ . 2 1 Dey and Morrison, 1979 2 Pratt et al., 1998, Pidlisecky et al., 2007

Simulated ER Model and Data

Inverse Electrical Resistivity Model � � dc − M dc ϕ ( σ ) � 2 2 + � L dc ϕ ( σ ) − s dc � 2 � d o min σ 2 Initial Estimates - Regularization � RLσ � 2 2 R = diag ( r 1 , . . . , r n ) , r i = 0 or 1 . L - 1st or 2nd derivative operator

Prior information - Ground Penetrating Radar (GPR) • GPR survey at BHRS acquired across center of gridded ER survey. • GPR sampled line collinear with ER survey. • GPR derived boundary gives prior boundary knowledge in the ER dataset.

Boise Hydrogeophysical Research Site Results • ER data inverted • Regularization in the form of subsurface boundary information in- ferred from GPR data

Summary - Joint inversion as regularization • New data can be used to form a regularization operator – This relies on secondary data processing or practitioner inter- pretation of data. – Will always produce a well-posed problem. • Additional data as derivative information - only requires knowledge of where parameter values change, rather than parameter values.

Joint Inversion - additional data with physics

GPR Model u = 1 µ ∇ 2 u + s w ε ¨ u + σ ˙ u -electric field E y , ε -permittvity, µ -constant permeability, s w -source wavelet. u = L w s w d w = M w u Forward model: finite difference time domain on a Yee grid with PML. 3 Inverse model: full waveform inversion to solve for ε and σ . 4 3 Yee, 1966; Berenger, 1994 4 Ernst et al., 2007

Simulated GPR Model and Data

GPR inversion � � � d o w − M w u � 2 2 + � u − L w s w � 2 min σ,ǫ 2 • d o w - shot gather for all time and all receivers, function of ( σ, ǫ, s w ). • Steepest descent: compute gradient using adjoints ( g ǫ , g σ ); step size ( α ǫ , α σ ) by line search method. • Update for each source ( ∆ ǫ s , ∆ σ s ) and sum ( � ∆ ǫ s /n s , � ∆ σ s /n s ) for parameter update. • Invert for source wavelet with Wiener filter. • Transform 2d wavefield into 2.5d data by filtering wavefield in the frequency domain.

Complementary data GPR ER • High frequency • Low frequency • Conductivity through • Directly sensitive to conductivity attenuation and reflection

Inverting ER and GPR jointly - full physics

Inverted images - full physics

Inverted cross section - full physics

Combining updates

Data weights

Summary - Jointly inverting with full physics • We have developed a joint inversion algorithm to solve for both permittivity ǫ and conductivity σ using complementary GPR and ER data. • Features were recovered that neither GPR or ER can individually resolve. • Data weights must reflect the physics that can be captured by each particular data set during iterations of the inversion.

Assessing Effectiveness of Joint Inversion

Discrete Joint Inversion � � � G 1 m − d 1 � 2 2 + || G 2 m − d 2 || 2 min m 2 2 � �� � � � G 1 d 1 � � = min m [ m ] − � � � � G 2 d 2 � � 2 ≡ min m � G 12 m − d 12 � 2 2 Is G 12 better conditioned than G 1 or G 2 ?

Green’s function solutions L A u ( t ) = h ( t ) Forward problem: Given h , find u Inverse problem: Given u , find h Conditioning of the inverse problem depends on the forward operator A : H → H A � u ( t ) = Ah ( t ) ≡ K A ( t, s ) h ( s ) ds Ω

Continuous Least Squares - Singular Value Expansion � Ah − u � 2 H A has solution ∞ � ψ k , u � H A h = A † u = � φ k σ k k =1 Condition number: σ k → 0 , as k → ∞ Decay rate q : σ k ( A ) decays like k − q Decay rate of singular values allow us to classify model conditioning

Regularization with Tikhonov Operator min h � T λ h − ( u, 0) � 2 H A × H = min h � Ah − u � 2 H A + λ � h � 2 H √ λh ) with T λ : H → H A × H . 5 Optimal solution i.e T λ h = ( Ah, ∞ σ k h = A † � λ u = k + λ � ψ k , u � H A φ k σ 2 k =1 so that σ k k + λ → 0 , as k → ∞ σ 2 and λ restricts the solution space. 5 Gockenbach, 2015

Continuous Joint Inversion � Ah − u 1 � 2 H A + � Bh − u 2 � 2 H B with A : H → H A and B : H → H B . Joint operator: C : H → H A × H B = { ( h A , h B ) : h A ∈ H A , h B ∈ H B } so that Ch = ( Ah, Bh ) , h ∈ H Continuous analog to stacking matrices

Singular Values of Joint Operator σ 2 φ = C ∗ Cφ = A ∗ Aφ + B ∗ Bφ Galerkin method � A ( n ) � C ( n ) = B ( n ) Approximate A with A ( n ) , a ( n ) = � q i , Ap j � , where { q i ( s ) } n i =1 and { p j ( t ) } n j =1 ij are orthonormal bases. Then V ( n ) � T , A ( n ) = U ( n ) Σ ( n ) � Σ ( n ) = diag � σ ( n ) 1 , σ ( n ) � 2 , . . . σ ( n ) n

Special case: Self-Adjoint Operator Use singular functions in Galerkin method � σ i i = j a ( n ) = � φ j , Aφ i � = � φ j , σ i φ i � = ij 0 i � = j then A ( n ) = Σ ( n ) and B ( n ) = Σ ( n ) A B so that � 2 + � 2 C ( n ) � T � � C ( n ) � � Σ ( n ) � Σ ( n ) = A B and � � A ( n ) � 2 + σ i � C ( n ) � � B ( n ) � 2 σ i = σ i

Joint Singular Value Example L A u = − u ′′ , u (0) = u ( π ) = 0 , L B u = u ′′ + b 2 u, u (0) = u ( π ) = 0 , and b / ∈ Z Green’s functions:  1 b ( π − x ) y, 0 ≤ y ≤ x ≤ π,  K A = 1 b ( π − y ) x, 0 ≤ x ≤ y ≤ π.   − sin( by ) sin[ b ( π − x )] , 0 ≤ y ≤ x ≤ π,  b sin( bπ ) K B = − sin( bx ) sin[ b ( π − y )] , 0 ≤ x ≤ y ≤ π.  b sin( bπ )

Individual Singular Values - b = π σ k ( A (200) ) = 1 σ k ( B (200) ) = 1 k 2 k 2 + π 2

Joint Singular Values - b = π �� � 2 + � 2 � σ k ( C (200) ) = 1 1 k 2 k 2 + π 2

Joint Singular Values b = 10 . 1 b = 0 . 1 �� 1 �� 1 � 2 + � 2 � 2 + � 2 � � σ k ( C (200) ) = 1 σ k ( C (200) ) = 1 k 2 k 2 +10 . 1 2 k 2 k 2 +0 . 1 2

Summary - Assessing joint inversion with singular values of joint operators • Adding data may not completely eliminate ill-posedness in individual inversions. • Even in situations where regularization is necessary, joint inversion will reduce the amount of regularization. • Theoretical models and analysis can identify properties and or situa- tions where different data types effectively regularize each other.

Thank you! Jodi Mead Mathematics Department Boise State University jmead@boisestate.edu