G.A. Newman, New approach to joint geophysical imaging of

Report
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
M
j=1
j=1
j = a S {(djobs - djp)/ej}2 + b S {(Zjobs - Zjp)/pj}2
+ lh mh WT W mh + lv mv WT W mv
s.t. mv  mh
dobs and dp are N observed and predicted CSEM data
Zobs and Zp are M observed and predicted MT impedance data
e & p = CSEM and MT data weights
mh = horizontal conductivity parameters
mv = vertical conductivity parameters
W = 2 operator; constructs a smooth model
lh & lv = horizontal & vertical tradeoff parameters
a & b = 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
10 km
100 km
Optimal Grids
10 km
Ωm imaging grid
100 km
Ωs simulation grid
sail lines
GRID SEPARATION
EFFICIENCIES

Advantages
– Taylor an optimal simulation grid Ωs for each transmitterreceiver 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 & Ωs to Ωm
» Conductivity on Ωs edged based
» Conductivity on Ωm cell based
– 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
y=5 km cross section
-10
-5
0
x(km)
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) =  g (t )e  st dt
0
g (sˆ) and s are complex
Acoustic Wave Equation
Propagating Wave
Fourier/Frequency Domain
Time Domain
  2   2  2  2 
 1  2   2  2  2 
 p( x, y, z,  ) =  s( ).



 2 2   2  2  2  p( x, y, z, t ) =  s(t ).  v 2  x 2 y 2 z 2 



 v t  x y z 
Damped Diffusive Wave
Laplace/Fourier Domain
 s2
 2
 v
  2  2  2 
  2  2  2  pˆ ( x, y, z, s) =  sˆ( s).
 x y z 
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
Starting Velocity Model
Laplace Image
Standard FWI Image
s = 10.5 to 0.5
=0.5
s = 10.5 to 0.5
f = 6 to 0.5
=0.5
Laplace-Fourier Image
Taken from Shin & Cha, 2009
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 – 2nd order accuracy
– grid points per skin depth: 6 points – 2nd order accuracy
Accuracy Issues
– wavefield dynamic range
extreme ~ 70 orders
– iterative Krylov solvers
require tiny solution
tolerances
10
-10
Analytic
10
10
Amplitude

10
10
10
10
-30
10-50
-20
10-70
tol=
-30
10-90
-40
10
-110
10
-130
-50
-60
0
5
10
Offset (km)
15
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
LAPLACE-FOURIER IMAGE:
-0.9 km
Mahogany Prospect
0 km North
0 km North
13.5 km
13.5 km
2.5 km North
5.0 km North
-18.5 km
-0.9 km
2.5 km North
5.0 km North
7.5 km North
7.5 km North
10.0 km North
10.0 km North
12.0 km North
12.0 km North
18.5 km
-18.5 km
18.5 km
LAPLACE-FOURIER IMAGE:
Mahogany Prospect
1km below seabed
12.0 km North
12.0 km North
-5.0 km North
-5.0 km North
-18.5 km
West – East
18.5 km
-18.5 km
West – East
18.5 km
Joint EM-Seismic Imaging
Problem Formulation
N

j =a d
j =1
obs
em j
d
p
em j
/ e   b  dˆ
M
2
j
l =1
obs
sl
 
 dˆl l /  j
p
mc
 lemσ W Wσ  ls v W Wv    t i  t i
T
T
T
T
i =1
obs
p
dem
and d em are N observed and predicted EM data
dˆsobs and dˆsobs are M observed and predicted Laplace-Fourier seismic data
e and  = EM and seismic data weights
σ = m conductivity parameters
v = m acoustic velocity parameters
W = 2 operator; constructs a smooth model
lem and ls = conductivity & velocity tradeoff parameters
a and b = scaling factors for EM and seismic data types
t are mc cross gradient structural constraints; is a penalty parameter
2
Recipe for Auxiliary Parameters
Consider total objective functional : j = ajdem  bjds  lemjm  lsjmv jxg
First carry out separate inversions for seismic and EM
=> choose smoothing parameters (lem ls ) (cooling approach)
Next balance data funtionals for seismic and EM :j
=> set
jds
b = 1 ;a = em ; = 0
jd
=> rescale lem accordingly
Test out values for
=> selected 
s
d
; j
em
d
ajdem  bjds

 j xg balances aj ; bj
em
d
s
d
jxg
=> trail values tested out over a few inversion iterations
Initial Imaging Results
marine example
f = 0.25
s =5,4,3,2,1
f=0
conductivity image
correlated with velocity;
=1011
velocity image
correlated with conductivity;
=1011
conductivity image
no correlation to velocity;
=0
velocity image
no correlation to conductivity;
=0
Computational Requirements: 4250 cores – Franklin NERSC
Processing Time: ~22 hours
CSEM
16 km max. offsets
17 shots
161 detectors/shot
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 vx, y, z - velocity components,
s  vx = div  x   f x ,  x =  xx , xy , xz  ;
s  v y = div  y   f y ,  y =  xy , yy , yz  ;
 xx, xy, xz , yz , yy, zz
- stress components,
 - density,
l and  - Lame coefficients.
s  vz = div  z   f z ,  z =  xz , yz , zz  ,
s xy =    y vx   x v y  ;
s xz =    z vx   x vz  ;
s yz =    z v y   y vz  ;
s xx = l div  v   2 x vx ;
s yy = l div  v   2 y v y ;
s zz = l div  v   2 z vz .

Forces f x, y, z are defined via   Μ
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,
sv = b Dv τ  f .
 Dx

 Dy
D
Dτ =  z
 Dx
 0

D
 x
 b x

b = 0

 0

Dτ ; Dv - matrices of FD first derivative operators
Dz 
 Dx


0 
 Dy
D
Dx 
 ; Dv =  z
Dz 
 0

 0
Dy


Dz 
 0
Dy
Dx
0
Dy
Dz
Dy
0
b
0
y
0
Dx
0
Dy
Dz
0
T
0 

0 
Dx 
 ;
0 
Dy 

Dz 
 l 2 

xy



0

xz




0 ; λμ =


l
z


b


0

l

l

l
xy
0
l 2 

l
yz
System transformed : solve only for the velocity components
0
f x 
 vx 
 D11
ˆ   v  =  s f  , D
ˆ = D
D
y
y


 
 21
fz 
 v z 
 D31
 
D12
D22
D32
D13 
D23  .
D33 
200
400
600
800
1000
1200
D is complex non-symmetric
15 diagonals for 2nd order scheme
4th order scheme 51 diagonals
1400
1600
1800
0
500
1000
nz = 14798
1500
T


0


xz

 ;

l

yz



l  2  
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.
4000
500
1000
3500
1500
z, m
2000
3000
2500
3000
2500
3500
4000
2000
4500
y=4800 m
s=(3+18.85i) sec-1
5000
-2000
0
2000
4000 6000
x, m
8000 10000 12000
1500
Vertical cross section of velocity (p)
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
Computational Details

similar documents