System and method for analyzing model precision and updating model using delay time information

Abstract

The earthquake processing method comprises the following steps: receiving an earthquake waveform which is acquired from an earthquake receiver and contains underground area structure information; generating a back-propagating wavefield by back-propagating seismic data over time based on the initial model system; generating a forward propagation wavefield by forward propagation over time using a given wavelet or processed seismic data; generating a delay time common imaging point gather (DTCIGs) based on the back propagation wave field and the forward propagation wave field; extracting model precision information from the DTCIGs for analysis and / or quality control; determining a model error from the model precision information and / or DTCIGs for analysis and / or quality control; and adjusting the initial model system based on the model error to generate an adjusted / updated model.

Description

Cross-references to related applications

[0001] This application claims priority to U.S. Provisional Patent Application No. 63 / 492,827, filed March 29, 2023, the entire contents of which are incorporated herein by reference. Declaration of rights to inventions resulting from federally sponsored research and development

[0002] not applicable Technical Field

[0003] This disclosure generally relates to methods and systems for seismic exploration, and more specifically, to types of full-waveform inversion, which use time delay information generated by seismic imaging to assess the accuracy of seismic parameter models and improve models and seismic images of the explored subsurface strata. Background Technology

[0004] A widely used technique for finding energy resources is seismic exploration of subsurface geophysical structures. The seismic exploration process involves sending seismic waves (i.e., sound waves) from a seismic source into the subsurface region, collecting the seismic signals as recorded waves propagating through the Earth, including reflections, refractions, and transmissions. These signals carry subsurface geological information, and the data is analyzed to generate a profile (image) of the geophysical structure—the layers of the subsurface being explored. This final step in generating the subsurface image is called seismic imaging. This type of seismic exploration can be used both on land and on the seabed.

[0005] Figure 1 A simplified schematic diagram of the marine seismic acquisition process is shown when a reservoir 100 is present in the surveyed stratum 110. A vessel 130 tows a seismic source 140. Numerous receivers 150 are located either in the water (towed cable acquisition), on the seabed (seabed node / cable acquisition), or along well tracks within the Earth's interior (vertical seismic profile acquisition). The seismic excitation generated by the source penetrates the seabed and passes through the surveyed stratum 110, where it is then detected and recorded by the receivers, thereby generating seismic data. Although we use a marine data acquisition setup to illustrate seismic data acquisition, embodiments of this disclosure are applicable to all types of seismic data, including any type of marine and terrestrial seismic surveys acquired by various types of sources (air guns, vibrators, explosives, etc.) and receivers (seismographs, hydrophones, distributed acoustic sensors, etc.).

[0006] The recorded data is then processed to generate accurate subsurface models of the surveyed stratigraphic rock properties, or, in a suboptimal case, high-quality images of the surveyed stratigraphic structure. Therefore, using these subsurface rock properties or stratigraphic models or images, geophysicists can identify optimal locations of potential reservoirs, and petroleum engineers can design and develop methods for extracting these natural resources from reservoir sediments based on geomechanical analysis.

[0007] Subsurface models are a crucial component of seismic processing and imaging. As the name suggests, a model is a visual representation of the physical or mathematical properties of different locations within the subsurface. These properties may include, but are not limited to: P-wave velocity models, S-wave velocity models, anisotropic models, elastic parameters, Q-models (used to describe subsurface attenuation), density models, etc. Note that "subsurface" can refer to land areas or the area below the seabed (but can also include the different velocities of sound waves passing through different seawater layers). Using these models, the travel time, amplitude, phase, etc., of seismic data can be described or simulated at all subsurface locations within the model.

[0008] Figure 2 illustrates the use Figure 3a The Sigsbee2a formation P-wave velocity model (Paffenholz, J., B. McLain, J. Zaske, and P. Keliher, 2002, “Subsalt multiple attenuation and imaging: Observations from the Sigsbee2B synthetic dataset”: 72nd Annual International Meeting, SEG, Expanded Abstracts, 2122-2125) is a simplified synthetic towed seismic data acquired from survey strata. This seismic data was obtained through towed acquisition, and one co-source seismic record shows... Figure 2a In the middle, a shared channel seismic record shows Figure 2b In the middle. Ideally, if the data collected in Figure 2 could be used to recover... Figure 3a The formation velocity model allows for easy and direct interpretation of the central black main salt body (high velocity), the shallow gray gradient seafloor and scattered cavitation (contrasting with the surrounding velocities), and the bottom black horizontal mudstone strata (high velocity). Alternatively, at the end of the processing phase, it is possible to recover... Figure 3b The (reflectance) image. In this case, using this image and combining it with some complex rock property inversion processes, it is still possible to recover... Figure 3a The velocity model in the middle achieves the ideal benefits of seismic survey.

[0009] However, it is well known that many different geological factors can make creating accurate model descriptions difficult. For example, some subsurface areas have significant complex characteristics, such as strong variations in velocity or anisotropic parameters, or complex geological structures, such as salt bodies and basalt structures, heavily fractured zones, anisotropic environments caused by sedimentation or fracturing (anisotropic environment refers to a seismic wave propagation speed that depends on whether it propagates along or through rock strata), thrust faults, shallow gas, etc. Therefore, the treatment of reflected, refracted, and transmitted sound waves can become extremely complex.

[0010] Figure 4a The Sigsbee2a P-wave migration velocity model (Paffenholz, J., B. McLain, J. Zaske, and P. Keliher, 2002, “Subsalt multiple attenuation and imaging: Observations from the Sigsbee2B synthetic dataset”: 72nd Annual International Meeting, SEG, Expanded Abstracts, 2122-2125) is presented. This model represents the velocity profile used for imaging. Figure 3a The basic geological background of the accurate Sigsbee2a formation P-wave velocity model. Using the data acquired in Figure 2, the migration velocity model can be adjusted and established as shown in Figure 2. Figure 4b The image shown. If there is no image like... Figure 3a Such a precise speed model Figure 4b Image quality and Figure 3b The ratio will decrease significantly compared to the previous period.

[0011] Another simplified example of a model system is the 2007 BP Tilt Lateral Isotropic (TTI) baseline model created by Hemang Shah and BP Exploration Operations Ltd., which is illustrated in Figures 5, 6 and 7. Figure 5a and 5b The co-source and co-receiver seismic datasets synthesized using typical rock properties from the Gulf of Mexico are presented respectively. Figure 6a , 6b Figures 7a and 7b respectively show the P-wave velocity, anisotropic Thomsen parameter δ model, anisotropic Thomsen parameter ε model, and anisotropic polarization angle θ of the surveyed strata. Using the data recorded in Figure 5 and the TTI model system in Figures 6 and 7, we can obtain the following through migration: Figure 8 Underground images.

[0012] Full Waveform Inversion (FWI) has long been an important method for establishing velocity models for seismic imaging (see Tarantola, A., 1984, “Inversion of Seismic Reflection Data in the Acoustic Approximation”: Geophysics, 49, 1259-1266; and Virieux, J., and S. Opera, 2009, “An Overview of Full Waveform Inversion in Exploration Geophysics”, Geophysics, 74(6), WCC127-WCC152). Classical FWI inverts velocity models by iteratively minimizing the differences between recorded field data and simulated model data. However, two main problems can cause FWI to fail in production: period jumps and amplitude differences between field and model data, especially when there are sharp contrasts and large geological bodies in the model. When there are large phase or travel time differences between field and model data, the oscillation patterns of the seismic signal can introduce local minima. FWI may attempt to fit the simulated model data to erroneous periods in the observed field data. Furthermore, the amplitude distribution of field observation data is often very different from that of simulation data. If this difference is not properly handled, FWI will interpret it as a velocity error, leading to the inversion of an incorrect velocity model.

[0013] Since the linear relationship between travel time differences and velocity errors in seismic signals is stronger than that between waveform differences, and any inconsistencies in amplitude differences are not important, effectively extracting travel time information from field data is the core of velocity model building and FWI applications. To achieve this goal, researchers have made many influential advances, such as adaptive waveform inversion (Warner and Guasch, 2014, "Adaptive waveform inversion": 84th Annual International Meeting, SEG, Expanded Abstracts, 1089-1093), dynamic warping FWI (Ma and Hale, 2013, "Wave-equation reflection traveltime inversion with dynamic warping and hybrid waveform inversion": 83rd Annual International Meeting, SEG, Expanded Abstracts, 871-876), time-delay FWI (Luo and Schuster, 1991, "Wave-equation traveltime inversion": Geophysics, 56, 645-653), and optimal transport FWI (Yang et al., 2018, "Application of optimal transport and the quadratic Wasserstein metric to full-waveform inversion": Geophysics, 83, 43-62).

[0014] All the methods described above are applied only in the traditional data domain, and their success relies heavily on extracting the travel time mismatch between the recorded field data and the simulated data. However, extracting travel time mismatch in the data domain is not easy. When the geology is complex, the waveforms are also complex, and multiple reflection, refraction, or transmission events may overlap during recording. On the other hand, if the initial velocity is simple and far from the correct velocity, the simulated composite waveform will appear much simpler, with significant differences in travel time and amplitude. Therefore, it is not obvious how to correctly measure the travel time mismatch by comparing the two datasets. In conclusion, when the geological survey is complex or FWI starts from a poor initial velocity model, there is no guarantee that FWI based on data domain travel time mismatch can overcome the period jump problem.

[0015] Given the inherent limitations of existing methods, innovative systems and approaches are needed to provide a practical and effective framework. These systems and methods will ensure the inversion of subsurface velocity distributions from field data through measurements in the image domain, unaffected by the periodic skip problem, and will accurately and practically evaluate the precision of velocity models. This document introduces an enhanced system and approach that cleverly achieves both of these objectives. Summary of the Invention

[0016] One objective of this embodiment is to substantially address the period jumps and amplitude discrepancies that can cause classical FWI to fail to converge to a reasonably good geological model system. Here, the model system describing the geophysical properties of subsurface locations may include, but is not limited to: P-wave velocity models, S-wave velocity models, anisotropic models, elastic parameters, Q-models (used to describe subsurface attenuation), density models, etc.

[0017] Therefore, a general aspect of this embodiment is to provide a system and method for improving the application of FWI in seismic surveys to achieve better model resolution, thereby eliminating or minimizing problems of the type previously described.

[0018] Some implementations use enhanced FWI methods to determine the model of the subsurface strata being surveyed. One approach uses time delay information not only to measure model accuracy but also to provide feedback to the FWI for model updates. This approach addresses the period jumps and amplitude discrepancies observed when processing seismic data acquired for surveying strata using classical FWI.

[0019] According to one embodiment, a seismic exploration method is provided. The method includes acquiring an initial model of the subsurface strata and recorded field data, and generating Delayed Time Domain Common Imaging Gathers (DTCIGs) based on the initial model. The method also includes extracting delayed time migration values ​​at subsurface locations. An example of this extraction process can be defined by selecting imaging points with maximum amplitude along the delayed time direction. The method then uses the delayed time migration information to quantify, analyze, or quality control velocity accuracy. The method then updates the model using an optimization scheme (such as FWI), the cost function of which is minimizing the total delayed time migration of the subsurface locations. The updated model is used to acquire images of the subsurface strata that can be used to locate natural resources.

[0020] According to another embodiment, a seismic data processing apparatus is provided, comprising an interface and a data processing unit. The interface is configured to acquire an initial model and recorded data for surveying subsurface strata. The data processing unit is configured to generate DTCIGs based on the initial model, determine time delays on the DTCIGs at subsurface locations, quantify, analyze, or quality control the model accuracy using the time delay information, and update the model using an optimization scheme (such as FWI) whose cost function is to minimize the total time delay of the subsurface locations. Here, the updated model is also used to acquire surveyed subsurface strata images that can be used to locate natural resources.

[0021] According to yet another embodiment, a non-transitory computer-readable recording medium is provided, which stores executable code that, when executed by a computer, causes the computer to perform a seismic exploration method. The method includes acquiring an initial model and recorded data for surveying subsurface strata, and generating DTCIGs based on the initial model. The method also includes determining time delays on the DTCIGs at subsurface locations. The method then uses the time delay information to quantify, analyze, or quality control velocity accuracy. The method then updates the model using an optimization scheme (such as FWI) whose cost function is to minimize the total time delay at subsurface locations.

[0022] To locate natural resources, the acquired images can be further analyzed using geoscience techniques, such as seismic interpretation or reservoir characterization. These techniques can identify areas with favorable geological conditions, such as oil and gas reservoirs, mineral deposits, or other resources of interest. By combining these techniques with updated models, the location of natural resources in subsurface strata can be predicted more accurately. Detailed steps supporting these implementations are specified in the following literature: 1. Brown, AR, 2011, “Interpretation of Three-Dimensional Seismic Data”: Society of Exploration Geophysicists and American Association of Petroleum Geologists.

[92] 2. Chopra, S., & Marfurt, KJ, 2007, “Seismic Attributes for Prospect Identification and Reservoir Characterization”: Society of Exploration Geophysicists.

[93] 3. U.S. Patent No. 8,543,623, entitled “Method for generating a subsurface model and method for hydrocarbon prospecting”, filed May 12, 2010. These references are all incorporated herein by reference. Attached Figure Description

[0023] To gain a more complete understanding of the present invention, reference is now made to the following description taken in conjunction with the accompanying drawings, wherein:

[0024] Figure 1 This is a schematic diagram of traditional seismic data acquisition.

[0025] Figure 2 illustrates the use of Figure 3a The Sigsbee2a formation P-wave velocity model shown is based on synthetic seismic data obtained from explored formations. Typical common shot gather seismic records and typical common gather seismic records are shown in [the diagram / image / data]. Figure 2a and Figure 2b middle.

[0026] Figure 3a The P-wave velocity model of the Sigsbee2a formation is shown. Figure 3b The reflectance model of the Sigsbee2a formation is shown.

[0027] Figure 4a The Sigsbee2a offset P-wave velocity model is shown. Figure 4b This shows the use of seismic data from Figure 2 and Figure 3a The offset speed is obtained by creating a Sigsbee2a offset image using the reverse time offset tool.

[0028] Figure 5a and Figure 5b Seismic datasets acquired at a common source and a common receiver are shown separately. This dataset comes from the Tilt Lateral Anisotropy (TTI) benchmark established by Hemang Shah and BP Exploration Operations in 2007.

[0029] Figure 6a and Figure 6b The P-wave velocity model and the anisotropic Thomsen parameter δ model in the 2007 BP tilted transverse anisotropy (TTI) standard are shown respectively.

[0030] Figure 7a and Figure 7b The anisotropic Thomsen parameter ε model and the anisotropic polarization angle θ of the explored strata are shown in the 2007 BP tilted transverse anisotropy (TTI) benchmark.

[0031] Figure 8 This illustrates the use of seismic data from Figure 5 and... Figure 6a , 6b TTI velocity models 7a and 7b, offset images created using the reverse time offset tool.

[0032] Figure 9 The diagram illustrates a method for measuring, analyzing, or quality controlling model accuracy using DTCIGs according to one embodiment.

[0033] Figure 10a An accurate P-wave velocity model is shown. Figure 10b This is shown in system 902, utilizing Figure 10a Typical common shot set seismic data for velocity simulation.

[0034] Figure 11a An initial P-wave velocity model is shown for initiating velocity analysis or model update processes, as in system 904. Figure 11b Typical sub-wavelengths are shown, as in system 906.

[0035] Figure 12a This is shown in system 908. Figure 11a A snapshot of the reverse propagating wave field in the strata being explored, propagating backward over time. Figure 12b As shown in system 910, Figure 11a A snapshot of the forward propagating wavefield in the explored strata over time.

[0036] Figure 13a and Figure 13b The diagram shows the original and processed wavefields obtained from the initial model system at a typical spatial location in the explored formation, under System 912. The white line represents the maximum amplitude selected at all depths along the time delay direction. The dashed red line represents zero-delay time lag.

[0037] Figure 14a and Figure 14b The figures show the spatial locations of the explored strata, which are the same as those in Figure 13, through... Figure 10a The original wavefield and the processed wavefield of the accurate model system are shown in the DTCIG. The white line represents the maximum amplitude selected at all depths along the time delay direction. The dashed red line represents zero-delay time lag.

[0038] Figure 15 The velocity accuracy / error attribute extracted from all locations of the explored formation using formula (6) based on DTCIG in Figure 13 under system 914 is shown.

[0039] Figure 16a and Figure 16b The SEAM accurate P-wave velocity model is shown in slices at lateral line 875 and cross line 1001, respectively.

[0040] Figure 17a and Figure 17b The slices of the SEAM exact anisotropic Thomsen parameter δ model are shown at lateral line 875 and cross line 1001, respectively.

[0041] Figure 18a and Figure 18b The slices of the SEAM exact anisotropic Thomsen parameter ε model are shown at lateral line 875 and cross line 1001, respectively.

[0042] Figure 19a and Figure 19b The slices of the SEAM exact anisotropic Thomsen parameter dip model at lateral line 875 and cross line 1001 are shown respectively.

[0043] Figure 20a and Figure 20b The following are slices of the SEAM precise anisotropic Thomsen parameter azimuth model at lateral line 875 and cross line 1001.

[0044] Figure 21a and Figure 21b The SEAM accurate density model is shown in slices at lateral line 875 and cross line 1001, respectively.

[0045] Figure 22a The typical subwavelet is shown in system 906. Figure 22bThe following shows typical common receiver point seismic data obtained using velocity simulations in Figures 16, 17, 18, 19, 20, and 21 under system 902.

[0046] Figure 23a and Figure 23b The initial P-wave velocity models at lateral line 875 and cross line 1001 are shown in system 904, respectively, for initiating velocity analysis or model update processes.

[0047] Figure 24a and Figure 24b The differences between the initial P-wave velocity model in Figure 23 and the accurate P-wave velocity model in Figure 16 are shown at the lateral line 875 and the cross line 1001, respectively.

[0048] Figure 25a Figure 23 shows a snapshot of the backpropagating wave field in the explored formation as it propagates backward over time, under system 908. Figure 25b Figure 23 shows a snapshot of the forward propagating wave field in the explored formation over time, under system 910.

[0049] Figure 26a , Figure 26b , Figure 26c and Figure 26d The DTCIG obtained by the original wavefield of the initial model system in Figure 23 is shown below, with respect to the intersections 401, 801, 1201 and 1601 of the exploration stratum lateral direction 1001 in system 912.

[0050] Figure 27a , Figure 27b , Figure 27c and Figure 27d The DTCIG waveforms obtained by processing the initial model system in Figure 23 at intersections 401, 801, 1201, and 1601 of the exploration stratum lateral line 1001 under system 912 are shown respectively. The white line represents the maximum amplitude selected at all depths along the time delay direction. The dashed red line represents zero time delay.

[0051] Figure 28a and Figure 28b The velocity accuracy / error properties, based on the DTCIG in Figure 27 and extracted using Equation (6), are shown in System 914 for the exploratory strata at lateral line 875 and cross line 1001. The gray line represents the uncertainty boundary, below which the velocity accuracy / error properties are no longer accurate due to the lack of wave penetration and time-delay gather focusing.

[0052] Figure 29a and Figure 29bThe velocity accuracy / error attributes, based on DTCIG in Figure 27 and extracted using formula (6), are shown in System 914 at depths of 5000m and 6000m in the explored strata.

[0053] Figure 30a , Figure 30b , Figure 30c and Figure 30d The DTCIG values ​​obtained from the original wavefield of the precise model system in Figure 16 are shown in system 912 at the intersections 401, 801, 1201 and 1601 of the exploration stratum lateral line 1001.

[0054] Figure 31a , Figure 31b , Figure 31c and Figure 31d The DTCIG waveforms obtained by processing the initial model system in Figure 16 at intersections 401, 801, 1201, and 1601 of the exploration stratum lateral line 1001 under system 912 are shown respectively. The white line represents the maximum amplitude selected at all depths along the time delay direction. The dashed red line represents zero time delay.

[0055] Figure 32a and Figure 32b The velocity accuracy / error properties, based on the DTCIG in Figure 31 and extracted using formula (6), are shown in system 914 for the exploratory strata at lateral line 875 and cross line 1001. The gray line represents the uncertainty boundary, below which the velocity accuracy / error properties are no longer accurate due to the lack of wave penetration and time-delay gather focusing.

[0056] Figure 33a and Figure 33b The velocity accuracy / error attributes, based on DTCIG in Figure 31 and extracted using formula (6), are shown in System 914 at depths of 5000m and 6000m in the explored strata.

[0057] Figure 34 The diagram illustrates a flowchart for determining model error and model updating according to one embodiment.

[0058] Figure 35 The flowchart illustrates an iterative update of a model according to one embodiment.

[0059] Figure 36a Demonstrates the use of Figure 15 The recalibrated velocity gradient is obtained from the velocity accuracy / error attribute in the equation. This is the process required for equation (7) and for systems 3414, 3416, 3514 and 3516. Figure 36b Demonstrates the use of Figure 36a velocity gradient and Figure 11aThe update rate is obtained from the initial velocity model in the equation (7) and the process required for systems 3416, 3516 and 3518.

[0060] Figure 37 A schematic diagram of an example computing system is shown. Detailed Implementation

[0061] The following description of exemplary embodiments is made with reference to the accompanying drawings. In different drawings, the same reference numerals denote the same or similar elements. The following detailed description does not limit the invention. The scope of the invention is defined by the appended claims. The following embodiments are discussed based on terminology related to seismic data acquisition and processing.

[0062] Throughout this specification, references to "an embodiment" or "a particular embodiment" indicate that a feature, structure, or characteristic described in connection with that embodiment is included in at least one embodiment of the disclosed subject matter. Therefore, phrases such as "in one embodiment" or "in a particular embodiment" appearing in different locations throughout the specification do not necessarily refer to the same embodiment. Furthermore, the specific feature, structure, or characteristic may be combined in any suitable manner in one or more embodiments.

[0063] Classical FWI aims to obtain the optimal velocity model for explored formations by minimizing the following square amplitude mismatch objective function. Where t represents time. and These indicate the locations of the earthquake source and the receiving point, respectively. This represents a subsurface location within the explored strata. d represents measured seismic data, and F(v) represents a wave equation simulator, which simulates the data by solving the wave equation using numerical algorithms such as finite difference, finite element, or finite volume on a possible subsurface velocity model v.

[0064] Two main problems can cause classical FWI to fail in production: cycle-skipping and amplitude discrepancies between measured data (d) and simulated data (u = F(v)). First, if the initial model differs significantly from the actual model, the waveform differences between measured and simulated data are substantial, and the oscillation patterns of the seismic signal can introduce local minima. Therefore, FWI may attempt to fit the simulated data to incorrect cycles in the measured data, thus converging to an incorrect velocity. Second, the amplitude distribution of measured data often differs significantly from that of simulated data. If this difference is not properly handled, FWI may misinterpret it as a velocity error, leading to the inversion of an incorrect velocity model. Furthermore, the amplitude of measured data is not only determined by the velocity model but is also influenced by other subsurface properties (such as density). Therefore, even if the velocity model is correct, it cannot be guaranteed that the amplitude of the simulated data will be consistent with that of the measured data. To address the problems of cycle-skipping and amplitude discrepancies in FWI, many solutions have been proposed. However, most of these methods are built only in the traditional data domain, and their success depends on extracting the travel time error between measured and simulated data, which is not a simple task.

[0065] Figure 9 A flowchart is shown of a seismic processing method 900 that quantifies, analyzes, or performs quality control on the accuracy of a model system 904 using time delay offset information, according to one embodiment. Method 900 may include acquiring measured data 902 as input. For example, this measured data may be acquired by one or more physical seismic acquisition devices (such as seismographs), as described above. Furthermore, the acquired seismic data may be preprocessed, such as by removing noise, and / or subjected to other appropriate processing. Figure 10b Demonstrates the use of Figure 10a Simplified synthetic seismic data generated from the accurate P-wave velocity model is used as an example of measured data 902.

[0066] Method 900 may also include a model system 904 as input for describing or simulating waveforms or wavefields within the explored strata. This model system may include, but is not limited to: a P-wave velocity model, a S-wave velocity model, an anisotropy model, elastic parameters, a Q-model (for describing subsurface attenuation), a density model, etc. Furthermore, this model system may be constructed based on a subsurface strata model. For example, the model may be based on any prior information about the subsurface strata known before executing method 900. This information may include general characteristics of the subsurface strata and seismic waveforms. Figure 11a A simple P-wave velocity model is shown as an example of the initial model 904, which is similar to... Figure 10a The accurate P-wave velocity model is different.

[0067] Method 900 may further include acquiring a wavelet or data 906 as input. This wavelet may be a direct or processed acquisition measurement result, an analytical or numerical instrument response, an inversion of the instrument response, or the output of a data matching processing flow. Furthermore, the data may be acquired by one or more physical seismic acquisition devices (such as seismographs), as described above. The acquired seismic data may be preprocessed, such as noise removal, and / or subjected to other appropriate processing. Figure 11b A wavelet profile is shown as an example of wavelet 906.

[0068] Method 900 may include generating a reverse-propagating wavefield based on input data 902 and using model system 904, propagating the wavefield backward over time using a suitable wave propagation engine, and combining appropriate initial time conditions, spatial boundary conditions, and source terms, as shown in 908, as a simulation of the reverse-propagating wavefield. The reverse-propagating wavefield decreases over time. The simplest example generated is as follows: in, This indicates the longitudinal wave velocity. Figure 12a Demonstrates the use of Figure 11a velocity model for Figure 10b The seismic data was processed, and an example of a backpropagating wavefield snapshot was generated at position 908.

[0069] Method 900 may include generating a propagating wavefield based on wavelet or data 906 and utilizing a model system 904, propagating the wavefield forward over time using a suitable wave propagation engine, and combining appropriate initial time conditions, spatial boundary conditions, and source terms, as shown in 910, as a simulation of the propagating wavefield. The propagating wavefield increases over time. The simplest example generated is as follows: Figure 12b Demonstrates the use of Figure 11a velocity model for Figure 11b The wavelet is processed to generate a snapshot example of the propagating wave field at 910.

[0070] Method 900 may include forming delayed-time common imaging gathers (DTCIGs) (see Sava, P. and Fomel, S., 2006, “Time-shift imaging condition in seismic migration”: Geophysics, 71, 209-217), as shown in 912, as gather formation. It can be generated as follows: Where τ is the delay time, and Ψ and Φ are appropriate processing operators acting on the forward and reverse propagation wave fields, used to extract the required model information from DTCIGs. Figure 13a and Figure 13b The figures show, under system 912, a typical spatial location in the explored strata, using... Figure 11a The velocity model is an example of DTCIG generated from the original wavefield and the processed wavefield. The white line represents the maximum amplitude selected at all depths along the time delay direction. The dashed red line represents zero time delay.

[0071] Method 900 may include virtually, analytically, numerically, or statistically extracting model accuracy / error information (e.g., delay time picking, overlay, display, or attributes) from DTCIGs, as shown in 914, as a model accuracy analysis for measuring, analyzing, or quality controlling model accuracy.

[0072] Figure 14a and Figure 14b The figures show the spatial locations of the explored strata, which are the same as those in Figure 13, through... Figure 10a An example of a DTCIG generated from the original and processed wavefields of an accurate model system is shown. The white line represents the maximum amplitude selected at all depths along the delay time direction. The dashed red line represents the zero-delay time lag. For the correct velocity, the refracted waves in the forward and reverse propagation wavefields will meet at the correct time on the ray path; their energy is thus focused at the zero-delay time (τ = 0) on the DTCIG, as shown by the dashed red line in Figure 14. On the other hand, an incorrect initial model will cause the refracted waves to be out of focus and deviate from the correct time, as shown by the dashed red line in Figure 13 at the zero-delay time. Therefore, velocity accuracy / error can be quantified and analyzed by the focusing of the refracted waves on the DTCIG and the degree to which they deviate from the zero-delay time.

[0073] Therefore, a simple example of the extraction process at 914 is delayed time picking, as shown by the white lines in Figures 13 and 14. Figure 15 The time delay pick-up values ​​are displayed at all spatial locations of the explored strata, as shown in DTCIG in Figure 13. Figure 15 The delay time pick-up value in the model thus provides information on model accuracy / error.

[0074] For completeness, a flowchart of a seismic processing method 900 that quantifies, analyzes, or quality controls the accuracy of a model system 904 using delayed time offset information is shown using a 3D SEG Advanced Modeling (SEAM) model (Fehler, M. and K. Larner, 2008, “SEG Advanced Modeling (SEAM): Phase I first year update”: The Leading Edge 27:1006-1007).

[0075] Method 900 may include acquiring measured data 902 as input. For example, this measured data may be acquired by one or more physical seismic acquisition devices (such as seismographs), as described above. Furthermore, the acquired seismic data may be preprocessed, such as by removing noise, and / or subjected to other appropriate processing. Figure 22b Demonstrates the use of Figure 22a The source wavelet in Figures 902 and the simplified synthetic seismic data generated by the accurate P-wave velocity model, anisotropic Thomsen parameters δ, ε, dip, azimuth and density model in Figures 16, 17, 18, 19, 20 and 21 are used as an example of measured data 902.

[0076] Method 900 may also include a model system 904 as input to describe or simulate waveforms or wavefields within the explored strata. This model system may include, but is not limited to: a P-wave velocity model, a S-wave velocity model, an anisotropy model, elastic parameters, a Q-model (for describing subsurface attenuation), a density model, etc. Furthermore, this model system may be constructed based on a subsurface strata model. For example, the model may be based on any prior information about the subsurface strata known before executing method 900. This information may include general characteristics of the subsurface strata and seismic waveforms. Figure 23 shows a simple P-wave velocity model as an example of the initial model 904, which differs from the precise P-wave velocity model in Figure 16. The differences between them are shown in Figure 24.

[0077] Method 900 may further include acquiring a wavelet or data 906 as input. This wavelet may be a direct or processed acquisition measurement result, an analytical or numerical instrument response, an inversion of the instrument response, or the output of a data matching processing flow. Furthermore, the data may be acquired by one or more physical seismic acquisition devices (such as seismographs), as described above. The acquired seismic data may be preprocessed, such as noise removal, and / or subjected to other appropriate processing. Figure 22a A wavelet profile is shown as an example of wavelet 906.

[0078] Method 900 may include generating a reverse propagating wavefield based on input data 902 and using model system 904, propagating the wavefield backward over time using a suitable wave propagation engine, and combining appropriate initial time conditions, spatial boundary conditions, and source terms, as shown in 908, as a simulation of the reverse propagating wavefield. The simplest example of generating the reverse waveform is solving equation (2) and storing it. Figure 25a The velocity model shown in Figure 23 is used to illustrate the application of velocity model to the target velocity. Figure 22b The seismic data was processed, and an example of a backpropagating wavefield snapshot was generated at position 908.

[0079] Method 900 may include generating a forward propagating wavefield based on wavelet or data 906 and utilizing a model system 904, propagating the wavefield forward over time using a suitable wave propagation engine, and combining appropriate initial time conditions, spatial boundary conditions, and source terms, as shown in 910, as a simulation of the forward propagating wavefield. The simplest example of generating a forward propagating waveform is solving equation (3) and storing it. Figure 25b The velocity model shown in Figure 23 is used to illustrate the application of velocity model to the target velocity. Figure 22a The wavelet is processed to generate a snapshot example of the propagating wave field at 910.

[0080] Method 900 may include forming delayed-time common imaging gathers (DTCIGs) (see Sava, P. and Fomel, S., 2006, “Time-shift imaging condition in seismic migration”: Geophysics, 71, 209-217), as shown in 912, as gather formation. It can be generated using formula (4). Figure 26 shows examples of DTCIG generated from the original wavefield and the processed wavefield using the velocity model in Figure 23 at typical spatial locations of the explored strata under system 912. The white line represents the maximum amplitude selected at all depths along the time delay direction. The dashed red line represents zero time delay.

[0081] Method 900 may include virtually, analytically, numerically, or statistically extracting model accuracy / error information (e.g., delay time picking, overlay, display, or attributes) from DTCIGs, as shown in 914, as model accuracy analysis for measuring, analyzing, or quality controlling model accuracy. Except Figure 15In addition to the two-dimensional spatial illustration, Figures 28 and 29 summarize the model accuracy / error properties in three-dimensional space. At each three-dimensional spatial location in the explored strata, the delay time value picked up by the maximum amplitude in the DTCIG in Figure 27 is shown. Figures 28 and 29 illustrate cases where the delay time deviates from the ideal zero-delay time lag. Close observation reveals that the linear velocity perturbation from shallow to deep, as shown in Figure 24, matches the accuracy / error properties in Figures 28 and 29. An uncertainty boundary can also be extracted through the lack of wave penetration and delay time gather focusing; below this boundary, the velocity accuracy / error properties are no longer accurate, as shown in Figure 28.

[0082] In contrast, Figures 30 and 31 show the results at the same spatial locations in the explored strata as Figures 26 and 27, respectively, through... Figure 16a Examples of DTCIG for the original and processed wavefields generated by the accurate model system are shown. The white line represents the maximum amplitude selected at all depths along the time delay direction. The dashed red line represents zero time delay. For the correct velocity, the refracted waves in the forward and reverse propagation wavefields will meet at the correct time on the ray path; their energy is thus focused at the zero-delay time (τ=0) on the DTCIG, as shown by the dashed red line in Figure 31. On the other hand, an incorrect initial model will cause the refracted waves to be out of focus and deviate from the correct time, as shown by the dashed red line in Figure 27. Therefore, velocity accuracy / error can be quantified and analyzed by the focusing of the refracted waves on the DTCIG and the degree to which they deviate from the zero-delay time.

[0083] At each three-dimensional spatial location of the explored formation, the delay time value of the maximum amplitude picked up in the DTCIG in Figure 31 is shown, as in Figures 28 and 29. Figures 32 and 33 show the velocity accuracy / error attribute. Above the uncertainty boundary indicated by the gray line, the value of the velocity accuracy / error attribute approaches the ideal zero delay time lag, thus providing accurate quantitative measurements of velocity.

[0084] Therefore, a simple example of the extraction process at 914 is time delay picking, as shown by the white lines in Figures 27 and 31. Figures 28, 29, 32, and 33 show the time delay picking values ​​at all spatial locations of the explored strata, as indicated by the DTCIG in Figures 27 and 31. Therefore, the time delay picking values ​​in Figures 28, 29, 32, and 33 provide information on model accuracy / error.

[0085] Figure 34A flowchart is shown of a seismic processing method 3400 for updating a model system 3404 using time delay offset information, according to one embodiment. Method 3400 may include acquiring measured data 3402 as input. For example, this measured data may be acquired by one or more physical seismic acquisition devices (such as seismographs), as described above. Furthermore, the acquired seismic data may be preprocessed, such as by noise removal, and / or subjected to other appropriate processing. Figure 10b Demonstrates the use of Figure 10a Simplified synthetic seismic data generated from the accurate P-wave velocity model is used as an example of measured data 3402.

[0086] Method 3400 may also include an initial model system 3404 of the explored strata as input, used to describe or simulate waveforms or wave fields within the explored strata. This model system may include, but is not limited to: P-wave velocity models, S-wave velocity models, anisotropy models, elastic parameters, Q-models (for describing subsurface attenuation), density models, etc. Furthermore, this model system may be constructed based on a subsurface strata model. For example, the model may be obtained based on any prior information about the subsurface strata known before executing method 3400. This information may include general characteristics of the subsurface strata and seismic waveforms. Figure 11a A simple P-wave velocity model is shown as an example of the initial model 3404, which is similar to... Figure 10a The accurate P-wave velocity model is different.

[0087] Method 3400 may further include acquiring a wavelet or data 3406 as input. This wavelet may be a direct or processed acquisition measurement result, an analytical or numerical instrument response, an inversion result of an instrument response, or the output of a data matching processing flow. Furthermore, the data may be acquired by one or more physical seismic acquisition devices (such as seismographs), as described above. The acquired seismic data may be preprocessed, such as noise removal, and / or subjected to other appropriate processing. Figure 11b A wavelet profile is shown as an example of wavelet 3406.

[0088] Method 3400 may include generating a reverse propagating wavefield based on input data 3402 and using model system 3404, propagating the wavefield backward over time using a suitable wave propagation engine, and combining appropriate initial time conditions, spatial boundary conditions, and source terms, as shown in 3408, as a simulation of the reverse propagating wavefield. The simplest example of generating the reverse waveform is solving equation (2) and storing it. Figure 12a Demonstrates the use of Figure 11a velocity model for Figure 10b The seismic data was processed, and an example of a backpropagating wavefield snapshot was generated at position 3408.

[0089] Method 3400 may include generating a forward propagating wavefield based on wavelet or data 3406 and utilizing model system 3404, propagating the wavefield forward over time using a suitable wave propagation engine, and combining appropriate initial time conditions, spatial boundary conditions, and source terms, as shown in 3410, as a simulation of the forward propagating wavefield. The simplest example of generating a forward waveform is solving equation (3) and storing it. Figure 12b Demonstrates the use of Figure 11a velocity model for Figure 11b The wavelet is processed to generate a snapshot example of the propagating wave field at 3410.

[0090] Method 3400 may include forming delayed-time common imaging gathers (DTCIGs) (see Sava, P. and Fomel, S., 2006, “Time-shift imaging condition in seismic migration”: Geophysics, 71, 209-217), as shown in 3412, as gather formation. It can be generated using formula (4). Figure 13a and Figure 13b The figures at location 3412, representing a typical spatial location within the explored strata, are shown using... Figure 11a The velocity model is an example generated from the original wavefield and the processed wavefield using a DTCIG. The white line represents the maximum amplitude selected at all depths along the time delay direction. The dashed red line represents zero-delay time lag.

[0091] Method 3400 may include a model accuracy / error information extraction process T, such as delay time picking, overlaying, displaying, or attributeing, as follows: in, This represents the model gradient generated from the extracted time delay offset information, used to describe the model's accuracy / error. This process occurs at 3414, for example, as part of model gradient formation.

[0092] Method 3400 may include updating the model system using the model gradients in the initial model systems 3404 and 3414, as shown in 3434, as a model update, as follows: m+αδm(6) Here, α is a scalar that can be calculated manually, automatically, or by the system.

[0093] Figure 35A flowchart for seismic processing, such as method 3500 using a full-waveform inversion automatic model update system 3504, is shown according to one embodiment. Method 3500 may include acquiring measured data 3502 as input. For example, the measured data may be acquired by one or more physical seismic acquisition devices (such as seismographs), as described above. Furthermore, the acquired seismic data may be preprocessed, such as noise removal, and / or subjected to other appropriate processing. Figure 10b Demonstrates the use of Figure 10a Simplified synthetic seismic data generated from the accurate P-wave velocity model is used as an example of measured data 3502.

[0094] Method 3500 may also include an initial model system 3504 of the explored strata in the i-th iteration as input, used to describe or simulate waveforms or wave fields within the explored strata. This model system may include, but is not limited to: P-wave velocity models, S-wave velocity models, anisotropy models, elastic parameters, Q-models (used to describe subsurface attenuation), density models, etc. Furthermore, this model can be constructed based on a subsurface strata model. For example, this model system can be obtained based on any prior information about the subsurface strata known before executing method 3500. This information may include general characteristics of the subsurface strata and seismic waveforms. Figure 11a A simple P-wave velocity model is shown as an example of the initial model 3504, which is similar to... Figure 10a The accurate P-wave velocity model is different.

[0095] Method 3500 may also include acquiring a wavelet or data 3506 as input. This wavelet may be a direct or processed acquisition measurement result, an analytical or numerical instrument response, an inversion result of an instrument response, or the output of a data matching processing flow. Furthermore, the data may be acquired by one or more physical seismic acquisition devices (such as seismographs), as described above. The acquired seismic data may be preprocessed, such as noise removal, and / or subjected to other appropriate processing. Figure 11b A wavelet profile is shown as an example of wavelet 3506.

[0096] Method 3500 may include generating a reverse propagating wavefield based on input data 3502 and using model system 3504, propagating the wavefield backward over time using a suitable wave propagation engine, and combining appropriate initial time conditions, spatial boundary conditions, and source terms, as shown in 3508, as a simulation of the reverse propagating wavefield. The simplest example of generating the reverse waveform is solving equation (2) and storing it. Figure 12a Demonstrates the use of Figure 11a velocity model for Figure 10b The seismic data was processed, and an example of a backpropagating wavefield snapshot was generated at position 3508.

[0097] Method 3500 may include generating a forward propagating wavefield based on wavelet or data 3506 and utilizing model system 3504, propagating the wavefield forward over time using a suitable wave propagation engine, and combining appropriate initial time conditions, spatial boundary conditions, and source terms, as shown in 3510, as a simulation of the forward propagating wavefield. The simplest example of generating a forward waveform is solving equation (3) and storing it. Figure 12b Demonstrates the use of Figure 11a velocity model for Figure 11b The wavelet is processed to generate a snapshot example of the propagating wave field at 3510.

[0098] Method 3500 may include forming delayed-time common imaging gathers (DTCIGs) (see Sava, P. and Fomel, S., 2006, “Time-shift imaging condition in seismic migration”: Geophysics, 71, 209-235), as shown in 3512, as gather formation. It can be generated using formula (4). Figure 13a and Figure 13b The figures at location 3512, representing a typical spatial location within the explored strata, are shown using... Figure 11a The velocity model is an example generated from the original wavefield and the processed wavefield using a DTCIG. The white line represents the maximum amplitude selected at all depths along the time delay direction. The dashed red line represents zero-delay time lag.

[0099] Method 3500 may include the model accuracy / error information extraction process T in formula (5), such as delay time picking, overlay, display or attribute, as shown in 3514, as the model gradient formation process.

[0100] Method 3500 may include updating the model system using the model gradients in the initial model systems 3504 and 3514, as shown in 3534, as a model update, as follows: m i+1 = m i + αδm i (7) Where α is a scalar, which can be calculated manually, automatically, or by the system to minimize the following objective function: Where T is the model accuracy / error information extraction process using DTCIG, such as delay time picking, overlay, display, or attributes. Method 3500 may include the model accuracy / error information extraction process T in formula (5), as shown in 3514, for example as a velocity gradient formation process.

[0101] Method 3500 may include updating the model system m obtained in 3534. i+1 The initial model system in 3504 is reassigned for iterative updates in the next iteration, as shown in 3518, as an iterative model update.

[0102] Figure 36a Demonstrates the use of Figure 15 The inaccurate velocity model in the model is used to recalibrate the velocity gradient. This is the process required for formula (7) and 3414, 3434, 3514 and 3534. Figure 36b Demonstrates the use of Figure 18a velocity gradient and Figure 11a The update rate is obtained from the initial velocity model. This is the process required for formula (7) and 3434, 3534 and 3518.

[0103] The foregoing description has been made with reference to specific embodiments for ease of explanation. However, the illustrative discussion above is not intended to be exhaustive or to limit the invention to the precise forms disclosed. Many modifications and variations can be made based on the foregoing teachings. Furthermore, the order of the elements of the methods shown and / or described can be rearranged, and / or two or more elements can occur simultaneously. These embodiments were chosen and described to best illustrate the principles of the invention and its practical application, thereby enabling those skilled in the art to best utilize the invention and its embodiments in various modifications for a particular intended use.

[0104] The following disclosure relates to various operations that can be implemented by a computing device configured to execute program instructions instructing said operations. Operations may also be performed by circuitry specifically designed or configured for these operations. In some embodiments, a non-transitory computer-readable medium stores program instructions that enable it to perform the operations described herein. The terms “processor,” “processing unit,” or “processing element” refer to different elements or combinations thereof configured to execute program instructions. Processing elements may include, for example, application-specific integrated circuits (ASICs), custom processing circuitry or gate arrays, components or circuitry of a single processor core, an entire processor core, a single processor, programmable hardware devices (such as FPGAs or similar devices), and a large portion of a system comprising multiple processors, and any combination thereof.

[0105] See Figure 37The diagram illustrates a block diagram of an exemplary computing device 3710, according to certain embodiments. The computing device 3710 can be used to implement various parts of this disclosure. The computing device 3710 can be an instance of a mobile device, a server computing system, a client computing system, a distributed computing system, or any other computing system that implements parts of this disclosure. In different configurations, when programmed to execute a specific algorithm, the computing system 3710 can be configured to perform the functions corresponding to that specific algorithm.

[0106] The computing device 3710 may include, but is not limited to, a personal computer system, desktop computer, notebook or laptop computer, mobile phone, mainframe system, supercomputer, web server, workstation, or network computer. As shown, the computing device 3710 includes a processing unit 3750, a storage subsystem 3712, and an input / output (I / O) interface 3730 connected via an interconnect 3760 (e.g., a system bus). The I / O interface 3730 can connect to one or more I / O devices 3740. Furthermore, the computing device 3710 includes a network interface 3732 for connection to a network 3720 to communicate with other computing devices. Other bus architectures and subsystem configurations may also be employed.

[0107] As previously described, processing unit 3750 includes one or more processors. In some embodiments, processing unit 3750 includes one or more coprocessor units. In some embodiments, multiple instances of processing unit 3750 may be connected to interconnect 3760. Processing unit 3750 (or each processor therein) may include cache or other forms of on-chip memory. In some embodiments, processing unit 3750 may be implemented as a general-purpose processing unit, while in other embodiments it may be implemented as a special-purpose processing unit (e.g., ASIC). Computing device 3710 is not limited to any particular type of processing unit or processor subsystem.

[0108] Storage subsystem 3712 (which may include system memory and / or virtual memory) is used by processing unit 3750 (e.g., for storing executable instructions and data used by processing unit 3750). Storage subsystem 3712 may be implemented using any applicable type of physical storage medium, such as hard disk storage, floppy disk storage, removable disk storage, flash memory, random access memory (RAM—SRAM, EDO RAM, SDRAM, DDR SDRAM, RDRAM, etc.), read-only memory (ROM—PROM, EEPROM, etc.), etc. In some embodiments, storage subsystem 3712 may consist solely of volatile memory. Storage subsystem 3712 may store program instructions executed by computing device 3710 via processing unit 3750, including program instructions that cause computing device 3710 to implement the techniques described herein. In at least some embodiments, storage subsystem 3712 and / or medium 3714 may constitute instances of non-transitory computer-readable or machine-readable media for storing executable instructions.

[0109] In the illustrated embodiment, computing device 3710 also includes a non-transitory computer-readable medium 3714, which is distinguishable from storage subsystem 3712. As shown, computer-readable medium 3714 is configured as a peripheral or I / O device accessible via I / O interface 3730, but other interconnect configurations may also be used. In various embodiments, non-transitory medium 3714 may include persistent, tangible storage such as disk, non-volatile memory, magnetic tape, optical media, holographic media, or other suitable types of storage. In some embodiments, non-transitory medium 3714 can be used to store and transmit geophysical data and may be physically separate from computing device 3710 for portability. Therefore, in some embodiments, the aforementioned geophysical data products may be embodied in non-transitory medium 3714. Although shown as separate from storage subsystem 3712, in some embodiments, non-transitory medium 3714 may be integrated with storage subsystem 3712. Embodiments of the non-transient medium 3714 and / or storage subsystem 3712 may correspond to means for storing recorded seismic data, wavelets, model data, wavefields, DTCIGs, initial models, gradient models, and iteratively updated models, or to means for storing any I / O-related media in the middle of the process; these means may have different structures or may correspond to the same structure.

[0110] I / O interface 3730 may represent one or more interfaces and may be any of a variety of interfaces configured to couple and communicate with other devices. In some embodiments, I / O interface 3730 is a bridge chip from a front-side bus to one or more back-side buses. I / O interface 3730 may be coupled to one or more I / O devices 3740 via one or more corresponding buses or other interfaces. Examples of I / O devices include storage devices (hard disks, optical drives, removable flash drives, storage arrays, storage area networks (SANs) or their associated controllers), network interface devices, user interface devices, or other devices (e.g., graphics, sound, etc.). In some embodiments, the aforementioned geophysical data products may be embodied within one or more I / O devices 3740.

[0111] This specification contains references to "one embodiment," "some embodiments," or "a particular embodiment." The use of these phrases does not necessarily refer to the same embodiment. Specific features, structures, or characteristics may be combined in any suitable manner without departing from this disclosure.

Claims

1. A method for seismic data processing, comprising: Receive seismic data, including seismic waveform signals that characterize the geological structure of underground areas, collected by the seismic receiver; Based on the initial model system, backpropagation wavefields are generated using seismic data through backward propagation in time; Using a given source wavelet or processed seismic data, a forward propagation wavefield is generated by forward propagation in time. Delayed Time Common Imaging Point Gathers (DTCIGs) are generated based on the reverse propagation wavefield and the forward propagation wavefield. Extract model accuracy information from DTCIGs for analysis and / or quality control; Model errors are determined based on model accuracy information and / or DTCIGs and used for analysis and / or quality control; as well as The initial model system is adjusted based on model errors to generate an adjusted / updated model.

2. The method according to claim 1, wherein, Determining the model error involves identifying a delay-shift objective function in the model domain, which characterizes the delay-shift extracted from the DTCIGs.

3. The method according to claim 2, wherein, Determining model error also includes determining local time delay errors.

4. The method according to claim 3, wherein, Determining the local time delay error involves directly minimizing the time delay objective function as a local attribute.

5. The method according to claim 3, wherein, Determining model errors also includes converting time delay errors into amplitude and / or phase errors in the instantaneous model domain.

6. The method according to claim 5, wherein, Determining the model error also includes determining the objective function relationship of the instantaneous amplitude and / or phase based on the instantaneous amplitude and / or phase error.

7. The method according to claim 5, wherein, Determining model error also includes using the instantaneous amplitude and / or phase error to form the model update gradient.

8. The method according to claim 5, wherein, Determining model error also includes using the instantaneous amplitude and / or phase errors to calculate amplitude and / or phase objective functions.

9. The method of claim 1 further includes performing full waveform inversion using the adjusted model as a background model.

10. The method of claim 9, further comprising interpreting one or more features of the underground region based on the results of full waveform inversion.

11. The method of claim 1, further comprising performing seismic signal imaging on the underground area and locating natural resources in the underground area based on the image.

12. A computing system, comprising: processor; as well as The memory includes instructions stored on a non-transitory computer-readable medium, which, when executed by the processor, cause the computing system to perform the following operations: Receive seismic data, including seismic waveform signals that characterize the geological structure of underground areas, collected by the seismic receiver; Based on the initial model system, backpropagation wavefields are generated using seismic data through backward propagation in time; Using a given source wavelet or processed seismic data, a forward propagation wavefield is generated by forward propagation in time. Delayed Time Common Imaging Point Gathers (DTCIGs) are generated based on the reverse propagation wavefield and the forward propagation wavefield. Extract model accuracy information from DTCIGs for analysis and / or quality control; Model errors are determined based on model accuracy information and / or DTCIGs and used for analysis and / or quality control; and The initial model system is adjusted based on model errors to generate an adjusted / updated model.

13. The system according to claim 12, wherein, Determining the model error involves identifying a delay-shift objective function in the model domain, which characterizes the delay-shift extracted from the DTCIGs.

14. The system according to claim 13, wherein, Determining model error also includes determining local time delay errors.

15. The system according to claim 14, wherein, Determining the local time delay error involves directly minimizing the time delay objective function as a local attribute.

16. The system according to claim 14, wherein, Determining model errors also includes converting time delay errors into amplitude and / or phase errors in the instantaneous model domain.

17. The system according to claim 16, wherein, Determining the model error also includes determining the objective function relationship of the instantaneous amplitude and / or phase based on the instantaneous amplitude and / or phase error.

18. The system according to claim 16, wherein, Determining model error also includes using the instantaneous amplitude and / or phase error to form the model update gradient.

19. The system according to claim 16, wherein, Determining model error also includes using the instantaneous amplitude and / or phase errors to calculate amplitude and / or phase objective functions.

20. The system of claim 12 further includes performing full waveform inversion using the adjusted model as a background model.

21. The system of claim 20 further includes interpreting one or more features of the underground region based on the results of full waveform inversion.

22. A non-transitory computer-readable medium having instructions stored thereon that, when executed by a processor of a computing system, cause the computing system to perform the following operations: Receive seismic data, including seismic waveform signals that characterize the geological structure of underground areas, collected by the seismic receiver; Based on the initial model system, backpropagation wavefields are generated using seismic data through backward propagation in time; Using a given source wavelet or processed seismic data, a forward propagation wavefield is generated by forward propagation in time. Delayed Time Common Imaging Point Gathers (DTCIGs) are generated based on the reverse propagation wavefield and the forward propagation wavefield. Extract model accuracy information from DTCIGs for analysis and / or quality control; Model errors are determined based on model accuracy information and / or DTCIGs and used for analysis and / or quality control; as well as The initial model system is adjusted based on model errors to generate an adjusted / updated model.

23. The medium according to claim 22, wherein, Determining model errors includes: The time-shift delay objective function is directly minimized as a local property to determine the local time-shift delay error; Convert the time delay error into amplitude and / or phase error in the instantaneous model domain; Determine the objective function relationship of instantaneous amplitude and / or phase based on instantaneous amplitude and / or phase error; The gradient is generated using instantaneous amplitude and / or phase error to form the model update gradient; The objective function for amplitude and / or phase is calculated using instantaneous amplitude and / or phase error; Perform full waveform inversion using the adjusted model as the background model.

Images