The invention relates generally to the field of geophysical prospecting for hydrocarbons and, more particularly, to seismic data processing. Specifically, the invention relates to the technical fields of seismic simulation, reverse time depth migration, and full waveform inversion of seismic data to infer a subsurface model of velocity or other physical property. In addition, many other useful industrial simulators for doing electromagnetic propagation, reservoir simulation and heat flow may match the form needed for this invention to apply. As will be explained later in the invention disclosure, the necessary form corresponds to that of a stationary differential operator equation including mixed or non-mixed terms of spatial and temporal derivatives with coefficients that may vary with space but do not vary in time.
A seismic simulator used to compute either a forward simulation of a source or the adjoint simulation of a recorded wave field is typically implemented using a time stepping algorithm based upon a selected finite difference approximation to either a first or a second time derivative. Most historical implementations have been for 2nd order time stepping (error is proportional to (Δt)2, where Δt is the time step) because that is easy and efficient to implement and requires fewer resources. Using 2nd order time stepping provides a result with temporal dispersion artifacts. Higher-order approximations are better, because the error for approximation of order n is proportional to (Δt)n which→0 as n→∞ for Δt<1. But any finite-order approximation suffers from some degree of temporal dispersion. If the temporal dispersion is not corrected, the application of this type of simulator for forward simulation or to compute Reverse Time Depth Migration (RTM) images and Full Waveform Inversion (FWI) gradients and Hessians will have errors that degrade the value of these products for petroleum exploration and geophysical prospecting.
In one embodiment, the invention is a method for prospecting for hydrocarbons. Measured seismic data are obtained. Corresponding simulated seismic data are computer-generated using a finite-difference, time-stepping algorithm that approximates a time derivative operator to a selected order of approximation. The simulated seismic data are used to perform full-waveform inversion or reverse-time migration of the measured seismic data, wherein temporal numerical dispersion corresponding to the selected order of approximation is (i) removed from the simulated seismic data or (ii) introduced into the measured seismic data by steps including performing a Fourier transform in time on (i) the simulated or (ii) the measured seismic data, then resampling the transformed seismic data in frequency domain, and then performing an inverse Fourier transform from frequency domain back to time domain. The resampling utilizes a property of a class of stationary finite-difference operators whereby, in frequency domain, an aspect of the temporal numerical dispersion is that a desired numerical solution for a given frequency is computed at an incorrect frequency, and the resampling uses a mapping relationship that maps the incorrect frequency to the given frequency. The full-waveform-inverted seismic data or the reverse-time-migrated seismic data may then be used in known methods to prospect for hydrocarbons. The resampled frequency-domain seismic data may be scaled with a frequency-dependent scaling factor before performing the inverse Fourier transform back to time domain.
The present invention and its advantages will be better understood by referring to the following detailed description and the attached drawings, in which:
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.
To illustrate the problem that the present invention solves, consider an earth model consisting of a simple two-dimensional half-space with a free surface boundary on top as shown in
As stated above, a finite difference operator used for computing a temporal derivative has approximations that create numerical errors when used to solve partial differential equations. It is a realization of the present invention that for certain types of stationary differential operators that incorporate approximate temporal derivative operators, the correct solution is computed by the approximate equation, but at the wrong frequency. The invention uses this feature to correct an approximate solution into a more correct solution by a resampling operation in the frequency domain.
Virtually all of the time-domain forward and adjoint wave simulation algorithms used for seismic simulation, seismic migration and seismic full waveform inversion correspond to the type of stationary differential operators to which this invention applies. Electromagnetic equations and heat flow problems can also be formulated in a way such that the theory developed here applies. Thus, a wide range of time simulation processes may be able to use this method.
The form necessary for the invention to apply corresponds to that of a stationary differential operator equation including mixed or non-mixed terms of spatial and temporal derivatives with coefficients that may vary with space but do not vary in time. The differential operators for time stepping can be either explicit or implicit. For seismic simulation, RTM, and FWI applications, the basic import of this is that the earth model properties do not change during one seismic simulation. Alternatively, this assumption would be violated if the earth model were to change during the simulation. One example of a violation of this assumption would be if moving water waves on the air-water boundary were to make the earth model properties change with time during the simulation.
The present invention differs from Stork (2013) in many ways: (a) how the temporal numerical dispersion is corrected or applied (The invention uses Fourier-domain resampling instead of Stork's choice of time-domain filter banks to implement temporal numerical dispersion corrections.); (b) by extending the application to FWI objective function, gradient and Hessian computations; (c) by modifying dispersion from approximate derivative operators of any order of accuracy to any other order of accuracy and thereby enabling the match of temporal operator order to spatial operator, which enables (d) the application of temporal dispersion correction operators or their inverse to less perfect spatial operators. Shorter spatial operators enable more efficient halo exchanges for parallel domain-decomposed computations. The combination of all of these aspects can lead to improved efficiencies for a given level of accuracy for simulators, RTM applications and FWI applications.
The present invention differs from the patent application by Zhang et al. (2012) because the key step in the process is resampling in the frequency domain rather than filtering in the frequency domain. Filtering typically implies convolution but this invention is instead based upon resampling to change variables. The invention also has a more general range of applicability to a specific class of stationary differential operator equations. As a result the present invention applies to time-stepping differential equations in other fields, e.g. heat flow and reservoir simulation, in addition to seismic simulation and seismic migration. Zhang does not teach us to apply his method to FWI or other seismic applications.
Computing ωapprox (ω) for an example finite difference operator
A finite difference solution applies approximate derivative operators to solve a problem. The approximate derivative operators in the time domain can be Fourier transformed to the frequency domain and compared to the exact form of the derivative, which is iω.
The Fourier transform F(ω) of an explicit finite difference operator with coefficients fj at times tj is the Fourier transform of a digital filter. It is computed by the following equation.
On a regular grid, a second derivative finite difference operator is a symmetric digital filter.
Consider for example the explicit convolutional centered second derivative in time finite difference operator with norder+1 coefficents. Here, the coefficients aj represent the zero lag and positive lags of a symmetric filter used to approximate the second derivative operator. Those coefficients can be Taylor series coefficients or be optimized coefficients designed to fit a specified bandwidth with high accuracy. The exact operator in the Fourier domain would be −ω2. On a regular grid, a second derivative finite difference operator is a symmetric digital filter. The Fourier transform of a symmetric digital filter is a cosine transform
where norder is the order of the finite difference approximation and Δt is the finite time step. Comparing the exact operator in the Fourier domain to the approximate finite difference operator in the Fourier domain leads to the following relationship between approximate angular frequency and true angular frequency.
The error for this order of approximation will be proportional to (Δt)norder when equation (1.3) applies to explicit temporal derivative operators derived from Taylor series expansions.
The second order approximation for the angular frequency made while doing a finite difference temporal second derivative corresponds to the following in the frequency domain as a function of the true angular frequency ω and the time step increment Δt.
The inverse mapping from approximate angular frequency to true angular frequency can also be made.
The examples above have been for explicit convolutional-style temporal derivative operators. This type of mapping can also be done for implicit operators as described by Crank and Nicolson (1947) that could be designed to have properties of unconditional stability for time stepping with large time increments. The most common Crank-Nicolson approach as implemented by Claerbout (1985) for a wave equation would use the second-order bilinear Z transform. For that case, this mapping would apply:
Consider a differential operator L operating on a function u equal to a broad band source term s(x,t). Choose L to be a linear sum of terms, each scaled by the spatially varying coefficients and/or by spatial derivatives to any order and/or by mixed spatial derivatives to any order and/or by time derivatives to any order. The operator varies spatially with the kth operator coefficient term ck optionally a function of position x. However the operator is stationary with respect to time in that operator coefficients ck are not time dependent. The operator L may be dependent upon temporal derivatives of any order but not on time explicitly.
If the operator L contains time derivatives of any order but no coefficients that vary with time, then the frequency-domain equivalent operator {tilde over (L)} retains a similar form with each time derivative replaced by iω. Such an operator L is of the type of stationary differential operator to which the present inventive method applies.
Then if U(x,ω) is the temporal Fourier transform of the solution wavefield u(x,t) and S(x,ω) is the temporal Fourier transform of the source term s(x,t),
The Helmholtz equation (Morse and Feschbach, 1953) is an example of a differential equation with the form given in equation 1.8.
The following two operator equations have identical solutions for ω2=ω1. Mathematically, this is a trivial statement since the two equations are identical except that the variables have been renamed. The key point is that ω1 and ω2 can represent different temporal derivative operators. This is how one can recognize that the approximate solution contains within it the true solution at the wrong frequency.
In practice, the solution U1(x,ω1) to the operator ω2(ω1) is usually computed for the following equation.
If ω1=ωexact and ω2(ω1) is the Fourier transform of a temporal finite difference operator, then U1(x,ω1) would represent the solution to a specific temporal finite difference approximation. That solution can be mapped by a resampling and scaling operator in the frequency domain into a solution U2(x,ω2) with exact temporal derivatives consistent with equation (1.12) as follows.
Basic steps in this embodiment of the present inventive method are shown in the flow chart in
Equation 1.12 is quite general. However, one limitation is that if multiple seismic sources are being simultaneously simulated in a single simulation, as in simultaneous-source FWI or simultaneous source RTM (see, for example, U.S. Pat. No. 8,121,823 to Krebs, et al.), then the simultaneous sources need to have the same source time functions to within a scale factor. They can vary by scale factors of +1 and −1. They may be in multiple spatial locations.
In summary, the invention applies to correcting or modifying temporal numerical dispersion characteristics associated with solutions to stationary differential operator equations of the style discussed above. In one of its
(1) Begin with a solution u1(x,t1) to a stationary differential operator equation Lu1(x,t1)=s(x,t1) using a known approximate or true temporal derivative operator. This solution is the simulated seismic data 51 in the flow chart of
(2) In step 52, do a temporal Fourier transform from time t1 to angular frequency ω1 converting u1(x,t1) into U1(x,ω1).
(3) Take ω1(ω) to be a bijective function over some range of ω that relates the approximate or true temporal derivative operator associated with ω1 used for solution U1(x,ω1) to the true ω used for the exact solution U(x,ω). This can be readily computed by taking a Fourier transform of the approximate operator and comparing that to the Fourier transform of an exact operator as discussed in a prior section. Likewise, take ω2(ω) to be a bijective function over some range of ω that relates the approximate or true temporal derivative operator ω2 used for solution U2(x,ω2) to the true ω used for the exact solution U(x,ω). The invention converts an available solution U1(x,ω1) with the temporal numerical dispersion characteristics of temporal derivative operators associated with ω1(ω) to match a desired solution U2(x,ω2) with the temporal numerical dispersion characteristics of temporal derivative operators associated with ω2(ω) by finding the bijective relationship ω2(ω1) and resampling U1(x,ω1) to get the desired values of U2(x,ω2) via equation (1.12). This is step 53.
(4) In step 54, perform an inverse temporal Fourier transform of U2(x,ω2) to get the temporal dispersion corrected or modified result u2(x,t2) in the time domain (55).
Next, some aspects of the invention are described in more detail. The embodiment of the present invention implementing temporal numerical dispersion corrections, via resampling in the temporal Fourier domain, to simulated seismic data can be applied as a post-processing single-seismic-trace-at-a-time process applied either within the seismic simulator or as a separate seismic processing application. A flow chart showing basic steps in this embodiment of the present inventive method is given in
For a second embodiment of the present invention that applies temporal numerical dispersion via resampling in the temporal Fourier domain to field data to be input to RTM or FWI gradient or Hessian computation, this may be done as a preprocessing step. Basic steps in this embodiment are shown in
It may be noted that full waveform inversion can be corrected for temporal dispersion using either embodiment of the present inventive method, i.e. that of
The present invention can be used to either apply or remove all temporal numerical dispersion as shown in
The present invention can be applied to change the temporal numerical dispersion characteristics of field or synthetic seismic data from any operator order to any other operator order. Equation 1.3 gives the relationship between true frequency and approximate frequency for explicit centered temporal finite difference operators on a regular grid. This relationship is objective over the specific range of frequency and time increment parameters of interest. Therefore the inverse relationship can be found. Therefore one can map any operator to the true operator and then back to another approximate operator. A general expression for the inverse may be difficult to write for some operator choices, but a computer can easily tabulate these and then look up values in the table to solve the inverse relationship.
The effect of a second order temporal finite difference operator on phase velocity has been described by Fei (1994) and is shown in
A key aspect of the methods disclosed herein, as indicated by equation (1.12), is that the seismic data can be advantageously modified from one form of ω operator to another form by resampling in the frequency domain. The input and output ω operators can be exact or approximate, and if approximate, they can be explicit or implicit.
A test of the present invention is shown in
An ideal simulation result with no temporal dispersion would be three spikes, i.e. this would be the ideal result for
Low-frequency FWI is less affected by temporal numerical dispersion than broad-band FWI because temporal numerical dispersion is less important at low frequencies. High resolution images of the subsurface require broad band FWI, and temporal numerical dispersion corrections become more important. Temporal numerical dispersion corrections become very important when making accurate ties between inverted FWI earth model parameters and well logs.
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.
This application claims the benefit of U.S. Provisional Patent Application 62/009,593, filed Jun. 9, 2014, entitled A Method for Temporal Dispersion Correction for Seismic Simulation, RTM and FWI, the entirety of which is incorporated by reference herein.
Number | Date | Country | |
---|---|---|---|
62009593 | Jun 2014 | US |