UBC Faculty Research and Publications

Seismic reflector characterization by a multiscale detection-estimation method Maysami, Mohammad; 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-maysami07eage.pdf [ 702.25kB ]
JSON: 52383-1.0107405.json
JSON-LD: 52383-1.0107405-ld.json
RDF/XML (Pretty): 52383-1.0107405-rdf.xml
RDF/JSON: 52383-1.0107405-rdf.json
Turtle: 52383-1.0107405-turtle.txt
N-Triples: 52383-1.0107405-rdf-ntriples.txt
Original Record: 52383-1.0107405-source.json
Full Text

Full Text

ABSTRACT Seismic reflector characterization by a multiscale detection-estimation method Seismic transitions of the subsurface are typically considered as zero-order singularities (step functions). According to this model, the conventional deconvolution problem aims at recovering the seismic reflectivity as a sparse spike train. However, recent multiscale analysis on sedimentary records revealed the existence of accumulations of varying order singularities in the subsurface, which give rise to fractional-order discontinuities. This observation not only calls for a richer class of seismic reflection waveforms, but it also requires a different methodology to detect and characterize these reflection events. For instance, the assumptions underlying conventional deconvolution no longer hold. Because of the bandwidth limitation of seismic data, multiscale analysis methods based on the decay rate of wavelet coefficients may yield ambiguous results. We avoid this problem by formulating the estimation of the singularity orders by a parametric nonlinear inversion method.  EAGE 69th Conference & Exhibition — London, UK, 11 - 14 June 2007  Introduction The earth’s subsurface consists of layers of different materials separated by interfaces, also called transitions. Transitions are characteristic of regions where the acoustic properties of the earth vary rapidly compared to the length-scale of the seismic source wavelet. Extracting information on the locations and nature of the transitions from seismic data has recently received increasing interest, as it provides quantitative information that can be used in various applications, ranging from improving the geological interpretation of the subsurface to detecting changes in the lithology. Seismic deconvolution is aimed at finding the locations of seismic reflectors and is based on certain assumptions on the reflectivity. For instance, deconvolution methods have been developed for a reflectivity that behaves as a (colored) Gaussian random process (Saggaf & Robinson, 2000). Other approaches are based on the sparsity assumption, where the reflectivity is considered to be given by a sparse spike train. Multiscale analysis on well and seismic data have shown that neither assumption is rich enough to describe the different types of transitions present in sedimentary basins (Herrmann, 2001). These observations have led to the introduction of a new type of parametrization, where the observed reflectivity is written as a superposition of parametrized waveforms. This parametrization is designed to reflect the presence of transitions other than strictly zero-order transitions (blocked wells) and includes fractional-order transitions. The aim of this paper is two-fold, namely finding the locations of the reflectors, delineating the stratigraphy, and extracting information on the nature of the transitions. We present a new detection-estimation method, where the reflection events are first detected, then segmented, followed by an estimation based on a descent method (Boyd & Vandenberghe, 2004). The estimated parameters provide information on the transition sharpness that is related to the lithology (Herrmann, 2001; Liner et al., 2004), possibly through a critical point in the elastic moduli (Herrmann & Bernab´e, 2004). The Earth’s Model: We represent a vertical 1-D profile of the earth as a superposition of parametrized waveforms of the following type s(z) = cj Dαj ψ(z − zj ), (1) j  where z is depth, cj the amplitude for the j th transition, and Dα α-order operator for fractional differentiation (α > 0) or integration (α < 0). The ψ is some wavelet taken to be the Ricker wavelet. The parameters (attributes) of interest in this case are the location zj , the amplitude cj and the order αj .  The characterization problem Given the above signal representation, our task is to recover the different attributes from a seismic trace. We divide this task into a detection stage, where the main events in the data are detected, and an estimation stage, where the parameters of the individual waveforms are estimated. Since our problem does not fit into the classical deconvolution framework, we use a multiscale wavelet technique to locate the main events. After segmentation, the individual waveforms are submitted to a nonlinear inversion procedure to estimate the attributes. This procedure uses rough estimates for the location and scale from the detection stage. Event location by multiscale edge detection: The variety of different orders of transitions in the subsurface calls for a seismic-event-detection technique that does not make  EAGE 69th Conference & Exhibition — London, UK, 11 - 14 June 2007  any assumptions regarding the type of transitions. Edge detection based on the multiscale continuous (complex) wavelet transform modulus maxima (Mallat, 1999) offers an approach that is robust for different waveforms, reflecting different types of transitions. First, the method calculates the forward wavelet transform Ws(σ, t) = s ∗ ψ¯σ (t), (2) where ψ¯σ (t) = √1σ ψ ∗ ( −t ). The range of scales σ for the wavelet is adapted to the seisσ mic source function. After forming the modulus maxima lines (MML) from the wavelet coefficients (Mallat, 1999), the maximum points along these lines are calculated, yielding rough estimates for the scale (= band width) and position of the reflection events. The result of this stage is a set of locations and scales {τ (n) , σ (n) } with n = 1 · · · N , and N is the number of detected maxima, which corresponds to location and scale. These approximated values are subsequently used as initial guesses as part of the nonlinear inversion during the estimation stage. Seismic trace  8  6  4  Amplitude  2  0  -2  -4  -6  -8 0  200  400  Location  600  800  1000  (a) Synthetic seismic signal.  (b) Continuous wavelet transform of signal.  Figure 1: A typical example for the detection of a seismic trace (a) with seven reflection events. (b) The modulus of the continuous wavelet transform with warm colors corresponding to large magnitudes. The local maxima for the wavelet coefficients are used as preliminary estimates for the scale and location of the transitions. The plus signs show modulus maxima.  Partitioning: Given the estimates for the location and scale of the detected events, the trace is segmented into separate events (see Fig. 2). During segmentation, each individual waveform is calculated with s(n) (t) = W[τ (n) , σ (n) ]s(t) with n = 1 · · · N, (3) where τ (n) and σ (n) are the locations and scales of the nth detected waveform, W[.] is the windowing operator centered at τ (n) , and has a support proportional to σ (n) . The output of this procedure are N vectors with ’isolated’ events. Even though this segmentation procedure is somewhat arbitrary, e.g. it depends on a width parameter, we found this method to perform reasonably well for most cases. Sub-wavelength details are not extracted this way and are left to the ensuing estimation stage.  EAGE 69th Conference & Exhibition — London, UK, 11 - 14 June 2007  Segmentation of detected events 8  6  4  Amplitude  2  0  -2  -4  -6  -8 0  200  400  Location  600  800  1000  Figure 2: Segmentation of detected events. The solid lines show events and the dashed lines correspond to the action of windowing function.  Estimation: During the second stage of our method, the segmented events are subjected to a nonlinear descent-driven inversion procedure. To setup this inversion procedure, we first need to refine our mathematical model for the parametrized family of waveforms. Given our choice presented earlier, we define an individual waveform as a fractional derivative of a shifted and scaled Gaussian fθ (t) = Dα ( √  1 2πσ 2  e(t−τ )  2 /2σ 2  ),  (4)  where the scale is denoted by σ, the location by τ , and the fractional order by α. Motivated by the work of Wakin et al. (2005), we can parametrize the family of waveforms as M [θ] = {fθ : θ ∈ Θ} with θ = [σ, τ, α]. The nonlinear estimation procedure for each segmented waveform consists of minimizing e(n) (θ) = s(n) − M [θ] or  2 2  ,  (5)  θ˜(n) = arg min e(n) (θ),  (6)  θ  for θ ∈ Θ with Θ a feasible parameter range. To solve the above minimization problem, a descent method is employed that requires differentiability of the forward model with respect to its parameters (θ). Under that assumption, analytical expressions for these partial derivatives can be derived. With these partial derivatives w.r.t. θi , i = 1 · · · 3, the descent update can be calculated and is given by (n)  Ji with γθi =  ∂fθ , ∂θi  (n)  and Ji  =  ∂e(n) = 2 s(n) − M [θ], γθi , ∂θi  (7)  is the projected estimation error for each parameter. During (n)  each iteration of the method, e(n) (t) is formed, and J(n) = {Ji calculated, followed by the following update 1 θ˜(n),k+1 ←− θ˜(n),k − J(n) 2  : i = 1 · · · 3} is (8)  where k is the number of iterations of the descent method. Our experience using the above scheme have shown that this iterative optimization method provides acceptable solutions to the estimation problem (see Fig. 3).  EAGE 69th Conference & Exhibition — London, UK, 11 - 14 June 2007  Discussion In this paper, we have presented a new characterization method which allows for the estimation of fractional-order discontinuities. These scale attributes may lead in improvement of geological interpretation from the seismic trace. The examples we have presented, indicate that proposed characterization method is leads to accurate results. In addition, our method is well suited for the estimation of the scale exponent attributes from bandwidth limited data. As opposed to wavelet coefficient decay based methods, such as SPICE (Liner et al.,2004), our method does not lead to possibly ambiguous estimates since we do not rely on ‘infinite’ bandwidth. 0.3  Parameter estimation for 8th event, Initial iteration  0.3  0.2  Normalized amplitude  Normalized amplitude  0.2  0.1  0.0  -0.1  -0.2  -0.3  0.1  0.0  -0.1  -0.2  0  200  400  Location  600  800  1000  -0.3  0  Error of estimated events to seismic signal 8  4.5  Fractional Order  2  Amplitude  400  Location  600  800  1000  Error of Estimated attributes to actual values  3.5  4  0  -2  -4  -8 0  200  4.0  6  -6  Parameter estimation for 8th event, Final iteration Windowed event Estimation  Windowed Event Initial guess  3.0  2.5  2.0  1.5  1.0  0.5  Real Signal Estimated signal 200  400  Location  600  800  1000  0.0 0  200  400  Location  600  800  1000  Figure 3: Example of characterization with 10 reflectivity events.(top) Initial and final iteration of parameter estimation for one isolated event, where the actual values, initial guess and estimation are θ = (12.2, 667, 2.93), θinit = (7.81, 668, 0.7), and θ˜ = (12.72, 667, 3.01) respectively. The dashed line shows actual component and the solid line the estimation. (bottom left) Estimated seismic signal is formed by superposition of all characterized events and compared with the original seismic trace. (bottom right) The estimated attributes of events (τ, α) are compared to their actual values. The Blue diamonds show actual parameters whereas red the circles estimated values.  Acknowledgments: This work was in part financially 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 Felix J. Herrmann and was carried out as part of the ChaRM project supported by Chevron.  References [1] Boyd S., and Vandenberghe L. [2004] Convex Optimization, Cambridge University Press [2] Herrmann, F.J., Lyons W.J., and Stark, C. [2001] Seismic facies characterization by monoscale analysis, Geophysical Research Letter, 28(19), 3781-3784. [3] Herrmann, F.J., and Bernab´e Y. [2004] Seismic singularities at upper-mantle phase transitions: a site percolation model. Geophysical Journal International, 159, 949-960 [4] Liner, C., Li, C.-F., Gersztenkorn, A., and Smythe, J. [2004] SPICE: A new general seismic attribute, 72nd Annual International Meeting, SEG, 433436 (Expanded Abstracts). [5] Mallat, S. [1999] A wavelet tour of signal processing. Academic Press. [6] Saggaf, M.M., and Robinson, E. A. [2000] A unified framework for the deconvolution of traces of nonwhite reflectivity. Geophysics, 65(5), 16601676. [7] Wakin, M. B., Donoho, D. L., Choi, H., and Baraniuk R. G. [2005] The Multiscale Structure of NonDifferentiable Image Manifolds. SPIE Wavelets XI, 5914, 413-429, San Diego, (Expanded abstract).  EAGE 69th Conference & Exhibition — London, UK, 11 - 14 June 2007  


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