Method and system for reflection-based traveltime inversion using piecewise dynamic image regularization

Abstract

A computer-implemented method can include obtaining seismic data acquired in a time domain for a subsurface region of interest (300). The method can also include obtaining a property model for the subsurface region of interest (305). The method can also include determining one or more time shifts using a piecewise dynamic image warping function based on the seismic data and the property model (320). The method can also include determining a ghost operator using the derived time shifts and a one-way wave equation (325). The method can also include updating the property model using a gradient solver in a data domain reflection traveltime inversion (335). The method can also include outputting the updated property model for the subsurface region of interest (345). The method can also include generating a seismic image for the subsurface region of interest using the updated property model (345).

Description

Background Technology

[0001] Accurate velocity model construction is crucial for subsurface imaging and reservoir characterization. Full migration wavefield inversion algorithms can transform time-domain seismic data into a subsurface depth representation. Specifically, full waveform inversion can be performed to construct high-resolution models for seismic imaging and reservoir characterization. This model can represent the velocity values ​​of different subsurface particles. However, full waveform inversion may require a high-quality background model as an initial model, which then needs to be updated in multiple iterations at a significant computational cost. Summary of the Invention

[0002] This summary is provided to introduce a series of concepts that will be further described in the following detailed description. This summary is not intended to identify key or essential features of the claimed subject matter, nor is it intended to help limit the scope of the claimed subject matter.

[0003] In general, in one aspect, embodiments relate to a method comprising acquiring seismic data for a subsurface region of interest, acquired in the time domain. The method further comprises acquiring a property model for the subsurface region of interest. The method further comprises determining one or more time shifts based on the seismic data and the property model using a segment dynamic image warping function. The method further comprises determining an adjoint source operator using derived time shifts and one-way wave equations. The method further comprises updating the property model using a gradient solver in a data-domain reflection travel-time inversion. The method further comprises outputting the updated property model for the subsurface region of interest. The method further comprises generating a seismic image for the subsurface region of interest using the updated property model.

[0004] In general, in one aspect, embodiments relate to a system including a seismic survey system. The system also includes a seismic source and multiple seismic receivers. The system further includes a seismic interpreter comprising a computer processor. The seismic interpreter is coupled to the seismic survey system. The seismic interpreter acquires seismic data for a subsurface region of interest, acquired in the time domain. The seismic interpreter obtains a property model for the subsurface region of interest. Based on the seismic data and the property model, the seismic interpreter uses a piecewise dynamic image warping function to determine one or more time shifts. The seismic interpreter uses derived time shifts and one-way wave equations to determine an adjoint source operator. The seismic interpreter updates the property model using a gradient solver in a data-domain reflection travel-time inversion. The seismic interpreter outputs the updated property model for the subsurface region of interest. The seismic interpreter uses the updated property model to generate a seismic image for the subsurface region of interest.

[0005] In general, in one aspect, embodiments relate to a non-transitory computer-readable medium storing instructions executable by a computer processor. The instructions include obtaining seismic data acquired in the time domain for a subsurface region of interest. The instructions also include obtaining a property model for the subsurface region of interest using a piecewise dynamic image warping function. The instructions further include determining one or more time shifts based on the seismic data and the property model. The instructions further include determining an adjoint source operator using derived time shift and one-way wave equations. The instructions further include updating the property model using a gradient solver in a data-domain reflection travel-time inversion. The instructions further include outputting the updated property model for the subsurface region of interest. The instructions further include generating a seismic image for the subsurface region of interest using the updated property model.

[0006] Other aspects of this disclosure will become apparent from the following description and the appended claims. Attached Figure Description

[0007] Specific embodiments of the disclosed technology will now be described in detail with reference to the accompanying drawings. For consistency, similar elements are indicated by similar reference numerals in the drawings.

[0008] Figure 1 and Figure 2 A system according to one or more embodiments is shown.

[0009] Figure 3 A flowchart according to one or more embodiments is shown.

[0010] Figure 4 Examples according to one or more embodiments are shown.

[0011] Figure 5 A, Figure 5 B and Figure 5 C illustrates an example according to one or more embodiments.

[0012] Figure 6A , Figure 6B , Figure 6C and Figure 6D Examples according to one or more embodiments are shown.

[0013] Figure 7A , Figure 7B , Figure 7C , Figure 7D , Figure 7E and Figure 7F Examples according to one or more embodiments are shown.

[0014] Figure 8A and Figure 8B Examples according to one or more embodiments are shown.

[0015] Figure 9A , Figure 9B , Figure 9C , Figure 9D , Figure 9E , Figure 9F , Figure 9G and Figure 9H Examples according to one or more embodiments are shown.

[0016] Figure 10A , Figure 10B , Figure 10C and Figure 10D Examples according to one or more embodiments are shown.

[0017] Figure 11A and Figure 11B Examples according to one or more embodiments are shown.

[0018] Figure 12 A computing system according to one or more embodiments is shown. Detailed Implementation

[0019] Numerous specific details are set forth in the following detailed description of embodiments of the present disclosure in order to provide a more thorough understanding of the disclosure. However, it will be apparent to those skilled in the art that the disclosure may be practiced without these specific details. In other instances, well-known features have not been described in detail to avoid unnecessarily complicating the description.

[0020] Throughout the application, ordinal numbers (e.g., first, second, third, etc.) may be used as adjectives for elements (i.e., any noun in this application). Unless explicitly disclosed, such as by using the terms “before,” “after,” “single,” and other such terms, the use of ordinal numbers does not imply or create any particular order of elements, nor does it limit any element to a single element. Rather, the use of ordinal numbers is intended to distinguish between elements. As an example, a first element is distinct from a second element, and a first element may contain more than one element and be placed after (or before) the second element in the order of elements.

[0021] In general, embodiments of this disclosure include systems and methods for performing data domain reflection travel-time inversion (DRTI) using a piecewise dynamic image warping (SDIW) function. DRTI is applied to generate kinematically accurate property models (e.g., velocity models) based on reflection energy and one-way wave equations for seismic imaging and velocity model construction. For example, DRTI generates high-quality background property models to correctly image subsurface geological structures by matching travel-time information between predicted data (e.g., inverse migration data) and acquired seismic data. As another example, DRTI can generate kinematically accurate property models as initial models for full waveform inversion (FWI) to further refine models with short-wavelength components.

[0022] Furthermore, data domain reflection traveltime inversion uses a misfit function based on travel time difference to solve a least-squares optimization problem. This travel time difference is measured by applying a piecewise dynamic image warping function between the predicted data obtained using the migration function and the acquired seismic data. However, obtaining a high-quality property model presents several challenges. For example, the property model (e.g., velocity model) of the geological region of interest may be affected by strong noise, migration artifacts, and various unwanted signals due to subsurface complexity and seismic acquisition limitations (e.g., low-velocity zones). As another example, certain data domain reflection traveltime inversions may exhibit slow convergence. Thus, some embodiments address these issues by applying a piecewise dynamic image warping function within the data domain reflection traveltime inversion algorithm. In particular, piecewise dynamic image warping involves point-by-point segmentation to a piecewise matching function to measure the cumulative distance between two signals (e.g., Euclidean distance, negative cosine of the angle, sum of absolute differences, etc.) for the superimposed windowed polynomial fit of each segment. Therefore, the piecewise dynamic image warping function enhances the robustness of signal alignment and the reliable time shift estimation between the predicted data and the acquired seismic data.

[0023] Turning Figure 1 , Figure 1 A schematic diagram according to one or more embodiments is shown. Figure 1 The diagram illustrates a seismic survey system 100 and various formation paths of pressure waves (also referred to as seismic waves). The seismic survey system 100 includes a source 122 that includes functionality for generating pressure waves that penetrate the subsurface layer 124, such as reflected waves 136, latent waves A 142, or latent waves B 146. The pressure waves generated by the source 122 can penetrate the subsurface layer 124 along several paths according to a particle velocity V1, for detection at multiple seismic receivers 126 along a profile line. Similarly, particle velocity can refer to various velocity types, such as the two types of particle motion caused by seismic waves: the velocity of the first arrival wave (P-wave) and the different velocities of the second arrival wave (S-wave) through a specific medium. The source 122 can be a seismic vibrator, such as a seismic vibrator using controlled-source technology, an air gun in marine seismic surveys, explosives, etc. The seismic receivers 126 can include detectors, hydrophones, accelerometers, and other sensing devices. Similarly, the seismic receiver 126 may include single-component and / or multi-component sensors for measuring pressure waves on multiple spatial axes.

[0024] like Figure 1 As shown, the seismic source 122 generates an air wave 128, formed from a portion of the emitted seismic energy, which travels over the Earth's surface 130 to the seismic receiver 126. The seismic source 122 can also emit surface waves 132 traveling along the Earth's surface 130. The velocity of the surface waves 132 (also known as Rayleigh waves or roll waves) can correspond to a particle velocity that is generally slower than that of secondary waves. Although Figure 1 The seismic survey shown is a two-dimensional survey along a seismic profile in the longitudinal direction, but other embodiments, such as three-dimensional surveys, are also conceivable.

[0025] Furthermore, the subsurface layer 124 has a particle velocity V1, while the subsurface layer 140 has a particle velocity V2. In other words, different subsurface layers can correspond to different particle velocity values. Specifically, particle velocity can refer to the speed at which a pressure wave travels through a medium; for example, a latent wave B 146 travels through the subsurface layer 124, forming a curved ray path 148. Particle velocity may depend on the density and elasticity of a particular medium, as well as various wave properties, such as the frequency at which the pressure wave is emitted. When the particle velocities differ between the two subsurface layers, this seismic wave impedance mismatch can lead to seismic reflections of the pressure wave. For example, Figure 1 The diagram illustrates a pressure wave propagating downwards from the source 122 to the lower ground interface 138, which, in response to seismic reflection, becomes an upward-propagating reflected wave 136. The source 122 can also generate a direct wave 144 that travels directly from the source 122 through the lower ground layer 124 to the seismic receiver 126 at a particle velocity V1.

[0026] The source 122 can also generate a refracted wave (i.e., a latent wave A 142) by redirecting the refracted pressure wave. This refracted wave is refracted at the subsurface interface 138 and travels a certain distance along the subsurface interface 138 (e.g., ...). Figure 1 As shown), the refracted pressure wave travels upwards to the seismic receiver 126. Thus, the refracted pressure wave may include latent waves (e.g., latent wave A142, latent wave B146) that can be analyzed to map the subsurface layers 124, 140. For example, a latent wave can be a refracted wave that continuously refracts at various points beneath the Earth's surface. Therefore, latent waves can be generated where the particle velocity gradually increases with depth. Similarly, unlike reflected seismic energy, the vertex of a latent wave can be offset from the common midpoint (CMP). However, for analytical purposes, the vertex of the latent wave can be considered as the common midpoint of the refracted energy. Thus, the vertex can be used as a basis for organizing and classifying seismic survey datasets.

[0027] Furthermore, when analyzing seismic data acquired using the seismic survey system 100, rays can be used to approximate seismic wave propagation. For example, reflected waves (e.g., reflected wave 136) and latent waves (e.g., latent waves 142, 146) can be scattered at the subsurface interface 138. For example, in Figure 1 In this context, the latent wave B 146 can present a wide-angle ray path similar to a reflected wave, enabling mapping of the subsurface. For example, using latent waves, a velocity model can be generated for the subsurface, describing the particle velocities in different regions of different subsurface layers. An initial velocity model can be generated by simulating the velocity structure of the subsurface medium using seismic data inversion (often referred to as seismic inversion). In seismic inversion, the velocity model is iteratively updated until the velocity model and seismic data have a minimum mismatch; for example, the solution of the velocity model converges to a globally optimal value that meets predetermined criteria.

[0028] Regarding velocity models, they can map various subsurface layers based on particle velocities in different sub-regions (e.g., P-wave velocities, S-wave velocities, and various anisotropic effects within the sub-regions). For example, a velocity model can be used, along with the arrival times and directions of P-waves and S-waves, to locate seismic events. Anisotropic effects can correspond to subsurface properties that cause pressure waves to exhibit direction dependence. Thus, seismic anisotropy can correspond to multiple parameters in geophysics involving variations in wave velocity based on propagation direction. One or more anisotropic algorithms can be performed to determine anisotropic effects, such as anisotropic ray tracing localization algorithms or algorithms using deviated well acoustic logging, vertical seismic profiling (VSP), and core measurements. Similarly, a velocity model can include multiple velocity boundaries that define regions of rock type variation, such as interfaces between different subsurface layers. In some embodiments, the velocity model is updated using one or more fault photography updates to adjust the velocity boundaries within the velocity model.

[0029] Turning Figure 2 , Figure 2 A system according to one or more embodiments is shown. Figure 2 As shown, a seismic volume 290 comprises multiple seismic traces (e.g., seismic trace 250) acquired by multiple seismic receivers (e.g., seismic receiver 226) positioned on the Earth's surface 230. More specifically, the seismic volume 290 can be a three-dimensional cubic dataset of seismic traces in a 2D context. Individual cubic units within the seismic volume 290 may be referred to as voxels or volume pixels (e.g., voxel 260). In particular, different portions of the seismic traces may correspond to multiple depth points within the Earth's volume. To generate the seismic volume 290, a three-dimensional array of seismic receivers 226 is positioned along the Earth's surface 230, and seismic data is acquired in response to multiple pressure waves emitted from a source. Within voxel 260, statistics can be determined for the first arrival data assigned to a specific voxel to determine the multimodal distribution of wave travel times and to derive travel time estimates associated with azimuth sectors (e.g., based on mean, median, mode, standard deviation, kurtosis, and other appropriate statistical accuracy metrics). First arrival data can describe the initial arrival of refracted or latent waves generated by a specific source signal at seismic receiver 226.

[0030] Seismic data can refer to raw time-domain data acquired from seismic surveys (e.g., the acquired seismic data may generate seismic bodies 290). However, seismic data can also refer to data acquired at different time intervals, such as in the case of repeated seismic surveys to obtain time-shifted data. Seismic data can also refer to various seismic properties derived in response to processing the acquired seismic data. Furthermore, in some contexts, seismic data can also refer to depth data or image data. Similarly, seismic data can also refer to processed data (e.g., using seismic inversion operations to generate velocity models of subsurface strata) or migrated seismic images of rock layers within the Earth's surface. Seismic data can also be preprocessed data, such as time-domain data arranged within a two-dimensional shot gather.

[0031] Furthermore, seismic data can include multiple spatial coordinates, such as (x, y) coordinates for each shot point and (x, y) coordinates for each receiver. Thus, seismic data can be grouped into common shot point gathers or common receiver point gathers. In some embodiments, seismic data is grouped based on a common domain, such as a common midpoint (i.e., Xmidpoint = (Xshot + Xrec) / 2, where Xshot corresponds to the shot point location and Xrec corresponds to the seismic receiver location) and a common shot-receiver offset (i.e., Xoffset = Xshot - Xrec).

[0032] In some embodiments, seismic data is processed to generate one or more seismic images. For example, a process called migration can be used to perform seismic imaging. In some embodiments, migration can transform a pre-processed shot gather from the data domain to the image domain corresponding to depth data. In the data domain, seismic events in the shot gather can represent subsurface seismic events recorded in field surveys. In the image domain, the seismic events in the migrated shot gather can represent subsurface geological interfaces. Similarly, various types of migration algorithms can be used for seismic imaging. For example, one type of migration algorithm corresponds to wave equation migration (e.g., one-way wave equation migration, reverse time migration, etc.). In wave equation migration, seismic gathers can be analyzed by: 1) performing forward modeling of the seismic wavefield via mathematical modeling, starting from a synthetic source wavelet and velocity model; 2) performing backpropagation of the seismic data via mathematical modeling using the same velocity model; 3) performing cross-correlation of the seismic wavefield based on the results of the forward modeling and backpropagation; and 4) applying imaging conditions during the cross-correlation to generate seismic images at each time step. Under the fundamental assumption that the source wavefield represents the downlink wavefield and the receiver wavefield represents the uplink wavefield, the imaging conditions can be determined by estimating the cross-correlation between the source and receiver wavefields to determine how the actual image is formed. For example, in Kirchhoff and beamforming methods, the imaging conditions may include the sum of the contributions made by the input data channels after they have been partially unfolded along various isochronous lines (e.g., using the principles of constructive and destructive interference to form the image).

[0033] Furthermore, in some embodiments, seismic data is processed to generate one or more seismic images for model building. Seismic imaging can be performed using an iterative process known as inversion. Inversion transforms a pre-processed shot gather from the data domain to the image domain corresponding to the depth data. Similarly, various types of inversion algorithms can be used for seismic imaging. For example, data domain reflection travel-time inversion is performed to determine a kinematically accurate background model. As another example, full waveform inversion is performed to determine a high-resolution model. Because full waveform inversion suffers from the period jump problem, many applications are used to avoid the period jump effect during full waveform inversion. For example, full waveform inversion uses low-frequency data (e.g., the envelope or intensity of the data, artificial low-frequency data) obtained from the acquired seismic data using mathematical operations (e.g., Laplace-Fourier transform).

[0034] Furthermore, as another example, generating high-quality background property models relaxes the requirements for low-frequency data. Various methods exist for determining kinematically accurate background property models. Specifically, image-domain algorithms (e.g., migration velocity analysis (MVA)) determine background property models by maximizing the stacking capability of flat or inclined common imaging point gathers (CIGs) in the source-receiver offset domain or angular domain. Similarly, data-domain algorithms (e.g., reflection travel time inversion, reflection waveform inversion) determine background property models by matching travel time and / or waveform information between inverse migration data and observed seismic data.

[0035] Furthermore, image domain inversion algorithms are generally more robust than data domain inversion algorithms because they are less sensitive to poor initial models and low-frequency deficiencies. However, compared to data domain inversion algorithms, image domain inversion algorithms typically achieve lower-resolution models with greater computational cost and memory requirements. Therefore, a data domain inversion algorithm is needed to construct high-quality background models for seismic imaging and model building.

[0036] Continuing with seismic imaging, it can be near the end of the seismic data workflow before the seismic interpreter performs analysis. The seismic interpreter can then derive an understanding of the subsurface geology from one or more final migration images. To confirm whether a particular seismic data workflow accurately simulates the subsurface, normal time-of-flight (NMO) stacks can be generated, which consist of multiple NMO gathers with amplitudes sampled from a common midpoint (CMP). In particular, NMO correction can be based on a seismic imaging approximation that determines reflection travel times. However, in cases of complex subsurface geology with large inhomogeneities in particle velocities, or when seismic surveys are not acquired on a horizontal plane, NMO stacking results may fail to accurately indicate the subsurface geology. Ocean-bottom-node surveys and seismic surveys on rough terrain are examples of how NMO stacking results may fail to describe the subsurface geology.

[0037] Although Figure 2The diagram generally shows seismic traces with zero source-receiver offset, but these traces can be stacked, migrated, and / or used to generate attribute volumes derived from the underlying traces. For example, an attribute volume could be a dataset of seismic volumes subjected to one or more processing techniques, such as amplitude-versus-offse (AVO) processing. In AVO processing, seismic data can be classified based on variations in reflection amplitude caused by the presence of hydrocarbon accumulation in the subsurface strata. Using AVO, seismic properties of the subsurface interface can be determined based on the dependence of detected seismic reflection amplitude on the seismic energy incident angle. This AVO processing can determine the normal incidence coefficient and / or gradient components of seismic reflections. Similarly, seismic data can also be processed based on the vertices of pressure waves. In particular, vertices can be used as data gather points to classify the first arrival of seismic data records or traces into multiple source-receiver offset cells based on survey dimensional data (e.g., the xy position of seismic receiver 226 on the Earth's surface 230). These elements can include different numbers of channels and / or different coordinate dimensions.

[0038] Turning to seismic interpreter 261, seismic interpreter 261 may include hardware and / or software having functionality for storing seismic body 290, well logs, core sample data, and other data for seismic data processing, well data processing, training operations, and other corresponding data processes. In some embodiments, seismic interpreter 261 may include similar features to those described below. Figure 12 The computer system 1202 described in the corresponding description. While a seismic interpreter may refer to one or more computer systems used to perform seismic data processing, it may also refer to a human analyst who performs seismic data processing in conjunction with a computer. Although seismic interpreter 261 is shown at a seismic survey site, in some embodiments, seismic interpreter 261 may be located away from the seismic survey site.

[0039] Continuing with the discussion of seismic interpreter 261, seismic interpreter 261 may include hardware and / or software having the capability to perform one or more simulations using one or more components (e.g., forward migration operator 271, adjoint source operator 272, gradient solver 273, piecewise dynamic image warping function 274) through data domain reflection traveltime inversion to analyze seismic data and one or more subsurface strata. For example, seismic interpreter 261 may use the one-way wave equation to generate a sensitivity kernel to handle complex geological environments, such as low-velocity zones. The one-way wave equation is a partial differential equation whose solution includes only waves propagating in one direction. As another example, seismic interpreter 261 may use seismic data to generate a property model of interest (e.g., a velocity model) using data domain reflection traveltime inversion. Seismic interpreter 261 may iteratively update the property model during the inversion process using the forward migration operator 271 and the adjoint source operator 272. The forward migration operator 271 performs numerical simulations based on forward modeling of the one-way wave equation to generate synthetic data for the forward wave field and predictions used in the property model. The adjoint source operator 272 performs numerical simulations based on one-way wave equation modeling to generate the adjoint wave field for the property model. The adjoint source operator can be constructed by analyzing the explicit matrix formulas of forward propagation. For example, the forward migration operator and the adjoint source operator can be implemented using a finite-difference time-domain (TDFD) scheme.

[0040] Furthermore, the gradient solver (e.g., gradient solver 273) iteratively updates the model of the property of interest by numerically solving partial differential equations or optimization problems in a data domain reflection travel-time inversion. For example, the gradient solver can integrate the cross-correlation between the derivative of the source-side wavefield (e.g., the forward wavefield) and the receiver-side wavefield (e.g., the adjoint wavefield) over time up to the maximum recording time to determine the gradient of the current residual as defined by the mismatch function. As another example, the gradient solver can use the gradient of the current residual and the previous search direction to determine the conjugate gradient as the search direction for the current iteration. In some embodiments, for a particular system of linear equations with positive definite, large, and sparse matrices, the conjugate gradient algorithm is a direct method for finding an exact numerical solution after a finite number of iterations. Similarly, the conjugate gradient algorithm can provide a unique solution for a quadratic function. For example, the conjugate gradient algorithm can be applied to numerically solve optimization problems in partial differential equations or least squares optimization problems. At each iteration, the conjugate gradient algorithm can determine the search direction (e.g., the conjugate gradient) to find the final solution of the property model, which is conjugate with the gradient of the current residual defined by the mismatch function and the previous search direction.

[0041] Furthermore, piecewise dynamic image warping functions (e.g., piecewise dynamic image warping function 274) can include point-by-point segmentation to piecewise matching to determine the time shift between the back-migration data and the acquired seismic data. Dynamic image warping functions are multidimensional applications based on dynamic time warping (DTW), which uses local squeezing or stretching suitable for extracting the time shift between two signals to align two one-dimensional (1D) time series. Because the alignment error is defined based on amplitude matching, dynamic time warping functions are amplitude-sensitive algorithms. Various dynamic time warping functions (e.g., differential dynamic time warping, smoothed dynamic time warping) are used to match signal trends and noise signals. However, due to the lack of horizontal constraints, the time shift derived from dynamic time warping functions still contains significant horizontal inconsistencies.

[0042] Furthermore, the piecewise dynamic image warping function calculates the cumulative distance between two signals for the superimposed windowed polynomial fitting results of each segment. Therefore, for data domain reflection traveltime inversion based on the one-way wave equation, the piecewise dynamic image warping function enhances the robustness of signal alignment and the reliability of time shift estimation. Compared to the mismatch results derived by conventional dynamic image warping (DIW), the piecewise dynamic image warping function can handle noisy images with large and rapidly changing shifts. In some embodiments, the seismic interpreter 261 can apply one or more piecewise dynamic image warping functions to determine the time shift between the back-migration data and the acquired seismic data (e.g., see [reference to previous documentation]). Figure 3 (the corresponding description) to satisfy a predetermined mismatch function or other predetermined criteria.

[0043] In some embodiments, the seismic interpreter uses a positive migration operator and / or an adjoint source operator to determine the positive wavefield and / or adjoint wavefield in order to update the property model. For example, the property model may correspond to a model describing property values ​​such as anisotropy, attenuation, density, P-wave velocity, and / or S-wave velocity. Similarly, the complexity of the property model may be associated with the computational cost of updating the property model using the positive wavefield and / or adjoint wavefield. In some embodiments, the seismic interpreter applies one or more simulation algorithms (e.g., finite difference simulation algorithms) to determine the migration operator based on the one-way wave equation.

[0044] Turning Figure 3 , Figure 3 A flowchart according to one or more embodiments is shown. Specifically, Figure 3 A general approach is described to generate a model of the property of interest based on data domain reflection travel time inversion using piecewise dynamic image warping functions and one-way wave equations. Figure 3 One or more boxes in the structure can be formed by, for example Figure 1 and Figure 2 This is performed by one or more components described herein (e.g., seismic interpreter 261). Although Figure 3 The boxes in the document are presented and described in sequence, but those skilled in the art will understand that some or all of these boxes may be executed in a different order, may be combined or omitted, and may be executed in parallel. Furthermore, these boxes may be executed actively or passively.

[0045] In box 300, seismic data is obtained for a geological (or subsurface) area of ​​interest, according to one or more embodiments. The seismic data may be similar to that described above. Figure 1 and Figure 2 The seismic data described. The geological region of interest may be a portion of a geological region or volume that is desired or selected for further analysis, for example, for the purpose of hydrocarbon exploration of the corresponding reservoir or to enhance future hydrocarbon production or reservoir development.

[0046] In box 305, according to one or more embodiments, a property model is obtained for the geological region of interest. The goal of data domain reflection travel-time inversion is to find the optimal model that kinematically matches the back-migration data and the acquired seismic data. For example, data domain reflection travel-time inversion determines the updated model by minimizing a predetermined mismatch function (Equation 1), which is based on the back-migration data p(x r , t; x s ) and the obtained seismic data d(x r , t; x s The time shift between (Equation 2).

[0047] E=∫∫Δτ(x r , t; x s ) 2 dx r Equation 1

[0048] Δτ(x r , t; x s )=F[p(x r , t; x s ), d(x r , t; x s Equation 2

[0049] Where E is the predetermined mismatch function, Δτ is the observed seismic data d(x) r , t; x s ) and inverse offset data p(x) r , t; x s The time shift between t and x is calculated by the dynamic warping operator F, where t is the travel time and x is the time shift between t ... r and x s These are the receiver location and the source location, respectively.

[0050] In some embodiments, it is assumed that the background S-wave velocity has no effect on the source-side kinematics of the P-wave velocity model. In some embodiments, the property model is a model of the property of interest, which is updated using data-domain reflection travel-time inversion. For example, the property model may describe the reflectivity of P-waves at different regions beneath the surface. Specifically, the initial property model may have "0" values ​​at different regions beneath the surface, or the values ​​may be known from previous seismic data processing. As another example, the property model describes the P-wave velocity at different regions beneath the surface. Specifically, the initial property model is a smoothed background velocity model, or the model may be known from previous seismic data processing. Data-domain reflection travel-time inversion can update the property model using the adjoint source operator during the velocity model construction process.

[0051] In block 310, a positive migration operator is determined based on a property model and a one-way wave equation, according to one or more embodiments. The positive migration operator generates positive wave field and reverse migration data based on the property model. The positive migration operator can be represented in matrix form as a linear system, where velocity perturbations are represented as a vector of model parameters in an optimization problem to be iteratively solved by a gradient solver.

[0052] In block 315, according to one or more embodiments, the predicted data is determined based on a property model and a positive offset operator. For example, the predicted data can be determined based on the inverse offset equation (Equation 3) using a depth domain offset portion (e.g., a reflectivity model) and a Green's function extrapolated from the one-way wave equation.

[0053] p(x r , t; x s )=∫G(x′,t;x s )*s(t;x s )*G(x r Equation 3: z = 0, t; x′)m(x′)dx′ Where G is the Green's function extrapolated from the one-way wave equation, the symbol * denotes time convolution, m(x′) is the reflection coefficient at position x′ and is represented by the associated offset profile, z is the depth, t is the travel time, and x r and x s These are the receiver location and the source location, respectively.

[0054] In box 320, according to one or more embodiments, a piecewise dynamic image warping function is used to determine the time shift between the predicted data and the acquired seismic data. The piecewise dynamic image warping function improves the accuracy of dynamic image warping functions in handling strong noise in the input data. The mismatch function of the piecewise dynamic image warping function applies point-by-point piecewise to piecewise matching (Equations 4 and 5) to calculate the processed inverse migration data p′(x r ,t) and processed observed seismic data d′(xr The alignment error between images is e[x, t, τ(t)]. Therefore, a piecewise dynamic warping approximation solution (Equations 6 and 7) can be obtained by applying predetermined constraints to the time shift. Even for noisy data, the piecewise dynamic image warping function can accurately estimate the time shift between images. It solves the following optimization problem:

[0055]

[0056]

[0057] |τ(x,t)-τ(x,t-1)|≤∈ t Equation 6

[0058] |τ(x,t)-τ(x-1,t)|≤∈ x Equation 7

[0059] Where p′(x r ,t) and d′(x r ,t) represents the result of the superimposed polynomial fitting for each segment, and t′ represents the time dimension of each segment. ∈ t and ∈ x These are the predetermined constraints for time mismatch along the time and signal position directions, respectively. δt represents the half-segment length. n x and n t These are the image dimensions in the horizontal and vertical directions, respectively.

[0060] For the 1D case, the piecewise dynamic image warping function is also a piecewise dynamic time function. The piecewise dynamic time warping function can be applied in discretized form in two steps. For example, in the first step, the piecewise dynamic time warping function uses polynomial fitting to approximate the input channels f(t) and g(t) piecewise to form new signals f′(t) and g′(t) (Equations 8 and 9). In the second step, the piecewise dynamic time warping function uses f′(t) and g′(t) as new input channels to determine the time shift Δl(t) in the optimization problem (Equation 10) based on the alignment error e1 (Equation 11) and predetermined constraints (Equation 12).

[0061]

[0062]

[0063]

[0064]

[0065] |l(t)-l(t-1)|≤∈ t Equation 12

[0066] Where n is the signal length, F i，t Indicates in [t i -δt1,t i The polynomial fitting operator for the i-th segment within the range of +δt1], t i It is the midpoint of the i-th segment with a width of 2δt1+1. Δl(t) is the desired time shift, l(t) is the integer lag, δt2 represents the half-segment length in the second step, t′ represents the time dimension of each segment, e1 is the alignment error, ∈ t It is a threshold that limits the shift in the time direction t so that it neither decreases nor increases too quickly.

[0067] In block 325, according to one or more embodiments, based on a property model and a one-way wave equation, the determined time shift is used to determine the adjoint source operator for data domain reflection travel time inversion. For example, the adjoint source operator used for inversion can be derived from a predetermined connection function (Equation 13) and a mismatch function E relative to the inverse offset data p(x). r , t; x s The partial derivative of ) (Equation 15) is used to determine this. For the correct shift Δτ, the predetermined connection function reaches its minimum (Equation 14). Therefore, by substituting Equation 14 into Equation 15, the partial derivative of ) can be used to determine the partial derivative of ) based on the inverse offset data p(x). r , t; x s ) and observed seismic data d(x r , t; x s The time shift between (Equation 13) determines the adjoint source operator (Equation 16) for the mismatch function.

[0068] c(x r ,τ;x s )=∫[p(x r , t; x s )-d(x r , t+τ(t); x s )] 2 Equation 13

[0069]

[0070]

[0071]

[0072] Where c(x) r ,τ;x s ) is from source x s to receiver x r The predefined connection function, at time t, reverses the data p(x) r , t; x s) and observed seismic data d(x r , t+τ(t); x s There is a shift τ(t) between them. Δτ is the first derivative of the predetermined connection function c with respect to time t. Δτ is the time shift between the two signals. t is the travel time. It is relative to the inverse offset data p(x) r , t; x s The first derivative of ). It is the first derivative with respect to the time shift Δτ. It is the observed seismic data d(x) r , t+τ(t); x s The first derivative with respect to time t. It is the observed seismic data d(x) r , t+τ(t); x s The second derivative with respect to time t.

[0073] Furthermore, data domain reflection traveltime inversion utilizes the adjoint source operator used for the one-way wave equation of acoustic waves. For example, the adjoint source operator can be used to determine the gradient, which is derived from a predetermined mismatch function in the data domain reflection traveltime inversion algorithm. For instance, the adjoint source is used in backpropagation to compute the gradient to update the model of the property of interest.

[0074] In block 330, according to one or more embodiments, a gradient solver is determined in data domain reflection traveltime inversion using a forward migration operator and an adjoint source operator. Based on the time shift determined by using a piecewise dynamic image warping function, the forward migration operator, and the adjoint source operator, the gradient solver can determine an update to the property model of interest by minimizing a predetermined mismatch function (Equation 1). In particular, the gradient solver can back-project the adjoint source along the reflection path. The property model update can be determined by a process consisting of the following steps: 1) backpropagating and constructing a new gradient using the time shift between the back-migration data and the observed seismic data, and the gradient solver based on the adjoint source operator; 2) performing a line search to derive the step size; and 3) summing the gradient scaled by the derived step size and the property model from the previous iteration.

[0075] In box 335, according to one or more embodiments, a gradient solver is used to update the property model based on the adjoint source operator. For example, the updated property model can be determined by adding an update to the property model derived from the gradient solver to the property model from previous iterations. The convergence of the inversion depends on gradient preprocessing, the quality of the initial model, and the iteration stopping criterion. Convergence is enhanced by preprocessing the gradient using an illumination term generated by projecting back-migrated data along the wavepath.

[0076] In block 340, according to one or more embodiments, it is determined whether the updated property model has converged to a predetermined criterion. For example, the predetermined criterion could be a mismatch function (e.g., expressed in Equation 1). For example, convergence can be determined when the value of the mismatch function is less than the predetermined criterion (e.g., 5% of the initial mismatch function value). As another example, the maximum number of iterations reaches the predetermined criterion (e.g., the value "25"). If the property model has not yet converged to the predetermined criterion, the process can proceed to block 315. If the property model has converged to the predetermined criterion, the process can proceed to block 345.

[0077] In box 345, according to one or more embodiments, a final property model is output for full waveform inversion, and a seismic image is generated based on the final property model for the geological area of ​​interest. For example, the final property model can be used as an initial model for full waveform inversion to derive a high-resolution model. As another example, the seismic image can be a P-wave image of the subsurface reflectance based on the final property model of interest after one or more iterations. One or more post-processing procedures can be applied to the seismic image to further enhance image resolution and / or geological structure continuity. Thus, the high-resolution model and seismic image provide spatial and depth mapping of subsurface strata for various practical applications, such as predicting hydrocarbon deposition, predicting well paths for geological guidance, etc.

[0078] Turn Figure 4 , Figure 4 Examples of iterative data domain reflection time-lapse inversion based on one or more embodiments are provided. Figure 4As shown, the data domain reflection traveltime inversion applies the forward migration operator X 410, the adjoint source operator Y 420, the gradient solver Z 430, and the piecewise dynamic image warping function G 421 to obtain the updated property model X 440. Using property model A 411 and the synthetic source wavelet C 413, the seismic interpreter applies the forward migration operator X 410 to determine the forward wavefield and the predicted data obtained using the migration function (e.g., the inverse migration data B 412). The seismic interpreter applies the piecewise dynamic image warping function 421 to determine the time shift between the inverse migration data B 412 and the seismic data E 422. Similarly, the seismic interpreter applies the adjoint source operator Y 420 to determine the adjoint wavefield using the time shift obtained by using the piecewise dynamic image warping function 421, the inverse migration data B 412, the seismic data E 422, and the property model A 411. The seismic interpreter employs a gradient solver Z430 to cross-correlate the forward and adjoint wavefields to obtain gradient F434. Specifically, gradient solver Z430 applies imaging conditions and summations during cross-correlation to generate gradients at each iteration of the data-driven reflection traveltime inversion. Additionally, gradient solver Z430 determines a mismatch function I433 for each iteration, which is stored for the convergence history of the inversion process. Furthermore, the seismic interpreter uses the determined gradients (e.g., gradient F434) to obtain an updated property model (e.g., an updated property model X440), which is used for the next iteration 450 in the data-domain reflection traveltime inversion.

[0079] Go to Figure 5 A, Figure 5 B and Figure 5 C generates two signals from actual data to compare and verify the piecewise dynamic time warping function in the presence of strong random noise. Figure 5 A shows the trace f(t) with a root mean square (RMS) amplitude of “0.24”, which was extracted from field land data. Figure 5 B shows the result of... Figure 5 The channels in A are generated by applying known shifts. Figure 5 C shows how to... Figure 5 The noise-normalized channel g(t) is generated by adding strong random noise (e.g., an amplitude range from -0.5 to 0.5) to the channel in B. Figure 5 C indicates the presence of severe pollution caused by noise and numerous events.

[0080] Go to Figure 6A , Figure 6B , Figure 6C and Figure 6D ,exist Figure 5 A and Figure 5A comparison is made between piecewise dynamic time warping functions and dynamic time warping functions on two tracks in C. The same parameter lags l∈(-200,200) and ∈ are used in four tests. t =1. Figure 6A The study shows a significant mismatch between the time shift determined by the dynamic time warping function 602 and the actual time shift 601. Figure 6B The stark difference between the time shift determined solely by the first step of the piecewise dynamic time warping function 603 and the true time shift 601 is shown. Figure 6C The stark difference between the time shift determined solely by the second step of the piecewise dynamic time warping function 604 and the true time shift 601 is shown. Figure 6D A good match is shown between the time shift determined by the piecewise dynamic time warping function 605 and the actual time shift 601.

[0081] Turning Figure 7A , Figure 7B , Figure 7C , Figure 7D , Figure 7E and Figure 7F Applying dynamic image normalization function to determine Figure 7A The shot gathers in the data are obtained by adding time shifts and random noise with various amplitudes. Figure 7A The time shift between shot gathers generated from the shot gathers in the example. Figure 7A The shot gathers are shown with amplitudes ranging from -1.5 to 1.5 and an RMS amplitude of 0.22. Figure 7B The ground condition with time shift to be determined using the dynamic image warping function is shown. Figure 7C It shows how to... Figure 7B The time shift and random noise with moderate amplitude (e.g., in the range of -0.25 to 0.25) are applied. Figure 7A The cannon point set is obtained by using the cannon point set in the middle. Figure 7D The determined time shift and Figure 7B The ground conditions are very consistent, although there are slight variations in a small area 701 with poor time-shift estimates. Figure 7E It shows how to... Figure 7B Time shifts and random noise with strong amplitudes (e.g., in the range of -0.5 to 0.5) are applied to... Figure 7A The cannon point set is obtained by using the cannon point set in the middle. Figure 7F It shows when with Figure 7B When compared with the actual ground conditions, the determined time shift is degraded by strong noise.

[0082] Go to Figure 8A and Figure 8B The piecewise dynamic image normalization function is applied to determine Figure 7C and Figure 7EThe time shift of the artillery point set in the middle. Figure 8A It shows the target such as Figure 7C The time shift determined by the shot point gather shown is Figure 7B The input time shift matches (even in region 701 where the time shift estimation is poor). Figure 8B It shows the target such as Figure 7E The time shift determined by the shot point gather shown is Figure 7B The input time shift is matched (even for strong random noise). Figure 8A and Figure 8B It is also shown that the piecewise dynamic image warping function can determine the time shift of noisy input data more accurately than the dynamic image warping function because the piecewise dynamic image warping function uses windowed polynomial fitting to suppress noise and further estimates the time shift through piecewise-to-piece matching, which can capture local structural information better than point-to-point matching in the dynamic image warping function.

[0083] Go to Figure 9A , Figure 9B , Figure 9C , Figure 9D , Figure 9E , Figure 9F , Figure 9G and Figure 9H The seismic interpreter uses another example to apply data domain reflection travel time inversion based on piecewise dynamic image warping functions and one-way wave equations to automatically invert the background velocity model on a synthetic dataset of a partial Sigsbee 2A model. Figure 9A A ground-based real-world velocity model (3.75 km × 9 km) with a syncline structure embedded with several shallow low-velocity layers is shown. Figure 9B The ground-based depth-domain migration image is shown. The seismic interpreter used a Ricker wavelet with a peak frequency of 12 Hz to generate 120 sonic shot gathers. Each shot gather had 600 receivers evenly distributed on the surface at 15 m intervals. The maximum recording time was 6.0 s, and the time sampling interval was 3 milliseconds (ms). Figure 9C A uniform initial velocity model with a constant value of 1500 m / s is shown for use in data domain reflection travel time inversion. Figure 9D The image shown is a depth-domain offset image using an initial uniform velocity model. When compared with... Figure 9B Compared to the actual ground conditions shown, the incorrect velocity model used in seismic imaging... Figure 9D The geological structure shown is incorrectly imaged in the wrong location.

[0084] also, Figure 9E A smooth background velocity model is shown after 22 iterations, determined using data domain reflection travel time inversion based on a piecewise dynamic image warping function. Figure 9FIt shows the same as Figure 9E The image shown is a depth-domain offset image associated with the smooth background velocity model. When compared with... Figure 9B When comparing with the actual ground conditions shown, for Figure 9E The smooth background velocity model shown. Figure 9F The geological structures shown are correctly imaged in the correct locations. Similarly, Figure 9F A high-quality background velocity model is determined using data-driven reflection travel time inversion based on a piecewise dynamic image warping function. Figure 9G It shows the use of Figure 9E The background velocity model determined in the original model was used as the initial model, and the high-resolution velocity model was determined by full waveform inversion of the synthetic shot gather with a frequency higher than 4Hz. Figure 9G It shows the relationship with Figure 9A The ground speeds shown exhibit good consistency. Figure 9G The slight mismatch at the bottom left corner may be due to poor under-surface lighting. Figure 9H It shows that based on, Figure 9G The final depth domain offset image of the velocity model shown is very close to that of... Figure 9B The image shown is a real image.

[0085] Go to Figure 10A , Figure 10B , Figure 10C and Figure 10D The seismic interpreter compares the travel time and amplitude matching between the inverse migration data and the observation data. Figure 10A and Figure 10B It shows how to achieve this by using, for example Figure 9C The initial velocity model shown uses a piecewise dynamic image warping function to determine the time shift for two representative shot gathers (e.g., shot gather 40 and shot gather 80 in the synthetic dataset). Figure 10C and Figure 10D It shows how to achieve this by using, for example Figure 9E The final smoothed background velocity model derived from data domain reflection travel-time inversion is shown, used to determine the time shift of the same shot gather by applying a piecewise dynamic image warping function. Comparative results show that data domain reflection travel-time inversion significantly reduces the travel-time difference between the observed data and the inverse migration data, indicating that it provides a high-quality initial model for subsequent full-waveform inversion.

[0086] Go to Figure 11A and Figure 11B , Figure 11A and Figure 11B As shown Figure 11A The observation data shown (e.g., the 60th shot point gather in the synthetic dataset) and as shown Figure 11B The comparison shown is between the inverse offset data, which uses, as... Figure 9GThe final velocity model derived using full waveform inversion is shown. The good match between the inverse migration data and the observed data indicates that the velocity model derived from full waveform inversion has high fidelity, generating near-real shot gathers. Therefore, it is verified that data-domain reflection travel-time inversion using a piecewise dynamic image warping function can provide a good initial velocity model for full waveform inversion, even in the absence of low-frequency data, to construct a high-quality final velocity model.

[0087] The embodiments can be implemented on a computer system. Figure 12 This is a block diagram of a computer system for providing computing functions associated with algorithms, methods, functions, procedures, flows, and programs as described in this disclosure. The computer 1202 shown is intended to encompass any computing device, such as a high-performance computing (HPC) device, server, desktop computer, laptop / notebook computer, wireless data port, smartphone, personal data assistant (PDA), tablet computing device, one or more processors within these devices, or any other suitable processing device, including physical or virtual instances (or both) of computing devices. Additionally, computer 1202 may include a computer comprising: an input device, such as a keypad, keyboard, touchscreen, or other device capable of accepting user information; and an output device that transmits information associated with the operation of computer 1202, including digital data, visual or audio information (or a combination of information); or a GUI.

[0088] Computer 1202 may function as a client, network component, server, database, or other persistent storage, or any other component (or combination of roles) in a computer system for performing the subject matter described in this disclosure. The illustrated computer 1202 is communicatively coupled to network 1230 or the cloud. In some specific implementations, one or more components of computer 1202 may be configured to operate within an environment including a cloud-based environment, a local environment, a global environment, or other environments (or combinations thereof).

[0089] At a higher level, computer 1202 is an electronic computing device capable of receiving, transmitting, processing, storing, or managing data and information associated with the described subject. Depending on some specific implementations, computer 1202 may also include, or be communicatively coupled to, application servers, email servers, web servers, cache servers, streaming media data servers, business intelligence (BI) servers, or other servers (or combinations thereof).

[0090] Computer 1202 may receive requests via network 1230 or cloud from client applications (e.g., executing on another computer 1202) and respond to the requests by processing the received requests in a suitable software application. Additionally, requests may also be sent to computer 1202 from internal users (e.g., from a command console or via other suitable access methods), external or third parties, other automated applications, and any other suitable entity, individual, system, or computer.

[0091] Each component of computer 1202 can communicate using system bus 1203. In some specific implementations, any or all components (hardware or software, or a combination of hardware and software) of computer 1202 can interact with each other or with interface 1204 (or a combination of both) on system bus 1203 using application programming interface (API) 1212 or service layer 1213 (or a combination of API 1212 and service layer 1213). API 1212 may include descriptions of routines, data structures, and object classes. API 1212 may be independent of or dependent on a computer language and refers to a complete interface, a single function, or even a set of APIs. Service layer 1213 provides software services to computer 1202 or other components (whether shown or not) communicatively coupled to computer 1202. The functionality of computer 1202 is accessible to all service consumers using the service layer. Software services (such as those provided by service layer 1213) provide reusable, defined business functions through defined interfaces. For example, the interface may be software written in JAVA, C++, or other suitable languages ​​that provide data in Extensible Markup Language (XML) format or other suitable formats. Although shown as an integrated component of computer 1202, alternative concrete implementations may be shown as API 1212 or service layer 1213 as a separate component relative to or communicatively coupled to other components of computer 1202 (whether shown or not). Furthermore, any part or all of API 1212 or service layer 1213 may be implemented as a submodule or sub-module of another software module, enterprise application, or hardware module without departing from the scope of this disclosure.

[0092] Computer 1202 includes interface 1204. Although in Figure 12While shown as a single interface 1204, two or more interfaces 1204 may be used depending on specific needs, expectations, or a particular implementation of computer 1202. Interface 1204 is used by computer 1202 to communicate with other systems in a distributed environment connected to network 1230. Generally, interface 1204 includes logic coded in software or hardware (or a combination of software and hardware) and operable to communicate with network 1230 or the cloud. More specifically, interface 1204 may include software supporting one or more communication protocols associated with the communication, enabling the hardware of network 1230 or the interface to transmit physical signals both inside and outside of computer 1202.

[0093] Computer 1202 includes at least one computer processor 1205. Although in Figure 12 The computer 1202 is shown as a single computer processor 1205, but two or more processors may be used depending on specific needs, expectations, or a particular implementation of the computer 1202. Generally, the computer processor 1205 executes instructions and manipulates data to perform the operations of the computer 1202 and any algorithms, methods, functions, procedures, flows, and programs as described in this disclosure.

[0094] Computer 1202 also includes memory 1206, which stores data for computer 1202 or other components (or a combination of both) that can be connected to network 1230. For example, memory 1206 may be a database storing data consistent with this disclosure. Although in Figure 12 The memory 1206 is shown as a single memory unit, but two or more memories may be used depending on specific needs, expectations, or a particular implementation of the computer 1202 and the functions described. Although the memory 1206 is shown as an integrated component of the computer 1202, in alternative implementations, the memory 1206 may be external to the computer 1202.

[0095] Furthermore, memory 206 can be a computer-readable recording medium and can be composed of at least one of, for example, ROM (Read-Only Memory), EPROM (Erasable Programmable ROM), EEPROM (Electrically Erasable Programmable ROM), and RAM (Random Access Memory). Memory 1206 can be referred to as a register, cache, main memory (main storage device), etc. Memory 1206 can store programs (program code), software modules, etc., that can be executed to implement the radio communication method according to embodiments of the present invention.

[0096] Application 1207 is an algorithmic software engine that provides functionality (particularly with respect to the functionality described in this disclosure) for a specific need, expectation, or specific implementation of computer 1202. For example, application 1207 can be used as one or more components, modules, applications, etc. Furthermore, although shown as a single application 1207, application 1207 can be implemented as multiple applications 1207 on computer 1202. Additionally, although shown as integrated with computer 1202, in alternative specific implementations, application 1207 may be located external to computer 1202.

[0097] Any number of computers 1202 may exist, associated with or outside the computer system containing computer 1202, wherein each computer 1202 communicates on network 1230. Furthermore, the terms "client," "user," and other suitable terms may be used interchangeably where appropriate without departing from the scope of this disclosure. Moreover, this disclosure envisions a plurality of users using one computer 1202, or a single user using multiple computers 1202.

[0098] In some embodiments, computer 1202 is implemented as part of a cloud computing system. For example, the cloud computing system may include one or more remote servers and various other cloud components, such as cloud storage units and edge servers. In particular, the cloud computing system can perform one or more computing operations without direct, active management by user devices or local computer systems. Thus, the cloud computing system can have different functions distributed across multiple locations from a central server, which can be executed using one or more Internet connections. More specifically, the cloud computing system can operate according to one or more service models, such as Infrastructure as a Service (IaaS), Platform as a Service (PaaS), Software as a Service (SaaS), Mobile Backend as a Service (MBaaS), Artificial Intelligence as a Service (AIaaS), Serverless Computing, and / or Function as a Service (FaaS).

[0099] Although only a few exemplary embodiments have been described in detail above, those skilled in the art will readily understand that many modifications can be made to the exemplary embodiments without substantially departing from the invention. Therefore, all such modifications are intended to be included within the scope of this disclosure as defined by the appended claims. In the claims, the means-plus-function clause is intended to cover structures described herein as performing the listed functions and equivalent structures thereof. Similarly, any means-plus-function clause in the claims is intended to cover actions described herein as performing the listed functions and equivalent actions thereof. The applicant’s explicit intent is not to invoke Section 112(f) of the U.S. Patent Act to limit any claim herein, except those claims that expressly use the phrases “means for” or “steps for” along with associated functions.

[0100] Although this disclosure has been described with respect to a limited number of embodiments, those skilled in the art will recognize, with benefit from this disclosure, that other embodiments can be devised without departing from the scope of the invention disclosed herein. Therefore, the scope of this disclosure should be defined only by the appended claims.

Claims

1. A computer-implemented method, comprising: Obtain seismic data for the subsurface region of interest in the time domain; Obtain a property model for the underground region of interest; Using a computer processor, a multidimensional algorithm is employed to determine the time shift between the predicted data and the acquired seismic data, based on the seismic data and the property model, in order to: Windowed polynomial fitting is applied as input using both predicted and acquired seismic data for signal preprocessing. For each segment, the preprocessed signal of interest is aligned based on point-by-point segmentation to segment matching, and The threshold for the time shift in each dimension; The computer processor is used to determine the adjoint source operator using the time shift and one-way wave equations; The property model is updated using a gradient solver in data domain reflection time-lapse inversion via the computer processor. The computer processor outputs an updated property model for the underground area of ​​interest; and The computer processor uses the updated property model to generate seismic images of the subsurface area of ​​interest.

2. The method according to claim 1, in, The updated property model is used as the initial model for full waveform inversion to update the short-wavelength components of the property model. The accompanying source operator is based on the one-way wave equation of sound waves. The accompanying source operator generates multiple accompanying wavefields based on the seismic data. Specifically, a positive migration operator based on the one-way wave equation of acoustic waves is applied to generate multiple positive wavefields based on the seismic data, and The positive migration operator generates prediction data based on the seismic data and the reflectivity model.

3. The method according to claim 1, Using the computer processor, the seismic images are used to determine the presence of hydrocarbons in the geological region of interest.

4. The method according to claim 1, in, The property model is iteratively updated until the mismatch function corresponding to the time shift difference between the predicted data and the acquired seismic data converges to a predetermined criterion. The gradient solver determines the residual value based on the output of the mismatch function, and The property model is updated based on the residual value.

5. The method according to claim 1, Obtain a velocity model for the geological region of interest. in, The property model is a reflection model, and The velocity model is used to update the property model.

6. The method according to claim 5, further comprising: Seismic data about the geological area of ​​interest is obtained using a seismic survey system; and The velocity model is generated using the seismic data and seismic inversion operations.

7. A system comprising: An earthquake survey system, which includes a seismic source and multiple seismic receivers; as well as A seismic interpreter, including a computer processor, wherein the seismic interpreter is coupled to the seismic survey system, and the seismic interpreter includes the following functions: Obtain seismic data for the subsurface region of interest in the time domain; Obtain a property model for the underground region of interest; Based on the seismic data and the property model, a multidimensional algorithm is used to determine the time shift between the predicted data and the acquired seismic data, so as to: Windowed polynomial fitting is applied as input using both predicted and acquired seismic data for signal preprocessing. For each segment, the preprocessed signal of interest is aligned based on point-by-point segmentation to segment matching, and The threshold for the time shift in each dimension; The time-shift and one-way wave equations are used to determine the adjoint source operator; The property model is updated using a gradient solver in the data domain reflection travel time inversion. Output the updated property model for the subsurface region of interest; and The updated property model is used to generate seismic images for the subsurface region of interest.

8. The system according to claim 7, in, The adjoint source operator is based on the one-way wave equation of sound waves. The accompanying source operator generates multiple accompanying wavefields based on the seismic data. Specifically, a positive migration operator based on the one-way wave equation of acoustic waves is applied to generate multiple positive wavefields based on the seismic data, and The positive migration operator generates prediction data based on the seismic data and the reflectivity model.

9. The system according to claim 7, in, The property model is iteratively updated until the mismatch function corresponding to the time shift difference between the predicted data and the acquired seismic data converges to a predetermined criterion. The gradient solver determines the residual value based on the output of the mismatch function, and The property model is updated based on the residual value.

10. A non-transitory computer-readable medium storing instructions executable by a computer processor, the instructions comprising the following functions: Obtain seismic data for the subsurface region of interest in the time domain; Obtain a property model for the underground region of interest; Based on the seismic data and the property model, a multidimensional algorithm is used to determine the time shift between the predicted data and the acquired seismic data, so as to: Windowed polynomial fitting is applied as input using both predicted and acquired seismic data for signal preprocessing. For each segment, the preprocessed signal of interest is aligned based on point-by-point segmentation to segment matching, and The threshold for the time shift in each dimension; The time-shift and one-way wave equations are used to determine the adjoint source operator; The property model is updated using a gradient solver in the data domain reflection travel time inversion. Output the updated property model for the subsurface region of interest; and The updated property model is used to generate seismic images for the subsurface region of interest.

Images