Full waveform inversion method and device based on adaptive cross-correlation

Abstract

Methods and apparatus are provided for performing travel time based seismic full waveform inversion using frequency dependent Hann windows and source independent schemes to generate a final velocity model of a subterranean formation of a surveyed region. The method comprises the following steps: arranging a seismic data recording sensor in a survey area; blasting at an incidence point of the surveying area to generate seismic waves; and sensing and recording seismic waves using the seismic data recording sensor. The recorded seismic waves are seismic data. The method further includes transmitting the seismic data to a computer system including one or more memories and storing the seismic data in the one or more memories; storing the seismic source wavelets in one or more memories; performing, by the computer system, a forward modeling operation using the seismic source wavelet; and generating, by the computer system, a final velocity model using seismic full waveform inversion and forward modeling operations.

Description

Technical Field

[0001] The present disclosure relates to a seismic exploration method and apparatus for constructing a high-resolution geological model by performing full waveform inversion (FWI) based on adaptive cross-correlation to enhance the characterization of complex underground structures in an exploration area. The method first automatically evaluates the travel time difference between modeled synthetic data and field observation data by adaptive cross-correlation calculation based on frequency-dependent sliding Hann window scanning and with the help of a source-independent scheme. The evaluated travel time difference is inverted to generate a final velocity model to improve the image of the complex underground structure in the survey area. Background Art

[0002] 1. Overview

[0003] In oil and gas exploration, high-resolution subsurface volume attributes such as seismic propagation velocity, anisotropy, absorptivity, porosity, and reflectivity models can be generated by acquiring and processing seismic data. Combined with these geophysical properties, the subsurface structure can be effectively revealed. Seismic data processing usually includes seismic inversion, which is used to establish a propagation velocity model for medium and long wave numbers, followed by seismic migration to obtain seismic reflectivity images for short wave numbers. These seismic reflectivity images are used to determine the storage location and scale of natural resources (such as hydrocarbons of oil and gas), thereby providing a basis for exploration / drilling plans. In order to obtain high-resolution seismic reflectivity images, a high-magnification 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 may be sufficient for relatively simple geological targets (such as shallow sedimentary environments). However, for complex geological environments such as salt domes, subbasalt targets, thrust belts, and terrestrial foothills, ray-based tomography is not very effective, so FWI becomes a necessary tool for establishing velocity models.

[0005] FWI directly solves the seismic wave equation and matches it to seismic data, which can generate more accurate seismic propagation velocity models for complex underground structures (such as salt-related models). This propagation velocity model can be used to generate accurate, high-resolution seismic reflectivity images through seismic migration, facilitate time-lapse monitoring of hydrocarbon reservoirs, and can even be directly converted to generate high-resolution seismic image volumes, which are called FWI images.

[0006] 2. Earthquake wave equation

[0007] Seismic waves propagating beneath the Earth's surface can be simulated by solving the seismic wave equations. The seismic wave equations describe the Earth using different physical models and assume that the Earth is isotropic or anisotropic, elastic or acoustic, and attenuated or unattenuated. In most FWI developments, seismic waves are assumed to be purely acoustic waves because the acoustic wave equations are relatively simple and can be solved efficiently. Regardless of the assumptions, all wave equations can be expressed mathematically 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 velocity, density, and other physical properties; x is the spatial position; t is the time; F(m; x) is the corresponding (forward) modeling operator; and w s (x, t) is the source wavelet f excited at position s S (t)’s forward wave field.

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

[0011] 3. Full waveform inversion (FWI)

[0012] FWI is a data-driven tool that automatically establishes subsurface parameters m, such as velocity and / or density, by iteratively minimizing the difference between recorded data and modeled synthetic data. Given an initial estimate of the subsurface velocity m0, the subsurface velocity can be estimated by using the source wavelet f S (t) Solve the forward seismic wave equation (1) to predict synthetic data. Take the residual between the data and the synthetic data as the source, solve the adjoint equation of the forward wave equation (1), obtain the solution of the adjoint wave equation, and obtain the gradient by cross-correlating the solutions of the forward and adjoint wave equations. Then use the gradient to update the model along a certain direction to reduce the misfit between the modeled synthetic data and the field observation data. Repeat the iterative solution until the data misfit is small enough.

[0013] Mathematically, FWI based on least squares (L2) waveform differences uses the following unfitting:

[0014]

[0015] Where C(m) is the unfitting function, which is used to measure the waveform difference between the synthetic data and the recorded data, and the analog operation is represented by the symbol |·|. In the unfitting function, N s is the number of earthquake sources, N r is the number of receivers for a given source, t max is the longest recorded 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 for the same source s at the same receiver r 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 Solve the wave s (t) is resampled 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. FWI based on travel time

[0017] Travel time based FWI is based on selecting the travel time difference between the modeled synthetic data and the field observation data. The travel time difference can be obtained by automatically selecting the maximum cross-correlation value between the synthetic data and the recorded data, that is,

[0018]

[0019] On this basis, the travel-time-based FWI can be expressed as the least squares (L2) minimum of the travel-time differences as follows:

[0020]

[0021] Compared to the least squares cost function of waveform differences 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 mitigate the cycle skipping problem without requiring a maximum half-wavelength distance between the model synthetic events and the true observed data events.

[0022] In addition, the amplitude of the modeled synthetic data can easily be inconsistent with the amplitude of the recorded data. This is because seismic forward modeling cannot take into account all the physical properties of the subsurface and may ignore elastic and / or attenuation characteristics. For efficient calculation, acoustic modeling engines are usually used to simulate the propagation of seismic waves. In addition, seismic amplitudes are easily contaminated during seismic processing. The unfit function based on travel time (difference) can dilute the amplitude difference, thereby reducing the negative impact of the amplitude difference in the seismic velocity inversion process.

[0023] However, travel-time-based FWI is affected by inaccurate and low-resolution travel-time (difference) selection. First, the travel-time selection may be incorrect due to crosstalk between multiple seismic events that are not correctly matched. Second, the travel-time selection may have a low resolution because the calculation in equation (3) is usually dominated by strong events, while relatively weak events are not taken into account.

[0024] In addition, the source wavelet f s The accuracy of (t) will also reduce the accuracy of automatic travel time selection based on cross-correlation. Especially in terrestrial applications, the source wavelet may vary between different shot points, so it is difficult to obtain an accurate source wavelet.

[0025] Therefore, new methods are still needed to overcome the shortcomings of traditional travel-time-based FWI. Summary of the invention

[0026] In one embodiment of the present disclosure, a method for performing seismic full waveform inversion to generate a velocity model of underground strata in a survey area includes multiple steps: positioning seismic data recording sensors at different locations in the survey area, and / or positioning a logging tool including a seismic data recording sensor in a wellbore in the survey area; transmitting at an incident point in the survey area to generate seismic waves, which propagate through underground strata; using a seismic data recording sensor to observe the seismic waves, and recording seismic data based on the seismic waves; transmitting the seismic data observed by the seismic data recording sensor to a computer system including one or more memories, 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; and generating a seismic wave according to the forward wave direction by the computer system. The method comprises the steps of: performing a forward modeling operation using a source wavelet and a current medium parameter model to obtain a forward wave field; performing a Hanno window operation on the synthetic data and the observed seismic data respectively to obtain local synthetic data and local seismic data; selecting a travel time difference according to the cross-correlation between the local synthetic data and the local seismic data; solving an adjoint equation of the forward wave equation using the travel time difference and the local seismic data to obtain an adjoint source; generating 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 the updated velocity as a final velocity model to a display after convergence; and displaying a high-resolution image of the final velocity model on a display of the computer system.

[0027] Other embodiments provide FWI methods based on adaptive cross-correlation. In one embodiment, the unfitted measurements are automatically selected based on the adaptive cross-correlation calculation of the local seismic window. The adaptive cross-correlation indicates that there is no need to iterate to determine the optimal local window size and estimate the exact source wavelet for the seismic modeling data. The seismic data is automatically localized by a frequency-dependent Hann sliding window, and the possible source wavelet f in the forward seismic wave equation (1) is s (t) The error is corrected by the Wiener filter.

[0028] Also disclosed is a system for performing seismic full waveform inversion to generate a velocity model of underground formations in a survey area. The system includes: a plurality of seismic data recording sensors disposed at different locations in the survey area, and / or a logging tool including seismic data recording sensors disposed in a wellbore in the survey area; an blasting device disposed at each incident point in the survey area to generate seismic waves that will propagate through underground formations; and a plurality of seismic data recording sensors for sensing the seismic waves and recording seismic data based on the seismic waves. The seismic data recording sensors transmit seismic data to a computer system, the system including one or more memories and at least one processor, the one or more memories storing the transmitted seismic data, source wavelets, and instructions, the one or more processors executing instructions stored in the one or more memories to implement one of the methods disclosed in the present disclosure.

[0029] Although the embodiments disclosed herein describe an FWI algorithm with a travel time (difference) misfit function that utilizes a finite difference numerical solution of the scalar acoustic wave equation for seismic propagation, it is apparent to one of ordinary skill in the art that the algorithm may alternatively be applied to vector wave equations and elastic equations in isotropic and anisotropic media without departing from the true scope of the invention as defined in the following claims. BRIEF DESCRIPTION OF THE DRAWINGS

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

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

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

[0033] Figure 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 is a schematic diagram showing a high performance computing system according to one embodiment.

[0035] Figure 5 is a flow chart showing an FWI process using local Hann window based cross-correlation FWI according to the present disclosure.

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

[0037] Figure 7 A sliding Hann window is shown, whose size and sliding step are automatically determined based on the data frequency of the localized seismic data.

[0038] Figure 8 The true velocity model with a 50-m grid spacing is shown, used to generate seismic data with a maximum offset of 25 km.

[0039] Fig. 9 The initial velocity model used for the FWI test is shown.

[0040] Fig.10 Shows from Figure 8 The true velocity model shown is subtracted Fig. 9 The velocity differences obtained after the initial velocity model are shown, ranging from approximately -1365 m / s to 1980 m / s.

[0041] Fig.11 The inverted velocity model with bandpass filtered 2.5 Hz data using the adaptive cross-correlation FWI method of the present invention is shown, where the effective lowest frequency is assumed to be 2.5 Hz, while there are no lower frequencies in the observed data.

[0042] Fig.12 The final inverted velocity model with successively bandpass filtered 2.5 Hz, 4 Hz, 6 Hz and 10 Hz data using the adaptive cross-correlation FWI method of the present invention is shown.

[0043] Fig.13 The final inverted velocity model with successively bandpass filtered 2.5 Hz, 4 Hz, 6 Hz and 10 Hz data using the conventional L2 waveform difference based FWI method is shown.

[0044] Fig.14 The velocity model of a real land dataset generated by ray tomography is shown, which will serve as the initial velocity model for FWI.

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

[0046] Fig.16 Earthquake migration images obtained using the initial ray-based tomography model are shown.

[0047] Fig.17A seismic migration image obtained using the inverted velocity model by the adaptive cross-correlation FWI method of the present invention is shown.

[0048] The meanings of the reference numerals in the drawings are as follows:

[0049]

[0050] DETAILED DESCRIPTION

[0051] Several embodiments of the present disclosure will now be described in detail, examples of which are shown in the accompanying drawings. It should be noted that, where applicable, similar or similar reference numerals in the accompanying drawings represent similar or similar functions. The embodiments of the present disclosure are shown in the accompanying drawings only in an illustrative manner. Those skilled in the art can easily recognize from the following description that alternative embodiments of the structures, systems, and methods described herein can be adopted without departing from the principles disclosed herein.

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

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

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

[0055] Figures 5 to 17 Exemplary embodiments of apparatus, methods, and media are shown that improve lithology identification, fluid discrimination, and reservoir characterization in the field of seismic exploration by using improved seismic full waveform inversion techniques to improve the quality of seismic full waveform inversion results, including the use of a computer-implemented seismic full waveform inversion method. Figures 5 to 17 Exemplary embodiments of apparatus, methods, and media for generating one or more high-resolution geological models for use in high-resolution images for lithology identification, fluid discrimination, and reservoir characterization of complex subsurface structures in a survey area are shown. Figures 5 to 17An exemplary embodiment of generating inversion model parameters utilizing a frequency-dependent sliding Hann window and independent of source wavelets for performing seismic full waveform inversion to generate a high-resolution geological model for a high-resolution image of a survey area including complex subsurface structures is shown. Figures 5 to 17 Exemplary embodiments are shown that utilize a frequency-dependent sliding Hann window and are independent of the form of the source wavelet, which can generate one or more high-resolution geological models for high-resolution imaging to perform lithology identification, fluid discrimination, and reservoir characterization of the complex subsurface structure of the survey area.

[0056] Figure 1 is a schematic diagram showing a top view of a survey area including multiple incident points of a seismic source according to one embodiment. More specifically, Figure 1 A seismic survey area (survey area) 101 is shown, which is a land-based area, and reference numeral 102 represents the top strata of the land-based area. Those skilled in the art will recognize that a seismic survey area can generate a detailed image of the local geology to determine the location and size of possible hydrocarbon (oil and gas) reservoirs, thereby determining the well location 103. In these survey areas, seismic waves bounce back in the underground rock formations when they are emitted from one or more seismic sources at different injection points 104. Explosions are an example of seismic sources generated by seismic equipment. The seismic waves reflected back to the surface are captured by seismic data recording sensors 105, transmitted from the seismic data recording sensors 105 by one or more data transmission systems (usually wireless transmission), and stored for later processing and analysis by a high-performance computing system. Although this example shows the top strata of a land area, it is understood that this is only an example and the method and system can also be applied to survey areas at the bottom of the ocean.

[0057] Figure 2 is a schematic diagram according to an embodiment showing Figure 1 A cross-sectional view of a seismic survey area 101 is shown, showing the incident point of the seismic source, the seismic data recording sensor (seismic receiver), the well location, the wellbore, various transmission rays, and various incident angles. More specifically, in Figure 2 In the figure, reference numeral 201 represents a cross-sectional view of a part of the stratum on the seismic survey area, and reference numerals 202, 203 and 204 represent different types of strata. Although the seismic survey area in this example is based on land, it can be understood that this is just an example, and the method and system can also be applied to the survey area on the bottom of the ocean. Figure 2 shows a common midpoint gather, in which seismic data are sorted by surface geometry to simulate a single reflection point on the Earth. Exploration seismic data may also be called traces, gathers, or image gathers. Figure 2In the example of , data from one or more shots or blasts 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 blasting points represent seismic sources of different entry points or stations 104 located at the surface of the earth, at which one or more seismic sources are activated. Seismic energy or seismic sources from multiple entry points 104 will be reflected from the interface 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 the well 103. Since all entry 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 as gathers or image gathers in the art) will be recorded at different incident angles 208. The entry points 104 generate downward transmission rays 205 in the earth, and their upward transmission reflections are captured by the seismic data recording sensors 105. In this example, the well location 103 is an existing wellbore, which is connected to a wellbore 209, and multiple measurements can be made along the wellbore using techniques known in the art. The wellbore 209 is used to obtain logging data, which may include P wave velocity, S wave velocity, density, etc. Figure 2 The seismic data can be used to examine the dependence of amplitude, signal-to-noise ratio, time difference, frequency content, phase and other seismic attributes on the incident angle 208, offset measurement 210, azimuth and other geometric attributes, which are very important for data processing and imaging of the seismic survey area.

[0059] Figure 3 1 is a schematic diagram according to one embodiment showing a cross-section of a seismic exploration area having a wellbore and a logging tool, 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 the acoustic wave generator generates or produces one or more acoustic waves, which are also referred to as 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 sensors as seismic data recording sensor 105. In an embodiment of the present invention, oil and / or gas production may be stopped to generate seismic waves and record seismic data, which includes reflections of seismic waves as they move underground in one or more formations in the seismic exploration area.

[0060] Figure 3An oil drilling system 300 on land 305 is shown, including a drilling rig 310. The drilling rig 310 supports the lowering of a logging tool 315 into a wellbore 320. The logging tool 315 may include one or more acoustic wave generators (sound sources) to generate one or more acoustic waves that are transmitted into one or more formations to generate reflections or reflected waves in the one or more formations. Although this example shows one or more formations in a land-based survey area, it is understood that this is only an example, and the method and system may also be applied to survey areas at the surface or bottom of a body of water such as an 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 reflection data received by one or more logging data recording sensors in response to acoustic waves emitted by one or more acoustic wave generators into one or more formations. Logging data is an example of seismic data. The logging data may include compression wave velocity or P wave velocity (Vp), shear wave velocity (Vs), and density, the latter being an indicator of porosity. This logging process of recording logging data may also be referred to as acoustic logging. The vehicle 325 may be connected to the logging tool 315 to assist in lowering and lifting the logging tool 315, and communicate with the logging tool 315 to obtain logging data. Alternatively, in the method and system for surveying an area at the surface or bottom of a body of water such as an ocean, another device or system may be used to assist in lowering or lifting the logging tool 315, and 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 an embodiment of the present invention. Figure 1 and Figure 2 The seismic data recording sensors 105 and / or Figure 3 Seismic data recording sensors in Figure 3 A well data recording sensor (also called a well logging sensor) receives (often wirelessly) seismic data about seismic waves. Figure 4 The high performance computer system in the present invention stores the seismic data in at least one memory for subsequent processing and analysis by the computer-implemented method and apparatus of one or more embodiments. The analyzed or processed seismic data can be accessed by the personal computer system. More specifically, Figure 4A data transmission system 400 is shown for wirelessly transmitting 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 transmission system can also wirelessly transmit seismic data from the seismic data recording sensors 405 directly to one or more storage devices 410 for storing the seismic data in a database that can be accessed by the system computer 405. Wireless transmission is represented 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 this embodiment. The system computer 405 can be coupled (e.g., wirelessly coupled) to one or more output storage devices 420, which can receive the results of the computer-implemented process or method performed by the system computer 405. A personal computer 425 can be coupled (e.g., wirelessly coupled) to 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 obtain the results of the computer-implemented processing method performed by the system computer 405. The one or more storage devices 420 may 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 may include one or more of a keyboard, a mouse, a joystick, a button, a switch, an electronic pen or stylus, a gesture recognition sensor (e.g., for recognizing a user's gestures, including movements of body parts), an input seismic device or a voice recognition sensor (e.g., a microphone for receiving voice commands), an output seismic device (e.g., a speaker), a trackball, a remote control, a portable phone (e.g., a cellular phone or a smart phone), a tablet computer, a pedal or foot switch, a virtual reality device, etc. The user interface may also include a tactile device for providing tactile feedback to the user. For example, the user interface may also include a touch screen. In addition, the personal computer 425 may be a desktop, a laptop, a tablet computer, a mobile phone, or any other personal computing system.

[0063] The processes, functions, methods and / or computer software instructions or programs in the devices and methods described in the embodiments herein may be recorded, stored or fixed in one or more non-temporary computer-readable media (computer-readable storage / recording media), which include program instructions (computer-readable instructions) that can be executed by a computer so that one or more processors execute (implement) the program instructions. The medium may also include data files, data structures, etc., alone or in combination with program instructions. The media and program instructions may be specially designed and constructed, or they may be well-known and available to technicians in the field of computer software. Examples of non-temporary computer-readable media include: magnetic media, such as hard disks, floppy disks, and tapes; optical media, such as CD ROMs and DVDs; magneto-optical media, such as optical disks; and hardware devices specially configured for storing and executing program instructions, such as read-only memory (ROM), random access memory (RAM), flash memory, etc. Examples of program instructions include machine code (such as code generated by a compiler) and files containing high-level code (which can be executed by a computer using an interpreter). Program instructions can be executed by one or more processors. The hardware device may be configured to act as one or more software modules, which are recorded, stored or fixed in one or more non-transitory computer-readable media to perform the above operations and methods, and vice versa. In addition, the non-transitory computer-readable medium may be distributed between computer systems connected via a network, and program instructions may be stored and executed in a decentralized manner. In addition, the computer-readable medium may also be embodied in at least one application-specific integrated circuit (ASIC) or field-programmable gate array (FPGA).

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

[0065] Figure 5 1 is a flow chart showing a seismic full waveform inversion method using a Hann window according to an embodiment of the present disclosure. The survey area may be an underground structure under land or an underground structure under the ocean.

[0066] refer to Figure 5In 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 Figure 3 The logging data recording sensor of the logging tool 315 can detect seismic data and transmit the seismic data to the Figure 4 The high performance computing system shown in FIG. Figures 1 to 4 As described above, seismic data detected in the survey area may be stored in one or more memories, such as one or more storage devices 410 and one or more output storage devices 420 .

[0067] In addition, if Figure 5 As shown, the initial earth model m0 can be input into the seismic full waveform inversion method. The initial earth model m0 can be a P-wave velocity model.

[0068] It is understood that the initial velocity model can be a set of parameters. The initial velocity model m0 can be a predetermined velocity model based on established scientific specifications. The initial velocity model m0 can be input by a user through a user interface (such as a user interface of a personal computer 425), and can also be stored in one or more memories (such as one or more storage devices 410 and one or more output storage devices 420).

[0069] like Figure 5 As shown in the seismic full waveform inversion method, the initial model (such as the initial velocity model m0) is input into the loop. In operation 500, the current velocity model m is determined. i Initially, the velocity model m i is the initial velocity model m0. The variable i represents the number of iterations of the loop. Therefore, i = 0 at the beginning. Figure 5 In each iteration of 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 is output in operation 575. i+1 , which corresponds to the latest velocity model in operation 510, and is output to a storage device or a display.

[0070] Reference Figure 5 In operation 510 of the seismic full waveform inversion method shown in FIG. 1 , a velocity model m is input. iForward modeling of seismic data is a technique for creating (generating) synthetic seismic data based on geological information, in this case, the velocity model m i The forward modeling operation requires a source wavelet, which may be a known source wavelet function that is different from the true source characteristics of the seismic survey being conducted or analyzed. Examples of source wavelet functions include Ormsby wavelet, Ricker wavelet, etc. The source wavelet is used in the forward modeling operation to obtain the velocity model m. i To generate synthetic data d syn . Synthetic data d syn Reference numeral 520 represents the synthesized data output by the forward modeling operation 510 .

[0071] exist Figure 5 In operation 525, the synthesized data d syn is input into the Hann window, and in operation 515 the real data d obs Performs a Hann window operation. The Hann function is a window function used to separate a complete / long signal into a series of local / short signals along the recording time. Figure 7 A sliding Hann window is shown, whose size and sliding step are automatically determined based on the data frequency, for automatically localizing the seismic data. The sliding Hann window can effectively transform the data into the local domain without energy loss.

[0072] In operation 515, formula (5) represents the use of Hann function W k (t) The kth localized data d processed obs,s,r (t):

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

[0074] Formula (6) represents the sum of localized data:

[0075]

[0076] It is important to note that the sum of the localized data is equal to the original data. Therefore, the data localization process using Hann windowing does not cause energy loss.

[0077] Formula (5) can be expanded into

[0078]

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

[0080] The same Hann windowing process is also applied to the synthetic data in operation 525, and produces localized synthetic data d obs,s,r,k (t). Therefore, the total number of windows per seismic trace is N k The calculation formula is:

[0081]

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

[0083] In operation 530, the time shift of the kth local data set may 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 non-steady-state changes of travel time differences of different seismic events.

[0086] In operation 535, the localized observation data d is obtained by using ΔT(s, r, k) in operation 530 and the Hann windowing in operation 515. obs,s,r (t) and the associated source is obtained as follows:

[0087]

[0088] Referring to the FWI gradient g in operation 550 i In operation 550, one or more companion sources generated by formula (10) are back-propagated to obtain a companion wave field, which is calculated as follows:

[0089]

[0090] in, is the adjoint operator of the forward model operator F(m;x), d adj,s (t) is the associated source, and u s (x,t) is the companion wave field.

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

[0092]

[0093] Where F(m;x) is the forward operator of the wave equation, w s (x, t) is the forward wave field, and u s(x, t) is the inverse companion wave field of equation (11). In operation 550, the gradients as described in equation (11) from different shot points (blasting or sonicator) are added and superimposed to obtain a single gradient.

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

[0095] More specifically, after the gradient is calculated and outputted in operation 550, the step length a is calculated using an optimization method in operation 555. i and search direction P i For example, starting with 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] Among them, P i is the search direction, a i is the step size. The simplest inversion method is called the steepest gradient algorithm, where the search direction P i It 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] Then, the step size is determined in a way that minimizes the cost function.

[0102] In operation 555, a gradient is received and a search direction P is determined. i (used to estimate the increase or decrease in the velocity model) and the step size a i (used to estimate the magnitude of the increase or decrease in velocity toward the true velocity model). Therefore, in operation 555, the FWI gradient g in operation 550 is iThe output of is scaled to provide the optimal amplitude of the directional increase of the estimated velocity model to avoid excessive increase / decrease of the estimated velocity model (which will greatly increase the number of iterations) or excessive increase / decrease of the estimated velocity model (which will greatly increase the number of iterations). Thereafter, in operation 560, the velocity model m is updated according to formula (15) i+1 .

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

[0104] In one embodiment, a maximum number of iterations may be set in operation 565. The maximum number of iterations may range from 10 to 40. Once the number of iterations represented 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 output the final velocity model m F , which is consistent with the last updated velocity model m determined in operation 575 i+1 same.

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

[0106] Figure 6 is a workflow of another embodiment of the present disclosure, describing a source wavelet-independent FWI method that uses a matched filter to correct the synthetic data before cross-correlation calculation and travel time extraction. Figure 6 Operations 600, 605, 610, 615, and 620 in Figure 5 Operations 500, 505, 510, 515, and 520 in correspond to and perform the same functions. Figure 5 The difference between the method shown is that operation 680 obtains the matched filter M s (t), its mathematical expression is:

[0107]

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

[0109] Referring to the matched filter shown in equation (16), the matched filter M s (t) The observed seismic data d detected in the exploration area in operation 605 are obs The synthesized data d is obtained in operation 620 syn Match. Observed earthquake data d obs The data may be data recorded by sensors (eg, seismic data recording sensors 105 in the survey area, and well logging data recording sensors of the well logging tool 315 disposed in the wellbore of the survey area).

[0110] Matched filter M s An example of (t) may 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 uses the calculated synthetic data to determine in each iteration operation 665 in order to ultimately update the velocity model m i . Using synthetic data d syn The matched filter is convolved to generate matched synthetic data having the same source characteristics as the observed seismic data detected in the survey area.

[0111] As mentioned above, errors in the source wavelet may lead to errors in the travel time difference picking, thereby causing errors in the inversion model. Therefore, a reasonable source wavelet estimation is still very important for the success of cross-correlation based FWI. The typical process of this estimation is to window the first arrival event in the data and superimpose the event as the source wavelet. This estimated source wavelet does not change during the iteration process. However, it is difficult to achieve accurate estimation of the seismic source wavelet in industrial field applications. This is because the repeatability of the sound source characteristics of different shot points (different blasting or different sound sources) is poor, the coupling between the sound source and the earth is uncertain, and the coupling between the receiver and the earth is also uncertain. Therefore, in order to overcome these problems, a lot of effort has been spent on the source-independent misfitting function method. The operation 680 of the matched filter is different from other operations (such as Figure 6 These problems can be solved by combining the accompanying seismic source obtained in operation 635.

[0112] The synthesized data from operation 620 is first compared with the matched filter M obtained in operation 680. s (t) performs convolution. The synthetic data M after convolution s (t)*d syn,s,r,k(t) (* denotes the convolution operator) is more similar to the observed data, and then localized using the Hann window in operation 625 to obtain d obs,s,r,k (t).

[0113] In operation 630, the automatic travel time difference picking in the local window is defined as:

[0114]

[0115] Here, the shot-independent matched filter M is calculated by matching the synthetic data with the observed seismic data. s (t).

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

[0117]

[0118] In operation 635, the chain rule can be used to obtain the accompanying source:

[0119]

[0120] in, Denotes the cross-correlation operator. In each iteration, the Wiener matched filter M is determined shot by shot. s (t) to correct the phase of the companion source. The final companion 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 Operations 540, 545, 550, 555, 560, 565, 570, and 575 in correspond to and are substantially similar.

[0122] Therefore, in Figure 6 In the method shown, the gradient of the cost function can be constructed using the companion source from operation 635. For example, using the forward wavefield and the companion wavefield, the FWI gradient g can be obtained in operation 650 by solving the companion wave equation (11) i , which is then 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. Operations 550 and 650 are substantially similar.

[0123] Nevertheless, it is clear from comparing the companion sources obtained in operation 635 according to equation (19) with the companion sources obtained in operation 535 according to equation (10) that there is an additional cross-correlation operation in operation 635, namely, s (t) is the cross-correlation operation.

[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 require the updating of 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 independence of the source wavelet reduces the model inversion error caused by the source wavelet error and reduces the workload of manually estimating the source wavelet shot by shot. This makes cross-correlation-based FWI very robust and can be used to build high-resolution geological models to improve the image of the complex underground structure of 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 FWI based on adaptive cross-correlation can robustly invert medium parameters. According to one embodiment, each iteration includes a series of steps, for example: Step 1, using a theoretical wavelet (such as Ormsby wavelet) and the medium parameter model m in the current i-th iteration i , using the wave equation (1) to perform forward modeling, and use it to obtain the underground forward wave field w s (m i ; x, t), and estimate the synthetic data d recorded at the receiver location syn,s,r (m i ; t); Step 2, using formula (16) to determine a matched filter, such as a Wiener filter; Step 3, using formula (10) to automatically select the travel time, and using the sliding Hann window described in formula (17) to automatically determine the local window; Step 4, using formula (19) to calculate the companion source; Step 5: Use the same model to back-propagate the companion source to obtain the companion wave field u s (m i ; x, t); Step 6, using the obtained forward wave field w s (m i ; x, t) and the accompanying 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 length α iStep 9: Update the medium parameter model with the step size. Repeat steps 1 to 9 until the unfitting function, i.e., equation (18), is minimized to a predetermined value.

[0127] Figures 8 to 17 It is shown that the method disclosed in this paper can invert high-resolution subsurface velocity details while avoiding cycle jumps. Figure 8 The true velocity model with a 50-m grid spacing is shown, used to generate seismic data with a maximum offset of 25 km. Fig. 9 The initial velocity model used for FWI, the Marmousi model, is shown. In this inversion, the medium parameter m has only one parameter, the propagation velocity V P The data used by the FWI is from Fig. 9 The real Masmousi model in and generated by Ormsby wavelet. Figure 8 and Fig. 9 The velocity difference in the spherical plane is in the same range between 1500 m / s and 4500 m / s. Fig.10 Shows from Figure 8 The true velocity model in Fig. 9 The velocity difference is obtained after the initial velocity model is used. Fig.10 The velocity difference ranges from -1365 m / s to 1980 m / s.

[0128] Fig.11 Shows the use Figure 6 The inversion velocity model of the adaptive cross-correlation FWI method of the present invention has 2.5 Hz data that has been band-pass filtered, assuming that the effective lowest frequency is 2.5 Hz and there is no lower frequency data in the observed data. Fig.12 Shows the use Figure 6 The final inverted velocity model of the adaptive cross-correlation FWI method of the present invention has data of 2.5 Hz, 4 Hz, 6 Hz and 10 Hz that have been continuously band-pass filtered.

[0129] Fig.13 The final inverted velocity model using the conventional L2 waveform difference based FWI method according to equation (4) is shown. Fig.13 The same data is successively bandpass filtered at 2.5 Hz, 4 Hz, 6 Hz and 10 Hz. Fig.13 and Fig.12 By comparison, it is easy to see the inversion artifacts in the area enclosed by the rectangle, and the wrong update direction due to cycle jumps in the area enclosed by the ellipse.

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

[0131] Fig.16 shows an earthquake migration image obtained using the initial ray-based tomography model, while Fig.17 Shown through Figure 6 The seismic migration image obtained by the adaptive cross-correlation FWI method of the present invention and using the inversion velocity model is shown. Compare the area indicated by the arrow or surrounded by the ellipse. Fig.17 The geological features in the image are smoother and clearer, which shows the effectiveness of the method of the present invention.

[0132] Therefore, despite Fig.16 and Fig.17 shows an example of land data with unknown true source wavelets derived from different shot records and forward modeling using the Ormsby wavelet, but Figure 6 The described method of the invention successfully updates the velocity model, thereby improving the seismic image and obtaining a better geological interpretation, which shows the robustness of the method of the invention, especially when the true wavelets of different shot points of the land data are not known.

[0133] Although the embodiments disclosed herein describe an FWI algorithm with a travel time (difference) misfit function that utilizes a finite difference numerical solution of the scalar acoustic wave equation for seismic propagation, it is apparent to one of ordinary skill in the art that the algorithm may alternatively be applied to vector wave equations and elastic equations in isotropic and anisotropic media without departing from the true scope of the invention as defined in the following claims.

Claims

1. A method for performing seismic full waveform inversion to generate a velocity model of underground strata in a survey area, comprising the following operations: (a) placing seismic data recording sensors at a plurality of locations within the survey area, and placing a logging tool having seismic data recording sensors in a wellbore within the survey area; (b) transmitting at an entry point in the survey area to generate seismic waves, the seismic waves penetrating the subsurface formations; (c) observing seismic waves using the seismic data recording sensor, and recording seismic data based on the seismic waves; (d) transmitting observed seismic data from the seismic data recording sensor to a computer system including one or more memories, storing the observed seismic data in the one or more memories, and storing source wavelets and current medium parameter models in the one or more memories; (e) the computer system performs a forward modeling operation according to a forward wave equation using a source wavelet and a current medium parameter model to obtain a forward wave field, and generates synthetic data according to the forward wave equation; (f) processing the synthetic data and the observed seismic data by the computer system using a Hann window operation to obtain localized synthetic data and localized seismic data; (g) selecting the travel time difference by the computer system based on the cross-correlation between the localized synthetic data and the localized seismic data; (h) solving the adjoint equation of the forward wave equation by the computer system using the travel time difference and the localized seismic data to obtain an adjoint seismic source; (i) generating, by the computer system, an updated medium parameter model for seismic full waveform inversion and synthetic data using a forward modeling operation; (j) perform operations (e) to (i) until convergence; (k) after convergence, outputting the updated velocity as the final velocity model to a display; as well as (l) Displaying an image of the final velocity model on a display of the computer system.

2. The method according to claim 1, wherein: Also included is constructing, by the computer system, a matched filter that matches the synthetic data to observed seismic data.

3. The method according to claim 2, wherein: The synthetic data is convolved with the matched filter before performing a Hann window operation.

4. The method according to claim 3, wherein: The operation (g) further comprises, by the computer system: Reversely propagating the companion source to obtain a companion wave field; Calculating a gradient using the forward wavefield and the adjoint wavefield; Determine the search direction and step size; and The medium parameter model is updated according to the step size.

5. The method according to claim 1, wherein: When the value of the misfit function based on travel time is smaller than a predetermined value, it is determined that convergence has occurred.

6. The method according to claim 1, wherein: When the number of iterations reaches a predetermined value, convergence is determined.

7. The method according to claim 1, wherein: The medium parameter model is selected from a velocity model.

8. The method according to claim 1, wherein: The source wavelet is Ormsby wavelet or Ricker wavelet.

9. The method according to claim 2, wherein: The matched filter is a Wiener filter.

10. The method according to claim 4, wherein: The inversion method used in the updating step is selected from the steepest gradient algorithm, the nonlinear conjugate gradient method, the Gauss-Newton method and the quasi-Newton method.

11. A system for performing seismic full waveform inversion to generate a velocity model of underground strata in a survey area, the system comprising: a plurality of seismic data recording sensors disposed at different locations within the survey area, and / or a well logging tool including seismic data recording sensors disposed in a wellbore within the survey area; an blasting device arranged at each entry point of the survey area for generating seismic waves, wherein the seismic waves penetrate underground formations; as well as a plurality of seismic data recording sensors for sensing seismic waves and recording seismic data based on the seismic waves; The seismic data recording sensor transmits seismic data to a computer system including one or more memories and at least one processor, wherein the one or more memories store the transmitted seismic data, source wavelets and instructions, and the one or more processors execute the instructions stored in the one or more memories to implement: Performing forward modeling operations according to the forward wave equation and using the source wavelet and the current medium parameter model to obtain the forward wave field; Performing Hann window operation on synthetic data and observed seismic data respectively to obtain localized synthetic data and localized seismic data; selecting the travel time difference based on a cross-correlation between the localized synthetic data and the localized seismic data; Solving the adjoint equation of the forward wave equation using the travel time difference and localized seismic data to obtain an adjoint seismic source; and Forward modeling operations are used to generate updated medium parameter models for seismic full waveform inversion and synthetic data.

Images