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. Ldcϕ = sdc, ddc = Mdcϕ, Forward model: finite volume method to solve for ddc.1 Inverse model: adjoint method to invert for σ.2

1Dey and Morrison, 1979 2Pratt et al., 1998, Pidlisecky et al., 2007

Simulated ER Model and Data

Inverse Electrical Resistivity Model

minσ

• do

dc − Mdcϕ(σ)2 2 + Ldcϕ(σ) − sdc2 2

• Initial Estimates - Regularization

RLσ2

2

R = diag(r1, . . . , rn), ri = 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
• f where parameter values change, rather than parameter values.

Joint Inversion - additional data with physics

GPR Model

ε¨ u + σ ˙ u = 1 µ∇2u + sw u-electric field Ey, ε-permittvity, µ-constant permeability, sw-source wavelet. u = Lwsw dw = Mwu Forward model: finite difference time domain on a Yee grid with PML.3 Inverse model: full waveform inversion to solve for ε and σ.4

3Yee, 1966; Berenger, 1994 4Ernst et al., 2007

Simulated GPR Model and Data

GPR inversion

minσ,ǫ

• do

w − Mwu2 2 + u − Lwsw2 2

• do

w - shot gather for all time and all receivers, function of (σ, ǫ, sw).

• Steepest descent: compute gradient using adjoints (gǫ, gσ); step

size (αǫ, ασ) by line search method.

• Update for each source (∆ǫs, ∆σs) and sum ( ∆ǫs/ns, ∆σs/ns)

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

• High frequency
• Conductivity through

attenuation and reflection

ER

• Low frequency
• Directly sensitive to conductivity

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

minm

• G1m − d12

2 + ||G2m − d2||2 2

• = minm
• G1

G2

• [m] −
• d1

d2

• 2

2

≡ minmG12m − d122

2

Is G12 better conditioned than G1 or G2?

Green’s function solutions

LAu(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 → HA u(t) = Ah(t) ≡

• Ω

KA(t, s)h(s)ds

Continuous Least Squares - Singular Value Expansion

Ah − u2

HA

has solution h = A†u =

∞

• k=1

ψk, uHA σk φk 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

minh Tλh − (u, 0)2

HA×H = minh Ah − u2 HA + λ h2 H

i.e Tλh = (Ah, √ λh) with Tλ : H → HA × H.5 Optimal solution h = A†

λu = ∞

• k=1

σk σ2

k + λψk, uHAφk

so that σk σ2

k + λ → 0, as k → ∞

and λ restricts the solution space.

5Gockenbach, 2015

Continuous Joint Inversion

Ah − u12

HA + Bh − u22 HB

with A : H → HA and B : H → HB. Joint operator: C : H → HA × HB = {(hA, hB) : hA ∈ HA, hB ∈ HB} 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 C(n) =

• A(n)

B(n)

• Approximate A with A(n), a(n)

ij

= qi, Apj, where {qi(s)}n

i=1 and {pj(t)}n j=1

are orthonormal bases. Then A(n) = U(n)Σ(n) V (n)T , Σ(n) = diag

• σ(n)

1 , σ(n) 2 , . . . σ(n) n

Special case: Self-Adjoint Operator

Use singular functions in Galerkin method a(n)

ij

= φj, Aφi = φj, σiφi =

• σi

i = j i = j then A(n) = Σ(n)

A

and B(n) = Σ(n)

B

so that

• C(n)T

C(n) =

• Σ(n)

A

2 +

• Σ(n)

B

2

and σi

• C(n)

=

• σi

A(n)2 + σi B(n)2

Joint Singular Value Example

LAu = −u′′, u(0) = u(π) = 0, LBu = u′′ + b2u, u(0) = u(π) = 0, and b / ∈ Z Green’s functions: KA =

  

1 b (π − x) y,

0 ≤ y ≤ x ≤ π,

1 b (π − y) x,

0 ≤ x ≤ y ≤ π. KB =

  

−sin(by) sin[b(π−x)]

b sin(bπ)

, 0 ≤ y ≤ x ≤ π, −sin(bx) sin[b(π−y)]

b sin(bπ)

, 0 ≤ x ≤ y ≤ π.

Individual Singular Values - b = π

σk(A(200)) =

1 k2

σk(B(200)) =

1 k2+π2

Joint Singular Values - b = π

σk(C(200)) =

• 1

k2

2 +

• 1

k2+π2

2

Joint Singular Values

b = 10.1 b = 0.1

σk(C(200)) = 1

k2

2 +

• 1

k2+10.12

2 σk(C(200)) = 1

k2

2 +

• 1

k2+0.12

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