A nonlinear differential-semblance algorithm for waveform inversion Dong Sun
The Rice Inversion Project Annual Review Meeting March 30, 2012
,
1/15
Agenda
,
1
Full Waveform Inversion
2
Extension and Differential Semblance Optimization (DSO)
3
Examples: Inversion in Layered Medium
2/15
Full Waveform Inversion (FWI) A usual set-up: M = a set of models ({m(x)}), D = a space of data F : M → D - forward map FWI (least-squares inversion): given do ∈ D, solve min JLS :=
m∈M
1 2
2
kF[m] − do k
[+ regularizing term(s)]
Note: FWI - globally minimize the objective using local methods Q observed data do = s do (xr , t; s) - highly redundant & band-limited – possibilities for acquisition parameter s include source position, offset, slowness (plane wave data),...
,
3/15
Full Waveform Inversion (FWI)
Studied since the 80’s, becoming feasible in last 10 yrs: + accommodates any modeling physics & provides quantitative inferences of subsurface (some spectacular successes) ± applicable to long-offset surface data (diving and refracted waves); most of successful inversions so far involve transmitted data − inversion with reflection data still a challenge fundamental obstacles: . spectral data incompleteness - missing low frequencies in data leads to missing long scale components in estimated model . strong nonlinearity, many false local minima - descent methods fail
,
3/15
Full Waveform Inversion (FWI) Remedies for reflection inversion: automatic migration velocity analysis (MVA) (e.g., DSO variants): decompose model into slowly & fast varying parts (background & reflectivity), then alternately a. b.
build reflectivity via migration or linearized LS-inversion update background to reduce inconsistency of reflectivity gathers
+ infer global changes in model (long-scale updates) − limited by linearization and scale separation assumptions
X nonlinear DSO: import MVA’s concepts (image gathering, semblance measuring) into FWI and drop MVA limitations
,
3/15
Agenda
,
1
Full Waveform Inversion
2
Extension and Differential Semblance Optimization (DSO)
3
Examples: Inversion in Layered Medium
4/15
Formulate FWI via DSO Concept: relax FWI via extension then demand coherence of extended model
FWI - attempt to match all data simultaneously with one model min 1 m∈M 2
kF[m] − dk2 :=
1 2
X
kF[m](s) − d(s)k2
s
Extended WI - fit subsets of data with non-physical extended models min m∈M ¯
1 2
2
F[m] ¯ − d :=
Note: extended models M = {m(x, ¯ s)} m(x, ¯ s) = m(x) for all s)
1 2
X
kF[m(., ¯ s)](s) − d(s)k2
s
(M ⊂ M, regarding m ¯ = m iff
extended modeling F : M → D by F[m](s) ¯ = F[m(·, ¯ s)](s)
Equivalent problem to FWI - find coherent solution to extended WI
,
5/15
Formulate FWI via DSO
Differential Semblance (with surface oriented extension): s finely sampled ⇒ coherence criterion is ∂ m/∂s ¯ =0 Differential Semblance Optimization:
2 1 ∂
min JDS [m] ¯ := m ¯ (coherence measurement) 2 ∂s m∈M ¯
F[m] ¯ − d ≈ 0 (data-fitting constraints) s.t. Key to success: need a proper control in order to navigate through the feasible model set
¯ − d ≈ 0 m ¯ ∈ M : F[m]
,
5/15
Using LF data as control
Concept: Cannot use independent long-scale model as control, as in MVA: “low spatial frequency” not well defined, depends on velocity. However, temporal passband is well-defined, and observed data lacks very low frequency energy (0-3, 0-5,... Hz) with good s/n Generally, the impulsive inverse problem is solvable: LS inversion leads to “unique” model, if data d is not band-limited (good s/n down to 0 Hz) So: use low-frequency data as control (analogous to long-scale model in MVA)
,
6/15
Using LF data as control Approach I: use low-frequency data as control - define low-frequency source complementary to missing passband, low-frequency modeling op Fl and its extension Fl - define complementary low-frequency data dl - define m[d ¯ l ] to be minimizer of
¯ 2 2 2 ∂m
kF[m] ¯ + Fl [m] ¯ − (d + dl )]k + σ ∂s - nDSO: determine dl by solving
2
∂
¯ l ] min m[d
dl ∂s
Initial exploration: D. Sun (2008), D. Sun & W. Symes (2009), for plane wave / layered medium modeling – slices of DS objective are convex, i.e., an enlargement of the domain of attraction of the global minimum is achieved Further thought: find a way to supply meaningful/consistent low frequency data ,
6/15
Using LF data as control Approach II: generate low-frequency data from model - Given low frequency control model ml ∈ M, define extended model m ¯ = m[m ¯ l ] by minimizing
¯ 2 2 2 ∂m
kF[m] ¯ + Fl [m] ¯ − (d + Fl [ml ])]k + σ ∂s - nDSO: determine ml by solving
2
∂
min m[m ¯ ] l
ml ∂s Advantage of this approach: ml plays same role as migration velocity model, but no linearization, scale separation assumptions required by formulation. Current project: 2D constant-density acoustics, plane-wave sources (to minimize edge artifacts)
,
6/15
DSO Algorithm Key step: compute ml gradient: ∇JDS [ml ] = −DFl [ml ]T DF l [m[m ¯ l ]]H[m[m ¯ l ]]−1 where
∂2 m[m ¯ l] ∂s2
T ¯ + DFl [m] ¯ DF [m] ¯ + DFl [m] ¯ H[m] ¯ = DF [m]
Computational procedure: 1 2 3 4 5
,
initialize ml , etc. solve sub-LS problem for m[m ¯ l] evaluate JDS [m]. ¯ If stopping criterion satisfied, stop; else, continue. compute updating direction g = −∇JDS [ml ] update ml , cycle
7/15
Agenda
,
1
Full Waveform Inversion
2
Extension and Differential Semblance Optimization (DSO)
3
Examples: Inversion in Layered Medium
8/15
Model x (km) 0
0
2
z (km)
0.2
0.4
x10 4 0.6
0.8 MPa
1.0
1.2
Three layer bulk modulus model. (acoustic velocity v = 1.5, 2.5, 2 km/s, density ρ = 1 g/cm3 ) Two sets of inversion exercises: absorbing surface free surface ,
9/15
Absorbing Surf: LS Inversion (standard FWI)
Time (s)
0
sign(p)p^2 (s^2/km^2) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08
0.5 1.0 1.5 2.0 -5
0
5
10
Band limited (5-12-36-45 Hz) plane wave source - 31 slowness panels
,
10/15
Absorbing Surf: LS Inversion (standard FWI)
Time (s)
0
sign(p)p^2 (s^2/km^2) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08
0.5 1.0
time (s)
-1000
0.2 0.4 0.6 0.8 1.0 1.2
0
1000
sign(p)*p^2 (s^2/km^2) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08
-500
0
500
Target data (Top), Data residual (Bottom) : after 40 LBFGS iterations - RMS ∼ 32%
,
10/15
Absorbing Surf: LS Inversion (standard FWI) true:blk, init:blue, fnl:red
x (km) 0
x10 4
2
0 1.2
1.0
z (km)
bulk modulus (Mpa)
0.2
0.8
0.6
0.4
0.4
0 3000
4000
5000 MPa
6000
7000
0.2
0.4 Depth (km)
Left: FWI estimate of bulk modulus after 40 LBFGS iterations; Right: final model slice at x = 1.5( km)
,
10/15
Absorbing Surf: DS Inversion with LF control 0
sign(p)p^2 (s^2/km^2) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08
z (km)
0.2
0.4
3000
4000
5000 MPa
6000
7000
Inverted model m[m ¯ l ], ml = homogeneous model
,
11/15
Absorbing Surf: DS Inversion with LF control sign(p)p^2 (s^2/km^2) -0.10 -0.07 -0.05 -0.02 0.01 0.03 0
0.06
0.08
z (km)
0.2
0.4
3000
4000
5000 MPa
6000
Inverted gather m[m ¯ l ], ml = homogeneous model, x = 1.5 km
,
11/15
Absorbing Surf: DS Inversion with LF control x (km) 0
0
2
z (km)
0.2
0.4
x10 4 0.2
0.4
0.6
0.8 MPa
1.0
1.2
Low frequency control model ml after 7 updates
,
11/15
Absorbing Surf: DS Inversion with LF control sign(p)p^2 (s^2/km^2) -0.10 -0.07 -0.05 -0.02 0.01 0.03 0
0.06
0.08
z (km)
0.2
0.4
x10 4 0.4
0.6
0.8 MPa
1.0
1.2
Inverted gather m[m ¯ l ], 7 updates of ml , x = 1.5 km
,
11/15
Absorbing Surf: DS Inversion with LF control x10 4
init-models (true:blk, 1:ylw, 2-6:grn, 7:red, 8:blue)
x10 4
fnl-models (true:blk, 1:ylw, 2-6:grn, 7:red, 8:blue)
1.2 1.2
1.0
bulk modulus (Mpa)
bulk modulus (Mpa)
1.0
0.8
0.8
0.6
0.6
0.4
0.4
0
0.2
0.4 Depth (km)
0
0.2
0.4 Depth (km)
Left: ml at x = 1.5( km); Right: final model at x = 1.5( km)
,
11/15
Free Surf: DS Inversion with LF control
0
sign(p)p^2 (s^2/km^2) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08
z (km)
0.2
0.4
3000
4000
5000 MPa
6000
7000
Inverted model m[m ¯ l ], ml = homogeneous model ! multiple reflections got suppressed during LS inversion
,
12/15
Free Surf: DS Inversion with LF control
0
sign(p)p^2 (s^2/km^2) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08
z (km)
0.2
0.4
-1000
-500
0 MPa
500
1000
LS gradient in 1st iteration (i.e., residual migration) ! lots of multiple reflections
,
12/15
Free Surf: DS Inversion with LF control
Time (s)
0
sign(p)p^2 (s^2/km^2) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08
0.5 1.0
-1000
Time (s)
0
0
1000
2000
sign(p)p^2 (s^2/km^2) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08
0.5 1.0
-1000
0
1000
2000
Target Data (Top), Predicted Data (Bottom) :after 60 LBFGS iterations ,
12/15
Free Surf: DS Inversion with LF control
Time (s)
0
sign(p)p^2 (s^2/km^2) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08
0.5 1.0
-100
-50
0
50
100
Data residual after 60 LBFGS iterations - RMS ∼ 8%
,
12/15
Free Surf: DS Inversion with LF control 0
-0.10
-0.07
sign(p)p^2 (s^2/km^2) -0.05 -0.02 0.01 0.03
0.06
0.08
0.05 0.10 0.15 0.20
z (km)
0.25 0.30 0.35 0.40 0.45 0.50 0.55
Inverted model gather m[m ¯ l ], ml = homogeneous model, x = 1.5 km ,
12/15
Free Surf: DS Inversion with LF control x (km) 0
0
2
z (km)
0.2
0.4
x10 4 0.4
0.6
0.8 MPa
1.0
1.2
Low frequency control model ml in the 3rd DS-iteration
,
12/15
Free Surf: DS Inversion with LF control 0
sign(p)p^2 (s^2/km^2) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08
z (km)
0.2
0.4
x10 4 0.2
0.4
0.6
0.8 MPa
1.0
1.2
Inverted model m[m ¯ l ] in the 3rd DS-iteration
,
12/15
Free Surf: DS Inversion with LF control 0
-0.10
-0.07
sign(p)p^2 (s^2/km^2) -0.05 -0.02 0.01 0.03
0.06
0.08
0.05 0.10 0.15 0.20
z (km)
0.25 0.30 0.35 0.40 0.45 0.50 0.55
Inverted model gather m[m ¯ l ] after 3 DS-iterations, x = 1.5 km DS-objective value reduced by 60% ,
12/15
Free Surf: DS Inversion with LF control
Time (s)
0
sign(p)p^2 (s^2/km^2) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08
0.5 1.0
-1000
Time (s)
0
0
1000
2000
sign(p)p^2 (s^2/km^2) -0.08 -0.06 -0.04 -0.02 0 0.02 0.04 0.06 0.08
0.5 1.0
-200
-100
0
100
Target Data (Top), Data Residual (Bottom) :after 60 LBFGS iterations - RMS ∼ 15% (in 3rd DS iteration) ,
12/15
Free Surf: DS Inversion with LF control init-mdl (true:blk, 1:grn, 2:blue, 3:red, 4:purp)
fnl-mdl (true:blk, 1:grn, 2:blue, 3:red) x10 4
1.2
1.2
1.0
1.0
bulk modulus (Mpa)
bulk modulus (Mpa)
x10 4
0.8
0.8
0.6 0.6
0.4 0.4 0
0.2
0.4 Depth (km)
0
0.2
0.4 Depth (km)
Left: ml at x = 1.5( km); Right: final model at x = 1.5( km)
,
12/15
Free Surf: DS Inversion with LF control x (km) 0
0
2
z (km)
0.2
0.4
x10 4 0.4
0.6
0.8 MPa
1.0
1.2
Inverted model stack (in 3rd DS iteration)
,
12/15
Free Surf: DS Inversion with LF control Standard LS inversion staring from the final model stack from nDSO x (km) 0
0
2
z (km)
0.2
0.4
x10 4 0.5
1.0 MPa
Inverted model after 153 LBFGS iterations (RMS residual ∼ 6%) (initial RMS res ∼ 77%; after 30 LBFGS RMS res ∼ 14%) ,
12/15
Free Surf: DS Inversion with LF control Standard LS inversion staring from the homogeneous model x (km) 0
0
2
z (km)
0.2
0.4
4000
6000
8000
MPa
Inverted model after 30 LBFGS iterations (RMS residual ∼ 27%)
,
12/15
Free Surf: DS Inversion with LF control x (km)
Time (s)
0
0
3
6
9
0.5 1.0
x10 4 -4
-2
0
2
4
Data fitting (for shot with slowness 0): Left, target data; Middle, residual for LS inv from DS-fnl-stack; Right, residual for LS inv from homogeneous model
,
12/15
Free Surf: DS Inversion with LF control true:blk, init:blue, fnl:red
true:blk, init:blue, fnl:red x10 4
x10 4 1.2
1.2
bulk modulus (Mpa)
bulk modulus (Mpa)
1.0
0.8
1.0
0.8
0.6 0.6
0.4 0.4 0
0.2
0.4 Depth (km)
0
0.2
0.4 Depth (km)
Left: LS inv fnl model at x = 1.5( km); Right: nDS inv fnl model at x = 1.5( km)
,
12/15
Remarks
Pros: - successfully infers global model updates for both absorbing surface and free surface - stack at least gives a good initial model for LS-FWI
,
13/15
Remarks Cost: 3-7 DS-iterations, mainly consisting of: 1 LS inversion (for m[m ¯ l ]) 2
∂ ¯ l] direction computation g = DFl [ml ]T DF l [m[m ¯ l ]]H[m[m ¯ l ]]−1 ∂s 2 m[m
solve H[m[m ¯ l ]]q =
∂2 m[m ¯ l] ∂s2
for q (∼ 30 CG iterations)
to reduce cost: (a) pre-conditioning; (b) scaling strategies; (c) replacing H with I (might work) ... compute g = DFl [ml ]T DF l [m[m ¯ l ]] q (1 Born Sim + 1 Adj Comp)
step-length computation (1 CG or Lin-LS + 1 Born Sim): ∂ hψ, ∂s m[m ¯ l ]i hψ, ψi DF [m[m ¯ l ]]† DF [ml ]g
α=−
, where ψ =
,
∂ ∂s
13/15
Summary
This approach may: combine best features of MVA and FWI (no linearization, scale separation assumptions required), address the spectral data incompleteness and local-minima issue, infer global changes in model and provide good initial model for FWI. Next to-do: explore this strategy with more tests, improve efficiency, ...
,
14/15
Acknowledgements Great thanks to my Ph.D. committee: William Symes, Matthias Heinkenschloss, Yin Zhang, Colin Zelt Present and former TRIP team members: Xin Wang, Marco Enriquez, Igor Terentyev, Tanya Vdovina, Rami Nammour ExxonMobil URC FWI team especially Dave Hinkley, Jerry Krebs WesternGeco FWI Management for support and permission to present here Sponsors of The Rice Inversion Project NSF DMS 0620821
Thank you! ,
15/15