Method and apparatus for full waveform inversion based on adaptive cross-correlation

Abstract

Methods and apparatus are provided for performing traveltime-based seismic full waveform inversion using a frequency-dependent Hanning window and a source-independent approach to generate a final velocity model of a subsurface formation of a survey region. The method includes deploying seismic data recording sensors at the survey region; performing a blast at a shot point of the survey region to generate seismic waves; and sensing and recording the seismic waves using the seismic data recording sensors. The recorded seismic waves are seismic data. The method also includes transferring the seismic data to a computer system including one or more memories and storing the seismic data in the one or more memories; storing a source wavelet in the one or more memories; performing, by the computer system, a forward modeling operation using the source wavelet; and generating, by the computer system, the final velocity model using seismic full waveform inversion and the forward modeling operation.

Description

TECHNICAL FIELD

[0001] The present disclosure relates to a seismic exploration method and apparatus for building a high resolution geologic model by performing an adaptive cross-correlation based full waveform inversion (FWI) to enhance characterization of complex subsurface structures within an exploration area. The method first automatically evaluates traveltime differences between modeled synthetic data and field observed data by means of a frequency dependent sliding Hanning window scan based adaptive cross-correlation computation and by means of a shot independent scheme. The evaluated traveltime differences are inverted to generate a final velocity model to improve the image of complex subsurface structures within the survey area. BACKGROUND

[0002] 1. SUMMARY

[0003] In oil and gas exploration, high resolution subsurface volume properties such as seismic propagation velocity, anisotropy, absorption, porosity, and reflectivity models can be generated by acquiring and processing seismic data. These geophysical properties, in combination, can effectively reveal subsurface structures. Seismic data processing typically includes seismic inversion to build a propagation velocity model at medium to long wavenumbers, followed by seismic migration to obtain seismic reflectivity images at short wavenumbers. These seismic reflectivity images are used to determine the location and size of natural resource (such as hydrocarbon of oil and gas) reservoirs, thus providing the basis for exploration / drilling planning. To obtain high resolution seismic reflectivity images, a high fold acquisition system needs to be designed to obtain sufficient seismic data before seismic migration, as well as a good seismic velocity model with correct kinematic information.

[0004] Existing seismic velocity inversion methods include ray-based seismic tomography and full waveform inversion (FWI) methods. Ray-based seismic tomography is efficient and can invert smooth models, which can be sufficient for relatively simple geologic targets such as shallow sedimentary environments. However, for complex geologic environments such as salt domes, sub-basalt targets, thrust belts, and land foothills, ray-based tomography is less effective, and thus FWI becomes a necessary tool for building velocity models.

[0005] FWI directly solves the seismic wave equation and matches seismic data, which can generate more accurate seismic propagation velocity models for complex subsurface structures such as salt-related models. Such propagation velocity models can be used to generate accurate, high resolution seismic reflectivity images by seismic migration, to facilitate time-lapse monitoring of hydrocarbon reservoirs, and can even be directly converted to generate a high resolution seismic image volume, which is called a FWI image.

[0006] 2. SEISMIC WAVE EQUATION

[0007] Seismic waves propagating below the Earth's surface can be simulated by solving the seismic wave equation. The seismic wave equation describes the Earth with different physical models and assumes the Earth to be either isotropic or anisotropic, elastic or acoustic, and attenuating or non-attenuating. In most FWI developments, seismic waves are assumed to be purely acoustic because the acoustic wave equation is relatively simple and can be solved efficiently. Regardless of the assumptions, all wave equations can be mathematically represented as:

[0008] F(m; x)w s (x, t) = f s (t), (1)

[0009] where the vector m is the Earth model, i.e., a representation of the Earth's subsurface properties such as the distribution of seismic wave velocities, densities, and other physical properties; x is the spatial location; t is the time; F(m; x) is the corresponding (forward) modeling operator; and w s (x, t) is the forward wavefield of the source wavelet f S (t) excited at location s.

[0010] The numerical solution of equation (1) discretizes both the spatial variable x and the time variable t, resulting in a solution w s (x, t) being discretized into w s (x i , t j ), where i = 1,..., N, which is the grid point, and there are N of them, and j = 1,..., Nt, which is the time step, and there are Nt of them. Common numerical methods include finite-difference, finite-element, and spectral-element methods. As shown in the present disclosure, the methods disclosed herein are independent of the form of the wave equation (1) and the numerical methods used to solve the equation.

[0011] 3. Full-waveform inversion (FWI)

[0012] FWI is a data-driven tool that automatically builds the subsurface parameters m, such as velocities and / or densities, by iteratively minimizing the difference between the recorded data and the modeled synthetic data. Given an initial estimate of the subsurface velocities m0, the synthetic data can be predicted by solving the forward seismic wave equation (1) with the source wavelet f S (t). With the residuals between the data and the synthetic data as the source, the adjoint wave equation is solved to obtain the adjoint wave equation's solution and the gradient is obtained by cross-correlating the solutions of the forward and adjoint wave equations, and the gradient is then used to update the model in a certain direction to reduce the misfit between the modeled synthetic data and the field observed data. This iterative solution is repeated until the data misfit is small enough.

[0013] Mathematically, FWI based on the least-squares (L2) waveform difference adopts the following misfit:

[0014]

[0015] where C(m) is the misfit function that measures the waveform difference between the synthetic data and the recorded data, and the modulus operation is denoted by the symbol | · |. In the misfit function, N s is the number of sources, N r is the number of receivers for a given source, t max is the maximum recording time starting from 0, d obs,s,r (t) is the recorded data at receiver r for a given source s at time t, d syn,s,r (m; t) is the synthetic data at the same receiver r for the same source s at time t, obtained by solving the forward seismic wave equation (1), where d syn,s,r (m; t) = P s,r w s (m; t), where P s,r resamples the wavelet w s (t) to the receiver location r. The misfit function (2) is called the least-squares or L2 norm misfit function of the waveform difference, and is the most common misfit function in the oil and gas industry.

[0016] 4. Traveltime-based FWI

[0017] Traveltime-based FWI is based on picking the traveltime difference between the modeled synthetic data and the field observed data. The traveltime difference can be obtained by automatically picking the maximum cross-correlation value between the synthetic data and the recorded data, i.e.,

[0018]

[0019] On this basis, the traveltime-based FWI can be formulated as the least-squares (L2) minimization of the traveltime difference as follows:

[0020]

[0021] Compared with the least-squares cost function of the waveform difference in equation (2), it is more linear in the difference between the initial velocity model and the true velocity model. It has the potential to alleviate the cycle skipping problem without requiring the model synthetic events to be at most half a wavelength away from the true observed data events.

[0022] Furthermore, the amplitudes of the modeled synthetic data can easily be inconsistent with the amplitudes of the recorded data because seismic forward modeling cannot take into account all the subsurface physical properties and can miss elastic and / or attenuation characteristics. For efficient computation, an acoustic modeling engine is typically used to simulate the propagation of seismic waves. Furthermore, seismic amplitudes are easily contaminated during seismic processing. Traveltime-based misfit functions can attenuate the amplitude differences, thus reducing the negative impact of the amplitude differences during seismic velocity inversion.

[0023] However, traveltime-based FWI can be affected by inaccurate and low-resolution traveltime picks. First, the traveltime picks can be inaccurate due to cross-talk between multiple seismic events that are not properly matched. Second, the resolution of the traveltime picks can be low because the calculations in equation (3) are typically dominated by strong events, while relatively weak events are not taken into account.

[0024] Furthermore, the accuracy of the source wavelet f s (t) can also reduce the accuracy of the cross-correlation-based automatic traveltime picks. In particular, in land applications, the source wavelet can vary from shot to shot, thus making it difficult to obtain an accurate source wavelet.

[0025] Therefore, there is still a need for new methods to overcome the shortcomings of conventional traveltime-based FWI. SUMMARY

[0026] In one embodiment of the present disclosure, a method for performing seismic full waveform inversion to generate a velocity model of a subsurface formation of a survey region includes a plurality of steps of positioning seismic data recording sensors at different locations within the survey region, and / or positioning a logging tool including the seismic data recording sensors in a borehole of the survey region; making a transmission at an incident point in the survey region to produce seismic waves that propagate through the subsurface formation; observing the seismic waves using the seismic data recording sensors and recording seismic data from the seismic waves; transmitting the seismic data observed by the seismic data recording sensors to a computer system including one or more memories and storing the observed seismic data in the one or more memories, and storing a source wavelet and a current medium parameter model in the one or more memories; performing, by the computer system, a forward modeling operation according to a forward wave equation and using the source wavelet and the current medium parameter model and obtaining a forward wavefield; performing, by the computer system, a Han window operation on the synthetic data and the observed seismic data, respectively, to obtain local synthetic data and local seismic data; selecting, by the computer system, a traveltime difference according to a cross-correlation between the local synthetic data and the local seismic data; solving, by the computer system, an adjoint equation of the forward wave equation using the traveltime difference and the local seismic data to obtain an adjoint source; generating, by the computer system, an updated medium parameter model for the seismic full waveform inversion and the synthetic data using the forward modeling operation; performing operations (e) to (h) until convergence; outputting, after the convergence, the updated velocity as a final velocity model to a display; and displaying, on a display of the computer system, a high resolution image of the final velocity model.

[0027] Other embodiments provide FWI methods based on adaptive cross-correlation. In one embodiment, the traveltime difference that is automatically picked from the adaptive cross-correlation of local seismic windows is not fitted. The adaptive cross-correlation indicates that there is no need to iterate to determine the best local window size, nor to estimate an accurate source wavelet for the seismic modeling data. The seismic data is automatically localized by a frequency dependent Han sliding window, while the source wavelet f s (t) the error is corrected by a Wiener filter.

[0028] Also disclosed is a system for performing seismic full waveform inversion to generate a velocity model of a subsurface formation of a survey region. The system includes a plurality of seismic data recording sensors disposed at different locations in the survey region, and / or a logging tool including seismic data recording sensors disposed within a borehole in the survey region; a shot device disposed at each shot point within the survey region to generate seismic waves that propagate through the subsurface formation; and a plurality of seismic data recording sensors to sense the seismic waves and record seismic data from the seismic waves. The seismic data recording sensors transmit the seismic data to a computer system including one or more memories storing the transmitted seismic data, a shot wavelet, and instructions, and at least one processor executing the instructions stored in the one or more memories to implement one of the methods disclosed in the present disclosure.

[0029] While the embodiments disclosed herein describe a FWI algorithm with a traveltime (misfit) function that utilizes a finite-difference numerical solution of the scalar acoustic wave equation for seismic propagation, it will be apparent to those of ordinary skill in the art that the algorithm can alternatively be applied to the vector wave equation and elastic equation in isotropic and anisotropic media without departing from the true scope of the invention as defined by the claims below. BRIEF DESCRIPTION OF DRAWINGS

[0030] The teachings of the present disclosure can be readily understood by considering the following detailed description in conjunction with the accompanying drawings.

[0031] Figure 1 FIG. 1 is a schematic diagram showing a top view of a survey region with different seismic source shot points, according to one embodiment of the present disclosure.

[0032] Figure 2 FIG. 2 is a schematic diagram showing a cross-sectional view of an environment with a seismic source shot point, a seismic data recording sensor, a well site, a wellbore, various transmission rays, and various angles of incidence, according to one embodiment.

[0033] Figure 3 FIG. 3 is a schematic diagram showing a cross-sectional view of an environment with a wellbore and a logging tool including one or more acoustic wave generators and one or more logging data recording sensors, according to one embodiment.

[0034] Figure 4 FIG. 4 is a schematic diagram showing a high performance computing system, according to one embodiment.

[0035] Figure 5 FIG. 5 is a flowchart showing a FWI process using cross-correlation FWI based on a local Hanning window, according to the present disclosure.

[0036] Figure 6is a flowchart showing a FWI flow using local Hanning window and source-independent cross-correlation FWI according to the present disclosure.

[0037] Figure 7 shows a sliding Hanning window, the size and sliding step of which are automatically determined according to the data frequency of the localized seismic data.

[0038] Figure 8 shows a true velocity model with a grid spacing of 50 meters, used to generate seismic data with a maximum offset of 25 kilometers.

[0039] Figure 9 shows an initial velocity model used for the FWI test.

[0040] Figure 10 shows the velocity difference obtained after subtracting the initial velocity model shown in Figure 8 from the true velocity model shown in Figure 9 ranging from about -1365 m / s to 1980 m / s.

[0041] Figure 11 shows the final inverted velocity model with band-pass filtered 2.5 Hz data using the adaptive cross-correlation FWI method of the present invention, where the effective lowest frequency is assumed to be 2.5 Hz, and there is no lower frequency in the observed data.

[0042] Figure 12 shows the final inverted velocity model with continuously band-pass filtered 2.5 Hz, 4 Hz, 6 Hz, and 10 Hz data using the adaptive cross-correlation FWI method of the present invention.

[0043] Figure 13 shows the final inverted velocity model with continuously band-pass filtered 2.5 Hz, 4 Hz, 6 Hz, and 10 Hz data using the conventional L2 waveform difference based FWI method.

[0044] Figure 14 shows a velocity model of a real land data set generated by a ray tomography technique, which will serve as the initial velocity model for FWI.

[0045] Figure 15 shows the inversion model perturbation (the difference between the inverted model of FWI and the initial model of tomography), which is inverted using the adaptive cross-correlation based FWI method of the present invention.

[0046] Figure 16 shows the seismic migration image obtained using the ray-based initial tomography model.

[0047] Figure 17Seismic migrated images obtained using the inverted velocity model by the adaptive cross-correlation FWI method of the present invention are shown.

[0048] The meaning of the reference signs in the drawings is as follows:

[0049]

[0050] DETAILED DESCRIPTION

[0051] Several embodiments of the present disclosure will now be described in detail with reference to the drawings, which are by way of illustration only. It is noted that like or similar features are referred to using like or similar reference signs in the drawings. The embodiments of the present disclosure are shown by way of illustration only in the drawings. Alternative embodiments of the structures, systems and methods illustrated herein will be readily appreciated by one skilled in the art from this description.

[0052] Throughout the specification, the terms "method" and "approach" are used interchangeably and have the same meaning. The definitions of the symbols of the equations throughout the specification are listed in the table below.

[0053] The present disclosure relates to establishing high resolution geologic models by performing improved seismic full waveform inversion to improve the image of complex subsurface structures (strata) in a survey area by applying methods, apparatuses and media that include one or more source-independent misfit functions.

[0054] Figures 1 to 4 Exemplary embodiments of methods, apparatuses and media for acquiring and storing seismic data that are processed to generate one or more high resolution geologic models for high resolution images for lithology identification, fluid discrimination and reservoir characterization of complex subsurface structures in a survey area are shown. The survey area can be subsurface structures under land or under the ocean.

[0055] Figures 5 to 17 Exemplary embodiments of apparatuses, methods and media that improve the quality of seismic full waveform inversion results by using improved seismic full waveform inversion techniques to improve lithology identification, fluid discrimination and reservoir characterization in the field of seismic exploration, including using computer-implemented seismic full waveform inversion methods, are shown. Figures 5 to 17 Exemplary embodiments of apparatuses, methods and media for generating one or more high resolution geologic models for high resolution images for lithology identification, fluid discrimination and reservoir characterization of complex subsurface structures in a survey area are shown. Figures 5 to 17Exemplary embodiments are shown that generate inverse model parameters that utilize frequency dependent sliding Hanning windows and are independent of source wavelet form for performing seismic full waveform inversion to generate high resolution geologic models for high resolution images of survey regions that include complex subsurface structures. Figures 5 to 17 Exemplary embodiments are shown that utilize frequency dependent sliding Hanning windows and are independent of source wavelet form that can generate one or more high resolution geologic models for high resolution imaging to perform lithology identification, fluid discrimination, and reservoir characterization of complex subsurface structures of survey regions.

[0056] Figure 1 is a schematic diagram showing a top view of a survey region that includes a plurality of shot points of seismic sources in accordance with one embodiment. More specifically, Figure 1 A seismic survey region (survey region) 101 is shown that is a land-based region with reference numeral 102 representing the top formation of the land-based region. Those of ordinary skill in the art will recognize that seismic survey regions can generate detailed images of the local geology to determine the location and extent of possible hydrocarbon (oil and gas) reservoirs to determine well sites 103. In these survey regions, seismic waves are bounced off subsurface formations when emitted from one or more sources at different shot points 104. Blasts are one example of sources produced by seismic equipment. Seismic waves that are reflected back to the surface are captured by seismic data recording sensors 105, transmitted (usually wirelessly) from the seismic data recording sensors 105 by one or more data transmission systems, and stored for later processing and analysis by high performance computing systems. Although this example shows the top formation of a land-based region, it is understood that this is only an example and that the present methods and systems can also be applied to survey regions of the ocean floor.

[0057] Figure 2 is a schematic diagram showing a cross-sectional view of a seismic survey region 101 in accordance with one embodiment. Figure 1 More specifically, in Figure 2 In, reference numeral 201 represents a cross-sectional view of a portion of the formation of the seismic survey region, and reference numerals 202, 203, and 204 represent different types of formations. Although the seismic survey region in this example is based on land, it is understood that this is only an example and that the present methods and systems can also be applied to survey regions of the ocean floor. Figure 2 A common common midpoint gather is shown in which seismic data is sorted by surface geometry to simulate a single reflection point on the earth. The exploration seismic data can also be referred to as traces, gathers, or image gathers. In Figure 2In the example, data from one or more shot points or blast points and receivers can be combined into a single image gather, or they can be used individually depending on the type of analysis to be performed.

[0058] like Figure 2 As shown, one or more shot points or blast points represent seismic sources at different incident points or stations 104 located on the Earth's surface, where one or more seismic sources are activated. Seismic energy or seismic sources from multiple incident points 104 will be reflected from interfaces between different strata. These reflections will be captured by multiple seismic data recording sensors 105, each sensor 105 having a different position offset 210 relative to each other and relative to well 103. Since all incident points 104 and all seismic data recording sensors 105 are placed at different position offsets 210, exploration seismic data or seismic traces (also referred to in the art as gathers or image gathers) will be recorded at different incident angles 208. Incident points 104 generate downward propagating rays 205 in the Earth, whose upward propagating reflections are captured by seismic data recording sensors 105. In this example, well 103 is an existing drilled well connected to wellbore 209, where multiple measurements can be taken along the wellbore using techniques known in the art. The wellbore 209 is used to acquire logging data, which may include P-wave velocity, S-wave velocity, density, etc. It can be deployed within the exploration area. Figure 2 Other sensors, not shown in the diagram, are used to acquire seismic data. Seismic data can be used to examine the dependence of amplitude, signal-to-noise ratio, time difference, frequency content, phase, and other seismic properties on the incident angle 208°, migration measurement 210°, azimuth, and other geometric properties that are crucial for data processing and imaging in seismic survey areas.

[0059] Figure 3 This is a schematic diagram according to one embodiment, showing a cross-section of a seismic exploration area having a wellbore and logging tools, the latter including one or more acoustic wave generators and one or more logging data recording sensors. An acoustic wave generator is an example of a device that generates one or more acoustic waves. An acoustic wave generator may be referred to as a sound source because it generates or produces one or more acoustic waves, which are also called seismic waves. One or more logging data recording sensors are examples of one or more seismic data recording sensors (seismic receivers or seismic data recorders), and may be the same seismic data recording sensor as seismic data recording sensor 105. In embodiments of the invention, oil and / or gas production may be stopped to generate seismic waves and record seismic data, including reflections of seismic waves as they move underground through one or more strata in the seismic exploration area.

[0060] Figure 3A petroleum drilling system 300 on land 305 is shown, including a drilling rig 310. The drilling rig 310 supports the running of a logging tool 315 into a wellbore 320. The logging tool 315 can include one or more acoustic generators (sound sources) to generate one or more acoustic waves that are transmitted into one or more formations to produce reflections or reflected waves in the one or more formations. While the present example shows one or more formations in a land-based survey area, it is understood that this is merely an example, and that the present methods and systems can also be applied to survey areas at the surface or bottom of a body of water, such as the ocean. The logging tool 315 also includes one or more logging data recording sensors. As described above, the one or more logging data recording sensors receive and record logging data, including reflected data received by the one or more logging data recording sensors in response to acoustic waves transmitted into the one or more formations by the one or more acoustic generators. Logging data is an example of seismic data. Logging data can include compressional wave velocity or P-wave velocity (Vp), shear wave velocity (Vs), and density, which is an indicator of porosity. This logging process of recording logging data can also be referred to as sonic logging. A vehicle 325 can be coupled to the logging tool 315 to assist in running and pulling the logging tool 315, and to communicate with the logging tool 315 to obtain logging data. Alternatively, in methods and systems for survey areas at the surface or bottom of a body of water, such as the ocean, another device or system can be used to assist in running or pulling the logging tool 315, and to communicate with the logging tool 315 to obtain logging data.

[0061] Figure 4 is a schematic diagram showing a high performance computer system according to one embodiment, which receives (often wirelessly) seismic data about seismic waves at seismic data recording sensors 105 in Figure 1 and / or seismic data recording sensors (also referred to as logging data recording sensors in Figure 2 Figure 3 Figure 3 Figure 4 The high performance computer system in Figure 4 ​​​A data transfer system 400 is shown for wirelessly transferring seismic data from seismic data recording sensors to a system computer 405 connected to one or more storage devices 410 for storing the seismic data in a database. The data transfer system can also wirelessly transfer seismic data directly from the seismic data recording sensors 405 to one or more storage devices 410 for storing the seismic data in a database that the system computer 405 can access. The wireless transfer is indicated by reference numeral 402. The one or more storage devices 410 can also store other computer software instructions or programs to implement the apparatus and methods described in the embodiments. The system computer 405 can be coupled (e.g., wirelessly coupled) to one or more output storage devices 420 that can receive the results of the computer-implemented processes or methods performed by the system computer 405. A personal computer 425 can be coupled (e.g., wirelessly coupled) to the one or more output storage devices 420 and / or the computer system 405 so that a user can use the user interface of the personal computer 425 to input information or to obtain the results of the computer-implemented processes performed by the system computer 405. The one or more storage devices 420 can also store other computer software instructions or programs to implement the apparatus and methods described in some embodiments.

[0062] For example, the user interface of the personal computer 425 can include one or more of a keyboard, a mouse, a joystick, buttons, switches, an electronic or stylus pen, a gesture recognition sensor (e.g., to recognize gestures of a user, including movements of body parts), an input seismic device, or a voice recognition sensor (e.g., a microphone to receive voice commands), an output seismic device (e.g., a speaker), a trackball, a remote control, a portable telephone (e.g., a cellular phone or a smart phone), a tablet computer, a pedal or footswitch, a virtual reality device, etc. The user interface can also include a haptic device to provide haptic feedback to the user. For example, the user interface can also include a touch screen. In addition, the personal computer 425 can be a desktop computer, a notebook computer, a tablet computer, a cell phone, or any other personal computing system.

[0063] The processes, functions, methods, and / or computer software instructions or programs in the apparatus and methods described in the embodiments herein can be recorded, stored, or fixed in one or more non-transitory computer-readable media (computer readable storage / media) including a program instructions (computer readable instructions) executable by a computer to make one or more processors execute (implement) the program instructions. The media can also include, alone or in combination with the program instructions, data files, data structures, and the like. The media and program instructions can be specially designed and constructed for the purposes of the present disclosure, or they can be of the kind well known and available to those having skill in the computer software arts. Examples of non-transitory computer-readable media include magnetic media, such as hard disks, floppy disks, and magnetic tape; optical media such as CD ROMs and DVDs; magneto-optical media, such as optical disks; and hardware devices that are specially configured to store and execute program instructions, such as read-only memory (ROM), random access memory (RAM), flash memory, and the like. Examples of program instructions include both machine code, such as produced by a compiler, and files containing a high-level code that can be executed by the computer using an interpreter. The program instructions can be executed by one or more processors. The hardware devices can be configured to act as one or more software modules in order to perform the operations and methods described above, or vice versa. Additionally, a non-transitory computer-readable medium can be distributed over a network, such as the Internet, for example, by way of a wireless communication to effectuate execution of the program instructions by a computer. Moreover, a computer-readable medium can also be embodied in at least one special-purpose machine or apparatus, such as in at least one Application Specific Integrated Circuit (ASIC) or Field-Programmable Gate Array (FPGA).

[0064] The one or more databases can include a collection of data and supporting data structures, which can be stored, for example, in one or more storage devices 410 and 420. For example, the one or more storage devices 410 and 420 can be embodied in one or more non-transitory computer-readable storage media, such as non-volatile storage devices (e.g., read-only memory (ROM), programmable read-only memory (PROM), erasable programmable read-only memory (EPROM), and flash memory), USB drives, volatile storage devices (e.g., random access memory (RAM)), hard disks, floppy disks, Blu-ray discs, or optical media (e.g., CD ROMs and DVDs), or a combination thereof. However, examples of the storage devices 410 and 420 are not limited to the above descriptions, and storage can be implemented by other various devices and structures as understood by those skilled in the art.

[0065] Figure 5 is a flowchart showing a seismic full waveform inversion method employing a Hanning window according to one embodiment of the present disclosure. The survey area can be a subsurface structure under land or a subsurface structure under sea.

[0066] Reference Figure 5the seismic full waveform inversion method, in operation 505, input data (observed seismic data or d obs ). For example, Figure 1 and Figure 2 , the seismic data recording sensors 105 and / or the logging data recording sensors of the logging tools 315 in Figure 3 may detect seismic data and transmit the seismic data to the high performance computing system shown in Figure 4 . As described above with respect to Figures 1 to 4 , the seismic data detected in the survey area can be stored in one or more memories, such as the one or more storage devices 410 and the one or more output storage devices 420.

[0067] Further, as shown in Figure 5 , an initial earth model m0may be input to the seismic full waveform inversion method. The initial earth model m0may be a P-wave velocity model.

[0068] It can be appreciated that the initial velocity model can be a set of parameters. The initial velocity model m0may be a predetermined velocity model based on a given scientific specification. The initial velocity model m0may be input by a user through a user interface, such as a user interface of the personal computer 425, or stored in one or more memories, such as the one or more storage devices 410 and the one or more output storage devices 420.

[0069] As shown in the seismic full waveform inversion method of Figure 5 , the initial model (such as the initial velocity model m0) is input to a loop. In operation 500, a current velocity model m i is determined. Initially, the velocity model m i is the initial velocity model m0. The variable i represents the number of iterations of the loop. Thus, initially i = 0. In each iteration in Figure 5 , the velocity model m i is updated in operation 500. For example, after the first iteration, the velocity model is m1, where i = 1. After the second iteration, the velocity model is m2, where i = 2. Iterations will continue until convergence is detected in operation 570. After convergence is detected in operation 570, the final velocity model m i+1 corresponding to the latest velocity model in operation 510 is output in operation 575 to a storage device or a display.

[0070] Referring to operation 510 in the seismic full waveform inversion method shown in Figure 5 , in operation 510, a velocity model m ito perform forward modeling. Forward modeling of seismic data is a technique to create (generate) synthetic seismic data from geologic information, in this case, the velocity model m i . The forward modeling operation requires a source wavelet, which can be a known source wavelet function, different from the true source characteristics of the seismic survey being performed or analyzed. Examples of source wavelet functions include Ormsby wavelet, Ricker wavelet, etc. The source wavelet is used in the forward modeling operation to generate synthetic data d i from the velocity model m syn . The synthetic data d syn is represented by reference numeral 520, output from the forward modeling operation 510.

[0071] In Figure 5 operation 525, the synthetic data d syn is input to a Hanning window, while in operation 515, the real data d obs is Hanning windowed. The Hanning function is a windowing function that separates a full / long signal into a series of local / short signals along the record time. Figure 7 A sliding Hanning window is shown, whose size and sliding step are automatically determined according to the data frequency, for automatically localizing the seismic data. The sliding Hanning window can effectively convert the data into local domain without energy loss.

[0072] In operation 515, equation (5) represents the kth localized data d k (t) processed using the Hanning function W obs,s,r (t):

[0073] d obs,s,r,k (t) = W k (t) d obs,s,r (t), (5)

[0074] And equation (6) represents the summation of the localized data:

[0075]

[0076] It is noted that the summation of the localized data is equal to the original data. Therefore, the data localization process using Hanning windowing does not cause energy loss.

[0077] Equation (5) can be expanded as

[0078]

[0079] which is frequency dependent. T(f) represents a time wavelength at frequency f; the local window step t k is one wavelength; and the local window size tw It consists of two wavelengths.

[0080] The same Hann windowing process is also applicable to manipulating synthetic data in 525, producing localized synthetic data d. obs,s,r,k (t). Therefore, the total number of windows N for each seismic trace k The calculation formula is:

[0081]

[0082] Among them, t max This is the maximum recording time starting from 0. The number of windows increases with the center frequency.

[0083] In operation 530, the time shift of the k-th local dataset can be obtained by automatically selecting the maximum value of the cross-correlation between the synthetic data and the observed data:

[0084]

[0085] It can capture the unsteady changes in the travel time of different seismic events.

[0086] In operation 535, the localized observation data d obtained through ΔT(s,r,k) in operation 530 and the Hann windowing in operation 515 are used. obs,s,r (t) yields the associated earthquake source, as shown below:

[0087]

[0088] Refer to operation 550 for the FWI gradient g i The calculation involves backpropagating one or more associated sources generated by formula (10) in operation 550 to obtain the associated wave field, the calculation formula of which is as follows:

[0089]

[0090] in, It is the adjoint operator F(m;x), d of the positive model operator. adj,s (t) is the associated seismic source, while u s (x,t) is the adjoint wave field.

[0091] Using the forward wavefield (operation 540) and the adjoint wavefield as described above, the gradient for each iteration can be obtained using the following formula:

[0092]

[0093] Where F(m;x) is the positive operator of the wave equation, w s (x,t) represents the forward wave field, while u s(x, t) is the adjoint wavefield of equation (11). The gradients from different shot points (explosions or acoustic generators) as described by equation (11) are added and stacked in operation 550 to obtain a single gradient.

[0094] Operation 555 determines: (1) the magnitude of the increase in the estimated velocity toward the true velocity model, or the magnitude of the decrease in the estimated velocity model; and (2) whether the estimated value of the velocity model should be increased to enhance the modeling of the true velocity model, or whether the estimated value of the velocity model should be decreased to enhance the modeling of the true velocity model.

[0095] More specifically, after the gradient is computed and output by operation 550, an optimization method is used in operation 555 to compute the step size a i and the search direction P i . For example, starting from an initial guess of the subsurface parameters m i＝0 , the model update at iteration i+1 is

[0096] m i+1 = m k + a i P i i = 0, 1,..., (13)

[0097] where P i is the search direction and a i is the step size. The simplest inversion method is called the steepest descent algorithm, where the search direction P i is simply given by the negative gradient of the objective at iteration k:

[0098] P i = -g i , (14)

[0099] Thus, equation (13) becomes

[0100] m i+1 = m i - a i g i (15)

[0101] The step size is then determined in a way that minimizes the cost function.

[0102] In operation 555, the gradient is received and the search direction P i (for the increase or decrease of the estimated velocity model) and the step size a i (for the magnitude of the increase or decrease of the estimated velocity toward the true velocity model) are determined. Thus, in operation 555, the FWI gradient g ithe estimated velocity model to avoid an excessively large increase / decrease in the estimated velocity model (which would greatly increase the number of iterations) or an excessively small increase / decrease in the estimated velocity model (which would greatly increase the number of iterations). Thereafter, in operation 560, the velocity model m i+1 .

[0103] Referring to the convergence operation 570, the seismic full waveform inversion method determines whether additional iterations are needed. If additional iterations are needed, the value of i is increased in operation 565, and the updated velocity model m i+1 becomes the velocity model m i .

[0104] In one embodiment, a maximum number of iterations can be set in operation 565. The maximum number of iterations can range from 10 to 40. Once the number of iterations, denoted by i, is equal to the predetermined maximum number of iterations, it is considered that Figure 5 the seismic full waveform inversion method has converged, and the final velocity model m F is outputted, which is the same as the last updated velocity model m i+1 determined in operation 575.

[0105] In another embodiment, an additional threshold value can also be set in operation 575 according to the difference between the current updated velocity output in operation 575 and the previous velocity output in operation 575. If the difference is greater than (or greater than or equal to) a threshold value, the next iteration 575 can continue, provided that the maximum number of iterations has not been exceeded. However, if the difference is less than (or less than or equal to) the threshold value, it is considered that Figure 5 the seismic full waveform inversion method has converged. The image of the final velocity model m F may be displayed on the personal computer system 425.

[0106] Figure 6 is a workflow of another embodiment of the present disclosure, which describes a source-independent wavelet FWI method that employs a matched filter to correct the synthetic data before cross-correlation calculation and traveltime extraction. Figure 6 Operations 600, 605, 610, 615, and 620 in Figure 5 correspond to and perform the same functions as operations 500, 505, 510, 515, and 520 in Figure 5 However, unlike the method shown in s , operation 680 obtains a matched filter M s (t), which has a mathematical expression as follows:

[0107]

[0108] where dobs,s,r (ω) and d syn,s,r (ω) are the Fourier transforms of the time-domain observed data and the synthetic data, respectively.

[0109] Referring to the matched filter shown in equation (16), the matched filter M s (t) is convolved with the observed seismic data d obs obtained in operation 620. The matched filter M syn is matched to the observed seismic data d obs may be data recorded by sensors, such as the seismic data recording sensors 105 within the survey area, and the logging data recording sensors of the logging tool 315 disposed within the wellbore of the survey area.

[0110] The matched filter M s (t) can be a Wiener filter. The matched filter matches the phase and amplitude of the synthetic data to the phase and amplitude of the seismic data. The matched filter is determined in each iteration operation 665 using the computed synthetic data d i (t) to update the velocity model m syn Convolution of the matched filter M

[0111] As mentioned above, errors in the source wavelet can cause errors in the traveltime picks, which in turn can cause errors in the inverted model. Therefore, a reasonable source wavelet estimation is still very important for the success of the cross-correlation based FWI. A typical procedure for such estimation is to window the first arrival event in the data and stack the event as the source wavelet. This estimated source wavelet does not change during the iterations. However, it is very difficult to achieve an accurate estimation of the source wavelet in industrial field applications because of the poor repeatability of the source signature for different shots (different shots or different sources), the uncertainty in the coupling of the source to the earth, and the uncertainty in the coupling of the receivers to the earth. Therefore, to overcome these problems, a lot of effort has been spent on shot-independent misfit function methods. The operation 680 of the matched filter in combination with other operations, such as Figure 6 obtained in operation 635, can solve these problems.

[0112] The synthetic data from operation 620 is first convolved with the matched filter M s (t) obtained in operation 680. The convolved synthetic data M s (t)*d syn,s,r,k(t) is more similar to the observed data, then it is localized with a Hanning window in operation 625 to get d obs,s,r,k (t).

[0113] In operation 630, the automatic traveltime difference pick in the local window is defined as:

[0114]

[0115] where the match filter M s (t) is calculated by matching the synthetic data to the observed seismic data.

[0116] Therefore, the traveltime-based misfit function is defined as:

[0117]

[0118] In operation 635, the adjoint source can be derived using the chain rule as:

[0119]

[0120] where is the cross-correlation operator. The Wiener match filter M s (t) is determined shot by shot in each iteration to correct the phase of the adjoint source. The final adjoint source is the weighted stack result of the local window.

[0121] Figure 6 Operations 640, 645, 650, 655, 660, 665, 670, and 675 in Figure 5 correspond to and are substantially similar to operations 540, 545, 550, 555, 560, 565, 570, and 575 in

[0122] Therefore, in the method shown in Figure 6 the gradient of the cost function can be constructed using the adjoint source from operation 635. For example, using the forward wavefield and the adjoint wavefield, the FWI gradient g i can be obtained by solving the adjoint wave equation (11) in operation 650, then it is used to update the velocity model m using any inversion method (such as the steepest gradient algorithm) to minimize the misfit function equation (18). There are other more advanced inversion methods, such as the nonlinear conjugate gradient method, the Gauss-Newton method, and the quasi-Newton method L-BFGS. Operation 550 and 650 are substantially similar.

[0123] However, comparing the adjoint source obtained in operation 635 according to equation (19) with the adjoint source obtained in operation 535 according to equation (10) can clearly show that there is one extra cross-correlation operation in operation 635, i.e. the cross-correlation operation of M s (t) in equation (19).

[0124] Therefore, the cross-correlation based FWI method becomes independent of the source wavelet. The matched filter operation 680 and the adjoint source operation 635 do not need to update the source wavelet in the forward modeling operation 610.

[0125] The inversion model parameters generated by this method are independent of the choice of the wavelet. This wavelet independence reduces the model inversion error caused by the source wavelet error and reduces the workload of manually estimating the source wavelet from shot to shot. This makes the cross-correlation based FWI very robust and can be used to build a high-resolution geological model to improve the image of the complex subsurface structure in the exploration area, thereby improving lithology identification, fluid discrimination and reservoir characterization in the field of seismic exploration.

[0126] Figure 5 or Figure 6 The adaptive cross-correlation based FWI shown can robustly invert medium parameters. According to an embodiment, each iteration includes a series of steps, for example: Step 1, use a theoretical wavelet (such as the Ormsby wavelet) and the medium parameter model m i in the current i-th iteration to perform forward modeling using the wave equation (1) to obtain the forward wave field w s (m i ; x, t) of the subsurface and estimate the synthetic data d syn,s,r (m i ; t) recorded at the receiver locations; Step 2, determine the matched filter, e.g. the Wiener filter, using equation (16); Step 3, automatically select the travel time using equation (10) and automatically determine the local window using the sliding Hanning window described in equation (17); Step 4, use equation (19) to calculate the adjoint source; Step 5: use the same model to back-propagate the adjoint source to obtain the adjoint wave field u s (m i ; x, t); Step 6, use the obtained forward wave field w s (m i ; x, t) and the adjoint wave field u s (m i ; x, t) to calculate the gradient g i (x); Step 7, use the selected optimization method to obtain the search direction p i (x); Step 8, use linear perturbation theory to calculate the step size a iStep 9: Update the medium parameter model with the step size. Repeat steps 1 to 9 until the unfit function, i.e., equation (18), is minimized to a predetermined value.

[0127] Figures 8 to 17 This demonstrates that the method disclosed in this paper can invert high-resolution subsurface velocity details while avoiding loop jumps. Figure 8 The actual velocity model with a grid space of 50 meters is shown, which is used to generate seismic data with a maximum offset of 25 kilometers. Figure 9 The initial velocity model used for FWI, namely the Marmousi model, is shown. In this inversion, the medium parameter m has only one parameter, namely the propagation velocity V. P The data used by this FWI is from Figure 9 The real Masmousi model and Ormsby wavelet generation in the model. Figure 8 and Figure 9 The speed difference is within the same range of 1500 m / s to 4500 m / s. Figure 10 Showing from Figure 8 Subtracting the actual velocity model in Figure 9 The velocity difference obtained after the initial velocity model in the model. Figure 10 The speed difference ranges from -1365 m / s to 1980 m / s.

[0128] Figure 11 Showing the use Figure 6 The inversion velocity model of the adaptive cross-correlation FWI method of the present invention has bandpass filtered 2.5Hz data, assuming that the effective lowest frequency is 2.5Hz and that there is no lower frequency data in the observation data. Figure 12 Showing the use Figure 6 The final inversion velocity model of the adaptive cross-correlation FWI method of the present invention has data at 2.5 Hz, 4 Hz, 6 Hz and 10 Hz after continuous bandpass filtering.

[0129] Figure 13 The final inversion velocity model is shown using the FWI method based on the traditional L2 waveform difference according to formula (4). Figure 13 It has the same data at 2.5Hz, 4Hz, 6Hz, and 10Hz after continuous bandpass filtering. Figure 13 and Figure 12 By comparison, it is easy to see the inversion artifacts in the region enclosed by the rectangle, as well as the incorrect update direction caused by cyclic jumps in the region enclosed by the ellipse.

[0130] Figure 14 An inversion velocity model derived from a real land dataset is shown using ray-based tomography. Figure 15The inversion model perturbation (difference between the inversion model by FWI and the velocity model inverted by tomography) is shown using the adaptive cross-correlation based FWI method of the present invention as shown in Figure 6

[0131] Figure 16 The seismic migrated image obtained using the ray-based initial tomography model is shown while Figure 17 The seismic migrated image obtained using the adaptive cross-correlation FWI method of the present invention as shown in Figure 6 Figure 17 The geological features are more smooth and clear in

[0132] Thus, although Figure 16 and Figure 17 An example of land data is shown with unknown true source wavelets recorded from different shot points and forward modeled using Ormsby wavelets, but Figure 6 The method of the present invention as described successfully updates the velocity model and hence improves the seismic image and results in better geological interpretation, which shows the robustness of the method of the present invention, especially in the absence of knowledge of the true wavelets for different shots of the land data.

[0133] While the embodiments disclosed herein describe a FWI algorithm with a traveltime (differences) misfit function that utilizes a finite-difference numerical solution of the scalar wave equation for seismic propagation, it will be apparent to those of ordinary skill in the art that the algorithm can alternatively be applied to the vector wave equation and elastic equations in isotropic and anisotropic media without departing from the true scope of the invention as defined by the claims below.​​

Claims

1. A method of performing seismic full waveform inversion to generate a velocity model of a subsurface formation of a survey area, comprising the operations of: (a) arranging seismic data recording sensors at a plurality of locations within the survey area and arranging a logging tool having seismic data recording sensors in a wellbore of the survey area; (b) performing a shot at an entry point of the survey area to generate a seismic wave that travels through the subsurface formation; (c) observing the seismic wave using the seismic data recording sensors and recording seismic data from the seismic wave; (d) transmitting the observed seismic data from the seismic data recording sensors to a computer system comprising one or more memories and storing the observed seismic data in the one or more memories, storing a source wavelet and a current medium parameter model in the one or more memories; (e) performing, by the computer system, a forward modeling operation using the source wavelet and the current medium parameter model to obtain a forward wavefield and generate synthetic data from the forward wave equation; (f) processing, by the computer system, the synthetic data and the observed seismic data using a Hanning window operation to obtain localized synthetic data and localized seismic data, respectively; (g) selecting, by the computer system, a traveltime difference from a cross-correlation between the localized synthetic data and the localized seismic data; (h) solving, by the computer system, an adjoint equation of the forward wave equation using the traveltime difference and the localized seismic data to obtain an adjoint source; (i) generating, by the computer system, an updated medium parameter model for seismic full waveform inversion and synthetic data using the forward modeling operation; (j) performing operations (e) through (i) until convergence is reached; (k) outputting the updated velocity as a final velocity model to a display upon convergence; and (l) displaying an image of the final velocity model on a display of the computer system. Further comprising constructing, by the computer system, a matched filter that matches the synthetic data to the observed seismic data.

2. The method of claim 1, wherein, Convoluting the synthetic data with the matched filter prior to performing the Hanning window operation.

3. The method of claim 2, wherein, The operation (g) further comprises, by the computer system:

4. The method of claim 3, wherein, backpropagating the adjoint source to obtain an adjoint wavefield; calculating a gradient using the forward wavefield and the adjoint wavefield; determining a search direction and a step size; and updating the medium parameter model according to the step size. Convergence is determined when a value of a traveltime-based misfit function is less than a predetermined value.

5. The method of claim 1, wherein, Convergence is determined when a number of iterations reaches a predetermined value.

6. The method of claim 1, wherein, The medium parameter model is selected from a velocity model.

7. The method of claim 1, wherein, The source wavelet is an Ormsby wavelet or a Ricker wavelet.

8. The method of claim 1, wherein, The matched filter is a Wiener filter.

9. The method of claim 2, wherein, The inversion method used in the updating step is selected from a steepest gradient algorithm, a nonlinear conjugate gradient method, a Gauss-Newton method, and a quasi-Newton method.

10. The method of claim 4, wherein, 11. A system for performing seismic full waveform inversion to generate a velocity model of a subsurface formation of a survey area, the system comprising: ​ a plurality of seismic data recording sensors disposed at different locations within the survey region, and / or a logging tool including a seismic data recording sensor disposed in a wellbore within the survey region; a shot device arranged at each shot point of the survey region for generating seismic waves, wherein the seismic waves travel through the subsurface formation; and a plurality of seismic data recording sensors for sensing the seismic waves and recording seismic data from the seismic waves; wherein the seismic data recording sensors transmit the seismic data to a computer system including one or more memories storing the transmitted seismic data, shot wavelets, and instructions, and at least one processor executing the instructions stored in the one or more memories to implement: performing a forward modeling operation according to a forward wave equation and using a shot wavelet and a current medium parameter model to obtain a forward wavefield; performing a Han window operation on the synthetic data and the observed seismic data, respectively, to obtain localized synthetic data and localized seismic data; selecting a traveltime difference according to a cross-correlation between the localized synthetic data and the localized seismic data; solving an adjoint equation of the forward wave equation using the traveltime difference and the localized seismic data to obtain an adjoint shot; and using the forward modeling operation to generate an updated medium parameter model for seismic full waveform inversion and synthetic data.

Images