Patent Application
20120073824
The invention relates generally to the field of geophysical prospecting, and more particularly to geophysical data processing. Specifically, the invention is a method for inversion of data acquired from multiple geophysical sources such as seismic sources, involving geophysical simulation that computes the data from many simultaneously-active geophysical sources in one execution of the simulation.
Geophysical inversion [1,2] attempts to find a model of subsurface properties that optimally explains observed data and satisfies geological and geophysical constraints. There are a large number of well known methods of geophysical inversion. These well known methods fall into one of two categories, iterative inversion and non-iterative inversion. The following are definitions of what is commonly meant by each of the two categories:
Non-iterative inversion—inversion that is accomplished by assuming some simple background model and updating the model based on the input data. This method does not use the updated model as input to another step of inversion. For the case of seismic data these methods are commonly referred to as imaging, migration, diffraction tomography or Born inversion.
Iterative inversion—inversion involving repetitious improvement of the subsurface properties model such that a model is found that satisfactorily explains the observed data. If the inversion converges, then the final model will better explain the observed data and will more closely approximate the actual subsurface properties. Iterative inversion usually produces a more accurate model than non-iterative inversion, but is much more expensive to compute.
Iterative inversion is generally preferred over non-iterative inversion, because it yields more accurate subsurface parameter models. Unfortunately, iterative inversion is so computationally expensive that it is impractical to apply it to many problems of interest. This high computational expense is the result of the fact that all inversion techniques require many compute intensive simulations. The compute time of any individual simulation is proportional to the number of sources to be inverted, and typically there are large numbers of sources in geophysical data, where the term source as used in the preceding refers to an activation location of a source apparatus. The problem is exacerbated for iterative inversion, because the number of simulations that must be computed is proportional to the number of iterations in the inversion, and the number of iterations required is typically on the order of hundreds to thousands.
The most commonly employed iterative inversion method employed in geophysics is cost function optimization. Cost function optimization involves iterative minimization or maximization of the value, with respect to the model M, of a cost function S(M) which is a measure of the misfit between the calculated and observed data (this is also sometimes referred to as the objective function), where the calculated data are simulated with a computer using the current geophysical properties model and the physics governing propagation of the source signal in a medium represented by a given geophysical properties model. The simulation computations may be done by any of several numerical methods including but not limited to finite difference, finite element or ray tracing. The simulation computations can be performed in either the frequency or time domain.
Cost function optimization methods are either local or global [3]. Global methods simply involve computing the cost function S(M) for a population of models {M1, M2, M3, . . . } and selecting a set of one or more models from that population that approximately minimize S(M). If further improvement is desired this new selected set of models can then be used as a basis to generate a new population of models that can be again tested relative to the cost function S(M). For global methods each model in the test population can be considered to be an iteration, or at a higher level each set of populations tested can be considered an iteration. Well known global inversion methods include Monte Carlo, simulated annealing, genetic and evolution algorithms.
Unfortunately global optimization methods typically converge extremely slowly and therefore most geophysical inversions are based on local cost function optimization. Algorithm 1 summarizes local cost function optimization.
This procedure is iterated by using the new updated model as the starting model for another gradient search. The process continues until an updated model is found that satisfactorily explains the observed data. Commonly used local cost function inversion methods include gradient search, conjugate gradients and Newton's method.
Local cost function optimization of seismic data in the acoustic approximation is a common geophysical inversion task, and is generally illustrative of other types of geophysical inversion. When inverting seismic data in the acoustic approximation the cost function can be written as:
where:
The gathers can be any type of gather that can be simulated in one run of a seismic forward modeling program. Usually the gathers correspond to a seismic shot, although the shots can be more general than point sources. For point sources the gather index g corresponds to the location of individual point sources. For plane wave sources g would correspond to different plane wave propagation directions. This generalized source data, ψobs, can either be acquired in the field or can be synthesized from data acquired using point sources. The calculated data ψcalc on the other hand can usually be computed directly by using a generalized source function when forward modeling. For many types of forward modeling, including finite difference modeling, the computation time needed for a generalized source is roughly equal to the computation time needed for a point source.
Equation (1) can be simplified to:
where the sum over receivers and time samples is now implied and,
δ(M,wg)=ψcalc(M,wg)−ψobs(wg) (Eqn. 3)
Inversion attempts to update the model M such that S(M) is a minimum. This can be accomplished by local cost function optimization which updates the given model M(k) as follows:
M
(k+1)
=M
(k)−α(k)∇MS(M) (Eqn. 4)
where k is the iteration number, α is the scalar size of the model update, and ∇MS(M) is the gradient of the misfit function, taken with respect to the model parameters. The model perturbations, or the values by which the model is updated, are calculated by multiplication of the gradient of the objective function with a step length α, which must be repeatedly calculated.
From equation (2), the following equation can be derived for the gradient of the cost function:
So to compute the gradient of the cost function one must separately compute the gradient of each gather's contribution to the cost function, then sum those contributions. Therefore, the computational effort required for computing ∇MS(M) is Ng times the compute effort required to determine the contribution of a single gather to the gradient. For geophysical problems, Ng usually corresponds to the number of geophysical sources and is on the order of 10,000 to 100,000, greatly magnifying the cost of computing ∇MS(M).
Note that computation of ∇MW(δ) requires computation of the derivative of W(δ) with respect to each of the N model parameters mi. Since for geophysical problems N is usually very large (usually more that one million), this computation can be extremely time consuming if it had to be performed for each individual model parameter. Fortunately, the adjoint method can be used to efficiently perform this computation for all model parameters at once [1]. The adjoint method for the least squares objective function and a gridded model parameterization is summarized by the following algorithm:
While computation of the gradients using the adjoint method is efficient relative to other methods, it is still very costly. In particular the adjoint methods requires two simulations, one forward in time and one backward in time, and for geophysical problems these simulations are usually very compute intensive. Also, as discussed above, this adjoint method computation must be performed for each measured data gather individually, increasing the compute cost by a factor of Ng.
The compute cost of all categories of inversion can be reduced by inverting data from combinations of the sources, rather than inverting the sources individually. This may be called simultaneous source inversion. Several types of source combination are known including: coherently sum closely spaced sources to produce an effective source that produces a wavefront of some desired shape (e.g. a plane wave), sum widely spaces sources, or fully or partially stacking the data before inversion.
The compute cost reduction gained by inverting combined sources is at least partly offset by the fact that inversion of the combined data usually produces a less accurate inverted model. This loss in accuracy is due to the fact that information is lost when the individual sources are summed, and therefore the summed data does not constrain the inverted model as strongly as the unsummed data. This loss of information during summation can be minimized by encoding each shot record before summing. Encoding before combination preserves significantly more information in the simultaneous source data, and therefore better constrains the inversion [4]. Encoding also allows combination of closely spaced sources, thus allowing more sources to be combined for a given computational region. Various encoding schemes can be used with this technique including time shift encoding and random phase encoding. The remainder of this Background section briefly reviews various published geophysical simultaneous source techniques, both encoded and non-encoded.
Van Manen [6] suggests using the seismic interferometry method to speed up forward simulation. Seismic interferometry works by placing sources everywhere on the boundary of the region of interest. These sources are modeled individually and the wavefield at all locations for which a Green's function is desired is recorded. The Green's function between any two recorded locations can then be computed by cross-correlating the traces acquired at the two recorded locations and summing over all the boundary sources. If the data to be inverted have a large number of sources and receivers that are within the region of interest (as opposed to having one or the other on the boundary), then this is a very efficient method for computing the desired Green's functions. However, for the seismic data case it is rare that both the source and receiver for the data to be inverted are within the region of interest. Therefore, this improvement has very limited applicability to the seismic inversion problem.
Berkhout [7] and Zhang [8] suggest that inversion in general can be improved by inverting non-encoded simultaneous sources that are summed coherently to produce some desired wave front within some region of the subsurface. For example, point source data could be summed with time shifts that are a linear function of the source location to produce a down-going plane wave at some particular angle with respect to the surface. This technique could be applied to all categories of inversion. A problem with this method is that coherent summation of the source gathers necessarily reduces the amount of information in the data. So for example, summation to produce a plane wave removes all the information in the seismic data related to travel time versus source-receiver offset. This information is critical for updating the slowly varying background velocity model, and therefore Berkhout's method is not well constrained. To overcome this problem many different coherent sums of the data (e.g. many plane waves with different propagation directions) could be inverted, but then efficiency is lost since the cost of inversion is proportional to the number of different sums inverted. Herein, such coherently summed sources are called generalized sources. Therefore, a generalized source can either be a point source or a sum of point sources that produces a wave front of some desired shape.
Van Riel [9] suggests inversion by non-encoded stacking or partial stacking (with respect to source-receiver offset) of the input seismic data, then defining a cost function with respect to this stacked data which will be optimized. Thus, this publication suggests improving cost function based inversion using non-encoded simultaneous sources. As was true of the Berkhout's [6] simultaneous source inversion method, the stacking suggested by this method reduces the amount of information in the data to be inverted and therefore the inversion is less well constrained than it would have been with the original data.
Mora [10] proposes inverting data that is the sum of widely spaced sources. Thus, this publication suggests improving the efficiency of inversion using non-encoded simultaneous source simulation. Summing widely spaced sources has the advantage of preserving much more information than the coherent sum proposed by Berkhout. However, summation of widely spaced sources implies that the aperture (model region inverted) that must be used in the inversion must be increased to accommodate all the widely spaced sources. Since the compute time is proportional to the area of this aperture, Mora's method does not produce as much efficiency gain as could be achieved if the summed sources were near each other.
Ober [11] suggests speeding up seismic migration, a special case of non-iterative inversion, by using simultaneous encoded sources. After testing various coding methods, Ober found that the resulting migrated images had significantly reduced signal-to-noise ratio due to the fact that broad band encoding functions are necessarily only approximately orthogonal. Thus, when summing more than 16 shots, the quality of the inversion was not satisfactory. Since non-iterative inversion is not very costly to begin with, and since high signal-to-noise ratio inversion is desired, this technique is not widely practiced in the geophysical industry.
Ikelle [12] suggests a method for fast forward simulation by simultaneously simulating point sources that are activated (in the simulation) at varying time intervals. A method is also discussed for decoding these time-shifted simultaneous-source simulated data back into the separate simulations that would have been obtained from the individual point sources. These decoded data could then be used as part of any conventional inversion procedure. A problem with Ikelle's method is that the proposed decoding method will produce separated data having noise levels proportional to the difference between data from adjacent sources. This noise will become significant for subsurface models that are not laterally constant, for example from models containing dipping reflectors. Furthermore, this noise will grow in proportion to the number of simultaneous sources. Due to these difficulties, Ikelle's simultaneous source approach may result in unacceptable levels of noise if used in inverting a subsurface that is not laterally constant.
Source encoding proposed by Krebs et al. in PCT Patent Application Publication No. WO 2008/042081, which is incorporated herein by reference in all jurisdictions that allow it, is a very cost effective method to invert full wave field data. (The same approach of simultaneous inversion of an encoded gather will work for receivers, either via source-receiver reciprocity or by encoding the actual receiver locations in common-source gathers of data.) For fixed receivers, the forward and adjoint computations only need to be performed for a single effective source; see PCT Patent Application Publication No. WO 2009/117174, which is incorporated herein by reference in all jurisdictions that allow it. Given the fact that hundreds of shots are recorded for typical 2D acquisition geometries, and thousands in the case of 3D surveys, computational savings from this method are quite significant. In practice, a fixed receiver assumption is not strictly valid for most common field data acquisition geometries. In the case of marine streamer data, both sources and receivers move for every new shot. In addition, due to logistical problems, it is difficult to record data close to the source, and this means that near-offset data are typically missing. This is true for both marine and land surveys. Both of these factors mean that for a simultaneous source gather, every receiver location will be missing data for some source shots. In WO 08/042081, some of the disclosed embodiments may work better than others when the fixed receiver assumption is not satisfied. Therefore, it would be advantageous to have an accommodation or adjustment to straightforward application of simultaneous encoded sources (and/or receivers) inversion that would enhance its performance when the fixed receiver assumption is compromised. The present invention provides ways of doing this. Haber et al. [15] also describe an approach to the problem of moving receivers in simultaneous encoded source inversion using a stochastic optimization method, and apply it to a direct current resistivity problem.
In one embodiment, the invention is a computer-implemented method for full-wavefield inversion of measured geophysical data to determine a physical properties model for a subsurface region, comprising: (a) using a computer to invert a selected shallow time window of arrivals from the measured geophysical data by simultaneous encoded sources and/or receivers inversion to obtain a first physical properties model for the subsurface region; (b) using a computer to invert the measured geophysical data, or a selected deep time window of arrivals from the measured geophysical data, by iterative sequential source inversion, which may use only a sparse sampling of the measured data, to obtain a second physical properties model for the subsurface region, wherein the first physical properties model is used as a starting model and a set of source locations is used to update the second physical properties model in the iterative sequential source inversion; and (c) outputting, displaying, or saving to data storage the second physical properties model of the subsurface region.
In another embodiment, the invention is a computer-implemented method for full-wavefield inversion of measured geophysical data to determine a physical properties model for a subsurface region, comprising: (a) using a computer to invert a selected shallow time window of arrivals from the measured geophysical data by simultaneous encoded sources and/or receivers inversion to obtain a first physical properties model for the subsurface region; (b) using the first physical properties model to simulate, using a computer, synthetic data for longer offsets corresponding to arrivals from deeper then said shallow time window; (c) using a computer to invert the measured geophysical data, wherein the data with longer offsets are augmented, said inversion being simultaneous encoded-sources and/or encoded-receivers inversion, to obtain a second physical properties model of the subsurface region, wherein said augmented data with longer offsets are the sum of the synthetic data for longer offsets and the measured data at the longer offsets; and (d) outputting, displaying, or saving to data storage the second physical properties model of the subsurface region.
The patent or application file contains at least one drawing executed in color. Copies of this patent or patent application publication with color drawings will be provided by the Office upon request and payment of the necessary fee.
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.
One embodiment of the present invention is a hybrid combination of simultaneous encoded source inversion with traditional sequential source inversion. This embodiment uses simultaneous source encoding for the shallower time-window of the data and uses sparse sequential sources to invert the deeper part of the data. Krebs et al. [5,16] show that the encoded simultaneous source cost function can be computed more efficiently than conventional cost functions while still providing accurate inversions. The simultaneous source cost function is defined here as (compare with equation (2) above):
where a summation over receivers and time samples is implied as in Eqn. (2), and:
defines a sum over gathers by sub groups of gathers,
The outer summation in Eqn. (6) is over groups of simultaneous generalized sources corresponding to the gather type (e.g. points sources for common shot gathers). The inner summation, over g, is over the gathers that are grouped for simultaneous computation. For some forward modeling methods, such as finite difference modeling, the computation of the forward model for summed generalized sources (the inner sum over g ∈G) can be performed in the same amount of time as the computation for a single source. Therefore, as shown in Krebs et. al. [5] δ(M,Σcgwg) can be computed very efficiently using Algorithm 3.
wg) using a single run of the simulator
wg as the source,
ψobs(wg)),
Again as shown in Krebs et al. [5] this algorithm can compute Ssim(M) a factor of Ng/NG times faster than S(M) from Eqn. (2).
In fixed receiver geometry the shot gathers have complete receiver coverage. For the shallower window this can be achieved if the near offsets can be populated. The shallower windows are defined as corresponding to times smaller than the arrival times for the faster modes at the longest offsets. The size of the shallow window, over which the encoded simultaneous source approach is applicable, depends on the near-surface velocity and the maximum available offsets.
The negative offsets in the shot gathers may be filled using reciprocity (201). The missing near-offsets may be filled with estimated values, for example from neighboring traces (201). Filling the near-offsets in the shallow window will provide complete receiver coverage to use simultaneous source encoding. Another alternative is to simply eliminate all receivers that were a near offset for some shot. There are various prior approaches to fill in the near-offset data such as reconstruction methods typically used in surface-related multiple attenuation. This makes the shallow windowed shot gathers (202) conform to the fixed-receiver assumption preferable for simultaneous source encoding (203). Therefore, for the shallow window, one can generate simultaneous-source encoded data (204) and invert (205) for a subsurface model (206) that fits these data. Since the computation time is significantly less by several orders of magnitude, a more sophisticated simulation algorithm (e.g. elastic full wavefield inversion) can potentially be used if necessary. Following the simultaneous source inversion with the shallow windowed data, the next step is use the shallow model as a starting model (207) to invert the deeper window data (208) since the fixed receiver assumption breaks down for deeper windows.
The present inventive method includes two embodiments for inverting the deeper window data. A flowchart of these hybrid approaches is presented in
In the second approach (Approach II in
The present invention's hybrid approach is not only applicable to streamer data, but also to a variety of other acquisition geometries where the fixed-receiver assumption breaks down. For example, ocean bottom cable (OBC) acquisitions are typically patch-based (fixed receivers for a subset of shots), and they do not conform to the idealized situation where all receivers are fixed and recording for all shots. Similarly for land acquisitions due to logistical problems the fixed receiver geometry is difficult to achieve.
Many types of encoding functions cg can be used in equation (6) including but not limited to:
Some of these encoding techniques will work better than others depending upon the application, and some can be combined. In particular, good results have been obtained using frequency dependent random phase encoding and also by combining frequency independent encoding of nearby sources with frequency dependent random phase encoding for more widely separated sources. An indication of the relative merits of different encodings can be obtained by running test inversions with each set of encoding functions to determine which converges faster.
It should be noted that the simultaneous encoded-source technique can be used for many types of inversion cost function. In particular it could be used for cost functions based on other norms than L2 discussed above. It could also be used on more sophisticated cost functions than the one presented in Equation 2, including regularized cost functions. Finally, the simultaneous encoded-source method could be used with any type of global or local cost function inversion method including Monte Carlo, simulated annealing, genetic algorithm, evolution algorithm, gradient line search, conjugate gradients and Newton's method.
The present inventive method can also be used in conjunction with various types of generalized source techniques, such as those suggested by Berkhout [7]. In this case, rather than encoding different point source gather signatures, one would encode the signatures for different synthesized plane waves.
Some variations on the embodiment described above include:
The starting model in
The model obtained from shallow window simultaneous source inversion (
In the example presented above, the measured near-offset data were used. In practice, the marine streamer data typically have missing near-offsets. In the next example, offset data up to 200 m are eliminated. The goal here is to show the impact of the missing near offsets in the simultaneous inversion and its importance to be included in the encoded data.
Instead of interpolating the near-offsets traces that is a common procedure for multiple-elimination, the nearest available trace at 200 m is used and multiple copies are made to populate the missing near-offset traces, with this trace shown in
The foregoing application 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 in the appended claims. Persons skilled in the art will readily recognize that in preferred embodiments of the invention, at least some of the steps in the present inventive method are performed on a computer, i.e. the invention is computer implemented. In such cases, the resulting updated physical properties model may either be downloaded, displayed, or saved to computer storage.
This application claims the benefit of U.S. Provisional Application 61/386,828 filed on Sep. 27, 2010, which is incorporated herein by reference. This application is related to that particular application entitled “Simultaneous Source Encoding and Source Separation as a Practical Solution for Full Wavefield Inversion” filed on the same date as the present application and claiming the benefit of U.S. Provisional Application 61/386,831 filed on Sep. 27, 2010. This related application is also incorporated by reference herein.
| Number | Date | Country | |
|---|---|---|---|
| 61386828 | Sep 2010 | US | |
| 61386831 | Sep 2010 | US |
