Automatic Method and Apparatus for Modeling In-Layer Tomography Velocity During Reflection Travel

By comprehensively utilizing seismic and well logging data, the travel time of reflected waves is automatically picked up, and a depth-domain layer velocity model is constructed. This solves the problem of difficulty in picking up the travel time of reflected waves in existing technologies, improves the accuracy of pre-stack depth migration imaging, and realizes high-precision imaging of subsurface structures.

CN122307672APending Publication Date: 2026-06-30CHINA PETROLEUM & CHEMICAL CORP +1

Patent Information

Authority / Receiving Office
CN · China
Patent Type
Applications(China)
Current Assignee / Owner
CHINA PETROLEUM & CHEMICAL CORP
Filing Date
2024-12-28
Publication Date
2026-06-30

Smart Images

  • Figure CN122307672A_ABST
    Figure CN122307672A_ABST
Patent Text Reader

Abstract

This disclosure relates to the field of seismic exploration technology, and particularly to a method and apparatus for automatically acquiring reflection travel time for layer-by-layer tomographic velocity modeling. The method includes: acquiring a seismic time-averaged velocity function constructed from raw seismic data and a well-logging time-averaged velocity function constructed from raw well-logging data; calibrating the seismic time-averaged velocity function using the initial seismic average velocity and well-logging average velocity to obtain a calibrated seismic average velocity model; performing time-depth conversion using the calibrated seismic average velocity model to obtain calibrated depth-domain interpretation horizons; performing stacking velocity inversion layer by layer on the calibrated depth-domain horizons to obtain optimized depth-domain layer velocities; obtaining observation travel time using the stacking velocity as a criterion; and performing tomographic inversion using the calibrated depth-domain interpretation horizons, optimized depth-domain layer velocities, observation travel time, and CMP observation system to construct a depth-domain layer velocity model.
Need to check novelty before this filing date? Find Prior Art

Description

Technical Field

[0001] This disclosure relates to the field of seismic exploration technology, and in particular to a method and apparatus for automatically picking up layer-by-layer tomographic velocity modeling during reflection travel. Background Technology

[0002] As oil exploration targets become increasingly complex and concealed, pre-stack depth migration is being applied more and more frequently in production. Pre-stack depth migration can achieve high-precision seismic imaging of complex subsurface structures, accurately reflecting subsurface structures in areas with complex structures and dramatic lateral velocity variations. However, establishing an accurate velocity model is crucial for pre-stack migration imaging. The correctness or accuracy of the velocity model directly affects the success or failure of the imaging process. Summary of the Invention

[0003] This disclosure provides a method and apparatus for automatically picking up the velocity modeling along the layer during reflection travel, in order to solve the problems existing in the related art.

[0004] In a first aspect, this disclosure provides a method for automatically picking up tomographic velocity modeling along layers during reflection travel, including:

[0005] Obtain the seismic time-averaged velocity function constructed from the original seismic data and the well logging time-averaged velocity function constructed from the original well logging data; wherein, the seismic time-averaged velocity function is used to obtain the initial seismic average velocity of vertical change within the layer at or near the well location, and the well logging time-averaged velocity function is used to obtain the well logging average velocity of vertical change within the layer at the well location.

[0006] Using the initial seismic average velocity and the well logging average velocity, the seismic time-averaged velocity function is calibrated to obtain a calibrated seismic average velocity model; the time-domain interpretation horizon is converted to the depth-domain interpretation horizon using the calibrated seismic average velocity model to obtain the calibrated depth-domain interpretation horizon.

[0007] The stacked velocity inversion is performed layer by layer on the calibrated depth domain to obtain the optimized depth domain layer velocity; and the observation travel time is obtained based on the stacking velocity, wherein the observation travel time is the reflection travel time of the CMP gather at each reflection point;

[0008] The depth domain interpretation horizon, the optimized depth domain layer velocity, the observation travel time, and the CMP observation system are used to perform tomographic inversion and construct a depth domain layer velocity model.

[0009] In some embodiments, the method further includes:

[0010] The raw seismic data is processed to obtain target seismic data; the target seismic data includes root mean square velocity volume, seismic reference surface, pre-stack time migration data volume or stacking data volume and CMP gather;

[0011] The root mean square velocity volume is transformed using the DIX formula to obtain the seismic layer velocity volume and mean velocity volume in the time domain and depth domain.

[0012] The time-depth relationship is obtained by performing time-depth conversion on the mean velocity volume, and the depth domain layer is obtained by performing time-depth conversion on the time domain structural interpretation layer.

[0013] The time-depth relationship, time-domain and depth-domain layer velocities, and average velocities are displayed and cross-analyzed to fit the seismic time-averaged velocity function.

[0014] In some embodiments, the method further includes:

[0015] The raw logging data is processed to obtain target logging data; the target logging data includes sonic transit time or sonic transit time curve, density curve and stratigraphic interpretation results;

[0016] The target seismic data is used to perform crystal seismic calibration on the target well logging data to obtain calibration results;

[0017] From the calibration results, the time-depth relationship, layer velocity, and average velocity at the well location are obtained;

[0018] The time-depth relationship, layer velocity, and average velocity at the well location are displayed and cross-analyzed to fit the logging time-averaged velocity function.

[0019] In some embodiments, the step of scaling the seismic time-averaged velocity function using the initial seismic average velocity and the well logging average velocity to obtain a scaled seismic average velocity model includes:

[0020] The seismic time-averaged velocity function and the well logging time-averaged velocity function are fitted to obtain a time-averaged velocity model at the well location; wherein, the time-averaged velocity model at the well location is used to calculate the time-averaged velocity at the well location;

[0021] The time-averaged velocity at the well location is used to calibrate the time-averaged velocity at the corresponding location or near the well to obtain a calibration factor.

[0022] Perform Bayesian skryz interpolation on the scale factor to obtain the scale factor at all grid points;

[0023] The average velocity is calibrated using the aforementioned scale factor to obtain the calibrated average velocity.

[0024] In some embodiments, the formula for calculating the observation travel time is:

[0025]

[0026] Where T(h) represents the reflection travel time of the CMP gather at each reflection point, and V rms The root mean square velocity is represented by h, which is half the distance between the shot and receiver, i.e., half the offset. T0 represents the time of the time domain interpretation layer corresponding to the depth domain interpretation layer.

[0027] In some embodiments, the method further includes:

[0028] Obtain the system of inversion equations for constructing tomography;

[0029] Solving the tomographic inversion equations yields the depth domain layer velocity update and the depth domain layer depth update.

[0030] The optimized depth domain layer velocity is determined using the depth domain layer velocity update amount and the depth domain layer depth update amount.

[0031] Secondly, this disclosure provides an apparatus for automatically picking up tomographic velocity modeling along layers during reflection travel, comprising:

[0032] The acquisition module is used to acquire the seismic time-averaged velocity function constructed from the original seismic data and the well logging time-averaged velocity function constructed from the original well logging data; wherein, the seismic time-averaged velocity function is used to acquire the initial seismic average velocity of vertical change within the layer at or near the well location, and the well logging time-averaged velocity function is used to acquire the well logging average velocity of vertical change within the layer at the well location.

[0033] The processing module is used to calibrate the seismic time-averaged velocity function using the initial seismic average velocity and the well logging average velocity to obtain a calibrated seismic average velocity model; and to use the calibrated seismic average velocity model to convert the time-domain interpretation horizon to the depth-domain interpretation horizon to obtain a calibrated depth-domain interpretation horizon.

[0034] The processing module is also used to perform stacking velocity inversion on the calibrated depth domain layers one by one to obtain the optimized depth domain layer velocity; and to obtain the observation travel time based on the stacking velocity, wherein the observation travel time is the reflection travel time of the CMP gather at each reflection point;

[0035] The processing module is also used to perform tomographic inversion using the calibrated depth domain interpretation layer, the optimized depth domain layer velocity, the observation travel time, and the CMP observation system to construct a depth domain layer velocity model.

[0036] Thirdly, this disclosure provides a computer device including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the method described in the foregoing aspects.

[0037] Fourthly, this disclosure provides a computer-readable storage medium having a computer program stored thereon that, when executed by a processor, implements the steps of the methods described in the above aspects.

[0038] Fifthly, this disclosure provides a computer program product, including a computer program / instructions that, when executed by a processor, implement the steps of the methods described in the foregoing aspects.

[0039] This disclosure provides a method and apparatus for automatically picking up reflection travel time and modeling tomographic velocity along layers. It acquires a seismic time-averaged velocity function constructed from raw seismic data and a well-logging time-averaged velocity function constructed from raw well-logging data. The seismic time-averaged velocity function is used to obtain the initial seismic average velocity of vertical variation within the layer at or near the well location, while the well-logging time-averaged velocity function is used to obtain the well-logging average velocity of vertical variation within the layer at the well location. The seismic time-averaged velocity function is calibrated using the initial seismic average velocity and the well-logging average velocity to obtain a calibrated seismic average velocity model. The time-domain solution is then obtained using the calibrated seismic average velocity model. The interpretation horizon is converted to a depth-domain interpretation horizon to obtain a calibrated depth-domain interpretation horizon. Stacking velocity inversion is performed layer by layer on the calibrated depth-domain horizon to obtain an optimized depth-domain layer velocity. Based on the stacking velocity, observation travel time is obtained, which is the reflection travel time of the CMP gather at each reflection point. Using the calibrated depth-domain interpretation horizon, optimized depth-domain layer velocity, observation travel time, and the CMP observation system, tomographic inversion is performed to construct a depth-domain layer velocity model. This model integrates seismic and well logging data, and utilizes seismic tectonic interpretation results without the need for reflection point picking, automatically picking the reflection wave travel time, thus achieving automatically picked layer-by-layer tomographic velocity modeling. Attached Figure Description

[0040] The present disclosure will be described in more detail below based on embodiments and with reference to the accompanying drawings:

[0041] Figure 1 A flowchart illustrating an automatic method for modeling tomographic velocity during reflection travel provided in an embodiment of this disclosure;

[0042] Figure 2 A main flowchart of an automatic picking-based method for modeling tomographic migration velocity along layers provided in this disclosure embodiment;

[0043] Figure 3This is a schematic diagram of the structure of an automatic tomography velocity modeling device for picking up reflections during travel, provided in an embodiment of the present disclosure. Detailed Implementation

[0044] To enable those skilled in the art to better understand the technical solutions of this disclosure, and to fully understand and implement the process of how this disclosure applies technical means to solve technical problems and achieve corresponding technical effects, the technical solutions in the embodiments of this disclosure will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this disclosure, not all embodiments. The embodiments of this disclosure and the various features within them can be combined with each other without conflict, and the resulting technical solutions are all within the protection scope of this disclosure. All other embodiments obtained by those skilled in the art based on the embodiments of this disclosure without creative effort should fall within the protection scope of this disclosure.

[0045] It should be noted that the terms "first," "second," etc., in the specification, claims, and accompanying drawings of this disclosure are used to distinguish similar objects and are not necessarily used to describe a specific order or sequence. It should be understood that such data can be interchanged where appropriate so that the embodiments of this disclosure described herein can be implemented in orders other than those illustrated or described herein. Furthermore, the terms "comprising" and "having," and any variations thereof, are intended to cover non-exclusive inclusion; for example, a process, method, system, product, or apparatus that comprises a series of steps or units is not necessarily limited to those steps or units explicitly listed, but may include other steps or units not explicitly listed or inherent to such processes, methods, products, or apparatus.

[0046] It should be noted that the steps shown in the flowchart in the accompanying drawings can be executed in a computer system such as a set of computer-executable instructions, and although a logical order is shown in the flowchart, in some cases the steps shown or described may be executed in a different order than that shown here.

[0047] As oil exploration targets become increasingly complex and concealed, pre-stack depth migration is being applied more and more frequently in production. Pre-stack depth migration can achieve high-precision seismic imaging of complex subsurface structures, accurately reflecting subsurface structures in areas with complex structures and dramatic lateral velocity variations. However, establishing an accurate velocity model is crucial for pre-stack migration imaging. The correctness or accuracy of the velocity model directly affects the success or failure of the imaging process.

[0048] Velocities can be obtained from well logging data, core measurement data, and seismic data. Velocity acquisition from seismic data is generally referred to as velocity analysis, which includes stacking velocity analysis, migration velocity analysis, layer velocity analysis, and tomographic inversion. As their names suggest, different velocity analysis methods yield different velocities, and their underlying principles, required data, velocity properties, and applications also differ. Stacking velocity analysis uses Common Middle Point (CMP) gathers as input, obtaining the time-domain velocity that best achieves the stacking effect. Migration velocity analysis, depending on its migration algorithm, can use different input data, such as common shot datasets and common offset datasets, and the resulting velocity model can be a time-domain velocity (pre-stack time migration) or a depth-domain velocity model (pre-stack depth migration). Layer velocity analysis generally requires CMP gathers and reflection interfaces as input, obtaining a depth-domain layer velocity model. Tomographic inversion can use reflection travel time (reflection tomography) to obtain a depth-domain velocity model, or it can simultaneously use travel time and amplitude information to perform waveform inversion to obtain a depth-domain velocity model. Among these methods, stacking velocity analysis is relatively simple and efficient, and therefore the most widely used, but its accuracy is relatively low, and it only yields velocities in the time domain. Migration velocity analysis produces velocity models with higher accuracy than those obtained from stacking velocity analysis, making it more suitable for complex geological conditions.

[0049] Early velocity models were constructed by performing velocity analysis through stacking or pre-stack time migration to obtain root-mean-square velocities. These velocities were then converted into layer velocities using the DIX formula, and the initial model velocity volume was obtained through DIX constraint inversion. Tomographic methods were then used, with interactive human-computer interaction, to extract remaining velocities for model optimization. However, stacking velocity analysis is no longer suitable for complex geological conditions such as lateral velocity variations and stratigraphic dips, and therefore cannot meet the accuracy requirements of velocity analysis.

[0050] Therefore, based on the initial velocity model provided by the above method, residual velocity analysis is performed using the sensitivity of pre-stack depth migration to migration velocity errors. The initial residual velocity analysis involves creating a gamma spectrum, followed by manual picking and interpretation of the gamma spectrum, similar to stacking velocity analysis. On one hand, picking and interpreting the gamma spectrum requires excessive manual operation, making it tedious and time-consuming. On the other hand, creating the gamma spectrum utilizes the common imaging point gather at each grid point (X, Y), which is similar to the CMP gather used in stacking velocity analysis. This method operates on a single-point basis and cannot adequately account for the influence of velocities at other points in the velocity model on the current point.

[0051] Tomographic inversion allows for the global adjustment of velocity models. Tomography was first applied in medicine and later in geophysics, but its primary application remains in transmitted waves. Mathematically, tomographic inversion involves solving a system of linear equations: L·ΔS=Δt. Here, L is a matrix whose elements represent the length of a ray within a grid, ΔS represents the change in slowness within each grid (slowness is the reciprocal of velocity), and Δt represents the observed travel time T. obs With calculation of travel time T cal The difference. Among them, T cal This is calculated using the initial velocity model. Although the tomographic inversion method is relatively mature, in practical applications, it is necessary to extract T from seismic observation data. obs In calculating T cal In addition to the initial velocity model, the locations of the excitation and receiving points are also required. This is because, in transmitted waves or first-arrival tomography, the T value is picked up... obs Relatively simple, tomographic inversion methods are mainly applied to transmission wave tomography and first-arrival tomography in seismic exploration. However, for surface seismic exploration, which uses reflected waves as the main effective signal, the target to be explored is the medium-deep layer. First-arrival tomography can only obtain the surface velocity model, while transmission wave tomography has high requirements for the number and spacing of wells. Therefore, these two tomographic inversion methods cannot be applied to the establishment of medium-deep velocity models in seismic exploration.

[0052] Reflection tomography utilizes reflected wave information from mid-to-deep layers, thus it can be used to establish mid-to-deep velocity models. Traditional reflection tomography, like first-arrival tomography and transmission tomography, requires travel time information. However, while the problem of rapidly and accurately obtaining the travel times of first-arrival and transmitted waves has been successfully solved in practical industrial applications, rapidly and accurately obtaining the travel times of reflected waves is extremely difficult. Therefore, reflection tomography has not been as widely used as first-arrival and transmission tomography.

[0053] Unlike traditional data-domain reflection tomography, tomographic migration velocity analysis utilizes common imaging point gathers generated by pre-stack depth migration, resulting in clearer and more visible reflected waves. Furthermore, the velocity model used in pre-stack depth migration can serve as the initial velocity model for tomographic inversion, leading to more stable inversions. Pre-stack depth migration can also generate depth-domain seismic profiles, making reflection point acquisition relatively easy. Moreover, tomographic migration velocity analysis leverages well-established inversion algorithms, and after considering factors such as implementation difficulty, computational complexity, and application effectiveness, it has begun to be widely adopted in industry. However, tomographic migration velocity analysis is based on imaging gathers provided by pre-stack depth migration, and the quality of these gatherings depends on the initial depth-domain velocity model; therefore, the construction of the initial depth-domain velocity model remains crucial.

[0054] Based on this, the embodiments of this disclosure utilize the traditional data domain reflection wave tomography concept, integrating seismic and well logging data. Using seismic tectonic interpretation results, reflection point picking is not required; instead, the travel time T of the reflected waves is automatically picked up by a software program. obs This enables automatic picking of layer-by-layer tomography velocity modeling.

[0055] Example 1

[0056] Figure 1 This is a flowchart illustrating a method for automatically picking up tomographic velocity modeling along layers during reflection travel, provided in an embodiment of this disclosure. Figure 1 As shown, an automatic method for modeling tomographic velocity along layers during reflection travel includes:

[0057] S101, obtain the seismic time-averaged velocity function constructed from the original seismic data and the well logging time-averaged velocity function constructed from the original well logging data; wherein, the seismic time-averaged velocity function is used to obtain the initial seismic average velocity of the vertical change within the layer at or near the well location, and the well logging time-averaged velocity function is used to obtain the well logging average velocity of the vertical change within the layer at the well location.

[0058] S102, using the initial seismic average velocity and well logging average velocity, the seismic time-averaged velocity function is calibrated to obtain the calibrated seismic average velocity model; the time-domain interpretation horizon is converted to the depth-domain interpretation horizon using the calibrated seismic average velocity model to obtain the calibrated depth-domain interpretation horizon.

[0059] S103, perform stacking velocity inversion layer by layer on the calibrated depth domain to obtain the optimized depth domain layer velocity; and obtain the observation travel time based on the stacking velocity, which is the reflection travel time of the CMP gather at each reflection point;

[0060] S104 utilizes the calibrated depth domain to interpret the stratigraphic level, the optimized depth domain layer velocity, the observation travel time, and the CMP observation system to perform tomographic inversion and construct a depth domain layer velocity model.

[0061] According to the technical solution provided in this disclosure, a seismic time-averaged velocity function constructed from raw seismic data and a well-logging time-averaged velocity function constructed from raw well-logging data are obtained. The seismic time-averaged velocity function is used to obtain the initial seismic average velocity of vertical variation within the layer at or near the well location, and the well-logging time-averaged velocity function is used to obtain the well-logging average velocity of vertical variation within the layer at the well location. The seismic time-averaged velocity function is calibrated using the initial seismic average velocity and the well-logging average velocity to obtain a calibrated seismic average velocity model. The calibrated seismic average velocity model is then used to convert the time-domain interpretation horizon to depth. The depth domain interpretation horizon is obtained by measuring the depth domain interpretation horizon. Layer-by-layer stacking velocity inversion is performed on the measured depth domain horizon to obtain optimized depth domain layer velocities. Based on the stacking velocity, observation travel times are obtained, which are the reflection travel times of CMP gathers at each reflection point. Using the measured depth domain interpretation horizon, optimized depth domain layer velocities, observation travel times, and the CMP observation system, tomographic inversion is performed to construct a depth domain layer velocity model. This model integrates seismic and well logging data, and utilizes seismic tectonic interpretation results without the need for reflection point picking, automatically picking reflection wave travel times, thus achieving automatically picked layer-by-layer tomographic velocity modeling.

[0062] Example 2

[0063] Based on the above embodiments, the method further includes:

[0064] The raw seismic data is processed to obtain the target seismic data; the target seismic data includes the root mean square velocity volume, seismic reference surface, pre-stack time migration data volume or stacking data volume and CMP gather;

[0065] The root mean square velocity volume is transformed using the DIX formula to obtain the seismic layer velocity volume and mean velocity volume in the time domain and depth domain.

[0066] The time-depth relationship is obtained by performing time-depth conversion on the mean velocity volume, and the depth domain layer is obtained by performing time-depth conversion on the time domain structural interpretation layer.

[0067] The time-depth relationship, time-domain and depth-domain layer velocities, and average velocities are displayed and cross-analyzed to fit the seismic time-averaged velocity function.

[0068] Specifically, the raw ground seismic data collected in the field is first processed using conventional methods to obtain the target seismic data.

[0069] Here, the conventional processing can be pre-stack time migration (PSTM) processing, stacking processing, or a combination of both.

[0070] The target seismic data may include the obtained root mean square velocity volume, seismic reference surface, and processed seismic data. The processed seismic data includes pre-stack time migration data (or stacking data) and CMP gathers; then, based on the pre-stack time migration data (or stacking data), a fine-grained time-domain structural interpretation is performed to obtain the time-domain structural interpretation results (also known as time-domain structural horizon results).

[0071] Under the constraints of the above time-domain construction interpretation results, the root mean square velocity is transformed using the DIX formula to obtain the time-domain layer velocity model:

[0072]

[0073] Where (x,y) represents the planar coordinates of the grid point, V int V represents the velocity of the time-domain interpretation layer. rms Let T represent the root mean square velocity and T represent the reflection travel time.

[0074] Assuming the seismic reference depth to be processed is D0 = D(x,y,0), then the seismic reference depth of the time-domain interpretation layer and the time of the time-domain interpretation layer can be used to convert the time-domain interpretation layer to the depth-domain interpretation layer layer by layer using the following formula:

[0075]

[0076] Among them, T ow When indicating a one-way trip.

[0077] Furthermore, by interpreting the time of the horizon in the time domain and the corresponding horizon in the depth domain, the initial seismic mean velocity model is calculated using the following formula:

[0078]

[0079] The initial depth domain layer velocity model can be calculated from the time and depth between layers:

[0080]

[0081] The aforementioned initial earthquake mean velocity model, time-domain layer velocity model, and initial depth-domain layer velocity model are all constants within a given seismic interpretation layer. A more refined approach is to assume that the velocity within a given seismic interpretation layer varies with depth or time; that is, the velocity within a given seismic interpretation layer is a function of depth or time. The seismic interpretation layer velocity function is assumed to have the following form:

[0082]

[0083] Where V0(x,y,i) represents the initial velocity (at time i at the top of the layer) of a certain layer, and V(x,y,j) represents the velocity at a depth of Z within a certain layer. j The velocity at time t, where α and β are undetermined coefficients. This velocity can be the layer velocity or the average velocity.

[0084] By fitting the aforementioned seismic interpretation layer velocity function between each seismic interpretation layer, the coefficients α and β within each seismic interpretation layer can be obtained, thus obtaining an initial seismic velocity model with vertical variation within the layer.

[0085] Example 3

[0086] Based on the above embodiments, the method further includes:

[0087] The raw logging data is processed to obtain the target logging data; the target logging data includes sonic transit time or sonic transit time curve, density curve and stratigraphic interpretation results;

[0088] Crystal seismic calibration of target well logging data is performed using target seismic data to obtain calibration results;

[0089] From the calibration results, the time-depth relationship, layer velocity, and average velocity at the well location are obtained;

[0090] The time-depth relationship, layer velocity, and average velocity at the well location are displayed and cross-analyzed to fit the logging time-averaged velocity function.

[0091] Specifically, in this embodiment, after obtaining the raw logging data, the raw logging data can be processed to obtain target logging data. The target logging data may include sonic transit time (AC) or sonic transit time curves, density curves, and stratigraphic interpretation results. Then, the target logging data is calibrated using target seismic data to obtain calibration results. From the calibration results, the time-depth relationship, stratigraphic velocity, and average velocity at the well location are obtained. Finally, the time-depth relationship, stratigraphic velocity, and average velocity at the well location are displayed and cross-analyzed to fit a logging time-averaged velocity function.

[0092] Here, the fitting of the well logging time-averaged velocity function is similar to that of the seismic time-averaged velocity function, and will not be elaborated further.

[0093] Example 4

[0094] Based on the above embodiments, the seismic time-averaged velocity function is calibrated using the initial seismic average velocity and well logging average velocity to obtain a calibrated seismic average velocity model, including:

[0095] The seismic time-averaged velocity function and the well logging time-averaged velocity function are fitted to obtain the time-averaged velocity model at the well location; the time-averaged velocity model at the well location is used to calculate the time-averaged velocity at the well location.

[0096] The time-averaged velocity at the well location is used to calibrate the time-averaged velocity at the corresponding location or near the well to obtain the calibration factor;

[0097] Perform Bayesian skryz interpolation on the scale factor to obtain the scale factor at all grid points;

[0098] The average velocity is calibrated using a scale factor to obtain the calibrated average velocity.

[0099] Specifically, the initial seismic average velocity of vertical change within the layer at or near the well site is calibrated using the well-logging average velocity of vertical change within the layer, thus obtaining a calibration factor:

[0100]

[0101] Among them, V ave_well V represents the average logging velocity. ave This represents the initial average velocity of the earthquake.

[0102] By performing Bayesian skryz interpolation using the average velocity scale factor at the aforementioned well locations, the average velocity scale factor at all grid points can be obtained.

[0103] The initial seismic mean velocity is calibrated using the average velocity scale factor at the grid points to obtain the calibrated seismic mean velocity model.

[0104] The time-domain interpretation horizons are converted to depth-domain interpretation horizons using the calibrated seismic mean velocity model, thus obtaining the calibrated depth-domain interpretation horizons.

[0105] The local dip field of the depth-domain interpretation layer after automatic calculation of the scale: φ x (x,y,z), φ y (x,y,z).

[0106] Example 5

[0107] Based on the above embodiments, the formula for calculating the observation travel time is as follows:

[0108]

[0109] Where T(h) represents the reflection travel time of the CMP gather at each reflection point, and V rmsThe root mean square velocity is represented by h, which is half the distance between the shot and receiver, i.e., half the offset. T0 represents the time of the time domain interpretation layer corresponding to the depth domain interpretation layer.

[0110] Specifically, the time of the time-domain interpretation layer corresponding to the depth-domain interpretation layer is denoted as T0(x,y). The shot-receiver position of each trace in the CMP gather corresponding to each grid point (x,y) is extracted. The travel time on each CMP gather is then calculated according to the following formula, and this travel time is taken as the observation travel time T at the reflection point R(x,y,z) in the depth-domain interpretation layer. obs :

[0111]

[0112] Where h is half the distance between the gun and receiver, i.e., half the offset distance.

[0113] This formula shows the reflection travel time T for a reflection point with a half offset of h on a CMP gather. obs The relationship between h and T is hyperbolic, therefore T can be obtained without picking up this reflection during travel. obs That is, T obs =T(h). However, in more practical situations, when the surface or reflecting interface is not horizontal, or when the layer velocity of a certain overlying stratum has a relatively large lateral variation, or when the layer velocity of a certain overlying stratum exhibits azimuthal anisotropy, or when dynamic correction stretching occurs due to excessive offset, the reflection travel time T of a reflecting point at an upper half offset of h in the CMP gather is... obs The relationship with h is no longer hyperbolic. At this point, the hyperbolic relationship described above can be used to perform NMO correction on the CMP gather to obtain the NMO-corrected CMP gather. Then, on the corrected CMP gather, a window of duration τ is opened centered at T0(x,y) of each reflection point. Next, following the order of increasing offset, the trace data of each trace in the corrected CMP gather with duration τ and starting times from T0(x,y)-τ to T0(x,y)+τ are cross-correlated with the trace of duration τ centered at T0(x,y) with the smallest offset in the CMP gather (or the superposition energy is calculated). Then, the time shift γ (i.e., the cross-correlation criterion or superposition energy criterion) when the cross-correlation value or superposition energy is maximum during the scanning process is recorded. Then, T... obs =T0(x,y)+γ+ΔT NMO .

[0114] For each reflection point R(x, y, z), ray tracing is performed from the reflection point to the ground based on the updated depth-domain horizon, dip angle, and initial depth-domain layer velocity model, according to the shot-receiver relationship in the CMP gather for each reflection point. First, the normal direction of the reflection interface at the reflection point is obtained from each reflection point R(x, y, z) and its corresponding local dip field φ(x, y, z). Then, ray tracing is performed on both sides of the normal at an angle θ equal to the normal. The exit points of these two rays on the ground are denoted as the excitation point and receiver point, respectively. Ray pairs whose distances from the excitation point and receiver point to the shot-receiver point are less than a certain threshold are selected as valid rays. This allows the calculation of the travel time and ray path of each trace in the CMP gather for each reflection point, forming the ray path L and the calculated travel time T for the current reflection point. cal During this process, the layer velocity of the current layer is scanned within a certain range, centered on the initial layer velocity in the depth domain of the current layer. During the scan, the ray path L and travel time T are calculated again in the manner described above. cal Choose the T obtained in step (9) obs The closest T cal And the current layer velocity, as the new layer velocity for the current layer; in this way, a new depth domain layer velocity model is obtained, along with the corresponding ray path L and the calculated travel time T. cal This serves as the initial velocity model for the next step of tomographic inversion.

[0115] Example 6

[0116] Based on the above embodiments, the method further includes:

[0117] Obtain the system of inversion equations for constructing tomography;

[0118] Solve the tomographic inversion equations to obtain the depth domain layer velocity update and depth domain layer depth update.

[0119] The optimized depth domain layer velocity is determined by using the depth domain layer velocity update amount and the depth domain layer depth update amount.

[0120] Specifically, a set of tomographic inversion equations is constructed based on the above calculation results:

[0121] ΔT=T obs -T cal =ΔS·L+(cosθ) r +cosθ s )·S·ΔZ

[0122] Where S is the reciprocal of the layer velocity, derived from the layer velocity model obtained in step (10); θ s and θ rThese are the angles between the ray reaching the gun receiver and the normal to the reflecting interface.

[0123] Solving the above system of equations yields the updated velocity ΔS and the updated reflection interface depth ΔZ, thus updating the initial velocity model. There are many methods for solving the system of equations, such as the conjugate gradient method and the least squares QR decomposition method.

[0124] Through the above steps, the depth domain layer velocity model can be finally obtained by layer-by-layer inversion. It can be used as the initial velocity model for pre-stack depth migration. However, during the tomographic inversion process, it is not necessary to pick up the reflection interface or manually pick up the reflection travel time on the pre-stack gather.

[0125] Example 7

[0126] Based on the above embodiments, this embodiment provides an application example.

[0127] Figure 2 A main flowchart of an automatic picking-based layer-by-layer tomographic migration velocity modeling method provided in this disclosure is shown below. Figure 2 As shown, it specifically includes the following:

[0128] (1) Perform routine processing on the raw ground seismic data collected in the field to obtain the target seismic data.

[0129] Here, the conventional processing can be pre-stack time migration (PSTM) processing, stacking processing, or a combination of both.

[0130] The target seismic data may include the obtained root mean square velocity volume, seismic reference surface, and processed seismic data. The processed seismic data includes pre-stack time migration data (or stacking data) and CMP gathers; then, based on the pre-stack time migration data (or stacking data), a fine-grained time-domain structural interpretation is performed to obtain the time-domain structural interpretation results (also known as time-domain structural horizon results).

[0131] (2) Under the constraints of the above time-domain construction interpretation results, the root mean square velocity is transformed using the DIX formula to obtain the time-domain layer velocity model:

[0132]

[0133] Where (x,y) represents the planar coordinates of the grid point, V int V represents the velocity of the time-domain interpretation layer. rms Let T represent the root mean square velocity and T represent the reflection travel time.

[0134] Assuming the seismic reference depth to be processed is D0 = D(x,y,0), then the seismic reference depth of the time-domain interpretation layer and the time of the time-domain interpretation layer can be used to convert the time-domain interpretation layer to the depth-domain interpretation layer layer by layer using the following formula:

[0135]

[0136] Among them, T ow When indicating a one-way trip.

[0137] Furthermore, by interpreting the time of the horizon in the time domain and the corresponding horizon in the depth domain, the initial seismic mean velocity model is calculated using the following formula:

[0138]

[0139] The initial depth domain layer velocity model can be calculated from the time and depth between layers:

[0140]

[0141] The aforementioned initial earthquake mean velocity model, time-domain layer velocity model, and initial depth-domain layer velocity model are all constants within a given seismic interpretation layer. A more refined approach is to assume that the velocity within a given seismic interpretation layer varies with depth or time; that is, the velocity within a given seismic interpretation layer is a function of depth or time. The seismic interpretation layer velocity function is assumed to have the following form:

[0142]

[0143] Where V0(x,y,i) represents the initial velocity (at time i at the top of the layer) of a certain layer, and V(x,y,j) represents the velocity at a depth of Z within a certain layer. j The velocity at time t, where α and β are undetermined coefficients. This velocity can be the layer velocity or the average velocity.

[0144] By fitting the aforementioned seismic interpretation layer velocity function between each seismic interpretation layer, the coefficients α and β within each seismic interpretation layer can be obtained, thus obtaining an initial seismic velocity model with vertical variation within the layer.

[0145] (3) Using logging velocity and logging depth, the logging depth is converted into time to obtain the logging average velocity V. ave_well Then, similar to step (2), the time-domain layer velocity and time-domain average velocity at the well site can be obtained. Simultaneously, similar to step (2), assuming the velocity function is of the form... In this case, the logging velocity at the well site is fitted to obtain the coefficients α and β within each layer at the well site, thereby obtaining the logging velocity of vertical variation within the layer.

[0146] (4) The initial seismic average velocity of vertical change within the layer at or near the well site is calibrated using the well logging average velocity of vertical change within the layer at the well site, and the calibration factor is obtained:

[0147]

[0148] (5) Use the average velocity scale factor at the above well location to perform Bayesian skrygian interpolation to obtain the average velocity scale factor at all grid points.

[0149] (6) The initial seismic average velocity is calibrated using the average velocity scale factor on the grid points to obtain the calibrated seismic average velocity model.

[0150] (7) Use the calibrated seismic mean velocity model to convert the time domain interpretation horizon to the depth domain interpretation horizon to obtain the calibrated depth domain interpretation horizon.

[0151] (8) Automatically calculate the local dip field of the depth domain interpretation layer after calibration: φ x (x,y,z), φ y (x,y,z).

[0152] (9) Denote the time of the time domain interpretation layer corresponding to the depth domain interpretation layer as T0(x,y). Extract the shot-receiver position of each trace in the CMP gather corresponding to each grid point (x,y). Calculate the travel time on each CMP gather according to the following formula, and use this travel time as the observation travel time T at the reflection point R(x,y,z) in the depth domain interpretation layer. obs :

[0153]

[0154] Where h is half the distance between the gun and receiver, i.e., half the offset distance.

[0155] This formula shows the reflection travel time T for a reflection point with a half offset of h on a CMP gather. obs The relationship between h and T is hyperbolic, therefore T can be obtained without picking up this reflection during travel. obs That is, T obs =T(h). However, in more practical situations, when the surface or reflecting interface is not horizontal, or when the layer velocity of a certain overlying stratum has a relatively large lateral variation, or when the layer velocity of a certain overlying stratum exhibits azimuthal anisotropy, or when dynamic correction stretching occurs due to excessive offset, the reflection travel time T of a reflecting point at an upper half offset of h in the CMP gather is... obsThe relationship with h is no longer hyperbolic. At this point, the hyperbolic relationship described above can be used to perform NMO correction on the CMP gather to obtain the NMO-corrected CMP gather. Then, on the corrected CMP gather, a window of duration τ is opened centered at T0(x,y) of each reflection point. Next, following the order of increasing offset, the trace data of each trace in the corrected CMP gather with duration τ and starting times from T0(x,y)-τ to T0(x,y)+τ are cross-correlated with the trace of duration τ centered at T0(x,y) with the smallest offset in the CMP gather (or the superposition energy is calculated). Then, the time shift γ (i.e., the cross-correlation criterion or superposition energy criterion) when the cross-correlation value or superposition energy is maximum during the scanning process is recorded. Then, T... obs =T0(x,y)+γ+ΔT NMO .

[0156] (10) For each reflection point R(x,y,z) in step (9), ray tracing is performed from the reflection point to the ground based on the updated depth domain horizon and dip angle and the initial depth domain layer velocity model, according to the shot-receiver relationship in the CMP gather of each reflection point. First, the normal direction of the reflection interface at the reflection point is obtained from each reflection point R(x,y,z) and its corresponding local dip field φ(x,y,z). Then, ray tracing is performed on both sides of the normal at the same angle θ as the normal. The emission points of these two rays on the ground are recorded as the excitation point and the receiver point, respectively. Ray pairs with a distance of less than a certain threshold between the excitation point and the receiver point and the shot-receiver point are selected as effective rays. Thus, the travel time and ray path of each trace in the CMP gather of each reflection point are calculated, forming the ray path L of the current reflection point and the calculated travel time T. cal During this process, the layer velocity of the current layer is scanned within a certain range, centered on the initial layer velocity in the depth domain of the current layer. During the scan, the ray path L and travel time T are calculated again in the manner described above. cal Choose the T obtained in step (9) obs The closest T cal And the current layer velocity, as the new layer velocity for the current layer; in this way, a new depth domain layer velocity model is obtained, along with the corresponding ray path L and the calculated travel time T. cal This serves as the initial velocity model for the next step of tomographic inversion.

[0157] (11) Construct a set of tomographic inversion equations based on the above calculation results:

[0158] ΔT=T obs -T cal =ΔS·L+(cosθ) r +cosθ s )·S·ΔZ

[0159] Where S is the reciprocal of the layer velocity, derived from the layer velocity model obtained in step (10); θ s and θ r These are the angles between the ray reaching the gun receiver and the normal to the reflecting interface.

[0160] Solving the above system of equations yields the updated velocity ΔS and the updated reflection interface depth ΔZ, thus updating the initial velocity model. There are many methods for solving the system of equations, such as the conjugate gradient method and the least squares QR decomposition method.

[0161] Through the above steps, the depth domain layer velocity model can be finally obtained by layer-by-layer inversion. It can be used as the initial velocity model for pre-stack depth migration. However, during the tomographic inversion process, it is not necessary to pick up the reflection interface or manually pick up the reflection travel time on the pre-stack gather.

[0162] The foregoing mainly describes the solutions provided by the embodiments of this disclosure. It is understood that, in order to achieve the above functions, the electronic device includes hardware structures and / or software modules corresponding to the execution of each function. Those skilled in the art should readily recognize that, based on the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein, this disclosure can be implemented in hardware or a combination of hardware and computer software. Whether a function is executed in hardware or by computer software driving hardware depends on the specific application and design constraints of the technical solution. Those skilled in the art can use different methods to implement the described functions for each specific application, but such implementation should not be considered beyond the scope of this disclosure.

[0163] This disclosure embodiment can divide the electronic device into functional units according to the above method example. For example, each function can be divided into a separate functional module, or two or more functions can be integrated into one processing module. The integrated module can be implemented in hardware or as a software functional module. It should be noted that the module division in this disclosure embodiment is illustrative and only represents one logical functional division; other division methods may be used in actual implementation.

[0164] By dividing each function into corresponding functional modules, this disclosure provides an apparatus for automatically picking up the velocity modeling along the layer during reflection travel. Figure 3 This is a schematic diagram of a device for automatically picking up and modeling the velocity along a layer during reflection travel, as provided in an embodiment of this disclosure. Figure 3 As shown, the device 300 includes:

[0165] The acquisition module 301 is used to acquire the seismic time-averaged velocity function constructed from the original seismic data and the well logging time-averaged velocity function constructed from the original well logging data; wherein, the seismic time-averaged velocity function is used to acquire the initial seismic average velocity of vertical change within the layer at or near the well location, and the well logging time-averaged velocity function is used to acquire the well logging average velocity of vertical change within the layer at the well location.

[0166] Processing module 302 is used to calibrate the seismic time-averaged velocity function using the initial seismic average velocity and the well logging average velocity to obtain a calibrated seismic average velocity model; and to use the calibrated seismic average velocity model to convert the time-domain interpretation horizon to the depth-domain interpretation horizon to obtain a calibrated depth-domain interpretation horizon.

[0167] The processing module 302 is further configured to perform stacking velocity inversion on the calibrated depth domain layers one by one to obtain the optimized depth domain layer velocity; and to obtain the observation travel time based on the stacking velocity, wherein the observation travel time is the reflection travel time of the CMP gather at each reflection point;

[0168] The processing module 302 is also used to perform tomographic inversion using the calibrated depth domain interpretation layer, the optimized depth domain layer velocity, the observation travel time, and the CMP observation system to construct a depth domain layer velocity model.

[0169] In some embodiments, the processing module 302 is further configured to process the original seismic data to obtain target seismic data; the target seismic data includes root mean square velocity volume, seismic reference surface, pre-stack time migration data volume or stacking data volume and CMP gather; the root mean square velocity volume is converted using the DIX formula to obtain time-domain and depth-domain seismic layer velocity volumes and average velocity volumes; the average velocity volume is used to perform time-depth conversion to obtain time-depth relationships, and time-depth conversion is performed on the time-domain structural interpretation horizons to obtain depth-domain horizons; the time-depth relationships, time-domain and depth-domain layer velocities, and average velocities are displayed and cross-analyzed to fit a seismic time-averaged velocity function.

[0170] In some embodiments, the processing module 302 is further configured to process the original logging data to obtain target logging data; the target logging data includes sonic transit time or sonic transit time curve, density curve and layer interpretation results; the target logging data is calibrated using the target seismic data to obtain calibration results; the time-depth relationship, layer velocity and average velocity at the well location are obtained from the calibration results; the time-depth relationship, layer velocity and average velocity at the well location are displayed and cross-analyzed to fit a logging time-averaged velocity function.

[0171] In some embodiments, the processing module 302 is further configured to fit the seismic time-averaged velocity function and the well logging time-averaged velocity function to obtain a time-averaged velocity model at the well location; wherein, the time-averaged velocity model at the well location is used to calculate the time-averaged velocity at the well location; the time-averaged velocity at the well location is used to calibrate the time-averaged velocity at the corresponding location or near the well to obtain a calibration factor; Bayesian skewkin interpolation is performed on the calibration factor to obtain the calibration factor at all grid points; the calibration factor is used to calibrate the corresponding average velocity to obtain the calibrated average velocity.

[0172] In some embodiments, the formula for calculating the observation travel time is:

[0173]

[0174] Where T(h) represents the reflection travel time of the CMP gather at each reflection point, and V rms The root mean square velocity is represented by h, which is half the distance between the shot and receiver, i.e., half the offset. T0 represents the time of the time domain interpretation layer corresponding to the depth domain interpretation layer.

[0175] In some embodiments, the acquisition module 301 is further configured to acquire the system of tomographic inversion equations;

[0176] The processing module 302 is further configured to solve the tomographic inversion equation set to obtain the depth domain layer velocity update amount and the depth domain layer depth update amount; and to determine the optimized depth domain layer velocity using the depth domain layer velocity update amount and the depth domain layer depth update amount.

[0177] Based on the above embodiments, this embodiment provides a computer device, including a memory, a processor, and a computer program stored in the memory, wherein the processor executes the computer program to implement the steps of the method described in the above embodiments.

[0178] In some embodiments of this example, a computer-readable storage medium is provided, on which a computer program is stored, which, when executed by a processor, implements the steps of the method described in the above embodiments.

[0179] In some embodiments of this example, a computer program product is provided, including a computer program / instructions, which, when executed by a processor, implements the steps of the method described in the above embodiments.

[0180] The processor may include, but is not limited to, one or more processors or microprocessors. Each processor may be implemented as an Application Specific Integrated Circuit (ASIC), Digital Signal Processor (DSP), Digital Signal Processing Device (DSPD), Programmable Logic Device (PLD), Field Programmable Gate Array (FPGA), controller, microcontroller, microprocessor, or other electronic component, for performing the methods in the above embodiments.

[0181] Computer-readable storage media can be implemented by any type of volatile or non-volatile storage device or a combination thereof. Computer-readable storage media may include, but are not limited to, random access memory (RAM), read-only memory (ROM), flash memory, EPROM memory, EEPROM memory, registers, and computer storage media (e.g., hard disks, floppy disks, solid-state drives, removable disks, CD-ROMs, DVD-ROMs, Blu-ray discs, etc.).

[0182] Computer-readable storage media may also store at least one computer-executable program / instruction, such as computer-readable instructions. Computer-readable storage media include, but are not limited to, volatile memory and / or non-volatile memory. Volatile memory may include, for example, random access memory (RAM) and / or cache memory. Computer-readable storage media may include, for example, read-only memory (ROM), hard disk, flash memory, etc. For example, a non-transitory computer-readable storage medium may be connected to a computing device such as a computer, and then, when the computing device executes the computer-readable instructions stored on the computer-readable storage medium, the various methods described above can be performed.

[0183] In addition, the computer device may include (but is not limited to) a data bus, an input / output (I / O) bus, a display, and input / output devices (e.g., keyboard, mouse, speakers, etc.).

[0184] The processor can communicate with external devices via the I / O bus through wired or wireless networks.

[0185] In one embodiment, the at least one computer-executable instruction may also be compiled into or comprise a software product / computer program product, wherein one or more computer-executable instructions are executed by a processor to perform the steps of the various functions and / or methods in the embodiments described herein.

[0186] In the embodiments provided in this disclosure, it should be understood that the disclosed apparatus and methods can also be implemented in other ways. The apparatus embodiments described above are merely illustrative; for example, the flowcharts and block diagrams in the accompanying drawings illustrate the architecture, functionality, and operation of possible implementations of apparatus, methods, and computer program products according to various embodiments of this disclosure. In this regard, each block in a flowchart or block diagram may represent a module, segment, or portion of code containing one or more executable instructions for implementing a specified logical function. It should also be noted that in some alternative implementations, the functions marked in the blocks may occur in a different order than those marked in the drawings. For example, two consecutive blocks may actually be executed substantially in parallel, and they may sometimes be executed in reverse order, depending on the functions involved. It should also be noted that each block in a block diagram and / or flowchart, and combinations of blocks in block diagrams and / or flowcharts, can be implemented using a dedicated hardware-based system that performs the specified function or action, or using a combination of dedicated hardware and computer instructions.

[0187] It should be noted that, in this disclosure, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Without further limitation, an element limited by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes that element.

[0188] While the embodiments disclosed herein are as described above, the foregoing content is merely for the purpose of facilitating understanding of this disclosure and is not intended to limit this disclosure. Any person skilled in the art to which this disclosure pertains may make any modifications and changes in form and detail of the implementation without departing from the spirit and scope of this disclosure; however, the scope of patent protection of this disclosure shall still be determined by the scope defined in the appended claims.

Claims

1. A method for automatic picking of reflect traveltime along-bed tomographic velocity modeling, characterized by, include: Obtain the seismic time-averaged velocity function constructed from the original seismic data and the well logging time-averaged velocity function constructed from the original well logging data; wherein, the seismic time-averaged velocity function is used to obtain the initial seismic average velocity of vertical change within the layer at or near the well location, and the well logging time-averaged velocity function is used to obtain the well logging average velocity of vertical change within the layer at the well location. Using the initial seismic average velocity and the well logging average velocity, the seismic time-averaged velocity function is calibrated to obtain a calibrated seismic average velocity model; the time-domain interpretation horizon is converted to the depth-domain interpretation horizon using the calibrated seismic average velocity model to obtain the calibrated depth-domain interpretation horizon. The stacked velocity inversion is performed layer by layer on the calibrated depth domain to obtain the optimized depth domain layer velocity; and the observation travel time is obtained based on the stacking velocity, wherein the observation travel time is the reflection travel time of the CMP gather at each reflection point; The depth domain interpretation horizon, the optimized depth domain layer velocity, the observation travel time, and the CMP observation system are used to perform tomographic inversion and construct a depth domain layer velocity model.

2. The method of claim 1, wherein, The method further includes: The raw seismic data is processed to obtain target seismic data; the target seismic data includes root mean square velocity volume, seismic reference surface, pre-stack time migration data volume or stacking data volume and CMP gather; The root mean square velocity volume is transformed using the DIX formula to obtain the seismic layer velocity volume and mean velocity volume in the time domain and depth domain. The time-depth relationship is obtained by performing time-depth conversion on the mean velocity volume, and the depth domain layer is obtained by performing time-depth conversion on the time domain structural interpretation layer. The time-depth relationship, time-domain and depth-domain layer velocities, and average velocities are displayed and cross-analyzed to fit the seismic time-averaged velocity function.

3. The method of claim 2, wherein, The method further includes: The raw logging data is processed to obtain target logging data; the target logging data includes sonic transit time or sonic transit time curve, density curve and stratigraphic interpretation results; The target seismic data is used to perform crystal seismic calibration on the target well logging data to obtain calibration results; From the calibration results, the time-depth relationship, layer velocity, and average velocity at the well location are obtained; The time-depth relationship, layer velocity, and average velocity at the well location are displayed and cross-analyzed to fit the logging time-averaged velocity function.

4. The method of claim 1, wherein, The step of using the initial seismic average velocity and the well logging average velocity to calibrate the seismic time-averaged velocity function to obtain a calibrated seismic average velocity model includes: The seismic time-averaged velocity function and the well logging time-averaged velocity function are fitted to obtain a time-averaged velocity model at the well location; wherein, the time-averaged velocity model at the well location is used to calculate the time-averaged velocity at the well location; The time-averaged velocity at the well location is used to calibrate the time-averaged velocity at the corresponding location or near the well to obtain a calibration factor. Perform Bayesian skryz interpolation on the scale factor to obtain the scale factor at all grid points; The average velocity is calibrated using the aforementioned scale factor to obtain the calibrated average velocity.

5. The method of claim 1, wherein, The formula for calculating the observation travel time is: where T(h) represents the reflection travel time of CMP gathers at each reflection point, V rms denotes the root mean square velocity, h is half of the offset distance, i.e. half offset; T0 represents the time of the time domain interpreted horizon corresponding to the depth domain interpreted horizon.

6. The method according to any one of claims 1 to 5, characterized in that, The method further includes: Obtain the system of inversion equations for constructing tomography; Solving the tomographic inversion equations yields the depth domain layer velocity update and the depth domain layer depth update. The optimized depth domain layer velocity is determined using the depth domain layer velocity update amount and the depth domain layer depth update amount.

7. An apparatus for automatically picking up reflection traveltime tomographic velocity modeling, characterized by, include: The acquisition module is used to acquire the seismic time-averaged velocity function constructed from the original seismic data and the well logging time-averaged velocity function constructed from the original well logging data; wherein, the seismic time-averaged velocity function is used to acquire the initial seismic average velocity of vertical change within the layer at or near the well location, and the well logging time-averaged velocity function is used to acquire the well logging average velocity of vertical change within the layer at the well location. The processing module is used to calibrate the seismic time-averaged velocity function using the initial seismic average velocity and the well logging average velocity to obtain a calibrated seismic average velocity model; and to use the calibrated seismic average velocity model to convert the time-domain interpretation horizon to the depth-domain interpretation horizon to obtain a calibrated depth-domain interpretation horizon. The processing module is also used to perform stacking velocity inversion on the calibrated depth domain layers one by one to obtain the optimized depth domain layer velocity; and to obtain the observation travel time based on the stacking velocity, wherein the observation travel time is the reflection travel time of the CMP gather at each reflection point; The processing module is also used to perform tomographic inversion using the calibrated depth domain interpretation layer, the optimized depth domain layer velocity, the observation travel time, and the CMP observation system to construct a depth domain layer velocity model.

8. A computer device comprising a memory, a processor, and a computer program stored on the memory, wherein the computer program comprises instructions that, when executed by the processor, cause the processor to perform the method of any one of claims 1-7. The processor executes the computer program to implement the steps of the method according to any one of claims 1 to 6.

9. A computer readable storage medium having stored thereon a computer program, characterized in that, When executed by a processor, the computer program implements the steps of the method according to any one of claims 1 to 6.

10. A computer program product comprising computer programs / instructions, characterized in that, When executed by a processor, the computer program implements the steps of the method according to any one of claims 1 to 6.