Reflection seismic data Q tomography

Abstract

Method for reconstructing subsurface Q depth profiles from common offset gathers (92) of reflection seismic data by performing migration (40), ray tracing (100), CDP-to-surface takeoff angle finding (96, 98), kernel matrix construction (110), depth-to-time conversion and wavelet stretching correction (80), source amplitude spectrum fitting, centroid frequency shift calculation (90), and box-constrained optimization (120).

Claims

1. A ray-based Q-tomography method of generating a subsurface model for prospecting for hydrocarbons using recorded reflection seismic data from a subsurface region, comprising: migrating the seismic data to depth domain to obtain common image gathers, said common image gathers being common offset gathers, where a kernel matrix is constructed using raypath information obtained through a ray tracing procedure, and selected seismic traces from a common image gather are converted to time domain, then a centroid frequency shift relative to a seismic source wavelet centroid frequency is determined, then an expression relating the centroid frequency shift and the kernel matrix to seismic attenuation, as measured by a frequency independent quality factor Q, is iteratively solved for Q, a representation of seismic attenuation, using a computer, generating a subsurface model that identifies an anomaly in the subsurface corresponding to the Q; and using the subsurface model to prospect for hydrocarbons; wherein the kernel matrix construction includes, (a) in each image corresponding to a common offset, selecting at least one horizon, and on the horizon at least one common depth point, (b) tracing two rays from each selected common depth point to the surface, and (c) measuring ray lengths in each penetrated subsurface grid, each subsurface grid corresponding to a column index of the kernel matrix, and using a subsurface velocity model to form elements of the kernel matrix; wherein the selected seismic traces comprise a vertical trace for each selected common depth point, extracted from the common image gather; wherein the conversion to time domain of the selected seismic traces is performed while applying a wavelet stretching correction and is expressed by $t\left(z\right)=\frac{z\left(z\right)}{v\left(z\right)}\left[\cos \left({}_i+{}_d\right)+\cos \left({}_i-{}_d\right)\right]$ where t is a time duration corresponding to z, which is proportional to vertical mesh size in a computational grid, as a function of depth z; v(z) is vertical seismic velocity profile as a function of the depth z; .sub.i is seismic ray incident angle at the selected common depth point; and .sub.d is dip angle of the selected horizon at the selected common depth point; and wherein the centroid frequency shift is determined for each selected trace after transforming the selected trace to frequency domain.

2. The method of claim 1, wherein (b) comprises: (i) estimating a dip angle of the selected horizon where the selected common depth point is located; (ii) finding two takeoff directions making equal angles with the normal direction to the horizon at the selected common depth point such that two rays traced using Snell's Law from the selected common depth point at the two takeoff directions reach surface locations that are separated by the common offset; (iii) tracing the two rays defined in (ii).

3. The method of claim 1, wherein (ii) comprises: deploying a plurality of virtual sources on the surface; performing wavefront construction based ray tracing from the virtual sources to the selected common depth point and building a table giving ray incident angle at the selected common depth point for each virtual source location; and selecting from the table two incident angles for the selected common depth point that match the common offset, and using these two selected incident angles as the takeoff directions.

4. The method of claim 1, wherein the determination of the centroid frequency shift comprises analyzing the source wavelet frequency distribution and fitting an asymmetric frequency distribution with a frequency weighted exponential function of frequency, and further comprising computing a component of a measurement vector from the centroid frequency shift and using the measurement vector to represent the centroid frequency shift in the expression that is solved by iterative optimization.

5. The method of claim 4, wherein the frequency-weighted exponential function of frequency can be expressed in the form $F\left(f\right)=A{f}^n\exp \left(-\frac{f}{f_0}\right)$ where f is frequency, A is a constant for amplitude scaling, f.sub.0 is a characteristic frequency, and n is a symmetry index.

6. The method of claim 1, wherein the extracted vertical traces are truncated to include only the events associated with the selected horizon.

7. The method of claim 1, wherein the iterative solving for Q comprises a linear optimization using box constraints to keep estimated Q values within position-dependent ranges specified by upper boundaries and lower boundaries.

8. The method of claim 1, wherein the iterative solving for Q comprises a constrained optimization that is solved by a multi-index active-set method that allows updates of the active set by multiple computational grid indices at a time, wherein a grid index denotes subsurface location.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The present invention and its advantages will be better understood by referring to the following detailed description and the attached drawings in which:

(2) FIG. 1A shows a velocity model and transmission raypaths;

(3) FIG. 1B shows a scenario where there is a shallow Q anomaly covered by the transmission seismic rays, and a deep Q anomaly which is not covered by the transmission seismic rays but is penetrated by reflection seismic rays;

(4) FIG. 2 shows a transmission seismic shot-to-receiver raypath in a gridded model;

(5) FIG. 3 shows multiple reflection seismic raypaths for a single shot-receiver pair;

(6) FIG. 4 shows a diagram illustrating the migration process of converting time domain reflection seismic data to depth domain seismic data;

(7) FIG. 5 schematically illustrates a common offset gather obtained via migration of reflection seismic data;

(8) FIG. 6 schematically illustrates a common angle gather via migration of reflection seismic data;

(9) FIG. 7 shows the ray tracing from a selected CDP point P on a selected horizon in an image from the common image gather;

(10) FIG. 8 shows the depth domain seismic traces extracted from the common image gather;

(11) FIG. 9 shows the time domain seismic traces obtained through depth-to-time conversion and the wavelet stretching correction;

(12) FIG. 10 is a flow chart showing basic steps in one embodiment of the present inventive method;

(13) FIG. 11 shows the layered velocity model used for the example;

(14) FIG. 12 shows the post-migration seismic data for a specific offset;

(15) FIG. 13 shows the amplitude spectra of the source, the extracted post-migration reflection seismic data, and the centroid frequency shift due to the Q anomaly;

(16) FIG. 14A shows the true Q model for the example;

(17) FIG. 14B shows the reconstructed Q model using the present inventive method;

(18) FIG. 15 shows wavefront-construction-based ray tracing from a virtual source to the subsurface; and

(19) FIG. 16 shows ray tracing from a selected CDP point P on a selected horizon in an image from the common offset gather.

(20) FIGS. 1A, 1B, 11, and 14A-14B are black-and-white reproductions of colored original data displays, due to patent law restrictions.

(21) The invention will be described in connection with example embodiments. However, to the extent that the following detailed description is specific to a particular embodiment or a particular use of the invention, this is intended to be illustrative only, and is not to be construed as limiting the scope of the invention. On the contrary, it is intended to cover all alternatives, modifications and equivalents that may be included within the scope of the invention, as defined by the appended claims.

DETAILED DESCRIPTION OF EXAMPLE EMBODIMENTS

(22) The present invention includes a method for reconstruction of 2D or 3D seismic quality factor (Q) profiles from reflection seismic data, a technical field known as reflection seismic data Q tomography.

(23) The main features of the present invention in at least some embodiments are as follows. The reflection seismic data are migrated by any migration method including, but not limited to Kirchhoff migration to obtain a common image gather (CIG), which can be a common angle gather or a common offset gather. At least one common depth point (CDP) is selected in the CIG. In other words, the same point is located in each image of the CIG. If the migration velocity model happens to be perfect, then this selected CDP will be located at the same (x,y,z) coordinates in each image; otherwise not. With estimated dip information, two seismic rays are shot from each selected CDP to the surface using ray tracing according to Snell's Law. For a common angle gather, the takeoff angles for this CDP-to-surface ray tracing are known. For common offset gathers, the takeoff angles for the CDP-to-surface ray tracing may be found through deploying virtual sources on the surface, performing wavefront-construction-based ray tracing from the virtual sources to the subsurface and saving the incident angles at the CDPs and building the angle table, and implementing an angle selection procedure to find the matched takeoff angles for the selected CDPs for the specific offsets. The raypath information of the CDP-to-surface ray tracing is collected for all the rays shot from the selected CDPs to construct the kernel matrix. For each selected CDP, the associated vertical depth domain seismic trace in the CIG is extracted. The extracted vertical traces are converted to the time domain and then to the frequency domain through depth-to-time conversion, a wavelet stretching correction, and fast Fourier transform (FFT). The source wavelet amplitude spectrum is analyzed and fitted by a frequency weighted exponential function and the centroid frequency of the source wavelet is calculated. The centroid frequency shifts (with respect to the source wavelet centroid frequency) of the extracted traces are calculated to construct the measurement vector. The kernel matrix and the measurement vector are input to an optimization algorithm with box constraints to reconstruct the Q profile, where the ranges of the Q values are predetermined by a priori information. This constrained optimization is solved by employing the multi-index active-set method (Morigi et al., 2007).

(24) Some underlying theory of the invention is explained next.

(25) First, with the assumption that an accurate velocity model has been provided, the time domain reflection seismic data are migrated to obtain the common image gather (either common offset gather or common angle gather), which are the post-migration pre-stack reflection seismic data. By doing this, the time domain seismic data are converted to the depth domain. FIG. 4 is the diagram to show the migration process. In FIG. 4, the data domain consists all the shots and receivers represented by the dots on the surface and the model domain are the subsurface rectangular area. The seismic traces are back-propagated downward to the model domain to obtain the depth domain seismic data. For example, the horizon shown in FIG. 4 is the actually the seismic energy in the depth domain mapped from the time domain reflection seismic data. Separating the contributions from different offsets or from different angles, a common offset gather or common angle gather, which is a series of images as shown in FIG. 5 and FIG. 6, can be obtained.

(26) In each image in the common image gather (common offset gather or common angle gather), at least one horizon is selected and on this selected horizon, at least one CDP point is selected, as shown in FIG. 7 (point P). In FIG. 7, the dashed line is the normal direction to the horizon; the dot-dashed line and the dotted line are the vertical and horizontal directions, respectively. Then, for common angle gathers, from this selected CDP point P, because the takeoff angle .sub.i is known, with the dip angle .sub.d estimated by a dip estimation tool (Marfurt, 2006), two rays 71 and 72 are shot to the surface. After that, the ray length in each grid is measured. Let l.sub.j.sup.lk be the ray length in the j.sup.th grid cell for the k.sup.th selected point in l.sup.th image in the CIG, the elements in the kernel matrix A can be determined by (Hu, 2011)

(27) $\begin{array}{cc}{A}_{ij}=\frac{l_j^{kl}}{v_j},& \left(1\right)\end{array}$
where i is the row index of the kernel matrix corresponding to k and l (each combination of k and l corresponds to a distinct row index i of the kernel matrix, where the relationship between i, k and l depends on the number of images in the common image gather), and j is the column index of the kernel matrix, and v.sub.j is the seismic velocity in the j.sup.th grid. (In this document, the term grid, which might alternatively be called cell, means a single element in the gridded velocity model.) The whole kernel matrix can be constructed by repeating this procedure for all the selected images, horizons, and CDP points.

(28) Unfortunately, common angle gathers are not always available. In many applications, one may have only common offset gathers as shown in FIG. 5. Therefore, to shoot two rays from the selected CDP point to the surface, it is first necessary to find the corresponding takeoff angle at the selected CDP point. A procedure for finding the takeoff angle is as follows:

(29) (1) Deploy a series of virtual sources (S.sub.1, S.sub.2, . . . , S.sub.P, . . . ) on the surface at each grid as shown in FIG. 15.

(30) (2) From each virtual source, perform wavefront construction based ray tracing (see, for example, Vinje et al., 1993); then, at each CDP point i, save the incident angle .sub.i.sup.p, which is the angle between the ray direction and the vertical direction, to build a lookup angle table, where p is the index of the associated virtual source.
(3) At each selected CDP point in an image extracted from the common offset gather with offset d.sub.off, find two angles .sub.i.sup.p1 and .sub.i.sup.p2 satisfying .sub.i=.sub.r=(.sub.i.sup.p2.sub.d) and d.sub.off=|x.sub.p1x.sub.p2|, where .sub.d is the dip angle at the selected CDP point, p1 and p2 are the indices of the virtual sources S.sub.p1 and S.sub.p2 as shown in FIG. 16, whose positions are represented by x.sub.i1 and x.sub.i2.
(4) With the takeoff angle .sub.i.sup.p1 or .sub.i.sup.p2 and the dip angle .sub.d, two rays are shot from this selected CDP point to the surface, as shown in FIG. 16.
After that, similar to the procedure implemented in the common angle gather cases, a kernel matrix can be constructed by using equation (1). All variations and equivalent or similar approaches to the above 4-step method are within the scope of the present invention.

(31) A key part of the present invention is the construction of the measurement vector, which is explained as follows.

(32) In a specific image from a CIG, as shown in FIG. 7 or FIG. 16, a vertical post-migration reflection seismic trace W(z) is extracted at the selected CDP P and this depth domain trace is then truncated to contain only the wavelet w(z) around the position P. These extracted wavelets w(z) are shown in FIG. 8. In FIG. 8, it is observed that the widths of the wavelets are offset-dependent, which is a well-known wavelet stretching phenomenon. This wavelet stretching issue needs to be addressed in order to extract the correct frequency content information to be used in Q tomography. Another phenomenon observed is that, when Q anomalies (within which Q values are finite) are present, the wavelet widths are larger as long as the traces penetrate one of these Q anomalies. This phenomenon is expected because Q anomaly attenuates high frequency more significantly, which leads to a downshift of the frequency content. This frequency shift information will eventually be used to estimate the Q values in the Q tomography algorithm. However, these extracted and then truncated seismic traces shown in FIG. 8 are not ready for frequency content analysis because they are still in depth domain. Therefore, the extracted depth domain seismic traces w(z) need to be converted to time domain and the wavelet stretching correction needs to be applied (a single step). This procedure is implemented by using the following formulation for conversion between time and depth:

(33) $\begin{array}{cc}t\left(z\right)=\frac{z\left(z\right)}{v\left(z\right)}\left[\cos \left({}_i+{}_d\right)+\cos \left({}_i-{}_d\right)\right],& \left(2\right)\end{array}$
where v(z) is the vertical velocity profile at CDP P, .sub.i is the incident angle and .sub.d is the dip angle as shown in FIG. 7. The depth domain seismic trace w(z) can be converted to time domain signal w(t) shown in FIG. 9 by using formulation (2). In FIG. 9, it may be observed that, after the wavelet stretching correction is applied, the width of the time domain wavelet is independent of offset if the Q value is infinity, which corresponds to zero attenuation. However, if Q values are finite, the wavelet width is still dependent on the offset because the Q profile along the ray path determines the frequency content modification of the trace.

(34) The time domain seismic traces in FIG. 9 are now ready for frequency content analysis. First, the traces w(t) are converted to frequency amplitude spectrum w(t) by using the Fourier transform, preferably the fast Fourier transform (FFT). After that, the centroid frequency of the signal f.sub.c.sup.w and the centroid frequency of the source f.sub.c.sup.s are calculated respectively using

(35) $\begin{array}{cc}{f}_c^w=\frac{{}_0^{+}fw\left(f\right)df}{{}_0^{+}w\left(f\right)df}\mathrm{and}& \left(3\right)\\ f_c^s=\frac{{}_0^{+}fs\left(f\right)df}{{}_0^{+}s\left(f\right)df},& \left(4\right)\end{array}$
where s(f) is the frequency amplitude spectrum of the source wavelet. The centroid frequency shift is
f.sub.c=f.sub.c.sup.sf.sub.c.sup.w.(5)
Assuming that the source frequency amplitude spectrum can be approximated by a frequency weighted exponential function (see Hu 2011, which is incorporated herein by reference)

(36) $\begin{array}{cc}s\left(f\right)=A{f}^n\exp \left(-\frac{f}{f_0}\right),& \left(6\right)\end{array}$
where A, n, f.sub.0 are real numbers, one is able to find the values of A, n, f.sub.0 by fitting the frequency weighted exponential function (6) to the source frequency amplitude spectrum. (The constant A is unrelated to the kernel matrix A.) The measurement vector for centroid frequency shift based reflection seismic Q tomography can now be constructed as

(37) $\begin{array}{cc}b=\frac{f_c}{f_0{f}_c^w}.& \left(7\right)\end{array}$

(38) With equations (1) and (7), the reflection seismic data Q tomography can be posed as an optimization problem with box constraints (See Eqn. (14) in Hu 2011)
minAxb subject to l<x<u,(8)
where l and u are the vectors storing the lower boundaries and the upper boundaries of the 1/Q values, and x is a vector of the unknowns, i.e.,
x.sub.j=1/Q.sub.j.

(39) A preferred type of algorithm that may be employed to solve the optimization problem of Eqn. (8) is a multi-index active-set method such as was disclosed by Morigi et al. (2007).

(40) In one of its embodiments, the present invention can be implemented according to the flow chart shown in FIG. 10. In step 40, the velocity model 10, the seismic data 20, and the source/receiver positions 30 are input into the migration code to obtain the common image gather. In step 50, the poststack migration result is input into a dip estimation algorithm to estimate the dip angles (.sub.d in Eqn. (2)). In step 60, at least one horizon is selected and at least one point on this horizon is selected. If the obtained common image gather is a common offset gather (92), then in step 96 virtual sources are deployed on the surface and the wavefront construction based ray tracing is performed from the virtual sources to the subsurface; the ray incident angle for each CDP for each virtual source is saved and used to build an angle table. After that, in step 98, for each selected CDP in each image from the common offset gather, the associated takeoff angle is found in the angle table. Then, in step 100, two rays are shot from each of these selected points to the surface by using a ray tracing algorithm. As explained earlier, these are not arbitrary rays. If at step 92 the CIG is a common angle gather as shown in FIG. 5, then the takeoff angle for the ray tracing from the selected CDP point is known; if the CIG is a common offset gather, for each selected CDP point and the specific offset, the takeoff angle already has been found in step 98. With this takeoff angle and the estimated dip angle, two rays are traced from the CDP point according to Snell's law.

(41) The raypath information obtained in step 100 is input to step 110 to construct the kernel matrix for Q tomography. In step 70, the vertical seismic trace is extracted from the images obtained in step 40 for each CDP selected in step 60 and these extracted traces are truncated to contain only the events associated with the selected CDPs. A vertical seismic trace here means the extracted vertical slice from the migration image as shown in FIG. 7 and FIG. 16. Because it is extracted from an image, this trace is in depth domain.

(42) In step 80, the depth-to-time conversion and the wavelet stretching correction are applied on the extracted depth domain seismic traces obtained in step 70. In step 90, the centroid frequency shifts of the extracted traces with respect to the source wavelet centroid frequency are calculated to construct the measurement vector. The measurement vector and the kernel matrix are input into step 120 to formulate the box-constrained tomographic optimization problem and this optimization problem is solved iteratively for 1/Q. If the reconstructed Q model is accepted after the judgment in step 130, the Q tomography process finishes. Otherwise, the user may reselect the horizons and the points on the horizons and repeat steps 70 to 120 until the reconstructed Q model is satisfactory.

EXAMPLES

(43) In this section, a synthetic example of reflection seismic data Q tomography is presented. FIG. 11 shows the velocity model, which is a layered model. The first layer (from the surface to the depth of 1500 m) has a constant velocity of 1500 m/s. The second layer (from the depth of 1500 m to the depth of 2500 m) has a constant velocity of 3000 m/s. The third layer (from the depth of 2500 m to the depth of 3000 m) has a constant velocity of 5000 m. In the first layer, there is an embedded Q anomaly with a constant Q value of 20, as shown in FIG. 14A. Although this Q anomaly is located in a shallow region, the refraction seismic data still cannot be used. The reason is that, in this case, the refractive seismic rays only exist at the surface (i.e., z=0) because the seismic velocity in this example is a constant in the first layer. For a constant seismic velocity model, the refractive rays travel only on the surface.

(44) FIG. 12 shows the migrated reflection seismic data for a specific offset. At a selected CDP point P, a vertical post-migration seismic trace is extracted and then only a short segment of the trace (highlighted line in FIG. 12) is kept to contain only the event associated with CDP P. This segment of the trace (depth domain) is converted to the time domain and the wavelet stretching correction is applied for nonzero offset. After that, the time domain trace is converted to frequency domain. The frequency amplitude spectrum is plotted in FIG. 13. As shown in FIG. 13, in comparison with the source amplitude spectrum, as expected, the frequency content of the post-migration seismic data is different from the source frequency content. The amount of the centroid frequency downshift can be measured as f.sub.c shown FIG. 13. This procedure is repeated until all the f.sub.c are measured to build the measurement vector using formulation (7).

(45) Ray tracing (the procedure shown in FIGS. 7 and 16) is implemented for each selected CDP to build the kernel matrix A. Eventually, the tomographic inversion problem (8) is solved and the Q anomaly is reconstructed as shown in FIG. 14B, which is close to the true Q model.

(46) The foregoing description is directed to particular embodiments of the present invention for the purpose of illustrating it. It will be apparent, however, to one skilled in the art, that many modifications and variations to the embodiments described herein are possible. All such modifications and variations are intended to be within the scope of the present invention, as defined by the appended claims.

REFERENCES

(47) He, Yang and Jun Cai, Q tomography towards true amplitude image and improve sub-karst image, SEG Annual Meeting, Las Vegas (2012). Hu, W., J. Liu, L. Bear, and C. Marcinkovich, A robust and accurate seismic attenuation tomography algorithm, SEG Expanded Abstracts 30, 2727 (2011). Hu, W., Q Tomography Method, PCT Patent Application Publication WO2011/139419 (2011). Hung B., K. F. Xin, S. Birdus, and J. Sun, 3-D tomographic amplitude inversion for compensating transmission losses in the overburden, 70th Meeting, EAGE, H004 (2008). Liao, Q. and G. A. McMechan, Multifrequency viscoacoustic modeling and inversion, Geophysics 61, 1371-1378 (1996). Marfurt, K. J., Robust estimates of 3D reflector dip and azimuth, Geophysics 71, P29-P40 (2006). Morigi, S., L. Reichel, F. Sgallari, and F. Zama, An iterative method for linear discrete ill-posed problems with box constraints, Journal of Computational and Applied Mathematics 198, 505-520 (2007). Pratt, R. G., K. Bauer, and M. Weber, Crosshole waveform tomography velocity and attenuation images of arctic gas hydrates, 73rd SEG, Expanded Abstracts, 2255-2258, Society of Exploration Geophysics (2003). Quan, Y. and J. M. Harris, Seismic attenuation tomography using the frequency shift method, Geophysics 62, 895-905 (1997). Rossi, G., D. Gei, G. Bohm, G. Madrussani, and J. M. Carcione, Attenuation tomography: An application to gas-hydrate and free-gas detection, Geophysical Prospecting 55, 655-669 (2007). Vinje, V., E. Iversen, and H. Gjoystdal, Traveltime and amplitude estimation using wavefront construction, Geophysics 58, 1157-1166 (1993).