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