COMPUTERIMPLEMENTED METHOD AND SYSTEM FOR SMALL CAVE RECOGNITION USING SEISMIC REFLECTION DATA
A computerimplemented method and system implementing the method, are disclosed for computing small cave recognition models, using seismic reflection data. User inputs and earthmodel data are obtained over points of incidence of a survey region, at various angles of incidence. Various models are then computed that serve for cave identification and take part in preliminary seismic exploration and reservoir characterization. Therefore, the attributes developed by the computerimplemented method and system serve as indicators of low velocity and density cave recognition which are capable of separating the cave events from the normal layer events; identifying caves with size larger than half to one wavelength of the dominant signal; and identifying cave diffractions from caves that contain a local maximal/minimal at around nine degrees in amplitude versus angle models.
Latest China Petroleum & Chemical Corporation Patents:
 Tool for jet packing and fracturing and tubular column comprising same
 Computerimplemented method and system for small cave recognition using seismic reflection data
 Modified Ytype molecular sieve, catalytic cracking catalyst comprising the same, their preparation and application thereof
 FOAM PRODUCING METHOD, FIRE EXTINGUISHING METHOD, AND APPLIANCE FOR FOAM EXTINGUISHING
 Independent highspeed sampling for an oil drilling system
The present disclosure relates generally to computerimplemented methods and systems for seismic exploration of subsurface interfaces with impedance contrasts, image and angle gather containing information about the amplitude variation with angle lithology.
BACKGROUND OF INVENTION 1. OverviewSeismic exploration, also called seismic survey, involves the study of subsurface formations of interest and geological structures. In general, the purpose of seismic exploration is to image the subsurface of a survey region to identify potential locations of hydrocarbon underneath the surface. Typically, seismic exploration and reservoir characterizations are performed over a region that is surveyed for its soil, and fluid potential properties. Depending upon the properties found in the survey region, one or various hydrocarbon reservoirs (i.e., oil and gas) may be established. Sometimes, these reservoirs may be found inside subsurface structures or formations, like caves, caverns, or rocks.
In seismic exploration, one or more sources of seismic energies are placed at various locations near the surface of the earth to generate a signal in the form of waves, which travel downward through the earth while entering subsurface formations, like rocks, and caves. Once the waves, generated as a result of the emitted seismic energy, enter the subsurface formation, they get reflected, refracted, or scattered throughout the subsurface, which are then captured by receiving sensors that records, samples or measures said waves. The recorded waves are commonly referred to in the art, as seismic data or seismic traces. These data, or traces may contain information regarding the geological structure and properties of the survey region being explored. They are then analyzed using computerimplemented methods to extract details of the structure and properties of the survey region of the earth being explored. Many of these computerimplemented methods use an inversion model which comprise of an analytical approach in which depthdomain signals corresponding to the reflection of acoustic energy from reflective interfaces between subsurface strata in the earth are converted into one or more traces representative of physical attributes of the strata. In contrast, additional computerimplemented methods are required to provide for more specific lithology identification and fluid discrimination, in the field of seismic exploration and reservoir characterization, mainly using amplitude variation with angle (AVA), intercept and gradient volumes to successfully, efficiently, and accurately identify the diffractions associated with small karst caves (size smaller than half seismic wavelength) that contain low velocity (ν) and density (ρ) from reflection events of overburden and carbonate background layers.
Nevertheless, most subsurface formations of interest (or geological structures) have developed fractures which were originated during karstification, and these karst systems were then buried underground in the form of caves. These caves then turned into the main storage space for hydrocarbon but contained fractured zones which, in a sense made the key content of karst characterization, to become cave identification instead (Tian F, X. B Lu, S. Q Zheng, H. F Zhang, Y. S Rong, D. B Yang, and N. G Liu. 2017, Structure and filling characteristics of paleokarst reservoirs in the northern Tarim basin, revealed by outcrop, core and borehole images; Open Geoscience vol. 9 pp. 266280; and Tian F, Z. X Wang, F. Q Cheng, W Xin, O Fayemi, W Zhang, and X. C Shan. 2019, Threedimensional geophysical characterization of deeply buried paleokarst system in the Tahe oilfield, Tarim basin, China; Water vol. 11(5), p. 1045.). As this is not an easy task, a multiprong approach of core sample description, well logging interpretation, 3D seismic modeling, and highresolution impedance dataset was originally proposed by Fei Tian, et al., supra, (2019) to delineate the 3D geometry of the paleocaves and other paleokarst oil fields. Yet, this has still to assist those skilled in the art with small cave identification, as it further requires using innovative computerimplemented methods, to solve the problem of small cave identification.
2. Analysis of WavesIt is known by persons of ordinary skill in the art, that depthdomain seismic images from the reversetimemigration (RTM) model, reveal the reflectivity from subsurface interfaces with impedance contrasts; and angle gathers contain information about the amplitude variation with angle (AVA). It is also known that fracturedcave systems with different scales and sizes, are common features in carbonate fractured reservoirs. These systems contribute a lot to the production of oil and gas because they provide both storage spaces and migration pathways for hydrocarbons. Using seismic data, like prestack gathers and seismic attributes to identify the characteristics of the fracturedcave system, is a key approach to better understand the carbonate pathway of fractured reservoirs. As such, several methods have been proposed.
For instance, the coherence algorithm (See Marfurt K. J, R. L Kirlin, S. L Farmer and M. S Bahorich 1998, 3D seismic attributes using a semblancebased coherence algorithm Geophysics vol. 63 pp. 11501165), the variance algorithm (See U.S. Pat. No. 6,151,555), and the curvature algorithm (See Roberts A. 2001, Curvature attributes and their application to 3D interpreted horizons, First Break vol. 19 pp. 85100, and Marfurt K. J., 2006, Robust estimates of 3D reflector dip and azimuth, Geophysics vol. 71 pp. 297140) to just name a few, are popular methods in the art, used to characterize the physical properties of fractures using seismic attributes. On the other hand, the amplitude versus azimuth inversion for velocity anisotropy (Rüger A. and Tsvankin I., 1997, Using AVO for fracture detection: analytic basis and practical solutions, The Leading Edge vol. 16 pp. 14291434, and Bachrach R., Sengupta M., Salama A. and Miller P., 2009, Reconstruction of the layer anisotropic elastic parameters and high resolution fracture characterization from Pwave data: a case study using seismic inversion and Bayesian rock physics parameter estimation, Geophysics Prospect, vol. 57 pp. 253262.), and for attenuation azimuth anisotropy (Clark R. A., Benson P. M., Carter A. J., and Guerrero Moreno C. A., 2009, Anisotropic Pwave attenuation measured from a multiazimuth surface seismic reflection survey, Geophysics Prospect vol. 57 pp. 835845; and Shekar B. and Tsvankin I., 2012, Anisotropic attenuation analysis of crosshole data generated during hydraulic fracturing, The Leading Edge vol. 31 pp. 588593), also both using prestack seismic azimuth gathers are implemented to characterize the fractured reservoir parameters.
Nonetheless, hydrocarbon predictions from seismic amplitude and amplitudeversusoffset (“AVO”) still remain a difficult task. An approach used by persons having ordinary skill in the art, is to use seismic reflections to closely relate them to subsurface rock properties. Yet, the strongest AVO in the seismic data is often caused by hydrocarbon saturation in the rocks. Advances on the use of prestack seismic inversion for extracting information in terms of subsurface elastic parameters for seismic data have tremendously helped in characterizing lithofacies and predicting reservoir properties with minimum error thereby reducing the number of dry wells and drilling risks in some basins of the world (See e.g. Russel, B., 2014, Prestack seismic amplitude analysis: an integrated overview: Interpretation, v. 2, no. 2, SC19SC36). Such seismic models have been routinely applied for lithology prediction and fluid detection to identify potential targets for oil and gas exploration. Most recently, it has been widely used for estimating sweet spots in unconventional shale gas applications.
The essence of the seismic models are that the shear modulus of a rock does not change when the fluid saturation is changed. The bulk modulus in turn, changes significantly when the fluid saturation is changed. Since the shear modulus of the rock skeleton is equivalent to the shear modulus of the rock with pore fluid, the shear impedance can then be considered as a seismic attribute that mainly connects to the rock skeleton, whereas the acoustic impedance is dominated by both pore fluid and rock skeleton. As a result, the Poisson impedance can optimally eliminate the effect of rock skeleton in shear modulus from the one in acoustic impedance, which in return can produce better resolution of the fluid content. With seismic data, Poisson impedance is conventionally computed indirectly from Pwave, Swave velocities and density which can be inverted from seismic data directly (See Goodway, supra). However, the indirect way of parameters estimation creates more uncertainties caused by the indirect calculation (See e.g. Zong, Z., X. Yin, and G. Wu, 2013, Elastic impedance parameterization and inversion with Young's modulus and Poisson's ratio, Geophysics, v. 78, no. 6, N35N42). Wang, B., X. Yin, and F. Zhang, 2006, in their research paper Lamé parameters inversion based on elastic impedance and its application: Applied Geophysics, v. 3, pp. 120123; compared the direct way to estimate the Lame parameters from prestack seismic data, to the indirect way of the Lame parameters estimation using the inverted Pwave and Swave, but concluding that it still caused much more bias.
Thus, as a result of all these failed attempts to reduce the accumulative error, to provide accurate and reasonable results; direct and indirect extraction of fluid factors from seismic reflection data have gained much attention in recent years, especially given the wider availability and access to high performance computing systems, required to solve large scientific and engineer problems like the foregoing. Yet, none have proposed a unified, and consistent approach to solve these problems.
3. Fundamental Basis of Prestack MigrationPrestack Migration is a seismic data processing technique to map seismic events onto their appropriate positions (Sheriff, R. E., & Gerdart, L. P., 1995, Exploration Seismology, Cambridge University Press.). Prestack Migration is done either in time domain or depth domain, depending on the complexity of lithology. Prestack Time migration yields an inaccurate image in the presence of strong lateral velocity variation associated with complex overburden structure. In such cases, earth imaging is done by depth migration instead. Strong lateral velocity variation causes significant ray bending at layer boundaries, it gives rise to nonhyperbolic behavior of reflection times on common midpoint (CMP) gathers. As a result, amplitudes and travel times associated with the reflection events with nonhyperbolic moveout are distorted during conventional CMP stacking which is based on the hyperbolic moveout assumption. This causes CMP stack to depart from an ideal zero offset wave field. Therefore, when depth migration is needed, in principle, it is done before stack and not after stack (Yilmaz, Oz, 2001, Seismic Data Analysis, Society of Exploration Geophysicist). In most of these cases, Reverse Time Migration (RTM) is performed, as it is a prestack twoway waveequation migration algorithmsolution suited for imaging in areas of complex structure, being able to handle steep structural dips and high velocity contrasts.
The first step in prestack depth migration is to choose an interval velocity depth model. The quality of the depth image depends heavily on the input data, the inversion algorithm, and a chosen class of models (number of reflection interfaces, parameterization for interfaces, geometry, and velocities within the layers, etc.). Both time and depth migration use a diffraction term for collapsing energy along a diffraction hyperbola to its apex, only the depth migration algorithms implement the additional thinlens term that explicitly account for lateral velocity variation. The general workflow used by one of ordinary skills in the art, for pre stack depth migration is as given below:

 1) Stacking velocity analysis along time horizons;
 2) Rootmeansquare (RMS) velocity analysis along time migrated horizons;
 3) Stacking velocity refinement along time horizons;
 4) RMS velocity refinement along time migrated horizons;
 5) Interval velocity and depth model creation (coherency inversion); and
 6) Interval velocity and depth model refinement and modelling (tomography).
Nevertheless, because the geometry of acquisition in complex geological situations with rough topography is seldom regular, a Kirchhoff prestack depth migration seems to be the method of choice. The particular Kirchhoff integral that is used for data extrapolation and migration (See for instance, Schneider, William A., 1978, Integral Formulation for Migration in two and three dimensions, Geophysics vol. 43 pp. 4976; and Belistein, N. and S. H. Gray, 2001, From the Hagedoorn imaging technique to Kirchhoff migration and inversion, Geophysical Prospecting, vol. 49, pp. 629645) is based on Green's theorem (Morese, Philip M. and Feshbach, Herman, 1953, Methods of theoretical physics, International Series in Puer and Applied Physics, pp. 791869) and an integral solution of the wave equation, with 3D prestack seismic data recorded on separate sets of surface points. In practice the approximate Kirchhoff integral is (Sava, Paul and Biondi, Biondo, 2004, Waveequation migration velocity analysis—I: Theory, Geophysical Prospecting, Stanford Exploration Project, Stanford University), given by:
I(r)≈Σ_{i∈A}_{r}W_{i}(r,m,h)D[t_{D}(r,m,h),m_{i},h_{i}] (1)

 where the image I(r), is defined in a threedimensional space r=(x,y,z), and is equal to the summation of data values D(r,m,h) evaluated at the time t_{D}(r,m,h) and weighted by an appropriate factor W (r,m,h) for amplitude. m_{i }and h_{i }are midpoint and offset positions.
Equation (1) is numerically computed as a loop over each image point I(r) and sums the contributions of all input traces that are within the migration aperture. Kirchhoff prestack depth migration uses the actual raypath from every pointsource to every receivingsensor. This raypath is used to construct the diffraction surface. The migration of a seismic section is achieved by collapsing each diffraction hyperbola to its origin (apex). In practice, each offset plane is first migrated separately and then, all offsets are summed together to generate the stacked migrated image. Therefore, the Kirchhoff prestack depth migration is performed in two steps. First summing data points with equal offsets, and then summing along offsets.
Because of its flexibility, the aforementioned Kirchhoff prestack depth migration (PSDM) is one of the most commonly used tools for processing seismic data. CPUbased PSDM have been used in the art, especially on CPUclusters. Some algorithms, like Reverse Time Migration (RTM), run preferentially on Graphics Processing Units (GPU) based clusters (Xinyi Sun and Sang Suh, 2011, Maximizing throughput for high performance TTIRTM: From CPURTM to GPURTM, SEG Technical Program Expanded Abstracts, pp. 31793183), and as such require a complex structure of computer systems typically comprising of input devices, memory resources, nontransitory program storage computerreadable memory encoded for computing an array of subprocess or subroutines, a system computer, a computer system, and output devices. Investing in GPUbased clusters, means more software needs to be efficiently converted in order to maximize hardware utilization and also decrease seismic processing costs. Kirchhoff PSDM seems to be a prime candidate for such a conversion.
Compared to current CPU, GPU have a lot more cores. Unfortunately, to leverage this parallelism, an existing CPUonlycode needs to be modified. These changes are not straightforward, and they require more knowledge about the architecture of GPU hardware than when implementing a CPUonly software. Some of the main requirements for good performance are memory alignment and access patterns like Compute Unified Device Architecture (CUDA) programming. Nonetheless, to date, most of the articles about porting Kirchhoff migration algorithms to GPU are in time domain (See Panetta, Jairo & Teixeira, Thiago & Filho, P. R. P. & Finho, C. A. & Sotelo, David & Roxo, Fernando & Pinheiro, S. S. & Pedrosa, I. & Rosa, A. L. R. & Monnerat, Luiz & Carneiro, L. T. & Albrecht, C. H. B., 2009, Accelerating Kirchhoff Migration by CPU and GPU Cooperation, pp. 2632, 10.1109/SBACPAD.2009.29.; and Brouwer W., Natoli, V., and Lamont, M., 2011, A novel GPGPU approach to Kirchhoff Time Migration, SEG Technical Program Expanded Abstracts, pp. 34653469), which embodiments of the present invention solve.
4. Tomography InversionTomography of depth migrated gathers is a method for refining the velocitydepth model. Tomography is based on the principle that if migration is carried out, with a correct velocity depth model, the image gathers should have an event depth equal at all receiving sensors (Tianwen Lo and Philip Inder Weisen, 1994, Fundamentals of seismic Tomography, SEG monograph series). Therefore, tomography attempts to correct errors in the velocity depth model by analyzing the residual delays after PSDM. When prestack depth migration is performed with an initial incorrect velocity model derived from inversion methods based on nonglobal approaches, the depth gathers will exhibit nonflatness. The degree of nonflatness is a measurement of the error in the model. Tomography uses this measurement of nonflatness (residual moveout analysis) as input and attempts to find an alternative model, which will minimize the errors. The tomographic principle attributes an error in time to an error in both, velocity and depth.
Yet, in the field of inverse theory, an object is described based on measurements or observations that are associated with that object. For industrialscale applications, there is insufficient data to determine a unique solution and the data obtained may be noisy and/or unreliable. Consequently, an entire branch of mathematics has evolved with attempts to estimate a solution based on the interpretation of inaccurate, insufficient, and inconsistent data. In the case of travel time measurements made in seismic experiment, the inverse problem that is trying to get solved, is the velocity structure of the earth.
When merging the concept of tomography to inversion for multidimensional structures, one can observe that it has been utilized in oil field explorations since the late 1970's. However, the resolution for such inversions is typically in kilometers or even tens of kilometers. Details pertaining to the complex systems of oil field explorations cannot be resolved due to thousands of nodes incorporated into the models, so a 3D imaging that is still based on centralized offline processing and is typically accomplished by multiple activesource recordings where variations over multiple year spans are the main goal. To achieve this, new schemes and methodologies are required to solve the realtime seismic tomography inversion problems, that cannot be achieved using single computational algorithms. This is one motivation for developing computerimplemented method that can relay an accurate imaging systems and models.
Because tomography is based on ray tracing, it can be formulated for reflection, transmission, and refraction. Several techniques for computing statics corrections in seismic reflection surveys make use of refraction tomography, whilst transmission tomography is used for crosswell applications where both the source and the receiver are inside the medium (within the boreholes); hence there is access to any transmitted arrival information. Exploiting amplitude information in addition to arrival times can further assist raybased tomography in estimating a reliable velocity model (See for instance, Semtchenok N M, Popov M M, Verdel A R, 2009, Gaussian beam tomography, Extended Abstracts, 71st EAGE Conference & Exhibition, Amsterdam, U032). In addition to velocity estimation, tomography can be used to estimate other earth parameters, such as absorption. Further, tomography tries to solve a set of simultaneous equations using many raypaths traversing the cells in the model. For a commonmidpoint (CMP) gather, there are many travel time measurements for a given subsurface reflector element so the travel time expression for these five raypaths can be written as:
T=D*S (2)

 where T is the total travel time along the raypath, D is the length of the raypath in the cell of the velocity, and S is the slowness (i.e., the reciprocal of velocity) also in the cell of the velocity.
Unfortunately, in most cases, the matrix D is not invertible. To be invertible, it needs to be squared (i.e., the number of travel time measurements must just happen to be the same as the number of velocity cells in the model) and also to fulfil some other criteria. So, instead, the algorithm is then inverted into the time delay and the slowness update resulting in expression:
Δt=D*Δs (3)
However, there is a bit of a circular argument in the above description of the method because to estimate the raypath segment lengths in the cells by ray tracing, a slowness update and the local dip estimates of the reflector segments in each cell are needed. Hence when a method is initiated with forward modeling by ray tracing an initial guess of the model and the associated travel times computed, a firstguess or estimate model is generated. The tomography must then be iterated to converge on the best estimate of the true model, by minimizing the differences between the observed travel times and those computed by ray tracing for the current guess of the model. Having said that, one skilled in the art would realize that there are other ways of solving for algorithm (2), such as: a direct solver, which is a onestep solution but only suitable for smaller scale problems; or an iterative solver, such as the conjugate gradient method (Scales, J. A., 1987, Tomographic inversion via the conjugate gradient method, Geophysics, vol. 52, pp. 17985), which works well on largescale systems.
The procedure described so far for tomography inversion, was cast in terms of observed and raytraced travel times from unmigrated data. Over the past decade however, one skilled in the art has dealt primarily with data and measurements in the depth migrated domain, so there had also been a need to modify the tomography algorithms to account for the change of measurement domain. As such, tomography inversion in the migrated domain can be performed after either prestack time migration (See as an example Pierre Hardy and JeanPaul Jeannot, 1999, 3D reflection tomography in timemigrated space, SEG Technical Program Expanded Abstracts, pp. 12871290) or prestack depth migration (e.g., Christof Stork, 1992, Reflection tomography in the postmigrated domain, Geophysics vol. 57, pp. 680692.), yielding a depth model in both cases, but the most widespread industrial practice is to invert using measurements from prestack depth migrated data. Which again, when either of these models are used independently, they account for less than useful lithology identification for one skilled in the art.
5. Optic Flow EstimationOptic flow measurement is an early processing step in computer vision systems, that is used in a wide variety of applications, ranging from threedimensional scene analysis to video compression. This processing step was first described by Horn and Schunck (See B. K. P. Horn and B. G. Schunck, Determining optical flow, Technical Report A. I. Memo 572, Massachusetts Institute of Technology, 1980), who devised a simple way to compute the optic flow based on regularization. This first work was then followed by a great number of contributions that proposed alternative methods like: (a) the spatiotemporal filtering methods (See E. H. Adelson and J. R. Bergen, Spatiotemporal energy models for the perception of vision, J. Opt. Soc. Amer., A2:284299, 1985); (b) the split in energy based method (See D. J. Heeger, Optical flow using spatiotemporal filters, International Journal for Computer Vision, 1:279302, 1988); (c) the phase based methods (See D. J. Fleet and A. D. Jepson, Computation of component image velocity from local phase information, International Journal of Computer Vision, 5:77104, 1990); and (d) the region matching methods (See as e.g. P. J. Burt, C. Yen, and X. Xu, Multiresolution flowthrough motion analysis, In Proc. Conference Computer Vision and Pattern Recognition, pages 246252, Washington, 1983).
Many authors noticed that a good way to enhance the reliability of optic flow estimation was to perform a multiscale computation, thereby using properties from one level, and models from different levels. The multiscale approach proved to be very powerful, as in matching methods, it greatly reduced the dimension of the search space. In filtering based methods on the other hand, it increased the range of measurable displacement magnitudes, and relaxed the need for an a priori tuned frequency or scale parameter. This work was motivated by the observation that wavelets are a very welldesigned tool for optic flow measurement. Because of their multiscale structure, and because large scale filtering can be performed efficiently with the fast wavelet transform, they are a very natural tool to measure optic flow: the projection of the optic flow equation onto the wavelets yields a very fast variant of massive filtering.
As an image sequence is a real function l(t; x_{1}, x_{2}) of three variables t; x_{1}, and x_{2 }that it is supposed to be continuous in this first section, the concise notations x for (x_{1}, x_{1}) and x(t) for [x_{1}(t),x_{2}(t)] are used. The standard mathematical model used to find the optic flow is mostly based on a constant brightness assumption: a real point [X_{1}(t), X_{2}(t), X_{3}(t)] in the scene is projected onto the image plane to an image point [x_{1}(t), x_{2 }(t)]
[X_{1}(t),X_{2}(t),X_{3}(t)] →[x_{1}(t),x_{2}(t)] at time t (4)
The optic flow at time t and location x(t) is defined as the velocity of the image point:
The brightness constancy assumption consists in saying that the brightness of [x_{1}(t),x_{2}(t)] at time t is time independent, hence:
I[t;x(t)]=I_{0} (6)
The optic flow ν is therefore constrained by the following equation:
The variant from equation (8) can be further defined to take the possible illumination changes into account using a Lambertian surface aspect model:
I(t;x1,x2)=R(t;x1,x2)*L(t;x1,x2) (9)
In equation (9), R is the reflectance, or a picture sequence fulfilling the brightness constancy assumption. L is the illumination factor that accounts for brightness changes. For a single light source at finite distance, L is a product:

 where i is the angle of incidence of the source light falling on the object, and d is the distance between the source and the object.
Although changes in illumination are caused by relative moves between the light source and the object, one skilled in the art would assume that these changes have slow variations in space. This consists in assuming that the spatial derivatives
are negligible. By differentiation, the equation ends up being:
and a modified optic flow equation taking illumination changes into account:
The parameters that can then be estimated from this constraint are the optic flow ν=(ν1, ν2) and the time derivative of log
Note mat L cannot be measured more accurately than up to a multiplicative constant, since a pattern α×P and an illumination α^{−1}×L gives the same sequence I and thus the same constraints. This is reflected in the fact that that the time derivative of log L is not affected by a multiplication of L with a constant factor.
The problem of optic flow measurement has infinitely many solutions and is by essence illposed. The only way to reduce the number of solutions is to compute additional assumptions and apply it directly to a particular model. This rule applies in all optic flow measurement techniques and no method makes an exception to it.
6. Semblance AnalysisBuilding subsurface velocity model is one of the most important issue in exploration geophysics. There are generally four ways for building the velocity model. One of those four ways is by employing a normalmoveout (NMO) based velocity analysis, which requires selecting peaks in the velocity spectra (See Taner, M. T., and F. Koehler, 1969, Velocity spectra—Digital computer derivation and applications of velocity functions: Geophysics, 34, 859881; Fomel, S., 2009, Velocity analysis using AB semblance: Geophysical Prospecting, 57, no. 3, 311321; Luo, S., and D. Hale, 2012, Velocity analysis using weighted semblance: Geophysics, 77, no. 2, U15U22), which in turn is obtained by applying several NMOs with different velocities and then calculating their corresponding semblances. As it can observed, using semblance measures is an indispensable step for building an initial subsurface velocity model. Nevertheless, there exist several types of semblance methods, depending upon the specific datasets that is used. A person having ordinary skills in the art will recognize that the most commonly used semblance method, uses a conventional approach.
The conventional semblance approach was originally defined by Neidell and Taner in 1971 (See Neidell, N. S., and M. T. Taner, 1971, Semblance and other coherency measures for multichannel data: Geophysics, vol. 36, pp. 482497,) as:

 where i and j are time sample indices, s(i) denotes the semblance for time index i, 2M+1 is the length of the smoothing window (in this formulation a boxcar filter is used), along time axis, and a(j,k) is the trace amplitude at time index j and trace number k of the NMOcorrected common midpoint (CMP) gather.
Others have developed variations off the conventional semblance approach, as observed in Sarkar, D., J. P. Castagna, and W. Lamb, 2001, AVO and velocity analysis: Geophysics, vol. 66, pp. 12841293; and Fomel, S.; and E. Landa, 2014, Structural uncertainty of timemigrated seismic images: Journal of Applied Geophysics, vol. 101, pp. 2730), that modified the conventional semblance formulation and proposed amplitudeversusoffset (AVO) adaptive semblances. Similarly, Luo and Hale, supra, proposed increasing the resolution of the semblance map in order to distinguish the peaks between primary and multiple reflections. They did this, among other things to increase the resolution of the semblance spectra for picking the true NMO velocity, due to the existence of multiples common midpoint (CMP) gathers (Luo and Hale, supra).
7. Trimming StaticsDatasmoothing statics methods assume that patterns of irregularity which events have in common, result from nearsurface variations and hence staticcorrection trace shifts should minimize such irregularities. Most automatic staticsdetermination programs employ statistical methods to achieve the minimization. Datasmoothing methods are generally applied to remove small residual errors after first applying other methods. These secondorder statics correction methods are often called trim statics (See for instance Cox, Mike, 1999, Static Corrections for Seismic Reflection Surveys: Society of Exploration Geophysicists).
Several geophysical software packages for performing trim statics are known in the art. Examples of such products are available from companies like HampsonRussell. According to these known computerimplemented techniques, the trim static process they provide, can fix migration moveout problems on prestack data by adjusting the time an event occurs through userdefined controlled parameters. The software thereby automatically flattens the far offsets, yet still providing geophysicists full control and understanding of the conditioning processes applied, therefore giving persons skilled in the art, greater confidence in the results of the subsequent reservoir characterization.
Regardless of the software package use, when trimming statics, the summing of traces in a common midpoint (CMP) gather assumes that the reflection data is aligned, and the amplitude will increase with the fold or number of traces in the gather. As such, noise will also sum, but its root mean square (RMS) amplitude will increase with the squareroot of the fold. The corresponding gain in the signaltonoise ratio (SNR) will be proportional to the square root of the fold. It may be assumed then, that a high fold will give a better SNR. However, when the reflection data is misaligned, the amplitudes will not increase proportional to the fold, and a lower SNR may be obtained. There the objective of estimating statics is to align the reflection data to improve the SNR of the stacked section. However, the application of trim statics directly to each seismic trace can be very dangerous as any seismic response can be generated. Nonetheless, corrections are applied to seismic data to compensate for the effects of variations in elevation, nearsurface low velocitylayer (weathering) thickness, weathering velocity and/or reference to a datum.
Statics are estimated from the time shift found by crosscorrelating input data with a model trace. This time shift may be applied to the source trace when the time shifts are very small, say less than 6 ms, and is referred to as trim statics. However, the initial estimate of the crosscorrelation static is typically used as input to a secondary surfaceconsistent process. This secondary process decomposes all the estimated statics into source, receiver, offset, and structure statics, with only the corresponding source and receiver statics applied to an input trace. This secondary process is designed to reduce the timing errors that will occur if the correlation static is applied directly to the input trace.
When trim statics are applied to noisy data, the noise may be aligned with the model trace. The sum of a number of traces in a common midpoint (CMP) gather, in which trim statics have been applied, will therefore produce a stacked trace that tends to be similar to the model trace. Consequently, the resulting stacked section may also appear similar to the model trace giving the illusion of reflection data.
8. Flattening AlgorithmTo flatten an image gather, each image gather is shifted vertically to match a chosen reference point. A reference point can be the intersection of the horizon and a well pick. In three dimensions, this reference point becomes a vertical line or seismic trace reference. To flatten 3D cubes, it has been the study of individuals skilled in the art the need to find a mapping field T(x,y,t) such that each time slice of this field contains the locations of all data points along the horizon that happens to intersect the reference trace at that time slice. This is achieved by summing time shifts or depth dips to get total time shifts which are then used to flatten the data.
Clearly then, the first step in achieving a flattened image gather is to calculate local dips everywhere in the 3D seismic cube. Dips can be calculated efficiently using a local planewave destructor filter as described by Claerbout (see Claerbout, J. F., 1992, Earth soundings analysis: Processing versus inversion; Blackwell Scientific Publications) or with an improved dip estimator described by Fomel (see Fomel, S., 2002, Applications of planewave destruction filters; Geophysics, vol. 67, pp. 19461960). For each point in the data cube, two components of dip, b and q, are estimated in the x and ydirections, respectively. These can be represented everywhere on the mesh as b(x,y,t) and q(x,y,t). If the data to be flattened are in depth, then dip is dimensionless, but when the data are in time, dip has units of time over distance. Therefore, the goal then turns towards finding a timeshift (or depthshift) field τ(x,y,t) such that its gradient approximates the dip p(x, y, τ). The dip is a function of T because for any given horizon, the appropriate dips to be summed are the dips along the horizon itself. Using the matrix representation of the gradient operator
and the estimated dip p=(b q)^{T}, the regression turns to ∇T=(x,y,t)=p(x,y,τ) but because it is a function of the unknown τ(x,y,t) and GaussNewton iterating approach is utilized to achieve the following equation for flattening image gathers:
r=w[∇τ_{k}(x,y,t)−p(x,y,T_{k}) (14)
Planewave decomposition of a wavefield, such as a commonshot gather, can be achieved by applying linear moveout and summing amplitudes over the offset axis. This procedure is called slant stacking. An underlying assumption of slant stacking is that of a horizontally layered earth model. Conventional processing is done primarily in midpointoffset coordinates. Slant stacking replaces the offset axis with the ray parameter paxis. The ray parameter is the inverse of the horizontal phase velocity. A group of traces with a range of p values is called a slantstack gather.
The slant stacking algorithm can be traced back to the Radon transformation concept (See Radon, Johann (1917), “Ober die Bestimmung von Funktionen durch ihre Integralwerte langs gewisser Mannigfaltigkeiten”, Berichte Ober die Verhandlungen der KoniglichSächsischen Akademie der Wissenschaften zu Leipzig, MathematischPhysische Klasse [/Reports on the proceedings of the Royal Saxonian Academy of S7ciences at Leipzig, mathematical and physical section], Leipzig: Teubner (69): 262277; Translation: Radon, J.; Parks, P. C. (translator) (1986), “On the determination of functions from their integral values along certain manifolds”, IEEE Transactions on Medical Imaging, 5 (4): 170176, doi:10.1109/TMI.1986.4307775, PMID 18244009). In general, a slant stacking algorithm comprises transforming the offset axis, as part of a migration method. It extracts a local stepout or ray parameter in the form of
(where x is the offset, and t is the twoway traveltime) and enables immediate downward continuation, even when mixed apparent velocities are present as with diffractions, and multiple reflections.
When looking at gathers for events of some particular stepout,
it involves scanning hyperbolic events to find the places that are tangent to a straight line of slope p. The search and analysis becomes easier when the data is replotted with linear moveout analysis, which involves a source of energy located at offset x=g−s at time t in the (x,t)plane “movedout” to time τ=t−px in the (x, τ)plane (wherein τ is the intercept time at p=0). The linear moveout analysis converts the task of identifying tangencies to constructed parallel lines to the task of locating the tops of convex events. As such, it converts all events stepping out at a rate p in (x,t)space, to horizontal’ events in (x,τ)space. The presence of horizontal timing lines facilitates the search for and the identification and measurement of the locations of the events.
After the linear moveout analysis has been performed, the components in the data that have Snell parameters near p are slowly variable along the x′axis. To extract them, a lowpass filter is applied on the x′axis, and done so for each value of T. The limiting case of lowfrequency filtering extracts the mean, which leads to the concept of slant stacking. The slant stacking algorithm, then processes the linear moveout, with the x′axis, converting the entire gather of P(x,t) to a single trace that is a function of T, represented in the form of:
S(p,τ)=Σ_{x}P(x,τ+px) (15)

 where, S(p,τ) represents a plane wave with ray parameter
By repeating the linear moveout correction for a range of p values and performing the summation in equation (15), a complete slantstack gather is constructed. A slantstack gather, in practice, alternatively is referred to as a τ−p gather; it consists of all the dip components within the specified range of p values in the original offset data.
To restore amplitudes properly, rho filtering is applied before inverse mapping. This is accomplished by multiplying the amplitude spectrum of each slantstack trace by the absolute value of the frequency. Rho filtering is equivalent to differentiating the wavefield before the summation that is involved in the integral formulation of migration (migration principles).
A conventional slant stack yields an exact planewave decomposition when dealing with line sources; a slant stack algorithm in contrast, yields an exact planewave decomposition when dealing with point sources. A proper slant stack is generated using the same steps as described for the slant stack algorithm, except that a convolution of the linear moveoutapplied wavefield by a filter operator is performed before summation. This operator corrects for 3D effects by converting a wavefield that was obtained from a point source into a wavefield that was obtained from a line source. As far as kinematics is concerned, the two types of slant stacking are equivalent, yet only compatible to specific circumstances. This is because the slantstacking can satisfy a onedimensional wave equation in a laterally homogeneous model, which will depend on reflection coefficients for one ray parameter, a twodimensional problem, or a threedimensional problem with a line source.
Hence the reason why the analysis of seismic reflection and refraction data using the process of slant stacking has recently received much attention, and requires a faster way of analyzing, processing and computing that can only be achieved with the use of computer systems.
10. Methodology for Cave RecognitionExplored caves are only a limited portion of those actually existing underground (See White, W. B., 1990; Surface and nearsurface karst landforms, in Higgins, C. G., and Coates, D. R., eds., Groundwater Geomorphology: The Role of Subsurface Water in EarthSurface Processes and Landforms; Geological Society of America Special Paper vol. 252, pp. 157175; and Ford, D. C., and Williams, P., 2007; Karst Hydrogeology and Geomorphology: Chichester, John Wiley & Sons Inc.; p. 576). To obtain information about those hidden caves, or unknown or inaccessible continuations of known caves, a person skilled in the art must use a combination of methods (Parise, M., and Lollino, P., 2011; A preliminary analysis of failure mechanisms in karst and manmade underground caves in Southern Italy; Geomorphology, v. 134, no. 12, p. 132143; Margiotta, S., Negri, S., Parise, M., and Valloni, R., 2012; Mapping the susceptibility to sinkholes in coastal areas, based on stratigraphy, geomorphology and geophysics; Natural Hazards, vol. 62, no. 2; and Pepe, P., Pentimone, N., Garziano, G., Martimucci, V., and Parise, M., 2013; Lessons learned from occurrence of sinkholes related to manmade cavities in a town of southern Italy, in Land; Proceedings of the 13th Multidisciplinary Conference on Sinkholes and the Engineering and Environmental Impacts of Karst, Carlsbad (New Mexico, USA), National Cave and Karst Research Institute; pp. 393401). As such, near surface geophysical methods have recently become an important tool in karstcaves research.
The idea behind most of said geophysical methods focuses on a material property of the void that is significantly different from the surrounding host rock and thus makes a material contrast. This material contrast can then be detected using a specific geophysical techniques (See for example Gibson, P. J., Lyle, P., and George, D. M., 2004; Application of resistivity and magnetometry geophysical techniques for nearsurface investigations in karstic terranes in Ireland; Journal of Cave and Karst Studies, v. 66, no. 2, p. 3538; Mochales, T., Casas, A. M., Pueyo, E. L., Pueyo, O., Roman, M. T., Pocovi, A., Soriano, M. A., and Anson, D.; 2008; Detection of underground cavities by combining gravity, magnetic and ground penetrating radar surveys: A case study from the Zaragoza area, NE Spain; Environmental Geology, v. 53, p. 10671077; or Margiotta, S., and Varola, A., 2007, Il paleosito di Cutrofiano (Salento), proposta per l'istituzione di un parco—museo; Atti della Societa' Toscana Di Scienze naturali—Memorie. Serie A 112,18). Among some of the most frequently used geophysical techniques, a person skilled in the art would recognize as: electricalresistivity tomography and microgravity can be mentioned, albeit additional methods can provide very useful information as well which none have been tried.
One of those methods seems to be seismic amplitudeversus angle (AVA) or amplitude versus offset (AVO) analysis. These techniques which geophysicists use to determine mainly the formation properties, such as the thickness, porosity, density, velocity, lithology and fluid content of rocks. As such, it has been a powerful geophysical method in aiding the direct detection the presence of hydrocarbons.
Persons of ordinary skill in the art, will recognize that the angle in AVA analysis refers to the incidence angle θ, or half of the opening angle between the source and receiver wavefields. For an isotropic elastic medium, the Pwave reflectivity varies approximately linearly with respect to sin^{2}θ, such relationship is called the 2term Shuey's approximation—an approximation to the reflectivity versus incidence angle by Aki and Richards, 1980, and valid for small or medium angles below 30°:
The first term R_{o }is a linearized version of the zero angle/offset reflection coefficient and is a function of only P wave velocity V_{p }and density ρ. The second term is a gradient multiplied by sin^{2 }θ, it is dependent on changes in P wave velocity V_{p}, Swave velocity V_{s}, and density ρ, and has the biggest effect on amplitude changes as a function of angle/offset. For a function with the variable of sin^{2 }θ, equation (16) represents a linear equation with R_{o }as intercept and G as gradient.
For an acoustic medium, the equation for G is
a function of Pwave velocity, V_{p}. Nevertheless, equations 16 and 17 are only valid for the reflection energy, when setting V_{s}=0, or when
Other approaches that have been tried for cave detection using multiple geophysical algorithms have been successful only for caves within specific regions containing unique characteristics. Such has been the case of Kaufmann, Romanov and Nielbock in 2010 (see Kaufmann, G., Romanov, D., and Nielbock, R. (2011); “Cave detection using multiple geophysical methods: Unicorn cave, Harz Mountains, Germany”; GEOPHYSICS vol. 76; pp. B71B77) where the focus of the research had been in the Unicorn cave in the Southern Harz Mountains of Germany. In particular, their research involved using the geophysics methods of gravimetry and electrical resistivity imaging over the cave area, merely to identify the subsurface voids and the extent of the sediment infill of the cave passages. Their choice of methods was based on several conditions unique to the cave which involved: (1) having the well surveyed, (2) a shallow overburden, (3) a large airfilled passage, and (4) a thick sediment cover that, concealed true passage size. Using the cave survey as an initial model for the subsurface structure, they successfully identified the airfilled cave with both methods, and then use it to infer the thickness of the sediment infill by forward modeling and thereby identifying a possible southward continuation beyond the currently explored passages.
Most recently another example of successful cave detection using complex geophysical methods were used by Ptiska, et. al., in 2014 (See Putiska, R., Kusnirak, D., Dostal, I., Lacny, A., Mojzes, A., Hok, R., Pasteka, M., Krajnak, M., and Bosansky, M.; 2014; Integrated geophysical and geological investigations of karst structures in Komberek, Slovakia; Journal of Cave and Karst Studies, v. 76, no. 3, p. 155163). Their investigation was used over a small karst area aimed, solely to confirm geological localization of known sinkholes, and to find possible continuations of caves and voids below the surface. Using radiometric mapping and dipole electromagnetic profiling, supported by geological mapping, they were able to refine the geological boundaries of the lithological units, only within a certain karst area. Then using resistivity tomography and microgravity methods, a final geological crosssection model of the area was constructed. Yet the methods and evaluation techniques they used, were not successfully applied to other karst areas to identify hidden voids that possibly constitute karst hazards.
Therefore, each of the individually discussed methods, as well as the combination of the aforementioned methods as demonstrated, suffer from a wide array of limitations which makes it difficult to conduct reliable and accurate readings, and thereby limiting a person skilled in the art to properly assess low velocity, and density of subterranean caves. As such, in view of the known art, it is therefore seen as one object of the invention to improve and enhance known methods and systems for characterizing and recognizing subterranean caves using advanced computerimplemented systems, that can quickly, and accurate compute an array of functions, whilst providing the user of said computerimplemented systems, typically a person skilled in the art, with full control and understanding of the conditioning processes applied, thereby instilling greater confidence and less uncertainty in the results of the recognized cave, for a subsequent reservoir characterization over a survey region having at least one hydrocarbon well or reservoir location.
SUMMARY OF THE INVENTIONTypically, exploration and reservoir characterizations are performed over a region that is surveyed for its soil, and fluid potential properties. Depending upon the properties found in the survey region, one or various hydrocarbon reservoirs (i.e., oil and gas) may be established. Once these reservoirs, or wells, are established within the survey region computerimplemented methods are used for automating the interpretation of seismic results. Many of these computerimplemented methods are used independently by one skilled in the art because of their limitations to process a vast amount of data, perform distinct computational algorithms at the same time, or comprise of a single computerimplemented analytical approach in which depthdomain signals corresponding to the reflection of acoustic energy from reflective interfaces between subsurface strata in the earth are converted into one or more traces representative of physical attributes of the strata. Furthermore, these attributes are then analyzed offline by one skilled in the art, like geophysicists which in turn create a summary document in another existing computerimplemented program like word processors or spreadsheets.
In contrast, it is therefore an object of the embodiments of the present invention, to provide a single computerimplemented method for lithology identification and fluid discrimination, in the field of seismic exploration and reservoir characterization using an array of seismic reflection data, to successfully, efficiently, and accurately; identify the diffractions associated with small karst caves (size smaller than half seismic wavelength) with low velocity and density from reflection events of overburden and carbonate background layers.
Based on innetwork tomography inversion, embodiments of methods and systems to improve the realtime performance of the distributed approach are described herein. The key challenge of realtime tomography inversion is to update the system model incrementally before all information of seismic events is available. In one aspect, distributed incremental leastsquares algorithms can be used, in which history data are exponentially weighted according to oldness. In another aspect, use of a row action matrix method is described, which does not require the full design matrix to be in memory at one time. In fact, the approach of consecutive back projection can incorporate new information (raypaths or rows), in real time. In this case, the model is constantly being updated by new information as fractures occur and are incorporated to the computerimplemented method via a high performance computer system apparatus that minimizes computational costs, a person skilled in the art subjectivity, and the need to further summarize any analysis. Therefore, this approach is ideal for real time inversion of an evolving structure where fractures can be dynamically included in the analysis.
Further details, examples and aspects of the invention will be described below referring to the drawings listed in the following.
The teachings of the present invention can be readily understood by considering the following description in conjunction with the accompanying drawings.
Reference will now be made in detail, to several embodiments of the present disclosures, examples of which, are illustrated in the accompanying figures. It is noted that wherever practicable similar or like reference symbols may be used in the figures and may indicate similar or like functionality. The figures depict embodiments of the present disclosure, for purposes of illustration only. One skilled in the art will readily recognize from the following description that alternative embodiments of the structures, systems, and methods illustrated therein may be employed without departing from the principles of the disclosure described herein.
In
As shown on
Because all points of incidences 104, and all seismic data recording sensors 105 are placed at different offsets 210, the survey seismic data or traces (213 incidence trace, and 214 receiver trace), also known in the art as gathers, will be recorded at various angles of incidence represented by 208. The points of incidence 104 generate downward transmission rays 213, in the earth that are captured by their upward transmission reflection 214, through the recording sensors 105. Well location 103, in this example, is illustrated with an existing drilled well attached to a wellbore, 209, along which multiple measurements are obtained using techniques known in the art. This wellbore 209, is used to obtain well log data, that includes Pwave velocity, Swave velocity, Density, among others. Other sensors, not depicted in
As shown in
System computer 405, retrieves in a parallel operation the seismic reflection data, 304 to compute a source wavelet, 306, which originated as a packet of energy from the source point, 104, having a specific origin in time, and returned to the receivers 105, as a series of events distributed in time and energy. An initial set of wavelets is extracted from the retrieved image gathers found in seismic reflection data 304, which acts as an initial estimate of what a final source wavelet should look like. The extraction is computed within the present invention, using statistical wavelet extraction procedures which uses the seismic traces alone, to extract the source wavelets 307. However, in order to extract each source wavelet, the corresponding trace in the set of image gathers is used, an analysis window is extracted, and the start and end of the extracted analysis window are tapered in lengths equal to the lesser of 10 samples, or ¼ of the analysis window. An autocorrelation process is computed in a data window, wherein the length of the autocorrelation is equal to ½ of the desire wavelength. The amplitude spectrum of the autocorrelation is computed, and the square root of the autocorrelation is computed, by system computer 405 in order to compute the source wavelet 306. This operation is performed to approximate the amplitude spectrum of the source wavelet to be generated at 307. Thereafter, the desired phase (e.g. 0°, 5°, 10°, 15°, etc.) is added, and the inverse Fast Fourier Transform (FFT) is computed by system computer 405, to generate or extract at 307 the initial source wavelet and store it at memory resource 404. The extraction procedure 307 verifies if there were other wavelets produced before, generates and then stores the initial source wavelet. If other wavelets were created, the new wavelet is generated or parsed at 307, next to the previously generated wavelet from other traces of the same image gathers. The extraction procedure continues until all traces from the seismic reflection image gather data 304 have been analyzed. Once all wavelets have been statistically extracted by system computer 405, an extracted source wavelet is generated 307, and stored in memory resource 404 for later retrieval 308 by either the system computer 405 or the nontransitory program computer readable memory device, 406.
The system computer 405, then sends a message hook to the nontransitory program computer readable memory device, 406, to initiate subroutine 310. Subroutine 310, is then initiated by the nontransitory program computer readable memory device, 406 by retrieving from the memory resource 404, the final migration velocity model 309, the seismic reflection data 304, and the source wavelet 308. This subroutine computes a reverse time migration (RTM) model 310 in depthdomain using system computer 405 and nontransitory program computer readable memory device, 406, hybrid GPU/CPU computer hardware, where both computational components take part in the computation of wave propagation and wavefield crosscorrelation. This approach provides more efficient computational cost while reducing the artifacts produced by backscattering at sharp velocity contrasts. Once a final reverse time migration is generated, 626, and stored in the memory resource 404, it is then retrieved at 311, by system computer 405, and converted to timedomain 312. At 313, the final reverse time migration model in timedomain is then generated, and stored by the system computer 405, in memory resource 404. The system computer 405, then retrieves the final reverse time migration model in timedomain at 314, from the memory resource, 404, and initiates computing a rootmeansquare (RMS) velocity model 315, assuming a sequence of parallel horizontal layers (like 102, 203, 204, 211, and 212) of interval velocity V_{INT i }in the ith layer for all the receiver traces, 214, at various travel time values of t. This is done in order to derive the interval velocity profiles of the seismic traces, 213 and 214. The rootmeansquare velocity is then generated at 316, by system computer 405, and stored at memory resource, 404.
System computer then retrieves the rootmeansquare (RMS) velocity model at 317, and sends a message hook to the nontransitory program computer readable memory device, 406, to initiate subroutine 319. Subroutine 319, computes a residual normal moveout (RNMO) corrected image gather from the corrected image gather, using the corrected image gather from the final reverse time migration image gather in timedomain, 319, to generate a residual normal moveout corrected image gather 710, and store it in the memory resource, 404. This subroutine uses the amplitude versus offset (AVO) technique in order to automatically perform detailed residual moveout analysis. Automatic residual moveout operates as conventional normal moveout (2nd order) and as 4th order correction to correct gathers at large offsets which exhibit 4th order moveout. After the nontransitory program computer readable memory device, 406, generates and stores a residual normal moveout corrected image gather, 710; the system computer 405 retrieves the residual normal moveout image gather at 320, and processes another message hook to the nontransitory program computer readable memory device, 406, indicating the need to initiate subroutine 321.
The nontransitory program computer readable memory device, 406, initiates subroutine 321, using the retrieved residual normal moveout corrected image gather at 320, to compute a final flattened image gather. Computing a flattened image gather 321, from the corrected residual normal moveout image gather 320, generates an isoproportional representation of the data based on specified horizons, times, and userdefined tolerance values. As such, subroutine 321, applies an automatic and continuous velocity picking to the retrieved residual normal moveout corrected image gather 320, as well as corrects any residual moveouts from 320, that remain for final processing. A final flattened image gather is generated at 810, and stored in the memory resource, 404.
At 323, after receiving a message from the nontransitory program computer readable memory device, 406, the system computer 405, retrieves the final flattened image gather from the memory resource 404, and sends a message hook to the nontransitory program computer readable memory device, 406, to initiate subroutine 323 of computing a final intercept volume and a final gradient volume from the retrieved final flattened image gather, 321.
Subroutine 323 computes an intercept and gradient volumes using equation (15), to generate an initial intercept and gradient volume 905, that is then stored by the nontransitory program computer readable memory device, 406, to the memory resource 404, and the signals the system computer, 405, to retrieve the flattened image gather 324, and compute the multiplication 325 of the generated final intercept and gradient volumes 905. Said multiplication 325 applies the condition of R_{o}*G<0 for velocity and density changes in the opposite direction, to generate a small cave identification model, 326 (illustrated by
on the mean time.
As it pertains to
The computer system device, 407, acts as a user interface to system computer 405, and the nontransitory program computer readable memory storage device, 406; to input, set, select, and perform the operations of retrieving, computing, generating, invoking, determining, converting, and correcting functions (the message hook procedures). Said computer system device, 407, is connected to (wired and/or wirelessly) to the system computer 405, and nontransitory program computer readable memory storage device 406. The computer system device, 407, further includes other devices like a central processing unit (CPU), 408, a display or monitor, 409, a keyboard, 410, a mouse, 411, and a printer, 412.
The system computer device, 405, has firmware and software providing for the connection and interoperability of the multiple connected devices, like the memory resources for storing data, 404, the nontransitory program computer readable memory device storage, 406, and the computer system device, 407. The system computer, 405, includes an operating system, a set of message hook procedures, and a system application.
The operating system embedded within the system computer 405, may be a Microsoft “WINDOWS” operating system, OS/2 from IBM Corporation, UNIX, LINUX, Sun Microsystems, or Apple operating systems, as well as myriad embedded application operating systems, such as are available from Wind River, Inc.
The message hook procedures of system computer 405 may, for example, represent an operation or command of the memory resources, 404, the computer system device, 407, the nontransitory program computer readable memory storage device, 406, which may be currently executing a certain step process or subroutine from the computerimplemented method for small cave recognition using seismic reflection data.
The set of message hook procedures may be first initiated by an input from: the user, like the entering of userdefined values or parameters; the manipulation of the computer system device, 407; the processing of operations in the nontransitory program computer readable memory device storage, 406; or automatically once certain data has been stored or retrieved by either the memory resources, 404, or the nontransitory program computer readable memory device storage, 406. Based on any of these inputs, processes or manipulation events, the memory resources, 404, the nontransitory program computer readable memory storage device, 406, or the computer system device, 407; generate a data packet that is passed to the system computer, 405, which are indicative of the event that has occurred as well as the event that needs to occur. When system computer, 405, receives the data packet, it converts it into a message based on the event, and executes the required step of the computerimplement method. The computerimplement method includes a set of message hook lists that identifies the series of message hook procedures. When the operating system receives the message, it examines the message hook list to determine if any message hook procedures have registered themselves with the operating system. If at least one message hook procedure has registered itself with the operating system, the operating system passes the message to the registered message hook procedure that appears first on the list. The called message hook executes and returns a value to the system computer, 405, that instructs the system computer, 405, to pass the message to the next registered message hook, and either 404, 406 or 407. The system computer, 405, continues executing the operations until all registered message hooks have passed, which indicates the completion of the computerimplemented method by the generation of a small cave identification 326.
After the last message hook procedure has passed or been issued, the computer system device displays through a userinterface, on display 410, a message that indicates to the end user that the system computer, 405, has completed executing the computerimplemented method, and is ready to display on 410, and print on 412, a surface seismic reflection data represented in timedomain, a well log data represented in timedomain, an initial migration velocity model, a final migration velocity model, an image gather, a source wavelet, a final reverse time migration image gather in depthdomain, a final reverse time migration image gather in timedomain, a corrected image gather, a residual normal moveout corrected image gather, a final flattened image gather, a final amplitude versus angle volume, a set of initial moveout curves, a first, second, third and final Kirchhoff prestack depth migration gathers, an initial reflector dipangle image, a final reflector dipangle image, an updated initial migration velocity model, a user input flat value, a final intercept volume, a final gradient volume, and a small cave identification model.
Subroutine 305 of computing a final migration velocity model, as illustrated in
The tomography inversion used by the present invention involves computing statics corrections in the seismic reflection survey region making use of refraction tomography, and transmission tomography. The computed tomography uses equation (3) but with multiple iterations, in order to converge on the best estimate of the true model, thereby minimizing the differences between the observed travel times and those computed by ray tracing for the current guess of the model. Once computing the tomographic inversion model, 519 is complete, the nontransitory computer readable memory device, 406, generated a tomographic inversion model, at 520, and stores on memory resource 404. The system computer, 405, receives the signal to retrieve the generated tomographic inversion model at 521, and signals, the nontransitory computer readable memory device, 406, to begin 522 by updating the generated initial migration velocity model 506. Once updated, the nontransitory computer readable memory device, 406, starts computing, at 523 the final Kirchhoff prestack depth migration image gather from the updated initial migration velocity model 506. The nontransitory computer readable memory device, 406, generates, then stores on memory resource 404, a final prestack depth migration image gather, 524. The system computer, 405, then receives the signal to retrieve the generated final prestack depth migration image gather, at 525, and provides it to the nontransitory computer readable memory device, 406, to generate, at 526, a flatness value of the retrieved final Kirchhoff prestack depth migration image gather. The system computer, 405, then sends the inputted userdefined flatness value 502, to the nontransitory computer readable memory device, 406, to begin computing, at 527, the inequality between the generated flatness value 526 against the inputted userdefined flatness value 502. If the nontransitory computer readable memory device, 406, determines that the generated flatness value 526 is equal to, or less than, the inputted userdefined flatness value 502, then it generates the final migration velocity model of 528. If the generated flatness value 526 is greater than, the inputted userdefined flatness value 502, the nontransitory computer readable memory device, 406, sends a message to system computer, 405, that it will reinitiate from step 507 of computing an initial Kirchhoff prestack depth migration image gather, until it determines that the generated flatness value 526 is equal to, or less than, the inputted userdefined flatness value 502. Nevertheless, a person having ordinary skills in the art would soon realize that, to save on computational cost the user of computer system, 407, may at any time stop iteration 527, determining through the use of the display 409, that the stored data, and the inequality are complete, thereby indicating through the use of keyboard 410, and mouse 411 to the CPU 408, that the nontransitory computer readable memory device, 406 should stop computing iteration 527.
Subroutine 310 of computing a final reverse time migration image gather, using the retrieved surface seismic reflection data, the retrieved final migration velocity model, and the retrieved source wavelet, employing a finite different method is illustrated by
Step 605 computes a forward wave propagation in the subsurface by inputting the wavefield recorded at the surface, 102, and by stepping backwards in time, thereby propagating the seismic events to the subsurface location where they were generated, 211 and 212. The nontransitory computer readable memory device 406, does so by reversing the forward modeling operation, inputting the retrieve final migration velocity model 309, selecting the source location, and computing the wave equation f (x±νt) (where x is the sourcewavefield, ν the velocity of the wave, and t is time); in incremental time steps from the timestep value of zero to the timestep maximum value, 603. Nevertheless, the computational method used by the nontransitory computer readable memory device 406, to solve the multiple differential equations, uses the finite difference method by approximating the differential equations with difference equations, in which finite differences, approximate the derivatives. Thereafter, the nontransitory computer readable memory device 406, generates, then stores in memory resource 404, the forward sourcewavefield 606 for the system computer 405, to retrieve at 609 and message back to the nontransitory computer readable memory device 406, that the forwardsource wavefields should be decomposed at 608, into a set of source propagated angles, using optic flow method. Then, the set of sourcepropagated angles 609 are generated, then stored in memory resource, 404. At such point, the nontransitory computer readable memory device 406, messages the system computer 405 to retrieve, at step 610, the generated set of source propagatedangles from the memory resource 404.
Thereafter, at 611, the nontransitory computer readable memory device 406, messages the computer system device 407 to display on 409, a user interface showing that a new timestep value has been set by the nontransitory computer readable memory device 406, which equals the inputted userdefined maximum timestep value 602. Upon confirmation of said setting 611 by the user through the use of keyboard 410, and mouse 411 the nontransitory computer readable memory device 406, stores the recently set current timestep value in memory resource 404. The nontransitory computer readable memory device 406, messages the system computer 405 to retrieve from the memory resource, 404, the seismic reflection data containing the receiverwavefields. Once retrieved, the system computer messages the nontransitory computer readable memory device 406, who instead begins, at 612, backward propagating employing the finite difference method, the retrieved receiverwavefield at the set current timestep value 611 plus one incremental timestep value. Upon completion of the computing step 612, the nontransitory computer readable memory device 406, generates and stores in memory resource 404, the backward propagated receiverwavefield 613 for the system computer 405, to retrieve at 609 and message back to the nontransitory computer readable memory device 406, that the backward propagated receiverwavefields should be decomposed at 615, into a set of receiverpropagated angles, using optic flow method at the current timestep value. Then, the set of receiverpropagated angles 616 are generated, then stored in memory resource, 404. The nontransitory computer readable memory device 406, messages the system computer, 405, to retrieve 617 the generated set of receiver propagatedangles at current time step value from the memory resource. At which point, the nontransitory computer readable memory device 406, combines at 618 the retrieved set of source propagatedangle, with the retrieved set of receiver propagatedangles at current timestep value. Thereafter, the nontransitory computer readable memory device 406, begins computing at 619, a zerolag crosscorrelation value image condition for every retrieved sourcewavefield 607, and receiverwavefield 614, as well as for every combined source 610 and receiver 617 propagatedangles, at current timestep value. During step 619, the nontransitory computer readable memory device 406, is searching for similarities between the different set of signals (607, 614, 610, 617) being computed. The zerolag used by the nontransitory computer readable memory device 406, refers to the offset of the aforementioned signals, which a person of ordinary skill in the art will realize it provides the dot product for the retrieved signals.
Then, nontransitory computer readable memory device 406, generates, then stores at 620, a current image gather from the computed zerolag crosscorrelation value. The nontransitory computer readable memory device 406, messages the system computer, 405, to retrieve the generated current image gather 621, as well as the initial reverse time migration image gather 622. At which point, 623, the nontransitory computer readable memory device 406, computes the summation of the retrieved current image gather 621 to the retrieved initial reverse time migration image gather 622, and then reduces at 624, the current timestep value by the one incremental 611; and verifies at 625, whether the current image gather 620 for every time step value from the input maximum timestep value, to the timestep zero value, 603, was computed and added at 623 to the initial reverse time migration image gather, 622. The nontransitory computer readable memory device 406, messages the system computer 405, the successful execution of the verification process 625, and generates and store a final reverse time migration image gather at 626.
If the nontransitory computer readable memory device 406, messages the system computer, 405, indicating that step 625 was not successfully executed, the system computer, 405, messages the nontransitory computer readable memory device 406, to repeat subroutine 310 but starting from the step of backward propagating receiverwavefield from the retrieved seismic reflection data at the current timestep value, 612 instead. Nevertheless, a person having ordinary skills in the art would soon realize that, to save on computational cost, the user of computer system 407, may at any time stop verification 625.
According to the preferred embodiment of the invention, the subroutine in the computerimplemented method 301, that computes a residual normal moveout corrected image gather 319, is illustrated in
With subroutine 319 concluded, the system computer 405 messages the nontransitory computer readable memory device 406, to retrieve the residual normal moveout corrected image gather 320 and begin subroutine 321. This subroutine computes a final flattened image gather from the retrieved residual normal moveout corrected image gather, employing trim statics correction method; and is illustrated by
In
In
In
Unless specifically stated otherwise, terms such as “defining”, “creating”, “including”, “representing”, “preanalyzing”, “predefining”, “choosing”, “building”, “assigning”, “creating”, “introducing”, “eliminating”, “remeshing”, “integrating”, “discovering”, “performing”, “predicting”, “determining”, “inputting”, “outputting”, “identifying”, “analyzing”, “using”, “assigning”, “disturbing”, “increasing”, “adjusting”, “incorporating”, “simulating”, “decreasing”, “distributing”, “specifying”, “extracting”, “displaying”, “executing”, “implementing”, and “managing”, or the like, may refer to the action and processes of a computer system, or other electronic device, that transforms data represented as physical (electronic, magnetic, or optical) quantities within some electrical device's storage, like memory resources, or nontransitory computer readable memory, into other data similarly represented as physical quantities within the storage, or in transmission or display devices.
According the preferred embodiment of the present invention, certain hardware, and software descriptions were detailed, merely as example embodiments and are not to limit the structure of implementation of the disclosed embodiments. For example, although many internal, and external components of the computer system apparatus of
As used herein the term “survey region” refers to an area or volume of geologic interest, and may be associated with the geometry, attitude and arrangement of the area or volume at any measurement scale. A region may have characteristics such as folding, faulting, cooling, unloading, and/or fracturing that has occurred therein.
As used herein, the term “computing” encompasses a wide variety of actions, including calculating, determining, processing, deriving, investigation, look ups (e.g. looking up in a table, a database or another data structure), ascertaining and the like. It may also include receiving (e.g. receiving information), accessing (e.g. accessing data in a memory) and the like. Also, “computing” may include resolving, selecting, choosing, establishing, and the like.
As used herein, the term “trim static”, “static trimming” is used to refer to any time shift that is estimated from the crosscorrelation process and applied directly to the input trace.
As used herein, “subsurface”, and “subterranean” means beneath the top surface of any mass of land at any elevation or over a range of elevations, whether above, below or at sea level, and/or beneath the floor surface of any mass of water, whether above, below or at sea level.
Embodiments disclosed herein also relate to computerimplemented system, and computer system for performing the operations herein. This system may be specially constructed for the required purposes, or it may comprise a generalpurpose computer selectively activated or reconfigured by a computer program or code stored in the memory resources, or nontransitory computer readable memory. As such, the computer program or code may be stored or encoded in a computer readable medium or implemented over some type of transmission medium. A computerreadable medium includes any medium or mechanism for storing or transmitting information in a form readable by a machine, such as a computer (‘machine’ and ‘computer’ may be used synonymously herein). As a nonlimiting example, a computerreadable medium may include a computerreadable storage medium (e.g., read only memory (“ROM”), random access memory (“RAM”), magnetic disk storage media, optical storage media, flash memory devices, etc.). A transmission medium may be twisted wire pairs, coaxial cable, optical fiber, or some other suitable wired or wireless transmission medium, for transmitting signals such as electrical, optical, acoustical or other form of propagated signals (e.g., carrier waves, infrared signals, digital signals, etc.)).
A computer system as used herein, typically includes at least hardware capable of executing machine readable instructions, as well as the software for executing acts (typically machinereadable instructions) that produce a desired result. In addition, a computer system may include hybrids of hardware and software, as well as computer subsystems.
Hardware generally includes at least processorcapable platforms, such as clientmachines (also known as servers), and handheld processing devices (for example smart phones, personal digital assistants (PDAs), or personal computing devices (PCDs)). Further, hardware may include any physical device that can store machinereadable instructions, such as memory or other data storage devices. Other forms of hardware include hardware subsystems, including transfer devices such as modems, modem cards, ports, and port cards, for example.
Software includes any machine code stored in any memory medium, such as RAM or ROM, and machine code stored on other devices (such as nontransitory computer readable media like external hard drives, or flash memory, for example). Software may include source or object code, encompassing any set of instructions capable of being executed in a client machine, server machine, remote desktop, or terminal.
Combinations of software and hardware could also be used for providing enhanced functionality and performance for certain embodiments of the disclosed invention. One example is to directly manufacture software functions into a silicon chip. Accordingly, it should be understood that combinations of hardware and software are also included within the definition of a computer system and are thus envisioned by the invention as possible equivalent structures and equivalent methods.
Computerreadable mediums or memory resources include passive data storage, such as a randomaccess memory (RAM) as well as semipermanent data storage such as external hard drives, and external databases, for example. In addition, an embodiment of the invention may be embodied in the RAM of a computer to transform a standard computer into a new specific computing machine.
Data structures are defined organizations of data that may enable an embodiment of the invention. For example, a data structure may provide an organization of data, or an organization of executable code. Data signals could be carried across nontransitory transmission mediums and stored and transported across various data structures, and, thus, may be used to transport an embodiment of the invention.
The system computer may be designed to work on any specific architecture. For example, the system may be executed on a highperformance computing system, which typically comprise the aggregation of multiple single computers, physically connected, or connected over local area networks, clientserver networks, wide area networks, internets, handheld and other portable and wireless devices and networks.
An “output device” includes the direct act that causes generating, as well as any indirect act that facilitates generation. Indirect acts include providing software to an user, maintaining a website through which a user is enabled to affect a display, hyperlinking to such a website, or cooperating or partnering with an entity who performs such direct or indirect acts. Thus, a user may operate alone or in cooperation with a thirdparty vendor to enable the reference signal to be generated on a display device. A display device may be included as an output device, and shall be suitable for displaying the required information, such as without limitation a CRT monitor, a LCD monitor, a plasma device, a flat panel device, or printer. The display device may include a device which has been calibrated through the use of any conventional software intended to be used in evaluating, correcting, and/or improving display results (e.g., a color monitor that has been adjusted using monitor calibration software). Rather than (or in addition to) displaying the reference image on a display device, a method, consistent with the invention, may include providing a reference image to a subject. “Providing a reference image” may include creating or distributing the reference image to the subject by physical, telephonic, or electronic delivery, providing access over a network to the reference, or creating or distributing software to the subject configured to run on the subject's workstation or computer including the reference image. In one example, providing of the reference image could involve enabling the subject to obtain the reference image in hard copy form via a printer. For example, information, software, and/or instructions could be transmitted (e.g., electronically or physically via a data storage device or hard copy) and/or otherwise made available (e.g., via a network) in order to facilitate the subject using a printer to print a hard copy form of reference image. In such an example, the printer may be a printer which has been calibrated through the use of any conventional software intended to be used in evaluating, correcting, and/or improving printing results (e.g., a color printer that has been adjusted using color correction software).
A database, or multiple databases may comprise any standard or proprietary database software, such as Oracle, Microsoft Access, SyBase, or DBase II, for example. The database may have fields, records, data, and other database elements that may be associated through database specific software. Additionally, data may be mapped. Mapping is the process of associating one data entry with another data entry. For example, the data contained in the location of a character file can be mapped to a field in a second table. The physical location of the database is not limiting, and the database may be distributed. For example, the database may exist remotely from the server, and run on a separate platform. Further, the database may be accessible across the a local network, a wireless network of the Internet.
Furthermore, modules, features, attributes, methodologies, and other aspects can be implemented as software, hardware, firmware or any combination thereof. Wherever a component of the invention is implemented as software, the component can be implemented as a standalone program, as part of a larger program, as a plurality of separate programs, as a statically or dynamically linked library, as a kernel loadable module, as a device driver, and/or in every and any other way known now or in the future to those of skill in the art of computer programming. Additionally, the invention is not limited to implementation in any specific operating system or environment.
Various terms as used herein are defined below. To the extent a term used in a claim is not defined below, it should be given the broadest possible definition persons in the pertinent art have given that term as reflected in at least one printed publication or issued patent.
As used herein, “and/or” placed between a first entity and a second entity means one of (1) the first entity, (2) the second entity, and (3) the first entity and the second entity. Multiple elements listed with “and/or” should be construed in the same fashion, i.e., “one or more” of the elements so conjoined
Additionally, the flowcharts and block diagrams in the Figures illustrate the architecture, functionality, and operation of possible implementations of systems, methods, and computer program products according to various embodiments of the present disclosure. It should also be noted that, in some alternative implementations, the functions noted in the block may occur out of the order noted in the Figures. For examples, two blocks shown in succession may, in fact, be executed substantially concurrently, or the blocks may sometimes be executed in the reverse order, depending upon the functionality involved. It will also be noted that each block of the block diagrams and/or flowcharts illustration, and combinations of blocks in the block diagrams and/or flowchart illustration, can be implemented by special purpose hardwarebased systems that perform the specified hardware functions or acts, or combinations of special purpose hardware and computer instructions.
While in the foregoing specification this disclosure has been described in relation to certain preferred embodiments thereof, and many details have been set forth for purpose of illustration, the invention is not to be unduly limited to the foregoing which has been set forth for illustrative purposes. On the contrary, a wide variety of modifications and alternative embodiments will be apparent to a person skilled in the art, without departing from the true scope of the invention, as defined in the claims set forth below. Additionally, it should be appreciated that structural features or method steps shown or described in any one embodiment herein can be used in other embodiments as well.
Claims
1. A computerimplemented method for small cave recognition using seismic reflection data in a survey region the method comprising:
 retrieving a surface seismic reflection data represented in timedomain, comprising data from seismic traces, acquired from a plurality of field receiving sensors;
 retrieving a well log data represented in timedomain, from a well hole containing a small cave;
 computing a final migration velocity model from the retrieved surface seismic reflection data and the retrieved well log data;
 retrieving a generated final migration velocity model from the computed final migration velocity model;
 computing a source wavelet from the retrieved surface seismic reflection data;
 generating the source wavelet from the computed source wavelet;
 retrieving the source wavelet from the generated source wavelet;
 computing a final reverse time migration image gather, using the retrieved surface seismic reflection data, the retrieved final migration velocity model, and the generated source wavelet, employing a finite difference method;
 retrieving the final reverse time migration image gather in depthdomain from the computed final reverse time migration image gather;
 converting the retrieved final reverse time migration image gather, from the depthdomain to the timedomain;
 generating a final reverse time migration image gather in timedomain from the converted final reverse time migration image gather;
 retrieving the final reverse time migration image gather in timedomain from the generated final reverse time migration image gather in timedomain;
 computing an initial rootmeansquare velocity model in timedomain from the retrieved final migration velocity model;
 generating a final rootmeansquare velocity model in timedomain from the computed initial rootmeansquare velocity model;
 retrieving the final rootmeansquare velocity model in timedomain from the generated initial rootmeansquare velocity model;
 computing a residual normal moveout corrected image gather from the retrieved final reverse time image gather;
 retrieving a residual normal moveout corrected image gather from the computed residual normal moveout corrected image gather;
 computing a final flattened image gather from the retrieved residual normal moveout corrected image gather, employing trim statics correction method;
 retrieving the generated final flattened image gather from the computed final flattened image gather;
 computing a final intercept volume and a final gradient volume using the retrieved final flattened image gather;
 retrieving a generated final intercept volume and a generated final gradient volume from the computed final intercept volume, and the computed final gradient volume;
 computing a multiplication of the final intercept volume and the final gradient volume; and
 generating a small cave identification model from the multiplied intercept and gradient volumes.
2. A computing system apparatus programmed to perform a set of operations of a computerimplemented method for small cave recognition using seismic reflection data, comprising:
 a computer system device for inputting, setting, selecting, outputting, and performing the operations of retrieving, computing, generating, invoking, determining, converting, and correcting;
 a memory resource, for storing data corresponding to the operations of inputting, generating and retrieving, a surface seismic reflection data in timedomain, a well log data represented in timedomain, a final migration velocity model, a source wavelet, a final reverse time migration image gather in depthdomain, a final reverse time migration image gather in timedomain, a final rootmeansquare velocity model in timedomain, a residual normal moveout corrected image gather, a residual normal moveout corrected image gather, a final flattened image gather, an initial intercept volume, an initial gradient volume, a final intercept volume, a final gradient volume, a small cave identification model, a userdefined flatness value, an initial migration velocity model, an initial Kirchhoff prestack depth migration image gather, a set of initial moveout curves, a Kirchhoff prestack depth migration stacking image, a final reflector dipangle image, a tomographic inversion model, a final Kirchhoff prestack depth migration gather, a generated flatness value, a current timestep value, a maximum timestep value, a source wavelet, a forward sourcewavefield, a set of source propagatedangles for every sourcewavefield, a backward receiverwavefield, a set of receiver propagatedangles for every receiverwavefield, a current image gather, a corrected image gather, a final rootmeansquare velocity, a set of positive and negative residual normal moveout curves, an initial image gather, a set of relative time shifts for all seismic traces, and a trim statics correction time shift method;
 a nontransitory program storage computerreadable memory for performing the operations of inputting, setting, selecting, outputting, retrieving, computing, generating, invoking, determining, converting, and correcting;
 a computer system output device for outputting the operations of retrieving, computing, generating, invoking, determining, converting, and correcting; and
 a system computer, coupled to a computer system input device, coupled to a memory resource, coupled to a nontransitory computer readable memory device, and couple to the computer system output device, for performing the operations of: retrieving a surface seismic reflection data represented in timedomain, comprising data from seismic traces, acquired from a plurality of field receiving sensors; retrieving a well log data represented in timedomain, from a well hole containing a small cave; computing a final migration velocity model from the retrieved surface seismic reflection data and the retrieved well log data; retrieving a generated final migration velocity model from the computed final migration velocity model from the computed initial migration velocity model; computing a source wavelet from the retrieved surface seismic reflection data; generating the source wavelet from the computed source wavelet; retrieving the source wavelet from the generated source wavelet; computing a final reverse time migration image gather, using the retrieved surface seismic reflection data, the retrieved final migration velocity model, and the generated source wavelet, employing a finite difference method; retrieving the final reverse time migration image gather in depthdomain from the computed final reverse time migration image gather; converting the retrieved final reverse time migration image gather, from the depthdomain to the timedomain; generating a final reverse time migration image gather in timedomain from the converted final reverse time migration image gather; retrieving the final reverse time migration image gather in timedomain from the generated final reverse time migration image gather in timedomain; computing an initial rootmeansquare velocity model in timedomain from the retrieved final migration velocity model; generating a final rootmeansquare velocity model in timedomain from the computed initial rootmeansquare velocity model; retrieving the final rootmeansquare velocity model in timedomain from the generated initial rootmeansquare velocity model; computing a residual normal moveout corrected image gather from the retrieved final reverse time image gather; retrieving a residual normal moveout corrected image gather from the computed residual normal moveout corrected image gather; computing a final flattened image gather from the retrieved residual normal moveout corrected image gather, employing trim statics correction method; retrieving the generated final flattened image gather from the computed final flattened image gather; computing a final intercept volume and a final gradient volume using the retrieved final flattened image gather; retrieving a generated final intercept volume and a generated final gradient volume from the computed final intercept volume, and the computed final gradient volume; computing a multiplication of the final intercept volume and the final gradient volume; and generating a small cave identification model from the multiplied intercept and gradient volumes.
3. The computing system apparatus of claim 2, wherein the nontransitory program storage computerreadable memory is further programmed to perform the operation of computing a final migration velocity model from the retrieved surface seismic reflection data and the retrieved well log data, further including:
 inputting a userdefined flatness value for an image gather;
 retrieving the userdefined flatness value for an image gather from the memory resource;
 computing an initial migration velocity model in depthdomain with the retrieved surface seismic reflection data and the retrieved well log data;
 generating an initial migration velocity model in depthdomain from the computed initial migration velocity model;
 retrieving the generated initial migration velocity model in depthdomain from the memory resource;
 computing an initial Kirchhoff prestack depth migration image gather, using the retrieved initial migration velocity model;
 generating an initial Kirchhoff prestack depth migration image gather from the computed initial Kirchhoff prestack depth migration image gather;
 retrieving the generated initial Kirchhoff prestack depth migration image gather from the memory resource;
 computing a set of baseline moveout curves using the retrieved initial Kirchhoff prestack depth migration gather;
 generating a set of initial moveout curves from the computed set of baseline moveout curves;
 retrieving the generated set of baseline moveout curves from the memory resource;
 computing a Kirchhoff prestack depth migration stacking image, from the retrieved initial Kirchhoff prestack depth migration gather, and the retrieved set of baseline moveout curves;
 generating a Kirchhoff prestack depth migration stacking image from the computed Kirchhoff prestack depth migration stacking image;
 retrieving the generated Kirchhoff prestack depth migration stacking image from the memory resource;
 computing an initial reflector dipangle image from the generated Kirchhoff prestack depth migration stacking image, employing local slantstacking algorithm;
 generating a final reflector dipangle image from the computed initial reflector dipangle image;
 retrieving the generated final reflector dipangle image from the memory resource;
 computing a tomographic inversion model with the retrieved set of initial moveout curves and the retrieved final reflector dipangle image;
 generating a tomographic inversion model from the computed tomographic inversion model;
 retrieving the generated tomographic inversion model from the memory resource;
 updating the generated initial migration velocity model with the retrieved tomographic inversion model;
 computing a final Kirchhoff prestack depth migration image gather from the updated initial migration velocity model;
 generating a final Kirchhoff prestack depth migration image gather, from the computed final Kirchhoff prestack depth migration image gather;
 retrieving the generated final Kirchhoff prestack depth migration image gather from the memory resource;
 generating a flatness value from the retrieved final Kirchhoff prestack depth migration image gather;
 repeating the steps of computing an initial Kirchhoff prestack depth migration image gather, generating an initial Kirchhoff prestack depth migration image gather, retrieving the initial Kirchhoff prestack depth migration image gather, computing a set of baseline moveout curves, generating a set of baseline moveout curves, retrieving the set of baseline moveout curves, computing a Kirchhoff prestack depth migration stacking image, generating a Kirchhoff prestack depth migration stacking image, retrieving the Kirchhoff prestack depth migration stacking image, computing an initial reflector dipangle, generating an initial reflector dipangle, retrieving the an initial reflector dipangle, computing a tomographic inversion model, generating a tomographic inversion model, retrieving the tomographic inversion model, computing a final Kirchhoff prestack depth migration image gather, generating a final Kirchhoff prestack depth migration image gather, retrieving the final Kirchhoff prestack depth migration image gather, and generating a flatness values; until the generated flatness values is equal to, or less than, the userdefined flatness value; and
 generating a final migration velocity model.
4. The computing system apparatus of claim 2, wherein the nontransitory program storage computerreadable memory is encoded for computing a final reverse time migration image gather, using the retrieved surface seismic reflection data, the retrieved final migration velocity model, and the computed source wavelet, employing a finite difference method, further comprising the steps of:
 inputting a timestep value of zero, and a userdefined maximum time step value;
 retrieving the inputted userdefined maximum time step value from the memory resource;
 setting an initial reverse time migration image gather equal to zero and at time step value equal to zero;
 retrieving a source wavelet from the memory resource;
 forward propagating the sourcewavefield from the retrieved seismic reflection data, with the retrieved source wavelet using a finite difference method, from the retrieved timestep value of zero to the retrieved input maximum timestep value;
 generating a forward sourcewavefield from the forward propagated sourcewavefield, from the timestep value of zero to the input maximum timestep value;
 retrieving the generated forward sourcewavefield from the memory resource;
 decomposing the retrieved forward sourcewavefield into a set of source propagatedangles, using an optic flow method;
 generating a set of source propagatedangles, from the decomposed set of forward sourcewavefield for every sourcewavefield retrieved from the surface seismic reflection data, from the retrieved timestep value of zero to the retrieved userdefined maximum timestep value;
 retrieving the generated set of source propagatedangles from the memory resource;
 setting a current timestep value, equal to the userdefined maximum time step value;
 backward propagating the receiverwavefield from the retrieved seismic reflection data, at the set current timestep value plus one incremental timestep value, and employing a finite difference method;
 generating a backward receiverwavefield from the backward propagated receiverwavefield at the current timestep value;
 retrieving the generated backward receiverwavefield at the current timestep value, from the memory resource;
 decomposing the retrieved backward receiverwavefield, at the current timestep value, into receiver propagatedangles, using an optic flow method;
 generating a set of receiver propagatedangles at the current timestep value, from the decomposed set of backward receivedwavefield at the current timestep value;
 retrieving the generated set of receiver propagatedangles at current timestep value, from the memory resource;
 combining the retrieved set of source propagatedangles at current timestep value, with the retrieved set of receiver propagatedangles, at current timestep value;
 computing a zerolag crosscorrelation value image conditioning for every retrieved sourcewavefield and receiverwavefield, as well as for every combined source and receiver propagatedangle, at current timestep values;
 generating a current image gather from the computed zerolag crosscorrelation value at current timestep value;
 retrieving the generated current image gather from the memory resource;
 retrieving the initial reverse time migration image gather from the memory resource;
 adding the retrieved current image gather, to the retrieved initial reverse time migration image gather;
 reducing the current timestep value by a one incremental step;
 repeating the steps of backward propagating receiverwavefield from the retrieved seismic reflection data at the current timestep value, generating a backward receiverwavefield at the current timestep value, retrieving the generated backward receiverwavefield at the current timestep value, decomposing the retrieved backwardreceiverwavefield into receiver propagatedangles at the current timestep value, generating a set of receiver propagatingangles from the decomposed backward receivedwavefield at the current timestep value, retrieving the generated set of receiver propagatingangles at the current timestep value, combining the retrieved set of source propagatingangles at the current timestep value, with the retrieved set of receiver propagatingangles, computing a zerolag crosscorrelation value image conditioning, generating a current image gather from the computed zerolag crosscorrelation value at current timestep value, retrieving the generated current image gather, retrieving the initial reverse time migration image gather, and adding the current angle domain common image gather, to the retrieved initial reverse time migration image gather; until the set of current image gather for every time step value from the input maximum timestep values, to the timestep zero value, was computed and added to the initial reverse time migration image gather; and
 generating a final reverse time migration image gather.
5. The computing system apparatus of claim 2, wherein the nontransitory program storage computerreadable memory is encoded for computing a residual normal moveout corrected image gather from the corrected image gather, further including the steps of:
 retrieving the final reverse time migration image gather in timedomain from the memory resource;
 inputting a set of userdefined rootmeansquare velocities;
 retrieving the set of userdefined rootmeansquare velocities from the memory resource;
 computing a set of positive and negative residual normal moveout curves from the retrieved final reverse time migration image gather in timedomain with the retrieved set of userdefined rootmeansquare velocities;
 generating a set of positive and negative residual normal moveout curves from the computed set of positive and negative residual normal moveout curves;
 retrieving the generated set of positive and negative residual normal moveout curves from the memory resource;
 computing a semblance spectrum for the retrieved set of positive and negative residual normal moveout curves;
 selecting a set of residual normal moveout curves corresponding to the computed semblance spectrum with peak values, from the generated set of positive and negative residual normal moveout curves;
 computing a residual normal moveout corrected image gather from the final reverse time migration image gather in timedomain with the selected set of residual normal moveout curves corresponding to the peak semblance spectrum; and
 generating a residual normal moveout corrected image gather, from the computed residual normal moveout corrected image gather.
6. The computing system apparatus of claim 2, wherein the nontransitory program storage computerreadable memory is encoded for computing a final flattened image gather from the retrieved residual normal moveout corrected image gather, employing trim statics correction method, further comprising the steps of:
 computing a set of relative time shifts between seismic traces from the retrieved residual normal moveout corrected image gather, employing a crosscorrelation method;
 generating a set of relative time shifts for all seismic traces from the computed set of relative time shifts;
 retrieving the generated set of relative time shifts for all seismic traces from the memory resource;
 computing a set of trim statics correction timeshifts for the retrieved set of relative time shifts for all seismic traces, using the leastsquares inversion;
 generating a trim statics correction timeshift model for each seismic trace from the computed set of trim statics timeshifts;
 retrieving the generated trim statics correction time shift model from the memory resource;
 computing a flattened image gather from the retrieved initial image gather, applying the retrieved trim statics correction time shift model to each seismic trace in the retrieved initial image gather; and
 generating a final flattened image gather from the computed set of flattened image traces.
7. The computing system apparatus of claim 2, wherein the nontransitory program storage computerreadable memory is further programmed for computing a final intercept volume and a final gradient volume using the retrieved final flattened image gather, includes the steps of:
 computing an amplitude versus angle volume from the retrieved final flattened image gather;
 generating a final amplitude versus angle volume from the computed initial amplitude versus angle volumes;
 retrieving the generated final amplitude versus angle volume from the memory resource;
 computing a final intercept volume and a final gradient volume from the retrieved final amplitude versus angle volumes, using least squares regression analysis; and
 generating a final intercept volume and a final gradient volume from the computed final intercept volume and the computed final gradient volume.
8. The computing system apparatus of claim 3, wherein the computer system device is further programmed to display a userinterface;
9. The computing system apparatus of claim 3, wherein the computer system device is further programmed to display and print a surface seismic reflection data represented in timedomain, a well log data represented in timedomain, an initial migration velocity model, a final migration velocity model, an image gather, a source wavelet, a final reverse time migration image gather in depthdomain, a final reverse time migration image gather in timedomain, a residual normal moveout corrected image gather, a final flattened image gather, a final amplitude versus angle volume, a set of initial moveout curves, an initial and a final Kirchhoff prestack depth migration gathers, a Kirchhoff prestack depth migration stacking image, an initial reflector dipangle image, a final reflector dipangle image, an updated initial velocity migration model, a user input flat value, a final intercept volume, a final gradient volume, and a small cave identification model.
Type: Application
Filed: Feb 17, 2020
Publication Date: Aug 19, 2021
Patent Grant number: 11143776
Applicants: China Petroleum & Chemical Corporation (Beijin), Sinopec Tech Houston (Houston, TX)
Inventors: Min Zhou (Houston, TX), Franklin Ruiz (Houston, TX)
Application Number: 16/792,698