Transcript G.A. Newman, New approach to joint geophysical imaging of
A New Approach to Joint Imaging of Electromagnetic and Seismic Wavefields
International Symposium on Geophysical Imaging with Localized Waves Sanya, Hainin Island China July 24-28, 2011 Gregory A. Newman Earth Sciences Division Lawrence Berkeley National Laboratory
PRESENTATION OVERVIEW
Controlled Source EM and Magnetotellric Data Acquisition – Motivation for Joint Imaging Formulation of the Joint Inverse Problem – Large Scale Modeling Considerations – Need for High Performance Computing Joint MT/CSEM Imaging Results – Gulf of Mexico (synthetic example) Joint Imaging - EM & Seismic Data – Issues & Proposed Methodology – Laplace-Fourier Wave Field Concepts » Similarities and Differences to the EM Wave Field » Imaging Across Multiple Scale Lengths » 3D Elastic Wave Fields » Preconditioners and solvers for 3D Seismic Wave Field Simulations Conclusions
Marine CSEM & MT Surveying
CSEM Deep-towed Electric Dipole transmitter ~ 100 Amps Water Depth 1 to 7 km Alternating current 0.01 to 3 Hz ‘Flies’ 50 m above the sea floor Profiles 10’s of km in length Excites vertical & horizontal currents Depth of interrogation ~ 3 to 4 km Sensitive to thin resistive beds MT Natural Source Fields Less than 0.1 Hz Measured with CSEM detectors Sensitive to horizontal currents Depth of interrogation 10’s km Resolution is frequency dependent Sensitive to larger scale geology
JOINT 3D INVERSE MODELING
Minimize: N
j = a S {(d
j obs j=1
- d
j p
)/ e
j
}
2 M
+ b S {(Z
j obs j=1
- Z
j p
)/ p
j
}
2
+ l
h
m
h
W
T
W m
h
+ l
v
m
v
W
T
W m
v
d
obs and
d
p
s.t.
m
v are N observed and predicted CSEM data
Z
e obs and
Z
p are M observed and predicted MT impedance data & p = CSEM and MT data weights
m
h
m
h
=
horizontal conductivity parameters
m
v
=
vertical conductivity parameters
W =
2 operator; constructs a smooth model l h & a & l v b
=
horizontal & vertical tradeoff parameters
=
scaling factors for CSEM and MT data types
LARGE-SCALE 3D MODELING CONSIDERATIONS
Require Large-Scale Modeling and Imaging Solutions – 10’s of million’s field unknowns (fwd problem) » Solved with finite difference approximations & iterative solvers – Imaging grids 400 nodes on a side » Exploit gradient optimization schemes, adjoint state methods Parallel Implementation – Two levels of parallelization » Model Space (simulation and inversion mesh) » Data Space (each transmitter/MT frequency - receiver set fwd calculation independent) » Installed & tested on multiple distributed computing systems; 10 – 30,000 Processors Above procedure satisfactory except for very largest problems – To treat such problems requires a higher level of efficiency Optimal Grids – Separate inversion grid from the simulation/modeling grid – Effect: A huge increase in computational efficiency ~ can be orders of magnitude
Optimal Grids
10 km Ω m imaging grid Ω s 100 km simulation grid sail lines
GRID SEPARATION
EFFICIENCIES
Advantages – Taylor an optimal simulation grid Ω s for each transmitter receiver set – Inversion grid Ω m covers basin-scale imaging volumes at fine resolution – Simulations grids much smaller, a subset of the imaging grid – Faster solution times follow from smaller simulation grids What’s Required – A mapping of conductivity from Ω m to Ω s » Conductivity on Ω s edged based » Conductivity on Ω m cell based & Ω s to Ω m – An appropriate mixing law for the conductivity mappings
Joint CSEM - MT Imaging
Mahogany Prospect, Gulf of Mexico Study: 3D Imaging of oil bearing horizons with complex salt structures present Simulated Example: 100 m thick reservoir, 1 km depth, salt below reservoir Model: 0.01 S/m salt, 2 S/m seabed, 0.05 S/m reservoir, 3 S/m seawater MT Data: 7,436 data points, 143 stations & 13 frequencies 0.0005 to 0.125 Hz CSEM Data: 12,396 data points, 126 stations & 2 frequencies 0.25 and 0.75 Hz Starting Model: Background Model without reservoir or salt Processing Times: 5 to 9 hours, 7,785 tasks, NERSC Franklin Cray XT4 System Survey Layout -10 -5 0 x(km) y=5 km cross section 5 10
JOINT CSEM-MT IMAGING:
The Benefits
Joint CSEM - MT Imaging Mahogany Prospect Gulf of Mexico
Joint Imaging of EM and Seismic Data
Issues
– Rock Physics Model » links attributes to underlying hydrological parameters » too simplistic » difficult or impossible to define robust/realistic model – Differing Resolution in the Data » EM data 10x lower resolution compared to seismic – RTM & FWI of Seismic Data » requires very good starting velocity model » velocity can be difficult or impossible to define » huge modeling cost due to very large data volumes (10,000’s of shots; 100,000’s traces per shot)
Joint Imaging of EM and Seismic Data
A way forward
– Abandon Rock Physics Model » assume conductivity and velocity structurally correlated » employ cross gradients: t = » t = 0 => ; = 0 and/or =0 – Equalize Resolution in the Data » treating seismic and EM data on equal terms » Laplace-Fourier transform seismic data – Shin & Cha 2009
g
ˆ (
s
) = 0
g
(
t
)
e
st dt g
(
s
ˆ )
and s are complex
Acoustic Wave Equation
Propagating Wave 1
v
2 Time Domain
t
2 2
x
2 2
y
2 2
z
2 2
p
(
x
,
y
,
z
,
t
) =
s
(
t
).
Fourier/Frequency Domain
v
2 2
x
2 2
y
2 2
z
2 2
p
(
x
,
y
,
z
, ) =
s
( ).
Damped Diffusive Wave
s
2
v
2 Laplace/Fourier Domain
x
2 2
y
2 2
z
2 2 ˆ
p
(
x
,
y
,
z
,
s
) = ˆ
s
(
s
).
At first glance similar physics & similar resolution with EM fields skin depth: =
s r
Seismic Imaging: Laplace-Fourier Domain
337 shot gathers 151 detectors/shot maximum offsets 15km BP Salt Model
s
= 10.5 to 0.5
=0.5
Laplace Image
s
= 10.5 to 0.5
f
= 6 to 0.5
=0.5
Taken from Shin & Cha, 2009 Laplace-Fourier Image Starting Velocity Model Standard FWI Image New FWI Image
Laplace-Fourier Wavefield Modeling
There are differences compared to EM fields
– wavelength and skin depth are decoupled
Meshing Issues to Consider
– grid points per wavelength:10 points – 2 nd order accuracy – grid points per skin depth: 6 points – 2 nd order accuracy
Accuracy Issues
– wavefield dynamic range extreme ~ 70 orders – iterative Krylov solvers require tiny solution tolerances 10 -10 10 -20 10 -30 10 -40 10 -50 10 -60 0 5 tol= 10 Offset (km) 15 Analytic 10 -30 10 -50 10 -70 10 -90 10 -110 10 -130 20
LAPLACE-FOURIER IMAGING:
Mahogany Prospect Study: 3D Imaging of oil bearing horizons with complex salt structures present Simulated Example: 100 m thick reservoir, 1 km depth, salt below reservoir Model: 6 km/s salt, 3 km/s seabed, 2 km/s reservoir, 1.5 km/s seawater Seismic Data: 24,577 data points, 287 stations at 1 frequency (σ=1, ω=2π) Starting Model: Background Model without reservoir or salt Processing Times: 10.3 hours, 6,250 tasks, NERSC Franklin Cray XT4 System Survey line N 10000 km * 287 sources, (σ=1, ω=2π) Source & Receiver & Spacing 1 km & 300 m Max. offsets 17 km 50 m below sea surface Misfit Survey line N 7500 km Survey line N 5000 km Survey line N 2500 km Survey line N 0 km Survey line N -2500 km Survey line N -5000 km -20 km West-East 20 km
-0.9 km
LAPLACE-FOURIER IMAGE:
Mahogany Prospect -0.9 km 0 km North 0 km North 13.5 km 13.5 km 2.5 km North 2.5 km North 5.0 km North 7.5 km North 5.0 km North 7.5 km North 10.0 km North 10.0 km North 12.0 km North 18.5 km 12.0 km North -18.5 km -18.5 km 18.5 km
12.0 km North
LAPLACE-FOURIER IMAGE:
Mahogany Prospect 1km below seabed 12.0 km North -5.0 km North -18.5 km West – East 18.5 km -5.0 km North -18.5 km West – East 18.5 km
j
Joint EM-Seismic Imaging
Problem Formulation = a
j N
= 1
d em obs j
p d em j
/ e
j
2 b
l M
= 1
d
ˆ
s obs l
d
ˆ
l p l
/
j
2 l
em
σ
T
W
T
W
σ
l
s
v
T
W
T
W
v
i m c
= 1
t i
t i
obs d em
and
p d em
are
N
observed and predicted EM data
d
ˆ
s obs
e
obs
and are
s M
and observed and predicted Laplace-Fourier seismic data = EM and seismic data weights
σ
= m conductivity parameters
v
= m acoustic velocity parameters
W
l
em
=
2 operator; constructs a smooth model and
=
conductivity & velocity tradeoff parameters
t
a and b
c
=
scaling factors for EM and seismic data types
Recipe for Auxiliary Parameters
Consider total objective functional : j = aj
d em
bj
d s
l
em
j
m
l
s
j
m v
j
xg
First carry out separate inversions for seismic and EM => choose smoothing parameters (cooling approach)
em s
Next balance data funtionals for seismic and EM : => set b l = 1 ; a = j j
d s em d
; = 0 aj
d em
bj
d s
j
d s
; j
d em
Test out values for => selected j
xg
balances aj
d em
; bj
d s
j
xg
Initial Imaging Results
marine example
f
= 0.25
conductivity image correlated with velocity;
=10 11 conductivity image no correlation to velocity;
=0
CSEM 16 km max. offsets 17 shots 161 detectors/shot s
=5,4,3,2,1
f
= 0
velocity image correlated with conductivity;
=10 11 velocity image no correlation to conductivity;
=0 Computational Requirements: 4250 cores – Franklin NERSC Processing Time: ~22 hours
seismic 12-16 km offsets 85 shots 121-161 detectors/shot
Elastic Wave Field Simulator First- order system for velocity –stress components Laplace-Fourier Domain
s
xy s
xz s
s
xx s
yz s
zz yy x y z
= = = = = = = = =
div div div
x
f x
,
f y f z
, ,
x y z
= = =
v y x
z v z x v y
x v y
;
v x z
v y z
; ; l l l 2
v x x
; 2
y v y
2
v z z
.
;
xz xx xy xy
,
yy
,
yz
,
xz zz yz
, ; ;
v
l - velocity components, , , , , - stress components, - density, and - Lame coefficients.
f x
,
y
,
z
Μ
Moment-Tensor components (R. Graves 1996)
Boundary and Initial Conditions
The ordinary initial conditions for all components are zeros.
The boundary conditions are – a) PML absorbing boundary conditions for velocity – b) free surface boundary for normal stress component
τ
Solution Realization
Iterative Krylov Methods Coupled System
s
τ = λμ D v
,
s
v
=
b D τ
v
f
.
v
- matrices of FD first derivative operators
D τ
=
D D y D D
0
x z x D x
b
=
b x
0 0
D y D x
0
D y D z D y
0
b y
0
D
0
D x D D y z z
;
D v
D z
=
D x D y D
0 0 0
z
0 0
z
;
b
λμ
= l 0
D x
0
D y D z
0 0 0
D x
0
D y D z
T
; 2
xy
l
xz
0 l l
xy
l 0 2 l
yz
System transformed : solve only for the velocity components 0
v
v y =
s
f
f f
, =
D D D 11 21 31 D 12 D 22 D 32 D 13 D 23 D 33
.
200 400 600 800 1000 D is complex non-symmetric 15 diagonals for 2 nd order scheme 4 th order scheme 51 diagonals 1200 1400 1600 1800 0 500 1000 nz = 14798 1500 l l 0 l
xz
2
yz
T
;
Solution Accuracy
Two Half Spaces Model Test
SEG/EAEG SALT MODEL TEST
real 3D salt body 200x200x130 nodes y=4800 m s=(3+18.85i) sec -1 imaginary y=4800 m.
500 1000 1500 2000 2500 3000 3500 4000 4500 5000 -2000 0 2000 4000 x, m 6000 8000 10000 12000 4000 3500 3000 2500 2000 1500 Vertical cross section of velocity ( p ) y=4800 m s=(3+18.85i) sec -1 Snap-shots of velocity field y-component Point source at r=(800,800,500), m ).
Solution Time 1355 sec, 576 Iterations Solution Tolerance 1e-7
Laplace-Fourier Transformation
Benefits
– Wave field simulations » excellent choice for a preconditioner (frequency domain ) On a class of preconditioners for solving the Helmholtzs equation: Erlangga et al., 2004, Applied Numer. Mathematics, 50 409-425.
– Imaging » possibilities to image at multiple scales and attributes A consistent Joint EM seismic imaging approach » known to produce robust macro-models of velocity Critical to successful RTM and FWI of seismic reflection data
Conclusions
Demonstrated Benefits Massively Parallel Joint Geophysical Imaging – Joint CSEM & MT – acoustic seismic (Laplace-Fourier Domain) – HPC essential Future Developments in LF elastic wavefield imaging – massively parallel (MP) LF elastic wavefield simulator – exploit simulator as a preconditioner » frequency domain wavefield modeling – exam MP direct solvers – gradient based 3D LF elastic imaging code Future Plans in Joint Imaging – joint EM and elastic wavefield 3D imaging capability
ACKNOWLEDGEMENTS
My Colleagues: M. Commer, P. Petrov and E. Um Research Funding US Department of Energy Office of Science Geothermal Technologies Program