Method and apparatus for seismic characterization of subsurface structures

Abstract

A method and apparatus for performing seismic full waveform inversion to generate a final velocity model of subsurface formations of a survey region independent of a seismic wavelet is provided. The method includes deploying a plurality of seismic data recording sensors (geophones) at a survey region; exciting seismic waves at a plurality of source points at the survey region; and receiving and recording seismic wave signals by the seismic data recording sensors. The recorded seismic wave signals are seismic data. The method further includes transferring the seismic data to a computer system including one or more storage devices and storing the seismic data in the one or more storage devices; storing seismic wavelet signals in the one or more storage devices; performing, by the computer system, forward modeling operations based on the seismic wavelet signals; and generating, by the computer system, a final velocity model of seismic full waveform inversion by the forward modeling operations.

Description

Technical Field

[0001] This disclosure generally relates to a seismic full-waveform inversion method and apparatus. More specifically, this disclosure enhances the ability to characterize complex subsurface structural features of the survey area by generating a high-resolution geological model through an improved seismic full-waveform inversion method. Background Technology

[0002] As is well known, seismic imaging technology, which analyzes seismic wave data received at the surface to characterize complex subsurface structures, plays a crucial role in seismic exploration and reservoir characterization for lithological differentiation and fluid identification. The various seismic imaging techniques (techniques or methods) widely used in industry can be mainly classified into two categories: migration imaging and inversion. Seismic migration imaging methods can obtain reflectivity images of subsurface structures. The earliest such method involved weighted stacking of pre-stack seismic data, proposed by Smith and Gidlow (Smith GC and PMG Gidlow, 1987, "Weighted stacking for rock property estimation and detection of gas", published in Geophysical Prospecting, Vol. 35, pp. 993-1014). The drawback of this method is that the resulting seismic image is not a true mapping of subsurface reflectivity, and the provided location of the reflecting layer is not its true location in the time or depth domain. For seismic media with smooth lateral variations, this deficiency can be overcome by the Kirchhoff migration method.

[0003] By the 1980s, reverse time migration (RTM) had quickly become a high-end seismic imaging method for complex geological structures (see Baysal E., D.D. Kosloff, and J.W. C. Herwood, 1983, “Reverse time migration,” Geophysics, Vol. 48, pp. 1514–1524; and Etgen, J., 1986, “Pre-stack reverse time migration of shot profiles,” SEP-Report, Vol. 50, pp. 151–170). To obtain high-resolution seismic images using RTM, two main conditions must be met: (1) designing a high-density and high-coverage seismic data acquisition system; and (2) using a macroscopic model of the seismic medium or a background velocity model with at least correct kinematic information during the migration imaging process. For relatively simple geological targets, such as shallow water environments, establishing an accurate background velocity model is relatively straightforward; however, for complex geological environments (such as salt domes, subbasal structures, thrust zones, and piedmont areas), seismic modeling becomes extremely challenging.This type of seismic modeling process belongs to the second category of imaging methods, which are achieved by obtaining rock density and pressure propagation velocity. Examples include: refracted wave / first-arrival tomography inversion methods, such as those proposed by Osypov (Osypov K., 2001, "Refraction Tomography: A Practical Overview of Emerging Technologies," published in CSEG Recorder, Vol. 26); seismic travel-time tomography inversion methods, such as those proposed by Bording et al. (Bording RP, A. Gersztenkorn, LRLines, JAScales, and S. Treitel, 1987, "Applications of seismic travel-time tomography," published in Geophysical Journal International, Vol. 90, pp. 285–303); and full-waveform inversion methods, such as the work of Lailly (Lailly P., 1983, "The seismic inverse problem as a sequence of before stack migrations," published in the conference proceedings of Inverse Scattering, Theory and Application, Society for Industrial and Applied Sciences). (Mathematics, Extended Abstract, pp. 206–220) and Tarantola (Tarantola A., 1984, “Inversion of seismic reflection data in the acoustic approximation”, Geophysics, Vol. 49, pp. 1259–1266). Velocity models established using these tomographic modeling imaging (Tomography) or full waveform inversion (FWI) methods can be used in seismic migration imaging, significantly improving its imaging quality and resolution. While travel-time tomography typically only inverts smooth velocity models with correct kinematic information, which are sufficient for seismic imaging, full waveform inversion (FWI) can acquire high-resolution seismic P-wave (pressure wave) velocities, S-wave (shear wave) velocities, anisotropy parameters, and density models, which are crucial for geological interpretation.

[0004] Due to the complexity of subsurface media models, accurate modeling of actual seismic waves propagating underground is often difficult. However, subsurface media models can be approximated to varying degrees, such as using elastic or acoustic wave models in isotropic or anisotropic media, with or without considering attenuation effects. Seismic data is acquired through seismic data detectors (recorders), and this data relates to the reflection and refraction of seismic waves beneath the surface in response to a seismic source. The seismic source can be a seismic device that triggers an explosion, or other source-generating devices that can generate seismic waves below the surface. Full waveform inversion (FWI) involves numerically simulating the propagation of seismic waves in an iteratively updated subsurface model. In the early stages of FWI development, seismic waves were often simplified to pure acoustic wave models because acoustic equations are relatively simple and easy to solve efficiently. With improved computing power, elastic full waveform inversion has become feasible and more attractive because elastic waves more realistically reflect the characteristics of actual subsurface seismic waves compared to acoustic waves. Regardless of the wave equation used, it can be mathematically expressed as:

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

[0006] Where m represents the underground model vector, x is the spatial location vector, t represents time, F(m; x) is the corresponding forward modeling operator, and w s (x,t) represents a source wavelet f excited at spatial location s. s The orthogonal wave field generated by (t).

[0007] However, the analytical solution of the wave equation (1) exists only in simple models. For the vast majority of cases, numerical methods must be used. The core of all numerical methods lies in discretizing the spatial variable x and the time variable t, thereby discretizing the solution w s (x,t) is approximately w s (x i ,t j ), where i = 1, ..., N represents grid points, with a total of N spatial grid points, and j = 1, ..., M represents time sampling, with a total of M time sampling points. Commonly used numerical methods include the finite difference method, the finite element method, and the spectral element method (Fichtner, A., 2011, "Full Seismic Waveform Modelling and Inversion", Springer).

[0008] As mentioned above, reverse time migration (RTM) is a well-known seismic imaging method in this field. RTM is mainly used to map subsurface reflectivity structures using recorded seismic waveforms. RTM typically involves three steps: (a) generating a forward-modeled wavefield using a suitable velocity model; (b) backpropagating the recorded seismic data in reverse time using the same velocity model; and (c) imaging and overlaying the forward and reverse wavefields using imaging conditions. In most cases, the RTM method aims to obtain an initial image of subsurface reflectivity by achieving an optimal match between the backpropagated waveform and the forward-modeled waveform obtained based on the estimated velocity model and source parameters in image space. Therefore, the quality of the RTM image can serve as an important criterion for evaluating the quality of the velocity model (also known as the FWI velocity model or inverted velocity model) generated by full waveform inversion (FWI).

[0009] Full Waveform Inversion (FWI) has been further discussed in the work of Tarantola (Tarantola A., 1984, “Inversion of seismic reflection data in the acoustic approximation”, Geophysics, Vol. 49, pp. 1259–1266) and Virieux et al. (J., A. Asnaashari, R. Brossier, L. Métivier, A. Ribodetti, and W. Zhou, 2014, “An introduction to full waveform inversion”, Geophysical References Series, R1–R40). Full Waveform Inversion (FWI) is a data-driven tool for automatically constructing subsurface medium parameters m, such as velocity, anisotropy parameters, and / or density models. This method achieves this by iteratively minimizing the discrepancies between seismic data and simulated synthetic data. Given initial subsurface parameters m0 (such as velocity models, anisotropy parameter models, and density models), it is achieved by solving for the source wavelet f. s The forward-modeled seismic wave equation (Equation 1) of (t) is used to calculate synthetic seismic data. Then, using the mismatch between the seismic data and the synthetic data as the source, the adjoint equation corresponding to the forward-modeled wave equation is solved to obtain the solution of the adjoint wave field. Next, the gradient information is calculated by cross-correlation between the forward-modeled wave field and the adjoint wave field, and the subsurface model is updated along the gradient direction to further reduce the difference between the seismic data and the predicted data. The above process is repeated iteratively until the data mismatch is small enough to obtain the optimized subsurface medium model.

[0010] From a mathematical perspective, full waveform inversion (FWI) can be formulated as an optimization problem:

[0011]

[0012] Where C(m) is the mismatch (error) function, representing the distance of error between the seismic data and the synthetic data. The modulus operation is represented by the symbol |·|, and 2 represents the L2 norm. In the mismatch function, N... s It is the number of earthquake focal points ("s" represents "focal point"), N r It is the number of detectors corresponding to a given seismic source ("r" represents "detector"), T is the maximum recording time starting from zero, and d is the number of detectors. obs,s,r (t) represents the seismic data recorded at time t and at the detector r of the given source s, d syn,s,r (m; t) is the corresponding synthetic seismic data r, which is derived from the solution w of the forward seismic equation (1). s (t) Sampled at detector r: d syn,s,r (m;t)=P s,r w s (m;t), where P s,r This represents the corresponding sampling function. The mismatch function under the L2 norm (Equation 2) is called the least squares mismatch function.

[0013] Full waveform inversion (FWI), as an optimization problem, is a highly nonlinear process. Its success in practical applications depends on many factors, such as the choice of the wave equation used to describe the subsurface wavefield, the initial model m0 used in the wavefield simulation, and the source wavelet f. s (t), etc., these are only a part of it. All of these factors are unknown. Among these factors, the unknown source wavelet is a key issue in FWI. A small error in the source wavelet can lead to large discrepancies in the inversion model, and this error accumulates with increasing depth.

[0014] To improve the imaging quality of complex subsurface structures within the exploration area, a new technique is needed that is independent of the accuracy of the source wavelet in the wave equation (1), which is referred to in this disclosure as "source wavelet-independent" or "source wavelet-independent". This new technique should also be applicable to full waveform inversion of acoustic and elastic waves. Summary of the Invention

[0015] The inversion model parameters generated by one or more embodiments are independent of the source wavelet and are used to perform full-wavelength seismic inversion to generate a high-resolution geological model, thereby enabling high-resolution imaging of exploration areas containing complex subsurface structures.

[0016] One or more embodiments do not depend on the source wavelet in the wave equation (1) and are applicable to full waveform inversion of acoustic waves and full waveform inversion of elastic waves.

[0017] One or more embodiments overcome the shortcomings of unknown source wavelets in full waveform inversion (FWI) by providing a novel method to modify the error function in equation (2) above. This method eliminates the need to estimate the source wavelet f in the forward seismic wave equation (1). s (t) allows for inversion, generating model parameters independent of the source wavelet selection. This method overcomes a major deficiency in full waveform inversion (FWI), enabling the construction of high-resolution geological models to improve the imaging quality of complex subsurface structures in exploration areas. This, in turn, enhances the capabilities for lithology identification, fluid identification, and hydrocarbon reservoir description, and is widely applied in the field of seismic exploration.

[0018] In one aspect, a method is provided for generating a final velocity model of subsurface strata in an exploration area by running a full-waveform seismic inversion. The method may include the following steps: (a) deploying seismic data recording sensors (detectors) at different locations in the exploration area and / or arranging logging tools containing seismic data recording sensors in wells within the area; (b) generating (exploding) incident seismic waves at the source point in the exploration area, these seismic waves penetrating the subsurface strata; (c) measuring the seismic wavefield using the seismic data recording sensors (detectors) and recording the observed seismic data generated by these seismic wavefields; (d) transferring the seismic data from the seismic data recording sensors to a computer system including one or more storage devices, and storing the seismic data in one or more storage devices; (e) in a (f) The source wavelet is stored in one or more storage devices; (g) The computer system performs forward modeling using the source wavelet and the current velocity model; (h) The computer system matches the seismic data and synthetic data, generates the accompanying source and determines the inversion gradient and velocity update step size based on the mismatch between the seismic observation data and the synthetic data, and calculates the seismic full waveform inversion to generate the updated velocity model; (i) Steps (f) and (g) are repeated until convergence; (j) After convergence, the updated velocity model is output as the final velocity model to the display device; (h) A high-resolution image of the final velocity model is displayed on the display device of the computer system.

[0019] In one aspect, the seismic full-waveform inversion can be either acoustic inversion or elastic wave inversion, and the subsurface strata in the exploration area can be isotropic or anisotropic.

[0020] In one aspect, one or more storage devices store an initial speed model.

[0021] In one aspect, forward modeling operations generate forward wavefields and forward modeling simulation data.

[0022] In one aspect, via a computer system, step (g) may further include matching actual observed seismic data and synthetic seismic data, and calculating the mismatch (or residual) between the seismic data and the synthetic seismic data based on a mismatch function. The mismatch function includes:

[0023]

[0024] Where * denotes the convolution operator, m is the Earth model vector, t is time, and M s (t) is the matched filter function, which matches the synthetic data with the seismic data using the least squares Wiener filter operator.

[0025]

[0026] Where, d obs,s,r (ω) and d syn,s,r (ω) represents the Fourier transform of seismic data and synthetic seismic data, respectively.

[0027] In one aspect, step (g) may further include generating an accompanying seismic source based on the mismatch between the seismic data and the synthetic seismic data.

[0028] In one respect, the accompanying seismic source can be generated according to the following equation:

[0029]

[0030] in, This represents the cross-correlation operator.

[0031] In one aspect, step (g) may also include determining the inversion gradient based on the accompanying seismic source.

[0032] In one aspect, step (g) may further include determining the step size and direction for updating the velocity model, and updating the velocity model.

[0033] In one respect, the source wavelet can be an Ormsby or Ricker wavelet.

[0034] In one aspect, a system is provided for generating a subsurface velocity model of an exploration area by running a full-waveform seismic inversion. The system may include: deploying seismic data recording sensors at different locations within the exploration area and / or arranging logging tools containing seismic data recording sensors within exploration wells in the area; generating incident seismic waves at the source point of the exploration area, the incident seismic waves propagating through the subsurface strata; and measuring the seismic wave field using the seismic data recording sensors and recording the observed seismic data generated by the seismic wave field. The observed seismic data recording sensors can transmit seismic data to a computer system, the computer system including one or more storage devices and at least one processor. The one or more storage devices may store the transmitted seismic data, source wavelet, and instruction set, and the one or more processors may execute the instruction set stored in the one or more storage devices to perform the following functions: (a) performing forward modeling operations using the source wavelet; (b) generating an updated velocity model using full waveform inversion and forward modeling operations; (c) performing steps (a) and (b); (d) outputting the updated velocity model as the final velocity model to a display device after convergence; and (e) displaying a high-resolution image of the final velocity model on a display device of the computer system.

[0035] In one aspect, step (b) may further include matching seismic data with synthetic data and evaluating the mismatch between the seismic data and the synthetic data using a mismatch function, wherein the mismatch function includes:

[0036]

[0037] Where * denotes the convolution operator, m is the Earth model vector, t is time, and M s (t) is the matched filter function, which matches the synthetic seismic data with the observed seismic data using the least squares Wiener filter operator.

[0038]

[0039] Where, d obs,s,r (ω) and d syn,s,r (ω) represents the Fourier transform of seismic data and synthetic seismic data, respectively.

[0040] In one aspect, step (b) may also include generating an accompanying seismic source based on the mismatch between the seismic data and the synthetic data.

[0041] In one respect, the accompanying seismic source can be generated according to the following equation:

[0042]

[0043] in, This represents the cross-correlation operator.

[0044] In one aspect, step (b) may also include determining the inversion gradient based on the accompanying seismic source.

[0045] In one aspect, step (b) may further include determining the step size and direction for updating the velocity model, and updating the velocity model. Attached Figure Description

[0046] This patent or patent application document contains at least one color drawing. Upon request and payment of the appropriate fee, the competent authority will provide a copy of the published patent or patent application with the color drawing. The technical content of this invention can be clearly understood by reading the following detailed description in conjunction with the accompanying drawings.

[0047] Figure 1 This is a schematic diagram showing a top view of an exploration area according to an embodiment of the present disclosure, wherein the locations of multiple seismic source incident points are marked.

[0048] Figure 2 It is a schematic diagram showing an environmental profile according to an embodiment of the present disclosure, including the earthquake source incident point, the earthquake data recording sensor, the well location, the well casing, the propagation rays, and the incident angles.

[0049] Figure 3 This is a schematic diagram showing an environmental profile according to an embodiment of the present disclosure, including an exploration wellbore and a logging tool comprising one or more acoustic exciters and one or more logging data recording sensors.

[0050] Figure 4 This is a schematic diagram illustrating a high-performance computing system according to an embodiment of the present disclosure.

[0051] Figure 5 This is a flowchart illustrating a seismic full-waveform inversion method according to an embodiment of the present disclosure for generating a final velocity model to improve the imaging of complex subsurface structures in an exploration area.

[0052] Figure 6 It is a color image showing a comparison between an example of a real P-wave velocity model and an example of a preset initial P-wave velocity model in an exploration area profile.

[0053] Figure 7 This is an example of a source wavelet used for forward modeling operations in one embodiment of this disclosure.

[0054] Figure 8 This is an example of a source wavelet used for forward modeling operations in one embodiment of this disclosure.

[0055] Figure 9It is a chart that shows the convergence comparison when using different source wavelets in forward modeling operations.

[0056] Figure 10 It is a color image showing an example of a P-wave true velocity model in a vertical profile of an exploration area, compared with an example of a P-wave velocity model obtained using the traditional seismic full waveform inversion method.

[0057] Figure 11 It is a color image showing an example of a P-wave true velocity model in a vertical profile of the exploration area, compared with an example of a P-wave velocity model obtained using the seismic full waveform inversion method and source wavelet of this embodiment.

[0058] Figure 12 It is a color image showing an example of a P-wave true velocity model in a vertical profile of the exploration area, compared with an example of a P-wave velocity model obtained using the seismic full waveform inversion method of this embodiment and another source wavelet.

[0059] Figure 13 It is a color image showing the updated velocity model obtained using the traditional seismic full waveform inversion method and the reverse time migration image obtained using the updated velocity model.

[0060] Figure 14 It is a color image showing the updated velocity model obtained by the source wavelet independent seismic full waveform inversion method according to this embodiment, and the reverse time migration image obtained using the updated velocity model. Detailed Implementation

[0061] Reference will now be made in detail to embodiments of this disclosure, examples of which are illustrated in the accompanying drawings. It should be noted that, where feasible, similar or identical reference numerals may be used in the drawings; these numerals may denote similar or identical elements.

[0062] These illustrations are for illustrative purposes only and show embodiments of the present disclosure. Those skilled in the art will readily recognize from the following description that alternative embodiments exist without departing from the general principles of the present disclosure.

[0063] In this disclosure, the terms "method" and "approach" are used interchangeably and have the same meaning.

[0064] This disclosure relates to the establishment of high-resolution geological models by performing improved seismic full waveform inversion to improve imaging of complex subsurface structures (strata) in exploration areas, and the methods, equipment, and media used in performing the improved seismic full waveform inversion, as well as the inclusion of one or more source-independent mismatch (target) functions.

[0065] Figures 1 to 4 Exemplary embodiments of methods, apparatus, and media for acquiring and storing seismic data are illustrated. This seismic data, after processing, can generate one or more high-resolution geological models for high-resolution imaging of complex subsurface structures in an exploration area, enabling lithological identification, fluid identification, and reservoir description. The exploration area can be a subsurface structure beneath land or beneath the seabed.

[0066] Figures 5 to 14 Exemplary embodiments of apparatus, methods, and media for improving the quality of seismic full waveform inversion results are shown, which enhance the effectiveness of lithology identification, fluid identification, and reservoir description in the field of seismic exploration by using improved seismic full waveform inversion techniques, including computer-based seismic full waveform inversion methods. Figures 5 to 14 Exemplary embodiments of apparatuses, methods, and media for generating one or more high-resolution geological models are shown to enable high-resolution imaging of complex subsurface structures in exploration areas for lithological identification, fluid identification, and reservoir characterization. Figures 5 to 14 An exemplary embodiment is shown that can generate source wavelet-independent inversion model parameters, which can be used to perform full-waveform seismic inversion to generate high-resolution geological models, thereby enabling high-resolution imaging of exploration areas including complex subsurface structures. Figures 5 to 14 Exemplary models are shown that, under conditions independent of the source wavelet, they are applicable to full waveform inversion of acoustic waves and full waveform inversion of elastic waves, for generating one or more high-resolution geological models to achieve lithological identification, fluid identification, and reservoir description of complex underground structures in exploration areas. Figures 5 to 14 Exemplary embodiments of apparatus, methods, and media for improving the acquisition speed of seismic full waveform inversion results and reducing computational resource consumption by using an improved seismic full waveform inversion technique independent of the source wavelet are also demonstrated.

[0067] Figure 1 This is a schematic diagram showing a top view of a seismic exploration area according to an embodiment of the present disclosure, wherein the incident points of multiple seismic sources are marked. More specifically, Figure 1A seismic exploration area (exploration area) 101 is shown, which is a land area, denoted by reference numeral 102. Reference numeral 102 indicates the surface strata of the land area. Those skilled in the art will recognize that detailed images of the local geological structure are generated in the seismic exploration area to determine the specific location and size of possible oil and gas reservoirs (hydrocarbon reservoirs), and therefore also include a drilling location 103. In these exploration areas, seismic waves are reflected in the subsurface rock strata after being emitted from multiple incident points 104 from one or more seismic source points. An explosion is an example of a seismic source generated by a seismic device. The seismic waves reflected back to the surface are captured by a seismic data recording sensor 105 and transmitted by the seismic data recording sensor 105 via one or more data transmission systems (typically wireless), and then stored for subsequent processing and analysis by a high-performance computing system. Although this embodiment shows the surface strata 102 of a land area, it should be understood that this is only an example, and the method and system are equally applicable to exploration areas located on the seabed.

[0068] Figure 2 It is a diagram that shows Figure 1 A cross-sectional view of seismic exploration area 101, according to one embodiment, shows the incident point of the earthquake source, the seismic data recording sensor (seismograph), the drilling location, the wellbore, various propagation rays, and different incident angles. More specifically, Figure 2 A cross-sectional view of the Earth portion above the seismic exploration area, indicated by reference numeral 201, shows different types of strata, indicated by reference numerals 102, 203, and 204, respectively. Although the seismic exploration area in this example is located on land, it should be understood that this is merely an example, and the methods and systems of the present invention are equally applicable to exploration areas located on the seabed. Figure 2 This demonstrates a common midpoint (CMP) type gather, where seismic data is ordered according to surface geometry, approximating a single reflection point within the Earth. These exploration seismic data are also referred to as traces, gathers, or image gathers. Figure 2 In the example, data from one or more seismic sources (explosion points) and detectors can be combined into a single image gather, or used in combination individually depending on the type of analysis required.

[0069] like Figure 2As shown, one or more source points (explosion points) represent seismic sources located at different incident points or stations (reference number 104) on the surface. One or more seismic sources are excited at these incident points 104. Seismic energy (or seismic sources) from multiple incident points 104 will be reflected at different stratigraphic interfaces. These reflected wave signals will be captured by multiple seismic data recording sensors 105, each positioned at a different offset location 210 and distributed relative to the well location 103. Because all incident points 104 and all seismic data recording sensors 105 are arranged at different offsets 210, the seismic data or seismic traces acquired during exploration, i.e., gathers or image gathers as referred to in the art, will be recorded at different incident angles 208. Incident points 104 generate downward-propagating rays 205 that propagate through the formation and are reflected back to the surface and captured by the seismic data recording sensors 105. In this case, well location 103 illustrates a drilled well connected to a wellbore 209. Various measurements can be performed along the wellbore 209 using techniques known in the art. The wellbore 209 can be used to acquire drilling logging data, including P-wave velocity, S-wave velocity, and density. Other measurements not previously taken in the exploration area can also be performed. Figure 2 The sensors shown are used to capture seismic data. This seismic data can be used to analyze amplitude, signal-to-noise ratio, move-out velocity, frequency components, phase, and other seismic properties related to geometric attributes such as incident angle 20°, offset 21°, and azimuth. These parameters are of great significance for the processing and imaging of seismic exploration data.

[0070] Figure 3 The diagram illustrates a cross-section of a seismic exploration area, including a wellbore and a set of downhole logging tools comprising one or more acoustic generators and one or more downhole logging data recording sensors, according to one embodiment of this disclosure. An acoustic generator is an example of a device for generating one or more acoustic waves (acoustic waves). An acoustic generator may also be referred to as a sound source because it generates or emits one or more acoustic waves, which are also known as seismic waves. The one or more downhole logging data recording sensors are examples of seismic data recording sensors (seismic detectors or seismic data recorders), and may be identical to seismic data recording sensor 105. In one embodiment of the invention, oil and gas extraction activities are suspended in order to generate seismic waves and record seismic wave reflection data passing through one or more formations in the seismic exploration area.

[0071] Figure 3A land-based oil drilling system 300 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 generators (acoustic sources) for generating one or more acoustic waves that are transmitted to one or more formations and generate reflected or reflected wave signals in said formations. Although this example shows one or more formations in a land-based exploration area, it should be understood that this is only an example, and the method and system can also be applied to exploration areas on the surface or bottom of water bodies 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, which includes reflected data generated in one or more formations by acoustic waves emitted by the one or more acoustic generators. The logging data is an example of seismic data. The logging data may include compressive wave velocity or P-wave velocity (Vp), shear wave velocity (Vs), and density, which is an indicator of porosity. The logging process that records well logging data can also be called sonic logging. Vehicle 325 can be connected to logging tool 315 to assist in the raising and lowering of logging tool 315 and to communicate with logging tool 315 to obtain logging data. Alternatively, in methods and systems used for exploration areas on the surface or bottom of water bodies such as the ocean, other equipment or systems can be used to assist in the raising and lowering of logging tool 315 and to communicate with it to obtain logging data.

[0072] Figure 4 This is a schematic diagram of a high-performance computer system according to an embodiment of the present disclosure, which receives (typically wirelessly) data from... Figure 1 and Figure 2 Earthquake data recording sensor 105 and / or Figure 3 Earthquake data recording sensors (in) Figure 3 Seismic wave data (also known as well logging data recording sensor). Figure 4 The high-performance computer system shown stores seismic data in at least one memory for subsequent processing and analysis by the computer-implemented methods and apparatus of one or more embodiments. The analyzed or processed seismic data can be accessed via a 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. The system computer 405 is connected to one or more storage devices 410 to store the seismic data in a database. The data transmission system can also directly wirelessly transmit seismic data to one or more storage devices 410 and store it in a database, which is accessed by the system computer 405. Wireless transmission is indicated by reference numeral 402. The one or more storage devices 410 may also store other computer software instructions or programs for performing the apparatus and methods described in the embodiments. The system computer 405 may be connected (e.g., wirelessly) to one or more output storage devices 420 to receive results obtained after computer processing flows or methods performed by the system computer 405. A personal computer 425 may be connected (e.g., wirelessly) to the output storage devices 420 and / or the system computer 405, enabling a user to input information or obtain results of computer processing methods performed by the system computer 405 through the user interface of the personal computer 425. The one or more output storage devices 420 may also store other computer software instructions or programs for implementing the methods and apparatus described in the embodiments.

[0073] The user interface of the personal computer 425 may include, for example, one or more of the following components: keyboard, mouse, joystick, button, switch, electronic pen or stylus, gesture recognition sensor (e.g., for recognizing user gestures, including movement of body parts), input device or voice recognition sensor (e.g., microphone for receiving voice commands), output device (e.g., speaker), trackball, remote control, portable device (e.g., cellular phone or smartphone), tablet computer, pedal or foot switch, virtual reality device, etc. The user interface may also include a haptic device for providing haptic feedback to the user. The user interface may also include, for example, a touchscreen. Additionally, the personal computer 425 may be a desktop computer, laptop computer, tablet computer, mobile phone, or any other personal computing system.

[0074] The processes, functions, methods, and / or computer software instructions or programs in the apparatus and methods described in this embodiment may be recorded, stored, or fixed in one or more non-transitory computer-readable media (computer-readable storage (recording) media) containing program instructions (computer-readable instructions) that are executed (performed or implemented by a computer) to cause one or more processors to execute the program instructions. The media may also contain data files, data structures, etc., together with or separately from the program instructions. The media and program instructions may be specially designed and constructed, or may be of conventional types well known and available to those skilled in the art of computer software. Examples of non-transitory computer-readable media include magnetic media such as hard disks, floppy disks, and magnetic tapes; optical media such as CD-ROMs and DVDs; magneto-optical media such as optical discs; and hardware devices specifically 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 (e.g., code generated by a compiler), and files containing high-level code that can be executed by a computer through an interpreter. These program instructions may be executed by one or more processors. The hardware device can be configured to act as one or more software modules, which are recorded, stored, or embedded in one or more non-transitory computer-readable media to perform the operations and methods described above, and vice versa. Furthermore, the non-transitory computer-readable media can be distributed across multiple computer systems connected via a network, and program instructions can be stored and executed in a decentralized manner. Additionally, the computer-readable media can also be embodied in at least one application-specific integrated circuit (ASIC) or field-programmable gate array (FPGA).

[0075] The one or more databases may include data sets and their supported data structures, which may be stored, for example, in one or more storage devices 410 and 420. For example, the storage devices 410 and 420 may be embodied as 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 combinations thereof. However, examples of the storage devices 410 and 420 are not limited to those described above, and the storage may also be implemented using various other devices and structures understood by those skilled in the art.

[0076] Figure 5 It is a flowchart illustrating, according to one embodiment, a seismic full-waveform inversion method independent of the source wavelet to generate a final velocity model, thereby improving the imaging of complex underground structures in the surveyed area. Figure 5 The seismic full waveform inversion method in the illustrated embodiments is applicable to the full waveform inversion of acoustic waves and elastic waves in isotropic or anisotropic media. Figure 5 The illustrated embodiments can be used to generate one or more high-resolution geological models, thereby enabling high-resolution imaging of lithological identification, fluid identification, and reservoir characterization of complex subsurface structures in the survey area. The survey area can be a subsurface structure beneath land or sea.

[0077] As mentioned above, full waveform inversion (FWI) is an optimization problem and a highly nonlinear process. Its success in practical applications depends on several factors, such as the choice of wave equation and the selection of source wavelet. Figure 5 The method presented aims to reduce reliance on the source wavelet during full waveform inversion in seismic exploration. FWI is a data-driven tool that automatically constructs velocity models by iteratively minimizing the discrepancies between observed and synthetic data. Generating synthetic data requires not only selecting appropriate wave equations and initial models but also the source wavelet, which is often unknown. Even small errors in the source wavelet can lead to significant biases in the inversion model, accumulating with depth. Therefore, accurate estimation of the source wavelet is one of the key factors for the success of FWI. Even with a wave equation that accurately describes the actual situation, inversion may still fail if the initial model and / or source wavelet are inaccurate. Figure 5 The method shown focuses on reducing or eliminating dependence on the selection of the source wavelet, provided that the chosen wave equation (acoustic or elastic) is sufficiently accurate and, with good initial model parameters, can accurately describe the propagation characteristics of seismic waves. This method significantly reduces dependence on the source wavelet by matching synthetic data with seismic data using a Wiener filter in each FWI iteration and by combining correlation and convolution operations to calculate the accompanying source.

[0078] refer to Figure 5 The earthquake full waveform inversion method shown uses earthquake data (observed earthquake data), denoted by reference number 500, which represents the earthquake data detected within the survey area. This data can be input into... Figure 5 The earthquake full waveform inversion method shown is illustrated. For example, as... Figure 2 The seismic data recording sensor 105 and / or shown Figure 3 The downhole logging data recording sensor in the logging tool 315 shown can detect seismic data 500 and transmit the seismic data 500 to... Figure 4 The high-performance computing system shown. For example... Figures 1 to 4As described above, the seismic data 500 detected in the survey area can be stored in one or more memories, such as one or more storage devices 410 and one or more output storage devices 420. Furthermore, as... Figure 5 As shown, the initial velocity model V0 (reference number 505) can also be input into the seismic full-waveform inversion method. This initial velocity model V0 505 can be a P-wave velocity model, an S-wave velocity model (elastic model), a density model, an anisotropic model (e.g., Thomsen anisotropic parameters epsilon (long migration effect) and delta (short migration effect)), or a combination of the above models. It should be understood that the initial velocity model can be a set of parameters. The initial velocity model V0 505 can be a velocity model pre-defined based on established scientific standards. The user can input the initial velocity model through a user interface such as a personal computer 425, or the model can be pre-stored in one or more memories, such as storage device 410 and / or output storage device 420. For simplicity, in an exemplary embodiment, the initial velocity model V0 505 is a P-wave velocity model for isotropic acoustic waves, or a combination of a P-wave velocity model for anisotropic acoustic waves and a Thomsen anisotropic parameter epsilon (long migration effect) and delta (short migration effect) model. However, Figure 5 The earthquake full waveform inversion method shown is also applicable to cases with other models, such as the S-wave velocity model as an elastic model.

[0079] like Figure 5 In the seismic full waveform inversion method shown, there is an initial model, such as the initial velocity model V0505, which is input into an FWI loop. In operation 510, the current velocity model is determined. Initially, the velocity model V... k This is the initial velocity model V0 505, which in this example is the initial P-wave velocity model. The variable k represents the number of iterations in the loop; therefore, initially k = 0. Figure 5 In each iteration shown (labeled as reference number 560), the velocity model V k The velocity model is updated in operation 510. For example, after the first iteration, the velocity model will become V1, where k = 1; after the second iteration, the velocity model will become V2, where k = 2. Iteration 560 continues until convergence is detected in operation 555. After convergence is detected in operation 555, the final velocity model (reference number 565) will be output, which corresponds to the most recent velocity model in operation 550. The final velocity model will be output and subsequently displayed in operation 570. Operations 550 and 555 will be described in more detail below.

[0080] refer to Figure 6 , Figure 6 This is a color image showing an example of a real P-wave velocity model in a seismic survey area profile. Figure 6 (Left side) and a preset initial P-wave velocity model example ( Figure 6 (On the right), this preset initial model corresponds to the input initial velocity model V0. As mentioned above, the input initial velocity model V0 is input to... Figure 5 The earthquake full-waveform inversion method shown includes an iterative loop for progressively approximating and improving the true velocity model. Although Figure 6 The example shown is a real speed model, but in reality, the real speed model cannot be known precisely. Figure 5 The full-waveform seismic inversion method in China generates an estimated (approximate) velocity model that is as close as possible to the actual velocity. Figure 6 The actual velocity model shown. Figure 6 In the diagram, the vertical axis represents the depth of the surveyed area, and the horizontal axis represents the offset distance (in kilometers, km), which is the horizontal distance from the origin (0 km). Figure 6 The color of the medium velocity model indicates the velocity value at a specific location in the survey area.

[0081] refer to Figure 5 Operation 515 in the earthquake full waveform inversion method shown, velocity model V k The data is input into operation 515 for forward modeling. Forward modeling of seismic data is a method based on geological information (in this case, the velocity model V). k This is a technique for generating synthetic seismic data. Forward modeling requires a source wavelet. This source wavelet can be a known source wavelet that is independent of the seismic wavelet actually used in current seismic surveys. Figure 7 An example of a source wavelet used in forward modeling operations according to this embodiment is shown, which may be referred to as the Ormsby wavelet. Figure 8 This demonstrates another example of a source wavelet, which can be called a Ricker wavelet. As shown in the figure, Figure 7 Ormsby wavelet in Figure 8 The Ricker wavelet in the model exhibits significant waveform differences. In the forward simulation, based on the velocity model V... k The source wavelet is used to generate synthetic seismic data. This synthetic seismic data is denoted by reference numeral 520, representing the synthetic data output from forward modeling operation 515.

[0082] The synthesized data 520 is input into the matched filter 525. For example, this synthesized data is generated based on the input initial velocity model V0. In addition, the forward simulation operation 515 generates and outputs a forward wavefield, denoted by reference numeral 535, which will be described in more detail below.

[0083] Referring to matched filter 525, the seismic data D detected in the exploration area obs Both the calculated synthetic data 500 and the calculated composite data 520 are input into the matched filter 525. An example of the matched filter 525 could be a Wiener filter, used to match the calculated synthetic data 520 with the seismic data D detected by the seismic data recording sensor. obs Matching is performed on the seismic data 500. The aforementioned seismic data recording sensor can be, for example, a seismic data recording sensor 105 deployed in the exploration area and a well logging data recording sensor deployed in the wellbore of the exploration area by the logging tool 315. Matching filter 525 matches the phase and amplitude of the synthetic data 520 with the phase and amplitude of the seismic data 500. Matching filter 525 is determined in each iteration 560 based on the calculated synthetic data 520, thereby ultimately updating the velocity model V. k 510. Convolve the seismic data 500 with the synthetic data 520 to produce matching synthetic data similar to the seismic data 500 detected in the exploration area.

[0084] As mentioned above, small errors in the source wavelet can eventually lead to large biases in the inversion model, and these biases accumulate with increasing depth below the Earth's surface. Therefore, in practical applications using the least squares mismatch function (Equation 2) for full waveform inversion (FWI), accurate source wavelet estimation is one of the key factors for successful FWI implementation. A typical estimation method is to extract the first arrival wave from the data using a time window and superimpose them to obtain the source wavelet. The source wavelet estimated in this way remains unchanged during FWI iterations. However, accurately estimating the source wavelet in practical industrial applications is very difficult due to several reasons, including poor repeatability of the source signal for each excitation (whether from an explosion or a sound source), uncertainties in the coupling between the source and the formation, and coupling problems between the detector and the formation. Therefore, much research has been dedicated to developing source-independent mismatch functions to overcome these problems. Figure 5 The matched filtering operation 525 shown, combined with other operations such as the accompanying source 530, can effectively solve these problems.

[0085] More specifically, the matched filter operation 525 transforms the least-squares mismatch function in equation (2) into equation (3) as shown below:

[0086]

[0087] Where * denotes the convolution operator, M s (t) is the matched filter function, which matches the synthetic data with the seismic data using the least squares Wiener filter operator.

[0088]

[0089] Where d obs,s,r (ω) and d syn,s,r (ω) represents the Fourier transform of the seismic data and the synthetic seismic data, respectively, and ∈ is a small-value regularization factor. Therefore, in each iteration, in order to minimize the mismatch function, in equation (3):

[0090] The formula for the matched filter 525 will be adjusted.

[0091] Referring to the accompanying source operation 530, the accompanying source operation 530 calculates one or more accompanying sources. The accompanying source operation 530 includes the following formula:

[0092]

[0093] in This represents the cross-correlation operator. In each iteration, the Wiener matched filter M... s (t) is used with the updated model parameters m k The corresponding time-domain synthesized data d syn,s,r (m k The accompanying source is determined by t). Operation 530 outputs an accompanying source, and the FWI gradient G can be calculated in operation 540. k The symbol * indicates convolution (e.g., M). s (t)*d syn,s,r (m;t) represents M s (t) and d syn,s,r (m;t) are convolved.

[0094] Since correlation and convolution are part of the accompanying source operation 530, as shown in equation (5), the FWI gradient G is calculated in operation 540. k At this time, the forward wavefield 535 generated by the forward simulation operation 515 and the accompanying source from the accompanying source operation 530 can be used. In the traditional seismic full-waveform inversion method, after the operation of matching the synthetic data 520 with the seismic data 500, an additional forward simulation operation must be performed using the updated source wavelet to provide the matched forward wavefield as input to the regular accompanying source operation 530 for FWI gradient G. k Operation 540. Therefore, the associated source operation 530 in equation (5) significantly improves computational efficiency and significantly reduces the use of computational resources in each iteration k. As mentioned above, in Figure 5The workflow eliminates the need for updated source wavelets, making the full-waveform seismic inversion method independent of the source wavelet. Matched filtering operation 525 and accompanying source operation 530 eliminate the need for updated source wavelets in the forward modeling operation 515 of the traditional FWI workflow.

[0095] The inversion model parameters generated by this method are independent of the selected wavelet, thus overcoming a major drawback of full waveform inversion (FWI). This method enables the construction of high-resolution geological models to improve imaging of complex subsurface structures in surveyed areas, thereby enhancing the accuracy of lithology identification, fluid identification, and reservoir characterization in seismic exploration.

[0096] Therefore, as described above, the forward simulation operation 515 outputs synthetic data 520. Furthermore, the forward simulation operation 515 also outputs a forward wavefield 535, which is a three-dimensional wavefield generated in each forward simulation time step. Since the forward wavefield 535 is generated only once by the forward simulation operation 515 in each iteration 560, instead of twice per iteration as in traditional seismic full-waveform inversion methods, therefore… Figure 5 The earthquake full-waveform inversion method shown significantly improves computational speed and reduces computational resources required for calculating the final velocity model. In traditional earthquake full-waveform inversion methods, the forward modeling operation needs to be performed twice, before and after the matched filtering operation, using two different source wavelets to generate an updated forward wavefield after the matched filtering operation. Therefore, traditional methods rely on the updated source wavelet. However... Figure 5 The embodiment shown illustrates a seismic full-waveform inversion method independent of the updated source wavelet of forward modeling operation 515. Figure 5 Another advantage of the earthquake full waveform inversion method shown is that the resolution of the generated offset image is significantly improved.

[0097] Based on the accompanying source output from operation 530 and the forward wavefield 535, the FWI gradient G is generated in operation 540. k G k It can be represented as G k (x), where k represents the number of iterations. Refer to the FWI gradient G in operation 540. k The calculation, operation 540, will backpropagate one or more associated seismic sources generated by formula (5), and use formula (6) to obtain the associated wave field, as shown below:

[0098]

[0099] in It is the adjoint operator of the forward operator F(m;x), f adj,s (t) is the associated seismic source, while u s(x,t) represents the adjoint wave field.

[0100] Using the forward wavefield 535 and the adjoint wavefield described above, a gradient can be obtained in each iteration using the following formula:

[0101]

[0102] Where F(m,x) represents the forward modeling operator of the wave equation; w s (x,t) represents the forward wave field; u s (x,t) represents the inverse adjoint wave field in equation (6). In operation 540, the gradients in equation (7) corresponding to different sources (explosions or acoustic generators) are summed and superimposed to obtain a unified gradient. Then, this gradient G... k This will be used to update the model parameters m via an inversion method to minimize the mismatch function, as described below regarding operation 545.

[0103] Operation 545 is used to determine: (1) by how much the update rate should be increased or decreased when moving closer to the true velocity model; and (2) whether the update rate should be increased or decreased in order to better approximate the true velocity model.

[0104] More specifically, once the gradient is calculated and output by operation 540, operation 545 calculates the step size A according to the inversion method. k and search direction P k For example, from the initial estimate of the subsurface parameter m k＝0 Initially, in the (k+1)th iteration, the model update formula is as follows:

[0105] m k+1 =m k +A k P k k = 0, 1, ..., (8)

[0106] Where P k This is referred to as the search direction, and A k This is the step size. The search direction P is obtained through the gradient (Equation 7), with the aim of minimizing the mismatch function; this process is called inversion. The simplest inversion method is called the steepest descent algorithm, where the search direction P... k The result is given by the negative gradient of the mismatch function at the k-th iteration:

[0107] P k =-G k (9)

[0108] Based on step size A k Operation 545, step size A k and search direction P kOperation 545 receives gradients and determines the magnitude and direction (increase or decrease) of velocity updates toward the true velocity model. Therefore, the step size A... k and search direction P k Operation 545 pairs of FWI gradients G k The output of operation 540 is scaled to provide the optimal increment for updating the velocity model, avoiding excessively large increases or decreases in velocity updates (which would significantly increase the number of iterations) or excessively small increases or decreases (which would also significantly increase the number of iterations). Then, operation 550 updates the velocity model V based on the following formula. k :

[0109] V k+1 =V k -A k G k (10)

[0110] Based on convergence operation 555, the seismic full-waveform inversion method determines whether additional iterations are needed. If additional iterations are needed, the value of k is increased in operation 560, and the updated velocity model V is updated in operation 520. k+1 Set as velocity model V k In one embodiment, operation 550 can set a maximum number of iterations k. The maximum number of iterations can be between 10 and 40. Once the number of iterations k equals the predetermined maximum number of iterations, the seismic full waveform inversion method is considered to have converged, and the final velocity model V is output. F This model is consistent with the last updated velocity model V determined in operation 550. k+1 Same. Final velocity model V F This can be represented as 565. In one embodiment, operation 550 can also set an additional threshold for the difference between the current update rate in operation 550 and the update rate of the previous iteration in operation 550. If the difference is greater than (or greater than or equal to) the threshold, then the next iteration 560 can continue if the maximum number of iterations has not been exceeded. However, if the difference is less than (or less than or equal to) the threshold, the seismic full waveform inversion method is considered to have converged, and the final velocity model V is output. F This model is consistent with the last updated velocity model V determined in operation 550. k+1 Same. Final velocity model V F This can be represented as 565. The final velocity model V F The image can be displayed in operation 570. For example, personal computer system 425 can display the final velocity model V. F The image.

[0111] The following is an example of synthetic data and field data inversion, which can be found by referring to... Figures 7 to 14 It showed Figure 5 The feasibility and stability of the method discussed herein are described, which is an improved source wavelet-independent full-waveform seismic inversion. Although the embodiments described below involve source wavelet-independent mismatch functions using the finite-difference approximation of the scalar acoustic equations, the embodiments also include source wavelet-independent mismatch functions using the vector wave equations and elastic wave equations in isotropic and anisotropic media.

[0112] For example, Figure 9 It is a diagram, reference number 900, showing the use of forward simulation operations. Figure 7 and Figure 8 The convergence comparison of different source wavelets is shown. Traditional seismic full-waveform inversion methods use a least-squares mismatch function. However, the seismic full-waveform inversion methods in these embodiments use a mismatch function independent of the source wavelet, thus enabling these embodiments to achieve excellent performance without depending on the source wavelet. Figure 7 The source wavelet shown is applied to the traditional seismic full-waveform inversion method, and its convergence corresponds to the curve shown in reference number 910. If... Figure 8 The source wavelet shown is applied to the traditional seismic full-waveform inversion method, and its convergence corresponds to the curve shown in reference number 920. If, in a certain embodiment of the seismic full-waveform inversion method, it is used in forward simulation operation 515... Figure 7 The convergence of the source wavelet shown is illustrated by reference number 930. If using... Figure 8 The convergence of the source wavelet shown is given by reference number 940.

[0113] like Figure 9 As shown, the number of iterations required to achieve convergence in traditional seismic full-waveform inversion methods varies considerably, as illustrated by curves in reference numbers 910 and 920. This indicates that convergence is highly dependent on the selected source wavelet, and the error rates differ significantly. However, in the embodiments of this disclosure, after employing a mismatch function independent of the source wavelet, the operations shown in reference numbers 930 and 940 both exhibit good convergence, are almost unaffected by the selected source wavelet, and have low error rates. This fully demonstrates that the seismic full-waveform inversion method of the embodiments of this disclosure is independent of the source wavelet. As mentioned above, the seismic full-waveform inversion method independent of the source wavelet greatly improves computational efficiency and reduces the computational resources required to perform seismic full-waveform inversion. Furthermore, when the selected source wavelet in the forward simulation operation is... Figure 8 When the source wavelet is shown, its error under the traditional full-waveform earthquake inversion method (curve shown in reference number 940) is much greater than the error when the same source wavelet is used in the embodiments of this disclosure (curve shown in reference number 920).

[0114] Figure 10This is a color image showing a comparison between the actual P-wave velocity model and the P-wave velocity model obtained using the traditional full-waveform inversion method in an exploration area profile. The left side of the image shows the actual P-wave velocity model, and the right side shows the model obtained using the traditional full-waveform inversion method (using...). Figure 7 The image shows the P-wave velocity model generated by the Ormsby subwavelet. Figure 10 It demonstrates that, under the ideal condition where the source wavelet is error-free, for Figure 6 Given an initial P-wave velocity model, the best final P-wave velocity result obtainable by the traditional full-waveform inversion method. Figure 10 In the graph, the vertical axis represents the depth of the exploration area, and the horizontal axis represents the offset (km), i.e., the horizontal distance from the origin (0 km). The colors in the graph represent the velocity values ​​at each location within the exploration area. As mentioned earlier, even a small error in the source wavelet can lead to significant deviations in the inversion model, and these deviations accumulate with increasing depth. Therefore, Figure 10 This only represents the "best case" of the traditional full waveform inversion method under ideal conditions; however, it is difficult to obtain a completely error-free source waveform in reality, so this method often cannot obtain an accurate underground velocity model.

[0115] Figure 11 This is a color image showing a comparison between the actual P-wave velocity model in a vertical structural profile of an exploration area and the P-wave velocity model obtained using the seismic full-waveform inversion method and source wavelet according to embodiments of the present invention. The left side of the image shows the actual P-wave velocity model, and the right side shows the P-wave velocity model obtained using the seismic full-waveform inversion method of this disclosure. Figure 5 The final velocity model V shown F (labeled 565) and using Ormsby wavelet ( Figure 7 The generated P-wave velocity model image is shown below. Figure 11 In the graph, the vertical axis represents the depth of the exploration area, and the horizontal axis represents the offset distance (km), which is the horizontal distance from the origin (0km). The colors in the graph represent the velocity values ​​at various locations within the exploration area.

[0116] Figure 12 This is a color image showing a comparison between the actual P-wave velocity model in a vertical structural profile of an exploration area and the P-wave velocity models obtained using the seismic full-waveform inversion method and different source wavelets according to embodiments of this disclosure. The left side of the image shows the actual P-wave velocity model, and the right side shows the P-wave velocity model obtained using the seismic full-waveform inversion method according to embodiments of this disclosure. Figure 5 The final velocity model V, marked as 565 in the middle. F ) and use Ricker wavelet ( Figure 8 The generated P-wave velocity model image is shown below. Figure 12In the graph, the vertical axis represents the depth of the exploration area, and the horizontal axis represents the offset (km), i.e., the horizontal distance from the origin (0 km). The colors in the graph represent the velocity values ​​at various locations. This graph shows that even when using... Figure 11 Despite different source wavelets (Ricker wavelets), the embodiments disclosed herein can still provide high-quality velocity inversion results. The inversion model is highly consistent with the real model, further verifying the stability and accuracy of the method under different source wavelet conditions.

[0117] like Figure 11 and Figure 12 As shown, the seismic full-waveform inversion method of this invention can generate a high-quality final velocity model V. f ( Figure 5 The value is marked as 565 and can be displayed on devices such as a personal computer system 425. Regardless of whether Ormsby wavelet ( Figure 7 ) or Ricker wavelet ( Figure 8 When the generated final velocity images are obtained, they all closely match the actual P-wave velocity model, indicating that the seismic full-waveform inversion method of this embodiment can generate excellent inversion images without relying on the source wavelet. Furthermore, Figure 11 and Figure 12 The final velocity model V in f Image quality can be compared with Figure 10 The velocity images obtained using traditional seismic full-waveform inversion methods (assuming no errors in the source wavelet) are comparable. This further verifies the superior performance of the embodiments of this disclosure in achieving high-quality inversion imaging even when there are uncertainties or errors in the source wavelet.

[0118] Figure 13 These are color images of the updated velocity model obtained by processing actual observation data using the traditional full-waveform inversion (FWI) method, and a reverse-time migration (RTM) image generated based on this updated velocity model. The updated velocity model generated based on the traditional full-waveform inversion method... Figure 13 The left side shows a reverse-time offset image, which exhibits a wavy structure. Similar to Figure 13, the vertical axis represents depth, and the horizontal axis represents the Common Midpoint (CMP) location. Figure 13 As shown, the velocity model perturbation (color image) obtained in the traditional full-waveform inversion method cannot effectively resolve geological structures. Figure 13 This is also reflected in the reverse-time migration images. Therefore, it is evident that traditional full-waveform seismic inversion methods have certain limitations in practical applications, making it difficult to accurately reconstruct complex underground geological features.

[0119] Figure 14The seismic full waveform inversion method (FWI) according to the embodiments of this disclosure is in conjunction with... Figure 13 Color plot of the updated velocity model generated under the same actual observation data and the same initial P-wave velocity model conditions, and a reverse time migration (RTM) image generated based on the updated velocity model. Figure 14 The vertical axis represents depth, and the horizontal axis represents the Common Midpoint (CMP). For example... Figure 14 As shown, from Figure 5 The FWI model perturbation of the embodiment shown (in) Figure 14 The color map (shown in the image) clearly reveals the underground geological structure. Figure 14 In the middle, reverse time migration (RTM) imaging shows that the subsurface layers are significantly flatter, which is more consistent with the geological interpretation of the exploration area, compared to Figure 13 The imaging results in the middle showed considerable fluctuations. This improvement is due to... Figure 5 The updated velocity model generated by the seismic waveform inversion method shown is illustrated. Figure 13 In comparison, this updated model significantly reduces structural fluctuations. The comparison uses the same source wavelet, demonstrating that the velocity model in this embodiment outperforms the velocity models used in conventional methods. This illustrates that the seismic full-waveform inversion method of this disclosure has higher imaging quality and stronger geological analysis capabilities in real-world data applications.

[0120] While embodiments of this disclosure have been shown and described, modifications can be made by those skilled in the art without departing from the spirit or teachings of the invention. The embodiments described herein are illustrative only and are not restrictive. Many variations and modifications of methods, systems, and apparatus are possible and all fall within the scope of this invention. Therefore, the scope of protection of this invention is not limited to the embodiments described herein, but is limited only by the claims. The scope of the claims should include all equivalents of its subject matter.

Claims

1. A method for performing seismic full-waveform inversion to generate a velocity model of subsurface strata in a survey area, the method comprising the following operations: (a) Deploy a first plurality of seismic data recording sensors at different locations in the survey area, and / or deploy logging tools containing a second plurality of seismic data recording sensors in the wellbore of the survey area; (b) Energy is emitted at the point of incidence in the survey area to generate seismic waves that propagate through the underground strata; (c) Using the first plurality of seismic data recording sensors and / or the second plurality of seismic data recording sensors to sense seismic waves and record seismic data based on the seismic waves; (d) Transmitting the seismic data from the first plurality of seismic data recording sensors and / or the second plurality of seismic data recording sensors to a computer system including one or more storage devices, and storing the seismic data in the one or more storage devices; (e) Storing the source wavelet in one or more of the storage devices; (f) Using the earthquake source wavelet and the current velocity model, the computer system performs forward modeling operations; (g) Using forward modeling operations, an updated velocity model is generated by a computer system for use in seismic full waveform inversion; (h) Perform operations (f) and (g) until convergence; (i) Upon convergence, the updated velocity is output as the final velocity model to the display; as well as (j) Display a high-resolution image of the final velocity model on the computer system's monitor. The forward modeling operation generates synthetic data and a forward wavefield. Operation (g) further includes matching seismic data with synthetic data and using a mismatch function to assess the mismatch between the seismic data and the synthetic data, wherein the mismatch function includes: ， Where * denotes the convolution operator, The vector represents the Earth model, and t represents time. As a matched filter function, it uses the least-squares Wiener filter operator to match the synthetic data with the seismic data. ， in, )and ) represent the Fourier transforms of seismic data and synthetic data, respectively.

2. The method according to claim 1, wherein, Earthquake full waveform inversion is one of the two types of acoustic wave inversion and elastic wave inversion, and the underground strata in the survey area are either isotropic or anisotropic.

3. The method according to claim 1, wherein, The one or more storage devices store the initial speed model.

4. The method according to claim 1, wherein, Operation (g) further includes generating associated seismic sources based on the mismatch between seismic data and synthetic data.

5. The method according to claim 4, wherein, The accompanying hypocenter is generated according to the following equation: ， in Represents the cross-correlation operator.

6. The method according to claim 4, wherein, Operation (g) further includes: determining the inversion gradient based on the associated seismic source.

7. The method according to claim 6, wherein, Operation (g) further includes: determining the step size and step direction for updating the velocity model; and updating the velocity model.

8. The method according to claim 1, wherein, The source wavelet is the Ormsby wavelet and the Ricker wavelet.

9. A system for performing seismic full-waveform inversion to generate a velocity model of subsurface strata in a surveyed area, the system comprising: The first plurality of seismic data recording sensors are arranged at different locations in the exploration area, and / or the logging tool including the second plurality of seismic data recording sensors is arranged in the wellbore in the exploration area; as well as Explosive devices were placed at each incident point in the survey area to generate seismic waves that propagated through the underground strata. Specifically, the first plurality of seismic data recording sensors and / or the second plurality of seismic data recording sensors sense seismic waves and record seismic data based on the seismic waves, and transmit the seismic data to a computer system including one or more storage devices and at least one processor, wherein the one or more storage devices store the transmitted seismic data, source wavelet, and instructions, and The one or more processors execute instructions stored in one or more storage devices to achieve the following: (f) Perform forward modeling using the source wavelet; (g) Use forward modeling to generate an updated velocity model for full-waveform seismic inversion; (h) Perform operations (f) and (g) until convergence; (i) Upon convergence, the updated velocity model is output as the final velocity model to the display; and (j) Display a high-resolution image of the final velocity model on the computer system's monitor. The forward modeling operation generates synthetic data and a forward wavefield. Operation (g) further includes matching seismic data with synthetic data and using a mismatch function to assess the mismatch between the seismic data and the synthetic data, wherein the mismatch function includes: ， Where * denotes the convolution operator, The vector represents the Earth model, and t represents time. As a matched filter function, it uses the least-squares Wiener filter operator to match the synthetic data with the seismic data. ， in, )and ) represent the Fourier transforms of seismic data and synthetic data, respectively.

10. The system according to claim 9, wherein, Earthquake full waveform inversion is one of the two types of acoustic wave inversion and elastic wave inversion, and the underground strata in the survey area are either isotropic or anisotropic.

11. The system according to claim 9, wherein, One or more storage devices store the initial speed model.

12. The system of claim 9, wherein operation (b) further comprises generating an accompanying seismic source based on a mismatch between the seismic data and the synthetic data.

13. The system according to claim 12, wherein, The accompanying hypocenter is generated according to the following equation: ， in Represents the cross-correlation operator.

14. The system according to claim 12, wherein, Operation (b) further includes: determining the inversion gradient based on the associated seismic source.

15. The system according to claim 14, wherein, Operation (b) further includes: determining the step size and step direction for updating the velocity model; and updating the velocity model.

16. The system of claim 9, wherein the source wavelet is an Ormsby wavelet or a Ricker wavelet.

Images