UBC Faculty Research and Publications

Phase transitions in explorations seismology : statistical mechanics meets information theory Herrmann, Felix J. 2007

Your browser doesn't seem to have a PDF viewer, please download the PDF to view this item.

Item Metadata


52383-herrmann07COIP.pdf [ 4.74MB ]
JSON: 52383-1.0107425.json
JSON-LD: 52383-1.0107425-ld.json
RDF/XML (Pretty): 52383-1.0107425-rdf.xml
RDF/JSON: 52383-1.0107425-rdf.json
Turtle: 52383-1.0107425-turtle.txt
N-Triples: 52383-1.0107425-rdf-ntriples.txt
Original Record: 52383-1.0107425-source.json
Full Text

Full Text

Phase transitions in exploration seismology: statistical mechanics meets information theory Felix J. Herrmann joint work with  Gilles Hennenfent*, Mohamma Maysami* and Yves Bernabe (MIT) *Seismic Laboratory for Imaging and Modeling slim.eos.ubc.ca Complexity in the oil industry meeting 2007, Natal, Aug 9 Research interests • Develop techniques to obtain higher quality images from (incomplete) data <=> seismic imaging of transitions • Characterization of reflectors <=> estimation of singularity orders of imaged reflectors • Understand physical processes that generate singular transitions in the earth <=> Percolation phenomena • Singularity-preserved upscaling Today’s topics Phase diagrams in the recovery of seismic data from incomplete measurements • old ideas in geophysics reincarnated in the new field of  “compressive sampling” • describes regions that favor recovery Phase diagrams in the description of seismic reflectors • mixture models with critical points <=> reflectors • a first step towards singularity-preserved upscaling Phase-transition behavior in compressive sampling joint work with Gilles Hennenfent “Non-parametric seismic data recovery with curvelet frames” in revision for GJI Seismic Laboratory for Imaging and Modeling Data nominal spatial sampling ~ 112.5m Seismic Laboratory for Imaging and Modeling CRSI spatial sampling ~ 12.5m Fourier example Consider n-random time samples from a signal with k- sparse Fourier spectrum, i.e. with                      the time restricted inverse Fourier transform. The signal with the k-non-zero spectrum can exactly be recovered. signal =y sparse representation of inverse extrapolated data x0 A A ∈ Cn×N f0 = FHx0 Solve Recovery for Gaussian matrices when [Donoho and Tanner ‘06] For arbitrary measurement sparsity bases [Candes, Romberg & Tao ‘06] for certain conditions on the matrix and sampling .... P1 :  x̃ = arg minx∈RN ‖x‖1 = ∑N i=1 |xi| s.t. y = Ax f̃ = FH x̃. data consistencysparsity enhancement Fourier example cont’d When a traveler reaches a fork in the road, the l1 -norm tells him to take either one way or the other, but the l2 -norm instructs him to head off into the bushes.  [Claerbout and Muir, 1973] n = k × 2 log(N/k) n = µ2k × log N mutual coherence Figure 5: Phase Diagram for !1 minimization. Shaded attribute is the number of coordinates of recon- struction which differ from optimally sparse solution by more than 10−4. The diagram displays a rapid transition from perfect reconstruction to perfect disagreement. Overlaid red curve is theoretical curve ρ!1 . 4.2 Phase Diagram A phase diagram depicts performance of an algorithm at a sequence of problem suites S(k, n,N). The average value of some performance measure as displayed as a function of ρ = k/n and δ = n/N . Both of these variables ρ, δ ∈ [0, 1], so the diagram occupies the unit square. To illustrate such a phase diagram, consider a well-studied case where something interesting happens. Let x1 solve the optimization problem: (P1) min ‖x‖1 subject to y = Φx. As mentioned earlier, if y = Φx0 where x0 has k nonzeros, we may find that x1 = x0 exactly when k is small enough. Figure 5 displays a grid of δ − ρ values, with δ ranging through 50 equispaced points in the interval [.05, .95] and ρ ranging through 50 equispaced points in [.05, .95]; here N = 800. Each point on the grid shows the mean number of coordinates at which original and reconstruction differ by more than 10−4, averaged over 100 independent realizations of the standard problem suite Sst(k, n,N). The experimental setting just described, i.e. the δ − ρ grid setup, the values of N , and the number of realizations, is used to generate phase diagrams later in this paper, although the problem suite being used may change. This diagram displays a phase transition. For small ρ, it seems that high-accuracy reconstruction is obtained, while for large ρ reconstruction fails. The transition from success to failure occurs at different ρ for different values of δ. This empirical observation is explained by a theory that accurately predicts the location of the observed phase transition and shows that, asymptotically for large n, this transition is perfectly sharp. Suppose that problem (y,Φ) is drawn at random from the standard problem suite, and consider the event Ek,n,N that x0 = x1 i.e. that !1 minimization exactly recovers x0. The paper [19] defines a function ρ!1(δ) (called there ρW ) with the following property. Consider sequences of (kn), (Nn) obeying kn/n→ ρ and n/Nn → δ. Suppose that ρ < ρ!1(δ). Then as n→∞ Prob(Ekn,n,Nn)→ 1. On the other hand, suppose that ρ > ρ!1(δ). Then as n→∞ Prob(Ekn,n,Nn)→ 0. The theoretical curve (δ, ρ!1(δ)) described there is overlaid on Figure 5, showing good agreement between asymptotic theory and experimental results. 9 Phase diagrams l1 solver [from Donoho et al ‘06] In the white region recovers exactly. measurement vector y ∈ Rn with y= Ax0, by solving the following optimization problem P1 :  x̂= argminx ‖x‖1 subject to Ax= y f̂= SH x̂. (9) The synthesis matrix A ∈ Cn×N is composed of three matrices, namely, A := RMSH with SH the sparsity matrix S := F with f the Fourier analysis or decomposition matri ; M := I the Dirac mea- surement basis with I the idendity matrix and R a restriction matrix. This restriction matrix extracts k rows from the N×N Fourier matrix since MSH = F. During the restriction, the columns of A are normalized such that ‖ai}‖= 1 for i= 1 · · ·N. The mathematical criteria for exact recovery, such as the !0-norm of the difference between the original and the recovered sparsity vectors ‖x̂−x0‖0, is are unfortunately impossible because of the finite precision in floating point arithmetic. Instead, we either call a recovery successful when an entry is to within a small constant equal to the corresponding entry in x0 or we call the recovery successful when the relative !2-error is smaller then some threshold, ‖x̂−x‖2/‖x0‖2 ≤ ε . The main results of compressed sensing is that it predicts the number of measurements n that are required to ’exacly’ recover an arbitrary sparsity vector x0 with k non-zero entries given a certain choice for the measurement and sparsity matrices. For individual realizations of the sparsity and measurement vector, these conditions are sharp. For instance, the recovery of a sinusiodal function of length N = 1024 with k =? non-zero entries in x0 is succesful for a measurement vector y con- sisting of n =??? elements and fails for n =???− 2 measurements. This behavior is numerically illustrated in Fig. 1 and has also been observed by (35) for the spiky deconvolution. Strong recovery conditions: An important result from the literature on compressed sensing states that exact recovery from incomplete measurements is possible as long as the synthesis matrix A 13 fully sampled sparse signal fully sampled rich signal undersampled sparse signal undersampled rich signal recovery failed recovery Phase transition Has a second-order phase transition at a oversampling of 5  transition becomes sharper for  conceptual but unexplored ‘link’ with percolation theory  k-neighborhoodness of polytopes undergoes a phase transition !"# !"#$ !"% !"%$ !"& !"&$ !"' !"'$ ! !"# !"% !"& !"' !"$ !"( !") !"* !"+ # ,-./0123.456-75869 ! : .8 ; - ; 5/5 70 18 <1 , - ./ 0 12 3 .4 56 - 75 8 6 1 1 61=1'!! 6=*!! 6=#(!! !"# !"#$ !"% !"%$ !"& !"&$ !"' !"'$ ! !"# !"% !"& !"' !"$ !"( !") *+,,-./0-1-21+34, ! 5 67 2 1+ 3 4 /8 69 :; /< 3 4 < 9 ::/ * +, , - . / / 4/=/'!! 4=>!! 4=#(!! Figure 8: Empirical Transition Behaviors, varying n. (a) Fraction of cases with termination before stage S. (b) Fraction of missed detections. Averages of 1000 trials with n = 400, 800, 1600 and k = !ρn", N = !n/δ", δ = 1/4 and ρ varying. Sharpness of each transition seems to be increasing with n. To make the comparison still more vivid, we point ahead to an imaging example from Section 9.1 below. There an image of dimensions d× d is viewed as a vector x of length N = d2. Again the system y = Φx where Φ is made from only n = δN rows of the Fourier matrix. One matrix-vector product costs V = 4N logN = 8d2 log d. How do the three algorithms compare in this setting? Plugging-in S = 10, ν = 10, and V as above, we see that the leading term in the complexity bound for StOMP is 960 · d2 log d. In contrast, for OMP the leading term in the worst-case bound becomes 4δ 3 3 d 6 + 16δd4 log d, and for $1 minimization the leading term is 16d4 log d. The computational gains from StOMPare indeed substantial. Moreover, to run OMP in this setting, we may need up to δ 2 2 d 4 memory elements to store the Cholesky factorization, which renders it unusable for anything but the smallest d. In Section 9.1, we present actual running times of the different algorithms. 6 The Large-System Limit Figures 6 and 7 suggest phase transitions in the behavior of StOMP , which would imply a certain well-defined asymptotic ‘system capacity’ below which StOMPsuccessfully finds a sparse solution, and above which it fails. In this section, we review the empirical evidence for a phase transition in the large-system limit and develop theory that rigorously establishes it. We consider the problem suite S(k, n,N ;USE,±1) defined by random Φ sampled from the USE, and with y generated as y = Φx0, where x0 has k nonzero coefficients in random positions having entries ±1. This ensemble generates a slightly ‘lower’ transition than the ensemble used for Figures 6 and 7 where the nonzeros in x0 had iid Gaussian N(0, 1) entries. 6.1 Evidence for Phase Transition Figure 8 presents results of simulations at fixed ratios δ = n/N but increasing n. Three different quantities are considered: in panel (a), the probability of early termination, i.e. termination before stage S = 10 because the residual has been driven nearly to zero; in panel (b) the missed detection rate, i.e. the fraction of nonzeros in x0 that are not supported in the reconstruction x̂S . Both quantities undergo transitions in behavior near ρ = .2. Significantly, the transitions become more sharply defined with increasing n. As n increases, the early termination probability behaves increasingly like a raw discontinuity 1{k/n≥ρFAR(n/N)} as n →∞, while the fraction of missed detections properties behave increasingly like a discontinuity in derivative (k/n−ρFAR(n/N))+. In statistical physics such limiting behaviors are called first-order and second-order phase transitions, respectively. 13 n→∞ 2-D curvelets curvelets are of rapid decay in space curvelets are strictly localized in frequency x-t f-k Oscillatory in one direction and smooth in the others! Curvelets live in wedges in the 3 D Fourier plane... 3-D curvelets Seismic Laboratory for Imaging and Modeling Model spatial sampling:  12.5 m Seismic Laboratory for Imaging and Modeling avg. spatial sampling:  62.5 m Data 20% traces remaining Seismic Laboratory for Imaging and Modeling spatial sampling:  12.5 m SNR = 9.26 dB Interpolated result using CRSI Seismic Laboratory for Imaging and Modeling Difference spatial sampling:  12.5 m (no minimum velocity constraint) SNR = 9.26 dB New paradigm Traditional data collection & compression paradigm  ‘over emphasis’ on data collection  extract essential features  throw away the rest .... New paradigm compression during sampling  project onto measurements that breaks aliases  recover with sparsity promotion Exploration seismology  ‘random’ sampling of seismic wavefields [Hennenfent & F.J.H ‘06]  compressive wavefield extrapolation where eigenfunctions of the Helmholtz operator are used as the measurement basis [Lin & F.J.H ‘07] Characterizing singularities joint work Mohammad Maysami “Seismic reflector characterization by a multiscale detection-estimation method” ‘07 Problem • Delineate the stratigraphy from seismic images. • Parameterize seismic transitions • beyond simplistic reflector models • consistent with observed intermittent behavior of sedimentary records • Estimate the parameters from seismic images: • location • singularity order • instantaneous phase Singularity characterization through waveforms [F.J.H ’98, ’00, ’03, ‘07] • generalization of zero- & first-order discontinuities • measures wigglyness / # oscilations / sharpness Singularity characterization through waveforms [F.J.H ’98, ’00, ’03, ‘07] 1st order • generalization of zero- & first-order discontinuities • measures wigglyness / # oscilations / sharpness Singularity characterization through waveforms [F.J.H ’98, ’00, ’03, ‘07] 1st order 0-order • generalization of zero- & first-order discontinuities • measures wigglyness / # oscilations / sharpness Approach [Wakin et al ‘05-’-07, M&H ‘07] • Use a detection-estimation technique • multiscale detection => segmentation • multiscale Newton technique to estimate the parameterization • Overlay the image with the parameterization CWT Seismic trace Singularity map 8 8.05 8.1 8.15 8.2 8.25 8.3 1.4 1.6 1.8 2 2.2 2.4 −1 0 1 2 3 4 0 5000 10000 15000 20000 Offset [m] 2.5 2.6 2.7 2.8 2.9 3.0 3.1 3.2 T im e  [ s]   Estimated alpha -6.4 -5.6 -4.8 -4.0 -3.2 -2.4 -1.6 -0.8 0.0 Courtesy CGG Veritas Observations • Stratigraphy is detected • Parameterization provides information on the lithology • evidence of changes in exponents along clinoforms • Method suffers from curvature in the imaged reflectors • Extension to higher dimensions necessary • Model that explains different types of transitions • A step beyond the zero-& first-order discontinuities Modeling seismic singularities Joint work with Yves Bernabe (MIT) “Seismic singularities at upper-mantle phase transitions: a site percolation model” GJI ‘04 Problem Earth subsurface is highly heterogeneous • sedimentary crust, upper-mantle transition zone & core-mantle boundary Smooth relation volume fractions and the transport properties. Homogenization/equivalent medium (EM) theory smoothes the singularities during upscaling • relatively easy for volumetric properties (density) • notoriously difficult for transport  properties (velocity) Q: How to preserve singularities in effective properties? Our approach Include connectivity in models for the effective properties of bi-compositional mixtures <=> SWITCH Start with binary mixtures, e.g. • sand-shale • gas-hydrate, Opal-Opal CT • upper-mantle mineralogy Studied two cases: • elastic properties upper mantle (H & B ‘04) • fluid-flow properties synthetic rock (B & H ‘04) Mixture laws for binary mixtures Elastic case: • Controlled by connectivity of stiffest phase Fluid-flow case: • Controlled by connectivity of high- conductivity phase Note: Stiff phase = low porosity, low conductivity phase No obvious link elastic-fluid flow properties ... Singularity modeling binary mixtures HP LP olivine β-spinel Site percolation volume fraction elastic properties random process Text Singularity modeling binary mixtures HP LP olivine β-spinel Site percolation volume fraction elastic properties random process Text Singularity modeling binary mixtures HP LP olivine β-spinel Site percolation volume fraction elastic properties random process Text Mixing model Homogeneous mixing (e.g., solid solution) of two phases (LP weak and HP strong) can only produce gradually varying elastic properties. If Heterogeneous (e.g. emerging random macroscopic inclusions) mixing, then a singularity in the elastic properties must arise at the depth where the strong, HP phase becomes connected (observed in binary alloys). Site-percolation model Assume volume fractions p and q =1–p, are linear functions of depth z. At a critical depth zc, which corresponds to the percolation threshold pc = p(zc), an "infinite", connected HP cluster is formed. For z ≥ zc • not all HP inclusions belong to the infinite cluster. • isolated HP inclusions can still be found, embedded in the remaining LP material and forming with it a mixture (M). Site-percolation model Above zc we have a weak LP matrix containing randomly distributed, non-percolating, strong HP inclusions. Below zc, a strong HP skeleton is intertwined with the weaker, mixed material M. Site-percolation model Volume fraction p* of HP material that belongs to the infinite cluster • is zero for p < pc (i.e., above zc) • has a power-law dependence on (p - pc) for 
 p ≥ pc. Hence, p* is given by:  p∗ = p ( p− pc 1− pc )β Site-percolation model Mixed M is given by q* = (1 - p*). For M, we need the volume fractions of its LP and HP parts,  qM = (1 - p)/{(1 - p) + (p - p*)}  pM = (1 - qM), yielding pM = 1− q 1− p ( p−pc 1−pc )β Site-percolation model Binary mixture: • Strong when its strong component is connected. • Weak otherwise. Assume locally isotropic • Bin. mixtures bounded by Hashin-Shtrikman (HS). • upper HS bound when strong component connects, the lower one applies otherwise. Bulk modulus K of the co-existence region above zc is given by the lower HS bound: Site-percolation model Below zc we must switch to the higher HS bound: Since the HP inclusions in M are isolated, KM is calculated using the lower HS bound: Site-percolation model Major consequence of this model is that it predicts: • a β-order, cusp-like singularity in the elastic moduli as the critical depth zc is approached from below (instead of a first-/ zero-order discontinuity). • Singularities that persist for vanishing contrasts. • Density that does not behave singularly. Site-percolation model Elastic contrasts between LP and HP are small: • Nearly coincident HS bounds. • Excessively small contrasts. Discard isotropy assumption: •Horizontally-oriented oblate ellipsoidal inclusions which coalesce below zc into long, vertical dendrites, leaving prolate M inclusions between them. •Transversely isotropic structure. Site-percolation model Near normal incidence, VP and VS approach limiting values as the aspect ratio is goes to zero. Same as replacing lower and higher HS bounds by Reuss and Voigt averages: Singularity model Singularity model upper-mantle transitions Modeled data vs seismic Singularity-preserved upscaling Joint work with Yves Bernabe (MIT) The problem • Equivalent medium based upscaling washes out the singularities • Reflection seismology lives by virtue of singularities in the elastic moduli (transport properties) • Propose a singularity preserving upscaling method: • upscales the lithology rather than the velocities • Singularities can be due to sharp changes in composition or due to the switch ... Volume fraction synthetic well 0 500 1000 1500 2000 2500 3000 3500 4000 4500 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 depth Vo lum e fra cti on  sh ale Synthetic well log Courtesy Chevron Switch vs no switch 0 500 1000 1500 2000 2500 3000 3500 4000 4500 1680 1700 1720 1740 1760 1780 1800 depth Ve loc ity Velocity from lithology Equivalent medium Percolation model 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1680 1700 1720 1740 1760 1780 1800 Volume fraction shale Ve loc ity Equivalent medium Percolation model Switch vs no switch singularity Lithological upscaling vo lu m e fr ac tio n 0 2 4 6 8 10 0 500 1000 1500 2000 2500 3000 3500 4000 Upscaling of the lithology (volume fractions) EM-upscaled reflectivity !1 0 1 2 3 4 5 6 7 8 9 10 0 500 1000 1500 2000 2500 3000 3500 4000 Reflectivity for the Equivalen medium model lost singularities Perco.-upscaled reflectivity !1 0 1 2 3 4 5 6 7 8 9 10 0 500 1000 1500 2000 2500 3000 3500 4000 Reflectivity for the Percolation model preserved singularities Relation to fluid flow & open problems joint work with Yves Bernabe (MIT) Sedimentary crust What does this model buy us? Some insight in • the complexity of transitions • the creation of a singularity for smooth varying composition, e.g. when clay lenses connect ... • the morphology at transitions • linking elastic and fluid properties remains a challenge .... Fluid flow • difficult to model • difficult to measure [B&H ‘04] Connectivity of the high conductive phase measured modeled numerical simulations kHP/kLP from ~1 to 106 Steady-state flow equation: Grid up to 50 X 50 X 50 Sand-clay model according HB model [HB ‘04] 0 0.2 0.4 0.6 0.8 1 1000 1200 1400 1600 1800 2000 sand volume fraction Vp  a nd  V s  (m /s ) Vp Vs 0 0.2 0.4 0.6 0.8 1 -15 -14 -13 -12 sand volume fraction lo g 1 0( pe rm ea bi li ty )  ( m2 ) HS+ HS- 1950 2000 2050 2100 2150 -15 -14 -13 -12 Vp   (m/s) lo g 1 0( pe rm ea bi li ty )  ( m2 ) HS+ HS- 840 860 880 900 920 940 960 980 -15 -14 -13 -12 lo g 1 0( pe rm ea bi li ty )  ( m2 ) HS+ HS- Vs   (m/s) sand clay sand clay velocities <=> permeability Elastic versus Fluid • singularities in elastic properties are small • singularities in fluid properties are large • waves are way more sensitive to singularities then diffusion driven fluid flow • models for fluid and elastic transport are not integrated • incorporate in Biot? 0.2 0.4 0.6 0.8 1 0.025 0.05 0.075 0.1 0.125 0.15 0.175 0.2 0.2 0.4 0.6 0.8 1 -0.175 -0.15 -0.125 -0.1 -0.075 -0.05 -0.025 0 dV p/ dp dk /d p dV s/ dp p p p Conclusions • Multiscale compressible signal representations that exploit higher-dim. geometry are indispensable for acquiring accurate information on the imaged waveforms. • Imaged waveforms carry information on the fine structure of the reflectors. • Percolation model provides an interesting perspective • on linking the micro-connectivity to singularities detected by waves • on providing an upscaling that preserves features = singularities that matter .... Acknowledgments The authors of CurveLab (Demanet, Ying, Candes, Donoho) Chevron, Statoil, CGG-Veritas for providing data. This work was mainly conducted as part of ChARM supported by Chevron. This work was also partly supported by the Natural Sciences and Engineering Research Council of Canada Discovery Grant (22R81254) and the Collaborative Research and Development Grant DNOISE (334810-05) of F. J. H matching ChARM and the SINBAD project with support, secured through ITF (the Industry Technology Facilitator), from the following organizations: BG Group, BP, Chevron, ExxonMobil and Shell.


Citation Scheme:


Citations by CSL (citeproc-js)

Usage Statistics



Customize your widget with the following options, then copy and paste the code below into the HTML of your page to embed this item in your website.
                            <div id="ubcOpenCollectionsWidgetDisplay">
                            <script id="ubcOpenCollectionsWidget"
                            async >
IIIF logo Our image viewer uses the IIIF 2.0 standard. To load this item in other compatible viewers, use this url:


Related Items