When performing surveying of a subterranean structure for identifying subterranean bodies of interest (e.g., hydrocarbon reservoirs, fresh water aquifers, gas injection zones, etc.), data points are collected by survey receivers in response to a stimulus. The survey receivers can be electromagnetic (EM) receivers or seismic receivers. The stimulus can be produced by a source, such as an EM source or a seismic source.
Imaging of a subterranean structure is based on building a forward modeling operator. The forward modeling operator relates physical properties of the subterranean structure to measurement data collected by a survey spread that includes survey receivers. Once formed, the forward modeling operator can be used to build an imaging operator for predicting a quantity of interest associated with a subterranean structure, by applying the imaging operator on input data, such as measured data and/or prior information.
Challenges associated with computing imaging operators include accuracies of forward modeling operators and complexities of deriving such imaging operators.
In general, according to an embodiment, a method is provided for deriving at least one weight to build an imaging operator for imaging a subterranean structure. The at least one weight is computed according to a criterion, where the criterion is selected from among: (1) the at least one weight is configured to provide an imaging operator that creates true amplitude image; (2) the at least one weight is configured to provide an imaging operator that minimize mean square error; and (3) the at least one weight is configured to provide an imaging operator that approximates a least squares inverse of a forward modeling operator (to image a subterranean structure). The computed at least one weight is then used to build the imaging operator which is used to produce an image of the subterranean structure.
Other or alternative features will become apparent from the following description, from the drawings, and from the claims.
Some embodiments are described with respect to the following figures:
In general, filtering techniques are provided for producing an imaging operator for use in imaging a subterranean structure. The imaging operator is defined at least in part by at least one weight. The at least one weight is computed according to a design criterion, which is selected from among the following: (1) the at least one weight is configured to provide an imaging operator that creates a true amplitude image; (2) the at least one weight is configured to provide an imaging operator that minimizes mean square error; and (3) the at least one weight is configured to provide an imaging operator that approximates a least squares inverse of a forward modeling operator. The at least one weight is used to build the imaging operator that is used to produce an image of the subterranean structure.
Use of an imaging operator is illustrated in a simplified diagram shown in
The techniques of producing an imaging operator according to some embodiments can be used in any one of several contexts, including the following: Gaussian beam migration, prestack depth migration (also referred to as Kirchhoff migration); and others.
In Gaussian beam migration, imaging of a subterranean structure is obtained by using Gaussian beams as a forward modeling operator. More precisely, the scattering potential, or reflectivity, of the subterranean structure is imaged using a Gaussian beams representation of the forward modeling operator. Gaussian beam migration enables imaging of relatively complex geologic structures without dip limitation, and with accuracy comparable to that of wave equation migration and efficient comparable to that of Kirchhoff migration. In some embodiments, a midpoint-offset domain Gaussian beam migration filtering technique is provided, which provides a filter that can produce a true amplitude image.
A true amplitude image of a quantity of interest is the image obtained from measurement data whose edge locations and strengths agree with an ideal/real quantity. The ideal/real quantity is the actual quantity that exists in the subterranean structure. Examples of a quantity of interest associated with a subterranean structure include formation reflectivity, scattering potential, resistivity, rock porosity, and so forth.
Designing weights according to some embodiments when applied in the prestack depth migration (Kirchhoff migration) context can produce a scattering potential, or similarly reflectivity image with true amplitude when composed with an adjoint of the forward modeling operator. The imaging operator that is produced by weight design techniques acts on measurement data (collected using survey receivers) based on local slant-stack transforms in source and receiver coordinates and based on convolution in time. The local slant-stack and convolution arrangement naturally provides a framework for parallel computations and for relatively efficient implementations. In the prestack depth migration context, a filter produced according to some embodiments can also provide a true amplitude image. For prestack depth migration, the filter can be a spatially varying space-wavenumber-time-frequency filter.
The survey source 106 and survey receiver 108 can be connected to a controller 110 which can be a central recording station. The controller 110 can be implemented with a computer system as well as with other components. The controller 110 has processing software 112 that can apply filtering techniques according to some embodiments. The processing software 112 is executable on one or plural processors 114. The one or plural processors 114 are connected to storage media 116, which can be implemented with disk-based storage media and/or integrated circuit or semiconductor storage media.
The storage media 116 stores measurement data 118 acquired by survey receivers, such as survey receiver 108. Moreover, the storage media 116 stores one or more filters 120 produced by filtering techniques according to some embodiments. The filter(s) 120 can be produced by the processing software 112, based on the measurement data 118.
The following discussion has two sections: a first section relates to Gaussian beam migration, while a second section relates to prestack depth migration (Kirchhoff migration).
Gaussian Beam Migration
In accordance with some embodiments, a weighted adjoint of the forward modeling operator is the imaging operator and provides true amplitude Gaussian beam migration of midpoint-offset sorted measurement data. A weighted adjoint operator designed according to some embodiments has the characteristic that the composition of the weighted adjoint operator with the forward modeling operator from right will approximate the identity operator. In other words, assume F represents the forward modeling operator, Fw* represents the weighted adjoint, and Iw represents the imaging operator. Then Iw=Fw* and Fw·F would produce the identity operator (I). The technique according to some embodiments is an analytic technique that can accommodate statistical information within its framework. The technique is optimal in the least square sense in the absence of noise in measurement data, and optimal in the minimum mean square error sense in the presence of noise in measurement data. Also, the notion of generalized semblance can be used in some embodiments. The MMSE (minimum mean square error) inverse of the forward modeling operator is computed using generalized semblance. When generalized semblance is used, the statistics relating to noise are no longer required to compute the MMSE inverse, which simplifies computations. Moreover, the generalized semblance can be treated as a depth oriented wavenumber gather which may be used as a focusing tool for tomography.
Gaussian beam migration involves expressing a Green's function using a discrete number of Gaussian beams. The following three approximations are used in some embodiments.
First, the Green's function G(r,r′,ω) is approximated in terms of the Gaussian beams uGB(r,r′,p′,ω):
where uGB(r,r′,p′,ω) is given by
uGB(r,r′,p,ω)=A(r,r′,p)exp(iωT(r,r′,p)). (Eq. 2)
Here A(r,r′,p) is the complex amplitude and T(r,r′,p) is the complex travel time. The vectors r and r′ are coordinates in space (r,r′εR3), ω is frequency, p and p′ are initial beam directions (p,p′εS2).
Second, for certain lattice points L, the following condition applies:
Third, for small r″,
uGB(r,r′+r″,p,ω)≈uGB(r,r′,p,ω)e−ip·r″. (Eq. 4)
The measurement data at receiver at location rd due to a source at location at rs is modeled by
D(rd,rs,t)=∫ω2P(ω)e−iωt∫dr′f(r′)G(r′,rs,ω)G(rd,r′,ω)dω, (Eq. 5)
where f(r) is referred to as the scattering potential, and P(w)=∫p(t)e−iωtdt is the Fourier transform of the source wavelet. Using the midpoint-offset representation
rm=(rd+rs)/2, h=(rd−rs)/2, (Eq. 6)
and the Gaussian beam representation of the Green's function, the measurement data can be expressed as:
Eq. 7 can be approximated by:
The forward modeling operator F is defined by the right hand side of Eq. 8, i.e., Ff(rm,h,t)≈D(rm+h,rm−h,t). For true amplitude imaging, techniques according to some embodiments design a left inverse of F.
The strategy for finding a left inverse of F and forming Gaussian beam migration images is to design a weighted version of the adjoint of the forward operator F. In this regard, the weighted adjoint operator Fw* (the imaging operator designed according to some embodiments) is defined by
where M(rd,rs,t) denotes measurement data. Note that for wL=1, Fw* is equal to the adjoint of the forward modeling operator F.
The weights wL are calculated as follows:
In the noise-free scenario, it is assumed that the measurement data is given by Eq. 8:
M(rm,h,t)=D(rm+h,rm−h,t). (Eq. 11)
Thus wL (weights of the imaging operator Iw) is computed such that Fw* is an approximate (least squares) inverse of F, i.e., Fw*F≈I. In some embodiments, with the choice of
wc(ζ,r)=∫└∫[{tilde over (b)}(r+r′,ζ′)]−1e−iζ·r′dr′, (Eq. 12)
Fw*Ff≈f is provided. Here
with Âh(ω,r,L′,qm)=Ah(r,L′,qm)exp(ωThℑ(r,L′,qm)). Also Ah(r′,L′,qm) and Th=+ are the complex amplitudes and travel times, respectively.
Thus, in the noise-free scenario in which measurement data is noise free, the weighted adjoint operator (as represented by Fw*) is a least squares inverse of the forward modeling operator (F).
In the presence of noise (noisy scenario), it is assumed that the measurement data is corrupted with additive noise and obey the following model:
M(rm,h,t)=D(rm+h,rm−h,t)+n(rm,h,t). (Eq. 14)
Similar to the noiseless case, the image is formed using weighted adjoint operator Fw* (the imaging operator Iw according to some embodiments), where this time the weights are designed such that the minimum mean square error J(w):
J(w)=∫E[|Fw*M(r)−f(r)|2]dr, (Eq. 15)
is minimized. Here ∥f∥22=∫|f(r)|2dr denotes the L2 norm and E denotes statistical expectation.
It is assumed that the noise is zero mean and is uncorrelated from D and can be modeled as a second order stationary process in time that is locally uncorrelated in space. Furthermore, it is assumed that the scattering potential is second order stationary with power spectrum Sf.
With the choice of
J(w) is minimized. Here {tilde over (b)}(r,ζ) is defined as in Eq. 13 and
where Sn is referred to as the localized power spectrum of n (noise). One can also think {tilde over (b)}(r,ζ) as the propagated noise spectrum observed at r.
Alternatively, one can use the notion of generalized semblance to compute Eq. 16 and eliminate the need for {tilde over (b)}n(r,ζ) and Sf(ζ) as follows:
where S(r,ζ) is the generalized semblance defined by
with w(r′,ζ′)(ξ,r)=exp[i(ξ−ξ′)·(r′+r)].
When generalized semblance is used, the statistics relating to noise are no longer required to compute the MMSE (minimum mean square error) inverse. This allows for more efficient computation of the weighted adjoint operator (Fw*).
The foregoing has discussed formation of the weighted adjoint operator Fw*. The following describes usage of the weighed adjoint operator once it has been derived.
Using Eq. 9, the image {circumflex over (f)}(r) to be produced is derived from input measurement data M(r) and the weighted adjoint operator Fw* according to the following:
which, by using the identity and change of variables rm→L+rm, pm=ps+pd and ph=pd−ps, becomes
Eq. 21 is used to form a target image {circumflex over (f)}h(r) representing a subterranean structure. The image {circumflex over (f)}h(r) is produced using the following procedure:
At this point one can either perform velocity analysis over the image volume {circumflex over (f)}h(r) by looking at its flatness with respect to surface offset h, or form the image by stacking:
{circumflex over (f)}(r)=∫{circumflex over (f)}h(r)dh.
According to the foregoing, a true amplitude Gaussian beam migration technique is provided, which relies on the design of weighted adjoint operators as imaging operators for the inversion of the forward modeling operator. The weighted adjoint operators are designed for both the noise-free and noisy scenarios. In the presence of noise, the weighted adjoint operator minimizes the mean square error, whereas in the noise-free scenario, the weighted adjoint operator simplifies to the least squares inverse of the forward modeling operator.
Prestack Depth Migration
In accordance with some embodiments, in the prestack depth migration context, a true-amplitude prestack depth migration technique is provided that includes a spatially varying space-wavenumber-time-frequency filtering technique. The imaging operator for prestack depth migration according to some embodiments is designed as follows.
Such an imaging operator is derived from a composition of a filtering operator with the adjoint of the forward modeling operator and has the characteristic that the composition of the imaging operator with the forward modeling operator from right will approximate the identity operator. In other words, assume that F represents forward modeling operator, Iw represents the imaging operator, F* represents the adjoint of F, and W represents the filtering operator. Then Iw=F*·W and Iw·F would produce the identity operator (I). The imaging operator formulation also can include the transmitted source wavelet in an approximation of the integral in computation of the adjoint of the forward modeling operator.
Let F be the forward modeling operator that maps scattering potential f to pressure data p:
Ff(xs,xr,t)=m(xs,xr,t)=∫ω2P(ω)A(xs,xr,y)ei2πω(t−φ(x
In the above, m(x,x′,t) represents ideal noise-free measurement data, xs represents a coordinate (position) of the survey source, xr represents a coordinate (position) of the survey receiver, t represents time, A(x,x′,y) represents the adjoint of the combined amplitude function A(x,x′,y), and P*(ω) represents the complex conjugate of P(ω), which is expressed as P(ω)=∫p(t)e−2πωtdt, where p(t) is the source wavelet. (The symbol * is used for complex conjugation as well as adjoint of operators—the meaning will be clear from the context.)
The filtering operator W is designed such that
f(y)≈F*Wm(y). (Eq. 24)
Eq. 25 defines the filtering operator W as a spatially dependent space-time-wavenumber-frequency filter. Due to the windowing functions the filtering is localized in source and receiver locations.
Furthermore, it is assumed that the filter (also referred to as a “weight”) w(xs,xr,t,ks,kr,ω,y) is a reparameterization of a standard-symbol of a pseudo-differential operator:
w(xs,xr,t,ks,kr,ω,y)=w(x
where w(x
For ease of manipulation and computation of the filter w(xs,xr,t,ks,kr,ω,y), a compound-symbol representation of wΨ(ξ,y) is used:
wΨ(ξ,y)=∫v1(y−y′,y′)ei2πξ·(y−y′)dy′.
With the choice of
v1(y−y′,y′)=∫[∫ei2π5·y″[∫e−i2πξ·y″a(xs,xr,ω,y)dxsdxrdω]dy″]−1e−i2πk·(y−y′)dk, (Eq. 28)
where
a(xs,xr,ω,y)=|A(xs,xr,y)|2ω4|P(ω)|2, (Eq. 29)
W satisfies Eq. 24 and gives a true-amplitude image.
Although the imaging operator is designed for adjoint imaging, the imaging operator can equivalenity be designed for matched-filtered backprojection imaging, i.e., BpWm(y)≈f(y), where
Bpm(y)=∫m(xs,xr,τ+t=φ(xs,xr,y))p*(τ)dτdxsdxr, (Eq. 30)
In this case, it is sufficient to set a(xs,xr,ω,y)=A(xs,xr,y)ω2|P(ω)|2 in the computation of the filter.
Once the filtering operator W is derived as discussed above, then the imaging operator can be computed by combining the filtering operator W with the forward modeling operator F, as Iw=F*·W.
To apply the imaging operator Iw designed above,
Instructions of machine-readable instructions described above (including processing software of
Data and instructions are stored in respective storage devices, which are implemented as one or more computer-readable or computer-usable storage media. The storage media include different forms of memory including semiconductor memory devices such as dynamic or static random access memories (DRAMs or SRAMs), erasable and programmable read-only memories (EPROMs), electrically erasable and programmable read-only memories (EEPROMs) and flash memories; magnetic disks such as fixed, floppy and removable disks; other magnetic media including tape; optical media such as compact disks (CDs) or digital video disks (DVDs); or other types of storage devices. Note that the instructions of the software discussed above can be provided on one computer-readable or computer-usable storage medium, or alternatively, can be provided on multiple computer-readable or computer-usable storage media distributed in a large system having possibly plural nodes. Such computer-readable or computer-usable storage medium or media is (are) considered to be part of an article (or article of manufacture). An article or article of manufacture can refer to any manufactured single component or multiple components.
In the foregoing description, numerous details are set forth to provide an understanding of the subject disclosed herein. However, implementations may be practiced without some or all of these details. Other implementations may include modifications and variations from the details discussed above. It is intended that the appended claims cover such modifications and variations.
Number | Name | Date | Kind |
---|---|---|---|
6027447 | Li | Feb 2000 | A |
6904368 | Reshef et al. | Jun 2005 | B2 |
6996470 | Kamps | Feb 2006 | B2 |
8060312 | Wang et al. | Nov 2011 | B2 |
20090184958 | Osypov et al. | Jul 2009 | A1 |
20110096627 | Hill | Apr 2011 | A1 |
Entry |
---|
Albertin, Uwe, et al., True-amplitude beam migration, SEG Int'l Exploration and 74th Annual Meeting, Oct. 10-15, 2004. |
Bleistein, N., et al. Mathematics of Multidimensional Seismic Imaging Migration, and Inversion, 2000, Springer-Verlan New York, Inc., New York, NY, United States of America. |
Cerveny, V., Seismic Ray Theory, 2001, Cambridge University Press, New York, NY, United States of America. |
Ferreira, Carlos A. S., et al., SEG/New Orleans 2006 Annual Meeting, 2006, pp. 2614-2618. |
Gray, Samuel H., Gaussian beam migration of common-shot records, Geophysics,Jul. 6, 2005, pp. S71-S77, vol. 70, No. 4. |
Haldorsen, jAKOB B.U., Multichannel Wiener deconvolution of vertical seismic profiles, Geophysics, Oct. 1994, pp. 1500-1511, vol. 59, No. 10. |
Hill, N. Ross, Gaussian beam migration, Geophysics, Nov. 1990, pp. 1416-1428, vol. 55, No. 11. |
Hill, N. Ross, Prestack Gaussian-beam depth migration, Geophysics, Jul. 8, 2001, pp. 1240-1250, vol. 66, No. 4. |
Liu, Jonathan, Multi-arrival Kirchhoff beam migration, SEG Las Vegas 2008 Annual Meeting, 2008, pp. 2311-2315. |
Nowack, Robert L., Gaussian beam migration for sparse common-shot and common receiver data, SEG, Expanded Abstracts, 2003, pp. 1114-1117, vol. 22(1). |
Prostasov, M. I., True-Amplitude Seismic Imaging, , Doklady Earth Sciences, 2006, pp. 441- 445, vol. 407A, No. 3. |
Symes, William W., et al., Migration velocity analysis and waveform inversion, 2008, Geophysical Prospecting, pp. 765-790, vol. 56. |
Virgilio, Massimo, Gaussian beam migration: prestack, common-shot, TTI, true amplitude in angular domain, 2008, SEG, Expanded Abstracts, SEG, Expanded Abstracts, pp. 2236-2240, 27(1). |
Wang, Bin, er al., Advances in velocity model-building technology for subsalt imaging, Geophysics, 2008, pp. VE173-VE181. |
Yilmaz, Oz, Seismic Data Analysis, SEG, 2000. vol. 1-vol. 2. |
Beylkin, G., Imaging of discontinuities i n the inverse scattering problem by inversion of a causal generalized Radon transform, J. Math. Phys. Jan. 1985, pp. 99-108, vol. 26. |
De Hoop, Maarten V., et al., GRT/AVA migration/inversion in anisotropic media, SPIE, 1994, pp. 15-27, vol. 2301. |
Miller, D. E., et al., Multiparameter inversion, dip-moveout, and the generalized radon transform, Geophysical Inversion, SIAM, 989/Sep. 27-29, 1990. |
Miller, D.E., et al., A new slant on seismic imaging, Migartion and Integral geometry, Geophysics,Jul. 1987 pp. 943-964, vol. 52, No. 7. |
Brandsberg-Dahl, Sverre, et al., Focusing in dip and AVA compensation on scattering-angle/azimuth common image gathers, Geophysics, Jan. 2, 2003, pp. 232-254, vol. 68, No. 1. |
Stein, Elias M., Harmonic Analysis: Real-Variable Methods, Orthogonality and Oscillatory Integrals, 1993, Princeton University Press, Princeion, NY, United States of America. |
Number | Date | Country | |
---|---|---|---|
20110317934 A1 | Dec 2011 | US |