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
Recommend
More recommend
