AUTOMATIC MULTICOMPONENT IMAGE REGISTRATION
by
A thesis submitted to the Faculty and the Board of Trustees of the Colorado School of Mines in partial fulfillment of the requirements for the degree of Master of Science (Geo-physics). Golden, Colorado Date Signed: Stefan Compton Signed: Dr. Dave Hale Thesis Advisor Golden, Colorado Date Signed: Dr. Terence K. Young Professor and Head Department of Geophysics
ABSTRACT
Multicomponent seismic images are composed of different combinations of downgoing and upgoing wavefields. Each wave mode has different propagation velocity and polariza-tion direcpolariza-tion and thus carries unique, direcpolariza-tion-dependent informapolariza-tion about the subsurface. Differences in propagation velocity cause events in converted- or shear-wave images to appear at later times than the compressional-wave image counterpart. Reflectivities are different for each wave mode and therefore, multicomponent images are not related simply by time shifts. These complications historically required that the alignment, also called registration of corresponding image features be done manually, a tedious process. This thesis devel-ops an approach to automatically register multicomponent images using a smooth dynamic warping algorithm that can be accurate with respect to problems unrelated to time shifts. Interval VP/VS ratios can be estimated from derivatives of time shifts that align reflections
in multicomponent images, and these VP/VS ratios may be used to assess the accuracy of the
automatic registration process. To improve the accuracy of estimated time shifts and VP/VS
ratios, we automatically construct a coarse lattice of points located on reflections with high amplitudes, and then estimate time shifts at only those points. By adjusting the coarseness of the lattice, we trade off resolution of changes in VP/VS with increased accuracy in VP/VS
estimates. The result is an efficient, robust, and automatic method for multicomponent image registration.
TABLE OF CONTENTS
ABSTRACT . . . iii
LIST OF FIGURES . . . vi
ACKNOWLEDGMENTS . . . x
CHAPTER 1 INTRODUCTION . . . 1
CHAPTER 2 SMOOTH DYNAMIC WARPING . . . 5
2.1 Introduction . . . 5
2.2 Dynamic time warping . . . 6
2.3 Smooth dynamic time warping . . . 12
2.4 Smooth dynamic image warping . . . 17
CHAPTER 3 MULTICOMPONENT IMAGE REGISTRATION . . . 20
3.1 Introduction . . . 20
3.2 Traveltime relations . . . 20
3.3 Geophysical constraints on time shifts . . . 22
3.4 SDIW versus DIW . . . 24
3.5 Coarse sample locations . . . 27
3.6 Interpolation of time shifts . . . 34
CHAPTER 4 SIMULTANEOUS MULTICOMPONENT REGISTRATION . . . 37
4.1 Introduction . . . 37
4.2 Traveltime relations . . . 38
4.4 Reducing computational cost . . . 43
4.5 Independent versus simultaneous registration . . . 46
CHAPTER 5 CONCLUSIONS . . . 52
LIST OF FIGURES
Figure 2.1 Synthetic PP trace in PP time and synthetic PS trace in PS time. . . 9 Figure 2.2 Alignment errors computed from synthetic traces in Figure 2.1. The
correct time-shift function is shown as the dashed white curves, and the correct VP/VS ratios are shown as the dashed red curves. The solid white
curve is the time-shift function computed using Algorithm 1, with the
corresponding VP/VS ratios shown in black. . . 10
Figure 2.3 Smoothed version of the computed time shifts in Figure 2.2, shown as the solid white curve in the top panel, and corresponding VP/VS ratios
shown in black in the bottom panel. . . 11 Figure 2.4 Zoomed views of alignment errors computed from synthetic traces in
Figure 2.1 and time shifts computed using Algorithm 1 (a) and Algorithm 2 (b) with a coarse sampling interval h = 10∆t. The computational stencils for each algorithm are shown in red. Correct
values are the dashed curves and computed values are the solid curves. . . 15 Figure 2.5 Time shifts computed from smooth dynamic warping (Algorithm 2), for
different coarse sampling intervals h. White dots show the coarse sample locations where time shifts are estimated. For h = ∆t and h = 10∆t dots are not shown because they would hide the time-shift curves completely. Dashed curves represent correct values and solid curves represent
computed values. The solid curves are obtained by piecewise-linear interpolation of computed values. As the coarse sample interval increases to h = 200∆t, time shifts and VP/VS ratios are more accurately
estimated. When the coarse sample interval is too large (h = 400∆t)
estimates are less accurate. . . 16 Figure 3.1 Time-migrated PP and PS images before image registration. Note that
reflectors are approximately aligned, but the PP and PS images are
Figure 3.2 Alignment errors and time shifts computed from a single pair of traces from the PP and PS images shown in Figure 3.1. Only the alignment errors for times 0 to 1 second are shown. Overlaid curves are the 1D time shifts extracted from the 3D time shifts u(x, y, tpp) at the
corresponding trace location. The solid white curve shows rough time shifts computed from DIW and the dashed blue curve shows these shifts after smoothing. The dashed yellow curve shows the time shifts
computed from SDIW with a minimum subsampling interval h = 100 ms. . 25 Figure 3.3 Slices through three 3D volumes with the same PP time window shown
in the alignment errors of Figure 3.2. The location for those alignment errors corresponds to the vertical black lines in the inline and crossline sections. The strong events at 0.7 s are well aligned in each warped image, but the time slice through the event at 0.38 s shows that time shifts computed using SDIW are more accurate. The differences in time
shifts can be seen clearly in Figure 3.2. . . 26 Figure 3.4 Interval VP/VS ratios estimated from smoothed DIW time shifts and
from SDIW time shifts. The VP/VS ratios from smoothed DIW time
shifts exhibit fine detail, but those smoothed shifts are less accurate, as
shown in Figure 3.2 and Figure 3.3, than shifts estimated using SDIW. . . 27 Figure 3.5 Alignment errors (a) computed from a single trace from the PP and PS
images shown in Figure 3.1, with time shifts (b) computed from 1D smooth dynamic time warping using two different subsampling grids. A nearly-uniform coarse grid was used for the blue curve, with subsampling locations represented by blue dots, and a reflector-aligned grid was used for the yellow curve with subsampling locations represented by yellow squares. The PP trace (c) overlaid with the envelope of the trace shows that strong PP events correlate with strong vertical features in the
alignment errors above. . . 29 Figure 3.6 The PP image (a) and time shifts (b) used to flatten the image so that a
constant time corresponds to a single seismic horizon. The coarse lattice of points on the flattened image (c) is shown by the white dots. This lattice can be mapped back to the original (unflattened) coordinate space to obtain a coarse lattice of points that is aligned with dipping reflectors (d). We use this lattice in SDIW to compute time shifts at
Figure 3.7 White dots indicate coarse lattice points used to compute time shifts with SDIW. Both nearly-uniform and reflector-aligned grids have minimum intervals of 500 m laterally. The nearly-uniform points are a minimum of 160 ms apart in time, while reflector-aligned points are a minimum of 100 ms apart but are aligned with strong reflectors in the image. Points in the nearly-uniform grid cross strong reflectors
horizontally, resulting in a noticeable stair step effect in the
corresponding VP/VS ratios. . . 33
Figure 3.8 Interval VP/VS estimates resulting from two different interpolation
methods for time shifts, piecewise cubic with monotonic splines and
piecewise linear. . . 34 Figure 3.9 Interpolated time shifts and interval VP/VS extracted from 3D SDIW
shift and VP/VS volumes. The three different curves represent three
different methods of interpolating time shifts: linear (black), cubic
monotonic spline (red), and cubic spline (blue). . . 36 Figure 4.1 3D alignment errors for simultaneous warping of synthetic traces shown
in Figure 4.3. The computed solution is shown as the white curve. . . 44 Figure 4.2 Alignment errors for synthetic PP and PS1 traces with no initial
squeezing of the PS trace (a). Yellow triangles indicate large areas where computation time is wasted. Alignment errors for the same traces, but after squeezing the PS1 trace by a constant factor c = 1.7 (b). The
dashed white curves show the known time-shift function. . . 45 Figure 4.3 Synthetic PP, PS1, and PS2 traces with a signal-to-noise (S/N) ratio of
1:1. . . 47 Figure 4.4 Time shifts estimated by individually warping three pairs of synthetic
traces shown in Figure 4.3, PP-PS1 (red), PP-PS2 (blue), PS1-PS2 (green). Correct time shifts are shown as dashed white curves. Results for noise-free synthetic traces are shown on the left and results for synthetic traces with S/N = 1:1 are shown on the right. The bottom panels show errors in γ(s) estimated from equations 4.15 (magenta), 4.16
Figure 4.5 Time shifts estimated by simultaneously warping the three synthetic traces shown in Figure 4.3. PP-PS1 (red) and PS1-PS2 (green) are directly computed from the simultaneous warping, PP-PS2 (blue) is computed from equation 4.7. Correct time shifts are shown as dashed white curves. Results for noise-free synthetic traces are shown on the left and results for synthetic traces with S/N = 1:1 are shown on the right. The bottom panels show errors in γ(s) estimated from equations 4.15
ACKNOWLEDGMENTS
I would like to thank my advisor, Dr. Dave Hale, for his guidance and support throughout my graduate studies. His valuable insights and enthusiasm were important factors in this research. Thanks to Jim Gaiser and Lee Bell for their helpful discussions and thanks to Geokinetics Inc. and Geophysical Pursuit, Inc. for providing the seismic data that was invaluable to this research. I would also like to acknowledge my colleagues in CWP for making my time here enjoyable.
Special thanks to my family, especially my wife Mara for her patience and support throughout this process.
CHAPTER 1 INTRODUCTION
Multicomponent seismic technology is a valuable tool in exploration geophysics. Histor-ically, much seismic analysis focused on compressional-wave (P-wave) seismic data which provide only partial information about the subsurface. Shear-wave (S-wave) data provide additional information about the same subsurface, and each wave mode contains unique in-formation about elastic constants, grain cementation, pore geometry, anisotropy axes, and lateral variations in rock and fluid types (Hardage et al., 2011).
Data acquisition options for multicomponent seismic surveys correspond to the types of wave modes recorded. Common acquisition options are referred to as 3C, 4C, 6C, or 9C, depending on the orientations of displacement generated by seismic sources and the orientations of displacement recorded by receivers. All possible wave modes are contained in 9C data, collected when a source generates three orthogonal displacements recorded by three orthogonal sensors.
In the simplest (3C) multicomponent survey the downgoing wavefield is generated by a P-wave source and the upgoing P-wavefield is recorded by a receiver containing three orthogonally oriented sensors. For an isotropic earth, 3C seismic surveys would record predominantly PP and PSV modes. The acronyms PP and PSV denote a downgoing P-wave, and an upgoing P- or SV-wave. The downgoing P-wave is reflected at a subsurface interface and travels upward as both a P-wave and mode converted SV-wave. The P- and SV-wave have orthogonal polarizations within the vertical plane containing the source and receiver. The SV mode is distinguished from the shear-wave SH mode, with polarization orthogonal to the vertical plane containing the source and receiver. The SH modes have significantly less energy than SV modes in 3C surveys. Therefore, the more complete mode description, PSV, is often abbreviated as simply PS. Because converted-wave surveys do not involve shear-wave
sources, converted-wave seismic data can also be acquired in the marine environment with the use of 4C (three orthogonal sensors plus a hydrophone) sensors on the seafloor.
The goal of converted-wave seismic data processing is to preserve and exploit the infor-mation provided by both PP and PS reflections. Joint PP and PS analysis has applications in structural imaging, lithology estimation, anisotropy analysis, fluid description, and reser-voir monitoring (Stewart et al., 2003). A typical first step is to generate PP and PS seismic images. Because traveltimes differ for seismic reflectors in PP images and corresponding reflectors in PS images, a crucial first step in this analysis is the alignment, or registration, of these reflectors.
The vertical shifts that align PP and PS images can be used to quantify the ratio (VP/VS)
of P- and S-wave velocities, an important subsurface property (Stewart et al., 2003). How-ever, PP and PS image registration is complicated because the PS image is not simply a time-shifted version of the PP image. PP and PS reflectivities differ, leading to amplitude and phase differences between the two images that can be significant. In fact, it is possible for a reflector in one image to not be apparent in the other image. For these reasons, image registration is often done by manually selecting corresponding events in each image. Never-theless, many techniques have been developed to automate the registration process and to estimate VP/VS ratios (Gaiser, 1996; Fomel and Backus, 2003; Nickel and Sonneland, 2004;
Yuan et al., 2008; Liang and Hale, 2012).
Liang and Hale (2012) used dynamic image warping (Hale, 2013) for this purpose. Dy-namic image warping computes vertical time shifts that align events in PP and PS images; these shifts are a globally optimal solution of a minimization problem with many possible local minima. The dynamic warping algorithm enforces strict bounds on both time shifts and derivatives of those time shifts. For these reasons, the dynamic warping algorithm may be preferable when compared to other methods.
For instance, Fomel and Backus (2003) describe a method that solves a similar minimiza-tion problem, and uses regularizaminimiza-tion terms that enhance, but do not strictly constrain, the
smoothness of estimated time shifts. This method also requires a good initial guess of the warping function, to avoid convergence to local minima. Constraints in dynamic warping are strict, hard bounds, and the solution is globally optimal, even though no starting model is required.
Yuan et al. (2008) use a global search method, simulated annealing, which solves problems with local maxima and allows for a poor initial model. However, the maximized cost function is the normalized crosscorrelation of PP and warped PS images. Hale (2013) shows that the success of local crosscorrelation methods depends on how rapidly shifts vary within windows used to compute correlation coefficients. If shifts vary rapidly (as for PP and PS images, especially near the surface) then a suitable correlation window may not exist. Dynamic warping does not require windows and is accurate even when shifts vary rapidly.
A common goal of the previously cited methods is to obtain, in conjunction with regis-tered PP and PS images, a high-resolution VP/VS volume. The term high-resolution implies
fine detail, but what level of detail can be extracted from traveltimes between PP and PS images, in which differences in noise and reflection waveforms are significant? To obtain high-resolution information, more expensive and thorough techniques are required, such as joint PP-PS inversion of seismic amplitudes. For PP and PS image registration, time shifts are related to time-variant averages of VP/VS ratios, and thus, are often smoothly varying.
The dynamic warping algorithm can be modified to exploit this inherent smoothness and, by sacrificing temporal resolution of time shifts, more accurately estimate the change in time shift.
Chapter 2 develops a smooth dynamic warping algorithm, that is quite general and has applications beyond multicomponent image registration. Chapter 3 discusses the practical application of smooth dynamic warping to registration of multicomponent images. Smooth dynamic warping is shown to be a robust and accurate method for automatic image regis-tration through an application to 3D PP and PS time-migrated images. Chapter 4 extends dynamic warping beyond the alignment of two images, for in the presence of anisotropy, the
PS mode can split into PS1 (fast-S) and PS2 (slow-S) modes. With three related seismograms (PP, PS1, and PS2), the algorithm and application developed in Chapter 2 and Chapter 3, respectively, could be applied independently to each pair of seismograms. However, time shifts estimated with that approach would likely be inconsistent. Therefore, in Chapter 4, the dynamic warping optimization problem is reformulated to simultaneously align PP, PS1, and PS2 seismograms with consistent time shifts.
CHAPTER 2
SMOOTH DYNAMIC WARPING
2.1 Introduction
Dynamic image warping is a robust method for computing shifts that align features in one image with corresponding features in another image. The computed shifts are a globally optimal solution to a non-linear error minimization problem that satisfies linear inequality constraints. For many problems, these linear inequality constraints can be directly related to physical parameters. Dynamic warping is a dynamic programming algorithm that solves a collection of nested subproblems to find the globally optimal solution. In the context of aligning features in seismic images, the error function may have many local minima, and thus the ability to find the globally optimal solution is a key advantage of this method.
A dynamic programming algorithm for spoken word recognition was described by Sakoe and Chiba (1978). They were the first to introduce strict constraints on the rate of change of time shifts. Their algorithm has become known as dynamic time warping (DTW) (e.g., M¨uller (2007) chapter 4). Anderson and Gaby (1983) propose applications of DTW to estimate time shifts in geophysical time series, and refer to this process as “dynamic waveform matching”. The use of DTW in geophysics is not new and has many potential applications. Time shifts that align seismic traces can increase or decrease rapidly with time, but the time derivative of shifts often varies slowly. We therefore introduce a modified version of the dynamic warping algorithm, called smooth dynamic warping (Hale and Compton, 2013), to exploit this inherent smoothness and estimate shifts at coarsely sampled times. The resulting decrease in temporal resolution yields more accurate estimates of the change in shift between the coarsely sampled times. This resolution and accuracy tradeoff is especially significant when differences between traces are not related to time shifts alone, but include differences in noise and reflection waveforms.
The application of DTW to 1D signals is useful for certain geophysical problems, but we must often estimate shifts between multidimensional images. DTW can produce highly in-accurate results when applied independently to pairs of corresponding traces in such images. Depending on the geometry of the geophysical problem, it is often reasonable to assume that neighboring traces should have similar time shifts, but there are no constraints for the rate at which shifts vary laterally if DTW is applied independently to pairs of traces. The exact extension of DTW to higher dimensions, with strict lateral constraints on shifts, has been shown to be computationally infeasible (Keysers and Unger, 2003). The time to compute an exact solution grows exponentially with image size. Therefore, an approximate solution to this multidimensional problem was proposed by Mottl et al. (2002) and extended by Hale (2013). These solutions are referred to as dynamic image warping (DIW) algorithms. The smooth dynamic time warping algorithm can be extended to multidimensional images in the same fashion.
2.2 Dynamic time warping
For two sampled sequences f and g, dynamic time warping will find a sequence of time shifts u[0 : ni− 1] ≡ {u[0], u[1], ..., u[ni− 1]} so that
f [i] ≈ g[i + u[i]], i = 0, 1, ..., ni− 1. (2.1)
The time shifts u[i] are the solution to the optimization problem:
u[0 : ni− 1] ≡ arg min l[0:ni−1] ni−1 X i=0 e[i, l[i]], (2.2) subject to constraints
ul ≤ u[i] ≤ uu, rl ≤ u[i] − u[i − 1] ≤ ru, (2.3)
e[i, l] ≡ (f [i] − g[i + l])2. (2.4)
From the inequality constraints 2.3, time shifts u[i] are bounded by specified ul and uu;
and the rate of change of those time shifts (time strains) are bounded by specified rl and
ru. For some problems these time strains can be directly related to geophysical parameters.
This is indeed the case for PP-PS registration, as discussed in Chapter 3.
For the 1D sequences f and g we first compute alignment errors e[i, l] for all sample indices i and lag indices l (equation 2.4). Alternative measures of alignment could be used, but here we use the squared difference between amplitudes in f and g. The lag l satisfies the constraint ul ≤ l ≤ uu. In practice, we use the bounds 0 ≤ l ≤ nl− 1, where nl= uu− ul+ 1
is the number of lags for which alignment errors are computed. This change simplifies array indexing, but the lower bound ul must now be added to l in equation 2.4.
The program to compute the time shifts is simple to implement, and the pseudocode for this program is shown below in Algorithm 1. The inputs to the procedure FindShifts are the bounds on strain rl and ru, a 2D array of alignment errors e[i, l], and an array u in which
to store the computed shifts.
In lines 2-3 a 2D array of accumulated errors d[i, l] is initialized. At sample index i = 0, the accumulated errors are equal to the alignment errors.
Lines 4-15 recursively compute accumulated errors at sample index i from accumulated errors at i − 1. Lines 7 and 8 compute lower and upper slopes that satisfy strain con-straints (the notation drle denotes the smallest integer not less than rl, and bruc denotes
the largest integer not greater than ru). The accumulation is nonlinear because all possible
slopes from ql to qu are queried to find the error minimizing move, ml (lines 9-13). The
error minimizing moves are stored in the 2D array m[i, l], to simplify the final step of the procedure. It is important to note that the slope q between samples i and i − 1 is an integer and therefore coarsely quantized, so that the computed time shifts u will be rough. This
shortcoming will be revisited later.
Algorithm 1 Find shifts u[i]
1: procedure FindShifts(rl, ru, e, u) 2: for l = 0 to nl− 1 . initialize 3: d[0, l] = e[0, l] 4: for i = 1 to ni− 1 . accumulate 5: for l = 0 to nl− 1 6: dl= ∞ 7: ql = max(drle, l − nl+ 1) 8: qu = min(bruc, l) 9: for q = ql to qu 10: dq = d[i − 1, l − q] + e[i, l] 11: if dq < dl 12: dl = dq 13: ml= q 14: d[i, l] = dl 15: m[i, l] = ml 16: i = ni− 1 . minimize d 17: di = ∞ 18: for l = 0 to nl− 1 19: if d[i, l] < di 20: di = d[i, l] 21: u[i] = l 22: while i > 0 . backtrack
23: u[i − 1] = u[i] − m[i, u[i]] 24: i = i − 1
When the accumulation step is complete, lines 16-21 find the lag l that minimizes d[ni−
1, l]. This lag is the optimal shift u[ni− 1], and the accumulated error d[ni− 1, u[ni− 1]] is
the minimized sum in equation 2.2.
Finally, lines 22-24 compute the optimal shifts for all indices i. After each optimal shift at index i is found and stored in u[i], we can use the error minimizing moves stored in array m[i, l] to compute the optimal shift u[i − 1] at index i − 1. This last step simply follows the optimal path backwards from index i = ni− 1 to index i = 0.
As an example, I apply this algorithm to a synthetic PP trace and a synthetic PS trace shown in Figure 2.1. These traces are plotted on two different time scales as seen in the
Figure 2.1: Synthetic PP trace in PP time and synthetic PS trace in PS time.
time axis labels of Figure 2.1. I denote the PP trace as f and the PS trace as g. The trace g was generated by convolving a random reflectivity series with a Ricker wavelet. The random reflectivity series used to construct g was warped by a time varying shift function and then convolved with a Ricker wavelet to form the trace f . All events in g appear at times not less than the corresponding times in f , because S-wave velocities cannot exceed P-wave velocities. Finally, different bandlimited random noise sequences were added to f and g.
The top panels of Figure 2.2 show alignment errors computed from the synthetic traces in Figure 2.1. The dashed white curve in each top panel is the correct time-shift function. The dashed red curve in each bottom panel is the corresponding correct VP/VS function:
VP VS (tpp) = 1 + 2 du(tpp) dtpp . (2.5)
This relationship will be derived in Chapter 3, but for this first synthetic example it is important only to observe that the VP/VS ratios are related to the derivative of the time
shifts u. This relationship is exhibited in Figure 2.2 by the decrease in VP/VS values as the
slope of the time-shift function decreases. Time shifts computed using Algorithm 1 are shown as the solid white curve. The VP/VS ratios computed from those time shifts (equation 2.5)
are shown in black.
Figure 2.2: Alignment errors computed from synthetic traces in Figure 2.1. The correct time-shift function is shown as the dashed white curves, and the correct VP/VS ratios are
shown as the dashed red curves. The solid white curve is the time-shift function computed using Algorithm 1, with the corresponding VP/VS ratios shown in black.
Recall that Algorithm 1 returns a sequence of integer shifts subject to inequality con-straints 2.3. The phrase “integer shifts” may be confusing here since Figure 2.2 shows time shifts with units of seconds that are clearly not integer values. However, remember that the 2D array of alignment errors, e[i, l], is a sampled function of PP time and time shift, with ni PP-time samples and nl time-shift samples. The uniform time sampling interval in this
example is ∆t = 2 ms, for both PP time and time shift, but Algorithm 1 operates with only integer sample indices. Thus the phrase “integer shifts” implies that computed time shifts can be only integer multiples of the time sampling interval ∆t. In this example, the bounds in the constraints 2.3 are the integers ul = 0, uu = 1331, rl = 0 and ru = 1. Therefore,
for each sample index i, the computed time shift can only increase by one or stay the same. In other words, there are only two possible slopes, 0 and 1. From equation 2.5, this coarse sampling of slope is the reason why estimated VP/VS values in Figure 2.2 oscillate wildly
between only two values, 1 and 3.
Figure 2.3: Smoothed version of the computed time shifts in Figure 2.2, shown as the solid white curve in the top panel, and corresponding VP/VS ratios shown in black in the bottom
panel.
Of course, in a practical application the estimated VP/VS ratios shown in Figure 2.2 are
worthless. Algorithm 1’s estimation of time shifts at every PP-time sample limits the possible values of slope and, hence, VP/VS ratios. To obtain more values, the integer shifts can be
smoothed, better approximating fractional changes in shifts. Figure 2.3 shows the result of smoothing the solid white time shifts from Figure 2.2. The smoothed estimated time shifts
better approximate the correct shifts at PP times greater than 1.2 seconds, so that the black curve in the bottom panel of Figure 2.3 more accurately approximates the correct VP/VS
ratios. However, the shifts estimated at earlier times remain significantly different from the correct shifts. This example demonstrates that we cannot obtain the correct shifts by simply smoothing incorrect shifts.
Another way to obtain more slope values, while still satisfying the slope constraints 0 ≤ u[i]−u[i−1] ≤ 1, would be to sample the lag axis more finely. That is, sample alignment errors at some fraction of the time sampling interval to obtain more possible changes in slope between each time sample. Hale (2013) makes this observation, but notes that this approach can greatly increase computation time and memory, especially when estimating shifts for multidimensional image warping.
In fact, there is another reason why finely sampling the lag axis is not a good solution. This approach would imply that between two consecutive time samples we could accurately estimate a change in time shift equal to a fraction of the sampling interval. When sequences f and g are not simply shifted versions of one another, it is unrealistic to think such minute changes in shift can be detected. The algorithm will find the globally optimal solution, but are the estimated changes in shift meaningful at that fine scale? For that matter, is it even reasonable to estimate time shifts between each consecutive time sample? For typical seismic data, the time sampling interval can be several times smaller than the dominant period of the data. This observation leads us to a modified version of the dynamic warping algorithm. 2.3 Smooth dynamic time warping
The smooth dynamic warping algorithm is a simple modification of Algorithm 1. Instead of estimating shifts at every time sample, the smooth warping algorithm estimates shifts at only coarsely sampled locations. This approach makes smooth dynamic warping more robust in the presence of noise and other differences between the sampled sequences. The coarse sampling sacrifices temporal resolution of time shifts to more accurately estimate changes in those time shifts. Such resolution and accuracy tradeoffs are common in signal processing.
Algorithm 2 shows how Algorithm 1 is modified to implement smooth dynamic warping. The algorithms are similar, but there are several key differences. The first difference is that Algorithm 2 computes shifts ui[j] ≡ u[i[j]] for only a subset of sample indices i[0 : nj− 1] ≡
i[0], i[1], ..., i[nj− 1]. Therefore, note that an array of nj subsampled indices i[0 : nj− 1] are
input to the procedure. Shifts u[i] for all sample indices i can be later interpolated from the subsampled shifts ui[j].
Algorithm 2 Find subsampled shifts ui[j] ≡ u[i[j]]
1: procedure FindShiftsI(rl, ru, e, i, ui) 2: for l = 0 to nl− 1 . initialize 3: d[0, l] = e[0, l] 4: for j = 1 to nj − 1 . accumulate 5: h = i[j] − i[j − 1] 6: for l = 0 to nl− 1 7: dl = ∞ 8: ql = max(dhrle, l − nl+ 1) 9: qu = min(bhruc, l) 10: for q = ql to qu 11: dq= d[j − 1, l − q] 12: for p = 0 to h − 1 13: dq = dq+ e[i[j] − p, l − pq/h] 14: if dq < dl 15: dl = dq 16: ml = q 17: d[j, l] = dl 18: m[j, l] = ml 19: j = nj − 1 . minimize d 20: dj = ∞ 21: for l = 0 to nl− 1 22: if d[j, l] < dj 23: dj = d[j, l] 24: ui[j] = l 25: while j > 0 . backtrack 26: ui[j − 1] = ui[j] − m[j, ui[j]] 27: j = j − 1
The next significant difference is in the computation of slope bounds on lines 8 and 9. The integers q and q depend on h, the difference between subsampled indices i. Thus, this
algorithm satisfies the same slope constraints as Algorithm 1, but with more slopes possible between the bounds ql and qu. Lines 11-13 replace line 10 in Algorithm 1. In Algorithm 2 we
find the q with minimum accumulated error, where q defines a linear trajectory from index i[j − 1] to i[j]. Finding the minimum accumulated error requires summing the alignment errors along each linear trajectory (lines 12 and 13) for each slope q between ql and qu.
Note that the term l − pq/h in line 13 may not have integer value, and thus interpolation of alignment errors e[i, l] at non-integer lags is required. Finally, note that if h = 1, then Algorithm 2 is equivalent to Algorithm 1.
Figure 2.4a shows a zoomed view of the correct and computed time shifts in Figure 2.2. Again, the dashed white curve is the correct time-shift function and the solid white curve is computed with Algorithm 1. The computational stencil is shown in red, and as previously discussed, has only two possible values of slope, 0 and 1. Figure 2.4b shows the smooth-dynamic-warping stencil for a coarse interval h = 10∆t, where ∆t is the time sampling interval (2 ms). Given this stencil, smooth dynamic warping estimates time shifts at only coarsely sampled locations approximately 10∆t apart. These locations are computed to be at least h samples apart, but some intervals may be larger to ensure that the first and last samples are included. The solid white curve was obtained by piecewise-linear interpolation of time shifts ui[j].
The smooth dynamic warping stencil shown in Figure 2.4b satisfies the same constraints as the stencil in Figure 2.4a, but has 11 possible values of slope. Therefore, smooth dynamic warping can more accurately estimate the change in time shift between coarse samples. The improvement is apparent in the bottom panel of Figure 2.4b, which shows estimated VP/VS
ratios in black. These values are still incorrect, but what happens if we increase the coarse sampling interval h even more?
Figure 2.5 shows a comparison for a range of possible values of h. In all Figures 2.5a-f, the dashed white curve shows the correct time shifts, the solid white curve shows the computed time shifts, the dashed red curve shows the correct VP/VS ratios, and the black
Figure 2.4: Zoomed views of alignment errors computed from synthetic traces in Figure 2.1 and time shifts computed using Algorithm 1 (a) and Algorithm 2 (b) with a coarse sampling interval h = 10∆t. The computational stencils for each algorithm are shown in red. Correct values are the dashed curves and computed values are the solid curves.
curve shows the computed VP/VS ratios. Figure 2.5a is the same as Figure 2.2b, but for
coarse sampling interval h equal to the time sampling interval ∆t, and Figure 2.5b shows the zoomed out version of Figure 2.4b that was already discussed. Figures 2.5c-f show results for larger intervals h. These figures also display white dots that indicate locations along the shift curve where time shift was actually computed by the smooth dynamic warping algorithm. The curves were obtained by piecewise-linear interpolation of those computed time shifts. In Figure 2.5a and Figure 2.5b these dots are not shown because they would be so dense that the time-shift curves would be invisible.
Figure 2.5 clearly shows that, as the interval h is increased, we are better able to recover the correct time shifts. Also, notice the improved accuracy of the estimated VP/VS ratios.
Shifts in Figure 2.5f are estimated at only 5 locations, and while this interval h is perhaps too coarse, the estimates are more accurate than those for intervals h < 200∆t. This example shows that, by increasing h, and thereby decreasing the temporal resolution of shifts, we are able to more accurately estimate the change in shift.
Figure 2.5: Time shifts computed from smooth dynamic warping (Algorithm 2), for different coarse sampling intervals h. White dots show the coarse sample locations where time shifts are estimated. For h = ∆t and h = 10∆t dots are not shown because they would hide the time-shift curves completely. Dashed curves represent correct values and solid curves represent computed values. The solid curves are obtained by piecewise-linear interpolation of computed values. As the coarse sample interval increases to h = 200∆t, time shifts and VP/VS ratios are more accurately estimated. When the coarse sample interval is too large
2.4 Smooth dynamic image warping
Hale (2013) describes how dynamic time warping can be extended to multidimensional image warping. This dynamic image warping (DIW) algorithm still estimates only vertical shifts. In fact, the last step of DIW is an application of DTW, via Algorithm 1, for each pair of traces in two multidimensional images. However, prior to using Algorithm 1, alignment errors are first smoothed vertically and horizontally so that strains (changes in shifts) are limited in each image dimension. Without this smoothing step, time-shift estimates can vary wildly from trace to trace where the two images are not related simply by time shifts.
Let f and g denote two 3D images with ni2× ni3 traces and ni1 samples per trace. With
DIW we can find vertical shifts u[i1, i2, i3] such that
f [i1, i2, i3] ≈ g[i1, i2, i3+ u[i1, i2, i3]]. (2.6)
The shifts u[i1, i2, i3] are constrained by
ul ≤ u[i1, i2, i3] ≤ uu, (2.7)
and
rl1 ≤ u[i1, i2, i3] − u[i1− 1, i2, i3] ≤ ru1, (2.8)
rl2 ≤ u[i1, i2, i3] − u[i1, i2− 1, i3] ≤ ru2, (2.9)
rl3 ≤ u[i1, i2, i3] − u[i1, i2, i3− 1] ≤ ru3. (2.10)
The bounds rl1 and ru1 are analogous to the slope constraints in equation 2.3, and these
control the change in vertical shift within each trace in the 3D image. The bounds rl2, ru2,
rl3, and ru3 constrain the rates at which the vertical shifts u[i1, i2, i3] vary laterally. These
slope constraints are used in the vertical and horizontal smoothing of all alignment errors computed from f and g.
The key to this smoothing step lies in recognizing that lines 4-14 in Algorithm 1 perform a recursive non-linear smoothing of alignment errors e[i, l]. For 3D images f and g,
align-ment errors are computed for all pairs of traces, resulting in a 4D array of alignalign-ment errors e[i1, i2, i3, l]. 2D alignment error subsets, e[ik, l], can be extracted from this 4D array in each
k dimension, where k = 1, 2, 3. To smooth alignment errors in the kth dimension, apply the
recursive non-linear smoothing in both forward and reverse directions. The forward direction is that in which errors are accumulated in order of increasing index ik, as is the case in line
4 of Algorithm 1. The reverse direction means that errors are accumulated in decreasing order of index ik, that is, from nik − 1 to 0. The sum
˜
e[ik, l] = ˜ef[ik, l] + ˜er[ik, l] − e[ik, l] (2.11)
is an array of alignment errors smoothed in the kth dimension. The 2D array ˜ef[ik, l] is the
result of smoothing e[ik, l] in the forward direction, and the 2D array ˜er[ik, l] is the result
of smoothing e[ik, l] in the reverse direction. Subtraction of the original alignment errors
e[ik, l] is necessary so that this error is not included twice in the sum. This smoothing is
first applied in the vertical direction (k = 1) and then in each horizontal direction (k = 2, 3). Smoothing in each dimension is constrained by the corresponding rlk and ruk. The final step
of DIW applies Algorithm 1 to ˜e[i1, i2, i3, l] for all traces indexed by i2 and i3.
A similar method applies to extend smooth dynamic time warping (Algorithm 2) to smooth dynamic image warping (SDIW). In Algorithm 2 it is lines 4-17 that perform a recursive non-linear smoothing as well as a subsampling of alignment errors. Alignment errors are summed along piecewise-linear paths (Figure 2.4b) and stored for subsampled indices ik[jk]. Therefore, the smoothed and subsampled alignment errors in dimension k can
be written
˜
e[jk, l] = ˜ef[jk, l] + ˜er[jk, l] − e[ik[jk], l]. (2.12)
The 2D array ˜ef[jk, l] is the result of smoothing and subsampling e[ik[jk], l] in the forward
direction, and the 2D array ˜er[jk, l] is the result of smoothing and subsampling e[ik[jk], l] in
(k = 1) and then in each horizontal direction (k = 2, 3). Smoothing in each dimension is constrained by the corresponding bounds on slope rlk and ruk, but the subsampling interval
hk = ik[jk] − ik[jk − 1] now allows slopes to be more accurately estimated between coarse
samples in the kthdimension. The final step applies a simple dynamic time warping (with no
further subsampling) to ˜e[j1, j2, j3, l] for all traces indexed by j2 and j3. This last dynamic
time warping is similar to Algorithm 1, but accounts for the distance between coarse samples, hk, when computing lower and upper slopes ql and qu (lines 7 and 8). The result is an array
of coarsely sampled vertical shifts u[j1, j2, j3] that can be interpolated to obtain the desired
shifts u[i1, i2, i3].
The subsampling for multidimensional images yields significant savings in computer mem-ory. The 4D array of alignment errors computed in DIW requires memory proportional to ni1ni2ni3nl; and, as image dimensions grow, this array can quickly become too large to fit in
computer memory. SDIW must also compute alignment errors for all pairs of traces in the two images f and g, but these alignment errors are, for each pair of traces, smoothed and subsampled vertically (k = 1).
Consider for example the alignment errors e[i1, l] in Figure 2.5e, where ni1 = 1941. With
subsampling interval h1 = 200∆t, nj1 = 10 and the array of smoothed alignment errors
˜
e[j1, l] is approximately 200 times smaller. For multidimensional images this smoothing and
subsampling is performed for all ni2×ni3 traces. Before smoothing horizontally (k = 2, 3), the
number of alignment errors is now proportional to nj1ni2ni3nl. Smoothing and subsampling in
horizontal directions reduces this number further to nj1nj2nj3nl, where nj2 and nj3 depend on
horizontal subsampling intervals h2 and h3, respectively. When hk = 1, SDIW is equivalent
to DIW. Increasing h1 allows more accurate estimation of changes in shifts between vertical
coarse samples, and the same is true for h2 and h3. Examples of dynamic image warping are
CHAPTER 3
MULTICOMPONENT IMAGE REGISTRATION
3.1 Introduction
We can use dynamic image warping to compute time shifts that align reflectors in a PP image with corresponding reflectors in a PS image. This alignment process, known as image registration, is complicated by the fact that a PS image is not simply a time shifted version of a PP image due to differences in noise and reflection waveforms. Because of these difference, we use the algorithm described in Chapter 2, smooth dynamic image warping (Hale and Compton, 2013), to estimate image-registration time shifts. This method can more accurately estimate changes in time shifts which are directly related to interval VP/VS
ratios, an important subsurface property. In this chapter I discuss the application of smooth dynamic image warping to 3D time-migrated PP and PS images shown in Figure 3.1. These data are licensed and provided courtesy of Geokinetics Inc. and Geophysical Pursuit, Inc. 3.2 Traveltime relations
Let tppdenote the PP traveltime for a reflector apparent in a PP image, and let tpsdenote
the PS traveltime for the same reflector in the corresponding PS image. For each tpp, we
have a corresponding tps defined by the function
tps(tpp) = tpp+ u(tpp), (3.1)
where u(tpp) denotes the time shift, which is a function of tpp.
The time shift u(tpp) is related to the ratio VP/VS of P- and S-wave velocities. For an
infinitesimal interval dz at depth z, the P- and S-wave velocities are VP = 2 dz dtpp , VS = 2 dz dtss , (3.2)
where tssdenotes the SS time corresponding to the PP time tpp. The factors of 2 are necessary
Figure 3.1: Time-migrated PP and PS images before image registration. Note that reflectors are approximately aligned, but the PP and PS images are plotted on two different time scales.
the times required for P- and S-waves to travel twice through the depth interval dz. From these definitions of interval velocities we define the interval VP/VS ratio by
γi(tpp) ≡ VP VS = dtss dtpp . (3.3)
The time tss is that at which we would observe a reflector in an SS image, if we had one.
If we instead have a converted-wave PS image, we must express the interval VP/VS ratio
γi(tpp) in terms of the time tps(tpp) or, equivalently, the time-shift function u(tpp). To do
this, we recall that PS time is the average of PP and SS times tps(tpp) =
tpp+ tss(tpp)
2 . (3.4)
Substituting this expression into equation 3.3, we obtain γi(tpp) = 2
dtps
dtpp
− 1. (3.5)
γi(tpp) = 1 + 2
du(tpp)
dtpp
. (3.6)
3.3 Geophysical constraints on time shifts
The linear inequality constraints 2.3 define lower and upper bounds, ul and uu,
respec-tively, on the time shifts u(tpp). To specify a lower bound ul, we may use the fact that
VP ≥ VS, which implies that tps ≥ tpp and that all time shifts u(tpp) must be non-negative,
so that the lower bound is ul = 0. Indeed, at the surface, where tpp = tps = 0, we should
theoretically find that u(0) = 0. In practice, we might relax this lower bound to permit negative shifts, because uncertainties in near-surface variations of VP and VS may lead to
in-consistencies in datums and statics corrections used to obtain PP and PS images. However, near the surface we are also likely to find ratios VP/VS 1, so that shifts u(tpp) will increase
rapidly with time tpp, and therefore quickly become positive. For all examples shown in this
chapter, we use the lower bound ul = 0.
We calculate the upper bound uu from an estimate of the average VP/VS ratio and a
maximum time of interest. Analogous to equation 3.5, the average VP/VS ratio is
γa(tpp) = 2
tps
tpp
− 1. (3.7)
Therefore, if Tps is the maximum time of interest in a PS image, then the corresponding time
in the PP image is
Tpp =
2Tps
1 + γa(Tpp)
. (3.8)
Note that Tpp is not the maximum recording time in the PP image. Rather, Tpp is the
maximum time in the PP image for which a corresponding event of interest appears in the PS image at the maximum PS time, Tps. We will ignore times tpp > Tppin the PP image for
which there are no corresponding events of interest in the PS image.
If we could solve equation 3.8 for Tpp, then we could simply compute the upper bound on
have estimated time shifts u(tpp), a precise upper bound uu is unnecessary. In practice, we
need only ensure that we do not set this bound too low, so that the correct shifts u(tpp) will
not exceed uu. This requirement means that we should generally overestimate the average
VP/VS ratio γa(Tpp) at the maximum PP time Tpp. For the images shown in Figure 3.1, we
find that γa(Tpp) ≈ 2, so that uu ≈ Tps/3 is a safe choice.
In addition to the bounds ul and uu on estimated time shifts u(tpp), the constraints 2.3
define lower and upper bounds on time strain, rl and ru, respectively. According to
equa-tion 3.6, these time-shift derivatives are directly related to interval VP/VS ratios γi(tpp).
Therefore, rland ru can be expressed in terms of lower and upper bounds on interval VP/VS,
γiM in ≤ γi ≤ γiM ax, from which we compute
rl= γiM in − 1 2 , (3.9) ru = γiM ax− 1 2 . (3.10)
Because γi ≥ 1, a safe lower bound is rl = 0. Since rl ≥ 0, the derivative of the time shifts
cannot be negative, and thus the time-shift curve must be monotonically increasing. The upper bound ru depends on the maximum interval VP/VS ratio γiM ax. Because this ratio is
seldom known in advance, an overestimate should generally be used.
For the PP and PS images shown in Figure 3.1 we must compute time shifts u(tpp) at
all x, y locations in the 3D image, thus the time-shift volume u(x, y, tpp) is the solution we
must compute with smooth dynamic image warping. The bounds 3.9 and 3.10 apply to only derivatives of time shifts with respect to time tpp. Although the time shifts u(x, y, tpp) can
never decrease with increasing time tpp, they may decrease or increase with increasing x or
y coordinates. The lateral constraints defined in section 2.4 control rates at which the time shifts u(x, y, tpp) may vary laterally. These lateral time-strain constraints are related to the
rates at which corresponding estimates of VP/VS ratios, γi(x, y, tpp), can vary laterally. In
practice, we do not know in advance what these rates will be, so we typically begin with loose bounds on lateral variations and tighten them as necessary to reduce suspicious oscillations
in estimated VP/VS ratios.
3.4 SDIW versus DIW
Chapter 2 showed how smooth dynamic time warping (Algorithm 2) could estimate changes in time shift more accurately than dynamic time warping (Algorithm 1), for a synthetic 1D example. In this section I compare those two algorithms applied to the real data shown in Figure 3.1. Algorithms 1 and 2, extended to multidimensional images, will be referred to as dynamic image warping (DIW) and smooth dynamic image warping (SDIW), respectively.
DIW was previously used by Liang and Hale (2012) to register 2D PP and PS images and to estimate VP/VSratios. As discussed in Chapter 2, DIW returns a sequence of integer shifts,
and the number of possible shift derivatives is small; thus from equation 3.6, the number of possible computed γi(tpp) values is correspondingly small. Therefore, an important step
in Liang and Hale’s (2012) application is to smooth the optimally computed sequences of integer shifts before estimating interval VP/VS. Since time shifts are related to the integral
of VP/VS, we expect u to be a smooth function, compared to γi; however, estimating VP/VS
by smoothing a rough sequence of integer shifts may be inaccurate. SDIW directly computes a smooth and more accurate time-shift function (Hale and Compton, 2013).
Figure 3.2 shows the alignment errors for a single pair of traces in the PP and PS images from Figure 3.1. Overlaid on the alignment errors are time-shift curves extracted from the 3D DIW and SDIW algorithms. For a more detailed comparison, these curves and the alignment errors are shown for only a 1-second time window.
The solid white curve in Figure 3.2 shows the time shifts computed from DIW. Note the roughness of the white curve, which results from the fact that DIW computes time shifts at every time sample. In fact, the white time-shift curve contains only two slopes, 0 and 1. These slopes correspond to the bounds rl and ru from equation 2.3, and thus γiM in = 1
and γiM ax = 3 from equations 3.9 and 3.10, respectively. Not only are γiM in and γiM ax the
Figure 3.2: Alignment errors and time shifts computed from a single pair of traces from the PP and PS images shown in Figure 3.1. Only the alignment errors for times 0 to 1 second are shown. Overlaid curves are the 1D time shifts extracted from the 3D time shifts u(x, y, tpp)
at the corresponding trace location. The solid white curve shows rough time shifts computed from DIW and the dashed blue curve shows these shifts after smoothing. The dashed yellow curve shows the time shifts computed from SDIW with a minimum subsampling interval h = 100 ms.
of γi(tpp) we can compute. Hence, the motivation to smooth these shifts, as Liang and Hale
(2012) do, is to obtain more estimates of interval VP/VS. The dashed blue curve in Figure 3.2
shows the white time-shift curve after smoothing.
The yellow curve in Figure 3.2 shows the interpolated time shifts computed from SDIW. With SDIW, time shifts are computed at only subsampled PP times and the time shifts u(x, y, tpp), are obtained by interpolation of those subsampled time-shift estimates. The
subsampled locations correspond to high-amplitude reflectors in the PP image that, in this example, are a minimum of 100 ms apart. This selection process is discussed in a section below.
Note the large difference between the yellow and blue curves between 0 and 0.6 s. To see that the yellow curve in Figure 3.2 is more accurate, look at the warped images in Figure 3.3. Compare the PP image with the PS image warped using smoothed DIW time shifts, and with the PS image warped using SDIW time shifts. Recall that the curves overlaid on alignment
Figure 3.3: Slices through three 3D volumes with the same PP time window shown in the alignment errors of Figure 3.2. The location for those alignment errors corresponds to the vertical black lines in the inline and crossline sections. The strong events at 0.7 s are well aligned in each warped image, but the time slice through the event at 0.38 s shows that time shifts computed using SDIW are more accurate. The differences in time shifts can be seen clearly in Figure 3.2.
errors in Figure 3.2 were extracted from 3D shift functions u(x, y, tpp). The location of those
alignment errors corresponds to the vertical black lines in the inline and crossline sections of Figure 3.3.
The reflectors at 0.7 s are well aligned in both warped PS images in Figure 3.3, and we can see in Figure 3.2 that the yellow and blue curves are consistent at this time. However, a large discrepancy between the two shift functions can be seen at 0.38 s in the alignment errors in Figure 3.2, and correspondingly in the constant-time (0.38 s) slices of the warped PS images in Figure 3.3. The warped PS image using the SDIW time shifts is a better match to the PP image, indicating that the yellow curve is more accurate.
The subsampling in SDIW greatly increases the number of possible changes in shift, and, as was shown in Figure 3.2 and Figure 3.3, increases the accuracy of estimated time shifts. The white curve in Figure 3.2 demonstrates that trying to resolve changes in time shifts at every image sample can result in misalignment of image features; furthermore, accurate time
shifts cannot always be obtained by smoothing inaccurate time shifts.
Figure 3.4: Interval VP/VS ratios estimated from smoothed DIW time shifts and from SDIW
time shifts. The VP/VS ratios from smoothed DIW time shifts exhibit fine detail, but those
smoothed shifts are less accurate, as shown in Figure 3.2 and Figure 3.3, than shifts estimated using SDIW.
Figure 3.4 shows the interval VP/VS ratios estimated from DIW time shifts and SDIW
time shifts. The two images exhibit some similar trends of high or low VP/VS; however,
the interval VP/VS ratios computed from DIW time shifts show a high level of detail that is
unwarranted given the resolution of seismic reflections and differences in noise and reflection waveforms between the PP and PS images.
3.5 Coarse sample locations
Coarse sampling improves the accuracy of computed shifts, but the placement of the coarse lattice of sample points affects the accuracy as well. Increasing the distance between coarse samples increases the smoothness of estimated time shifts. One simple approach
to choosing coarse sample locations is to use a nearly-uniform grid parameterized by a coarse sampling interval in each dimension. We examine this approach for a 1D example, where we compute time shifts u(tpp) using smooth dynamic time warping, and focus only on
subsampling the PP-time axis.
Figure 3.5a shows alignment errors computed from a single pair of traces from the 3D PP and PS images shown in Figure 3.1. In Figure 3.5b the blue curve represents the time shifts u(tpp) after interpolation and the blue dots represent the coarse sample locations where time
shifts were computed. In this example, we specified a coarse interval of 160 ms and computed sample locations that are almost uniformly spaced. The spacing is slightly variable so as to include the first and last time samples; however, all intervals are greater than or equal to the specified interval of 160 ms.
In Figure 3.5, note that strong events in the PP trace (Figure 3.5c) appear as pre-dominantly vertical features with high error (white) in the array of alignment errors. The diagonally crossing features with high error correspond to strong events in the PS trace. We can estimate shifts best at the locations where these strong events occur. However, by selecting coarse sample locations on a nearly-uniform grid, as shown by the blue dots in Figure 3.5b, strong events are sampled only when they happen to coincide with our coarse sampling grid. Rather than selecting coarse samples on a nearly-uniform grid, we should instead select coarse sample locations to coincide with high-amplitude events in the PP trace.
Let f (tpp) denote the PP trace and a(tpp) denote the trace envelope computed as
a(tpp) =
q f2(t
pp) + ˆf2(tpp), (3.11)
where ˆf (tpp) is the Hilbert transform of f (tpp). Figure 3.5c shows the PP trace in black and
the trace envelope in red. An easy way to improve the coarse sample locations is to use the envelope a(tpp) to preferentially select locations with the highest amplitudes, while also
satisfying a specified minimum interval. We search for coarse grid points in descending order of envelope amplitudes in a(tpp). A coarse sample is chosen when its location is at a distance
Figure 3.5: Alignment errors (a) computed from a single trace from the PP and PS images shown in Figure 3.1, with time shifts (b) computed from 1D smooth dynamic time warping using two different subsampling grids. A nearly-uniform coarse grid was used for the blue curve, with subsampling locations represented by blue dots, and a reflector-aligned grid was used for the yellow curve with subsampling locations represented by yellow squares. The PP trace (c) overlaid with the envelope of the trace shows that strong PP events correlate with strong vertical features in the alignment errors above.
samples. As with the nearly-uniform-grid scheme, we include the endpoints of the original sampling grid.
The yellow squares in Figure 3.5b represent the coarse sample locations where time shifts were computed along the yellow time-shift curve. These squares were preferentially selected to correspond with locations of strong events, while maintaining a minimum separation interval of 100 ms. We reduce the minimum interval used for the nearly-uniform grid (160 ms) in order to include some of the highest amplitude events that are within that interval of one another. In this example, we have selected the same number of subsamples for the nearly-uniform grid and the reflector-aligned grid.
The time-shift curves shown in Figure 3.5b for each grid are quite similar from 0 to 1.6 s except for a large discrepancy around 1.2 s. The blue dot indicates that a time shift was computed for that location, but there is no corresponding yellow square nearby. From both the alignment errors and the envelope of the PP trace in Figure 3.5c we see that there are weak events at, and around 1.2 s; therefore, we have little confidence in the accuracy of the computed time shift at the blue dot. After 1.6 s, the two curves are quite different, with a significant change in shift at approximately 1.6 s. Time shifts shown by the blue curve are inaccurate, based on 3D time shifts and warping results not shown here. This discrepancy in the blue curve is largely due to the fact that it is difficult to accurately register PP and PS events using a single pair of traces. However, the yellow curve corresponds to the correct shifts, even in the 1D case. Therefore, this example demonstrates that better placement of coarse sample locations can improve shift estimates.
When smoothing and subsampling alignment errors in 3D, we require the same number of coarse-grid time samples for all coarse inline and crossline sample locations. Therefore, we wish to select the highest amplitudes that are consistent throughout the image, i.e., the strongest reflectors. We cannot directly use the method described above because reflectors in seismic images are generally not horizontal, but when seismic reflectors are not structurally complex, we can make them horizontal using seismic image flattening (Lomask, 2006; Parks
et al., 2008; Parks, 2010).
Figure 3.6a shows the PP image, which is structurally simple, but contains dipping reflectors. Therefore, we compute a volume of shifts (Figure 3.6b) that flatten image features so that a constant time in the flattened image (Figure 3.6c) corresponds to a single seismic horizon. In this flattened image space we stack the envelopes of all image traces:
as(tpp) = ni3−1 X i3=0 ni2−1 X i2=0 ai2i3(tpp), (3.12)
where ni3 and ni2 are the number of inline and crossline samples, respectively. Now, we can
search for times corresponding to the strongest reflection events in the flattened 3D image and compute a coarse lattice of points in all three dimensions, indicated by the white dots in Figure 3.6c. We then map the coordinates of these lattice points from the flattened space to the original image space to obtain a reflector-aligned 3D lattice on which to compute time shifts using SDIW.
Figure 3.7 compares a simple nearly-uniform grid, defined by coarse-interval parameters in each dimension, with a reflector-aligned coarse grid. Both grids have the same nearly-uniform sampling in the inline and crossline directions. Time shifts are estimated at coarse sample locations, and we improve the accuracy of these estimates by aligning the coarse grid with strong reflectors. Figure 3.7 shows interval VP/VS ratios estimated from the
nearly-uniform grid and the reflector-aligned grid. Note the step in VP/VS at 1.7 and 1.8 s in
the crossline section in Figure 3.7. This step compensates for the nearly-uniform coarse grid that samples across a dipping layer in which VP/VS is consistently low throughout the
image. This layer is more apparent at the same location in the reflector-aligned VP/VS ratios
because coarse grid points are aligned with reflectors in the image.
Note that this automated approach to picking sample locations is independent of our smooth image warping algorithm. Thus, a more interpretive approach of selecting specific seismic horizons could be used instead, especially for a more complex image that might be more difficult to flatten.
Figure 3.6: The PP image (a) and time shifts (b) used to flatten the image so that a constant time corresponds to a single seismic horizon. The coarse lattice of points on the flattened image (c) is shown by the white dots. This lattice can be mapped back to the original (unflattened) coordinate space to obtain a coarse lattice of points that is aligned with dipping reflectors (d). We use this lattice in SDIW to compute time shifts at locations with strong reflectors in the PP image.
Figure 3.7: White dots indicate coarse lattice points used to compute time shifts with SDIW. Both nearly-uniform and reflector-aligned grids have minimum intervals of 500 m laterally. The nearly-uniform points are a minimum of 160 ms apart in time, while reflector-aligned points are a minimum of 100 ms apart but are aligned with strong reflectors in the image. Points in the nearly-uniform grid cross strong reflectors horizontally, resulting in a noticeable stair step effect in the corresponding VP/VS ratios.
3.6 Interpolation of time shifts
The smooth dynamic warping algorithm computes time shifts only at specified subsam-pled locations to improve accuracy, but we ultimately need time shifts for every image sample to align PP and PS images. Therefore, time shifts for all image samples not specified by the coarse grid are interpolated from the computed coarse lattice of time shifts.
For the 3D examples shown here, we do this interpolation in two steps. The interpolation is first done laterally using bilinear interpolation. With bilinear interpolation we ensure that the bounds on shifts ul and uu and our assumption of monotonically increasing time shifts
are satisfied. For vertical interpolation of our time shifts, we must consider the effect on the VP/VS ratios we are trying to estimate.
Figure 3.8: Interval VP/VS estimates resulting from two different interpolation methods for
time shifts, piecewise cubic with monotonic splines and piecewise linear.
As described in Chapter 2, the SDIW algorithm accumulates alignment errors at coarse sample locations by summing errors along piecewise-linear trajectories. Therefore, the most
consistent approach for interpolating time shifts is piecewise-linear interpolation. Figure 3.8 shows the estimated interval VP/VS ratios when time shifts are interpolated using
piecewise-cubic interpolation with monotonic splines and using piecewise-linear interpolation. Both results were obtained from time shifts computed at coarse sample locations shown in the reflector-aligned grid of Figure 3.7.
The VP/VS ratios from piecewise-linear interpolation of time shifts in Figure 3.8 are
piecewise-constant between subsampled time locations and are more consistent with the SDIW algorithm than VP/VS ratios from piecewise-cubic interpolation of time shifts. We
computed time shifts at only coarse sample locations, and thus estimated the change in shift between only these locations. When smoother estimates of VP/VS are required, we can
interpolate time shifts using a smoother method. However, smoother VP/VS ratios will be
less consistent with the SDIW method used to compute time shifts.
Figure 3.9 shows 1D time shifts u(tpp) and corresponding interval VP/VS ratios γi(tpp)
extracted from 3D functions u(x, y, tpp) and γi(x, y, tpp), respectively. The three curves in
each image show results from three methods of interpolation. The black curves show linearly-interpolated time shifts and the corresponding piecewise-constant VP/VS ratios. The other
two interpolation methods use piecewise-cubic polynomials. The red time-shift curve was interpolated using cubic polynomials that preserve monotonicity, while the blue curve was interpolated using the classic cubic spline with continuous first and second derivatives.
Linear interpolation is the only method that guarantees time shifts will be both monotonic and bounded by [rl, ru]. Monotonicity is still guaranteed for the red time-shift curve, but
the bounds [rl, ru] are no longer strictly satisfied. For example, at 1.7 s the red curve dips
slightly below the lower bound, rl = 1.5, used in this example. Cubic spline interpolation
is shown for completeness, but this method is worse still. Overshoot with this method will generally be more severe than with piecewise-cubic monotonic interpolation. For this method we can guarantee neither monotonicity of interpolated time shifts nor derivatives bounded by [rl, ru]. Therefore, it is possible to compute VP/VS < 1, but for the example shown in
Figure 3.9: Interpolated time shifts and interval VP/VS extracted from 3D SDIW shift and
VP/VS volumes. The three different curves represent three different methods of interpolating
time shifts: linear (black), cubic monotonic spline (red), and cubic spline (blue).
Figure 3.9 this does not occur.
If we want interval VP/VS estimates that are continuous, then cubic monotonic
interpo-lation might be satisfactory. Otherwise, the piecewise-constant VP/VS ratios are the most
CHAPTER 4
SIMULTANEOUS MULTICOMPONENT REGISTRATION
4.1 Introduction
The previous chapter discussed the application of smooth dynamic image warping to automatically register a pair of multicomponent PP and PS images. However, in multi-component surveys more image pairs can exist. For example, a 9C multimulti-component seismic survey can capture nine different wave modes, yielding nine different images. Even the sim-plest 3C multicomponent seismic survey can result in more than the PP and PS images previously considered.
In the presence of azimuthal anisotropy, which is associated with fractures or differen-tial stress, the mode-converted S-wave can split into fast and slow waves with orthogonal polarizations. This phenomenon is known as shear-wave splitting or seismic birefringence. The fast-S mode is commonly referred to as S1 and is polarized parallel to fracture strike or maximum stress. The slow-S mode, S2, is polarized parallel to fracture normal or min-imum stress. The S1 polarization angle and the time lag between S1 and S2 give valuable information about anisotropy and fracture characterization.
Thus the simple 3C multicomponent seismic survey can record three combinations of downgoing and upgoing waves: PP, PS1, and PS2. With multicomponent seismic data processing techniques, these different wave modes can be separated into three corresponding seismic images.
The automatic image registration method discussed in the previous chapter could be used to register independently the three image pairs: PP-PS1, PP-PS2, and PS1-PS2. This straightforward approach yields an estimated time-shift volume for each image pair and these three time-shift volumes should be related. However, independently warping image pairs will result in inconsistent time-shift estimates. Instead, the time-shift relations can be
incorporated into a single minimization problem involving all three images simultaneously. This chapter formulates the simultaneous dynamic warping problem and demonstrates its solution using synthetic PP, PS1, and PS2 seismograms.
4.2 Traveltime relations
Section 3.2 describes traveltime relations for PP and PS seismic traces, and shows how time shifts that align these traces can be used to estimate VP/VS ratios. This section derives
similar traveltime relations for PP, PS1, and PS2 traces. From these three traces we can estimate not only VP/VS ratios, but also anisotropy, which causes time shifts between PS1
and PS2 traces.
Let tpp denote the PP traveltime for an event in a PP trace, and let tps1 and tps2 denote
the traveltimes for the same event in PS1 and PS2 traces, respectively. For each tpp, we have
a corresponding tps1 and tps2 given by
tps1(tpp) = tpp+ u1(tpp), (4.1)
tps2(tpp) = tpp+ u2(tpp), (4.2)
where u1(tpp) denotes time shifts between PP and PS1 traces and u2(tpp) denotes time shifts
between PP and PS2 traces, and both time shifts are functions of tpp. Additionally, time
shifts can be found between PS1 and PS2 traces. Therefore, each tps1 has a corresponding
tps2 given by
tps2(tps1) = tps1 + ˜uS(tps1), (4.3)
where ˜uS(tps1) denotes time shifts between PS1 and PS2 traces and is a function of tps1.
To see the relationship between time shifts u1(tpp), u2(tpp), and ˜uS(tps1), we use
equa-tion 4.1 to rewrite equaequa-tion 4.3 in terms of tpp, as
Equations 4.4 and 4.2 give the relationship between the three time shifts
u2(tpp) = u1(tpp) + ˜uS(tpp+ u1(tpp)). (4.5)
The time shifts ˜uS(tpp+u1(tpp)) are specified at PS1 times given by tpp+u1(tpp), but ultimately
we want a function of PP time
uS(tpp) ≡ ˜uS(tpp+ u1(tpp)), (4.6)
which can be computed from the time shifts ˜uS(tps1) and u1(tpp). Of course, the accuracy of
uS(tpp) will depend partly on the accuracy of u1(tpp). Equation 4.6 simplifies equation 4.5 to
u2(tpp) = u1(tpp) + uS(tpp). (4.7)
This relationship is intuitive because of the relative traveltimes for PP, PS1, and PS2 traces. The longest traveltimes correspond to the PS2 trace and the shortest traveltimes correspond to the PP trace. Thus the largest time shifts are between PP and PS2 traces (u2(tpp)), and
are a combination of the smaller time shifts.
The time delay between two shear waves is conventionally described by the fractional difference
γ(s)≡ VS1− VS2 VS2
, (4.8)
where VS1 and VS2 represent the fast and slow shear mode velocities, respectively (e.g.,
Tsvankin (2012)). Equation 4.8 provides a measure of anisotropy, but requires the velocities VS1 and VS2. Alternatively, γ(s) can be computed from time shifts that align PP, PS1 and
PS2 traces. In fact, a combination of any two of the three time shifts u1(tpp), u2(tpp), and
uS(tpp), can be used to compute γ(s)(tpp).
