Speed model construction method, control device, and storage medium

Abstract

Embodiments of the present invention relate to the technical field of geological exploration, and provide a speed model construction method, a control device, and a storage medium. The method comprises: acquiring a three-dimensional seismic data volume of a preset area and well logging data of a pre-selected well; processing the three-dimensional seismic data volume and the well logging data to obtain a time-depth relationship and horizon interpretation data of the pre-selected well of a target layer; on the basis of the obtained time-depth relationship and horizon interpretation data of the pre-selected well of the target layer, carrying out frequency division inversion on the three-dimensional seismic data volume by means of a preset frequency division inversion strategy to obtain speed volumes of different frequency bands; and fusing the speed volumes of different frequency bands to obtain a speed model of the target layer. Well logging data and a three-dimensional seismic data volume are fully utilized, and frequency division inversion is carried out on the three-dimensional seismic data volume by means of a preset frequency division inversion strategy and a data fusion technique to obtain speed volumes of different frequency bands, so as to construct a speed model capable of finely characterizing the speed change features of the interior and boundary of a special geological body.

Description

Method, control device and storage medium for constructing velocity model

[0001] CROSS-REFERENCE TO RELATED APPLICATIONS

[0002] This application claims the benefit of Chinese patent application No. 202311033453.X filed on August 16, 2023, the contents of which are incorporated herein by reference. Technical Field

[0003] The present invention relates to the field of geological exploration technology, and in particular to a method for constructing a velocity model, a control device and a storage medium. Background Art

[0004] Velocity parameters are fundamental physical properties of underground strata and crucial parameters in seismic exploration. They influence every aspect of seismic exploration and play a crucial role in oil and gas exploration and development. Strata or geological bodies of different ages and lithologies typically exhibit varying interlayer velocities. Underground reservoir pores can also exhibit varying degrees of velocity differences due to differences in the fluids (gas, oil, water) they contain, or due to abnormalities or high-pressure closures in geological bodies. Therefore, an accurate and geologically consistent velocity model is not only crucial for prestack migration imaging and time-to-depth conversion, but also impacts numerous aspects of exploration and development, including structural accuracy, well depth determination, formation thickness, reservoir prediction, and pressure interface prediction.

[0005] As exploration and development continue, the geological conditions faced are becoming increasingly complex. We often encounter localized unique geological bodies, such as salt domes, ribbon reef formations, igneous intrusions, and carbonate anomalies. These unique geological bodies often differ significantly from the surrounding rock in lithology, velocity, and thickness, exhibiting strong heterogeneity. Current velocity modeling methods fail to reflect the local variations in velocity within these unique geological bodies, resulting in large errors in time-to-depth conversion, low structural accuracy, and the susceptibility to structural artifacts. These methods severely restrict the exploration and development of geological targets associated with these unique geological bodies. Therefore, complex geological conditions present new challenges for velocity modeling.

[0006] Currently, the main methods for establishing interval velocity models include: 1) using prestack time migration velocities or stacking velocities to calculate interval velocities using the Dix formula. This method, based on a homogeneous or horizontally layered medium model and estimated from seismic processing velocity spectra, can reveal velocity trends and relative relationships. However, this method suffers from poor accuracy for sedimentary stratigraphic sequences where spatial anisotropy is common. 2) Grid tomographic velocity modeling. This method uses a data-driven, computer-automated iterative solution approach to obtain more detailed velocity variations. However, this method is computationally intensive and, due to the lack of human intervention, suffers from high polynomial error. It is also difficult to control velocity errors in complex geological conditions. 3) Using seismic and acoustic logging data from existing wells, combined with interpretation of geological horizons, to establish interval velocity models. This method can generate high-precision interval velocity data that truly reflects formation lithology at locations with control wells. However, it is constrained by the planar distribution and drilling density of existing wells, resulting in poor velocity prediction accuracy, particularly in areas with few or no wells.

[0007] Summary of the Invention

[0008] The purpose of the embodiment of the present invention is to provide a method for constructing a velocity model. The velocity model constructed by the method for constructing a velocity model can accurately characterize the velocity variation characteristics inside and at the boundaries of a special geological body.

[0009] In order to achieve the above-mentioned purpose, an embodiment of the present invention provides a method for constructing a velocity model, which comprises: acquiring a three-dimensional seismic data volume and logging data of a preselected well in a preset area; processing the three-dimensional seismic data volume and the logging data to obtain the time-depth relationship and layer interpretation data of the preselected well of the target layer; based on the obtained time-depth relationship and layer interpretation data of the preselected well of the target layer, performing frequency-band inversion on the three-dimensional seismic data volume through a preset frequency-band inversion strategy to obtain velocity volumes of different frequency bands; and fusing the velocity volumes of the different frequency bands to obtain a velocity model of the target layer.

[0010] Optionally, before obtaining the three-dimensional seismic data body of the preset area and the logging data of the pre-selected wells, the method for constructing the velocity model also includes: determining the range of the target layer; obtaining the logging data of the wells in the preset area; and for the wells drilled through the target layer, screening the wells with acoustic wave curves, density curves and well diameter curves as the pre-selected wells according to the corresponding logging data.

[0011] Optionally, after acquiring the three-dimensional seismic data volume of the preset area and the logging data of the preselected well, the method for constructing the velocity model further includes: converting the acoustic wave curve of the preselected well into a velocity curve.

[0012] Optionally, before obtaining the time-depth relationship and stratigraphic interpretation data of the preselected wells of the target layer, the method for constructing a velocity model also includes: for the preselected wells of the target layer, analyzing the histogram of the single well through the logging data to determine the lithology and longitudinal variation characteristics of the velocity of the target layer; comparing and analyzing the velocity curves of connected wells through the logging data to determine the lateral variation characteristics of the velocity of the target layer; and determining the velocity variation characteristics and geological body type of the target layer through the longitudinal variation characteristics of the velocity of the target layer and the lateral variation characteristics of the velocity of the target layer.

[0013] Optionally, the logging data of the preselected well is processed to obtain the time-depth relationship of the preselected well in the target layer, including: obtaining the acoustic wave curve and density curve in the logging data of the preselected well; and making a synthetic seismic record through the acoustic wave curve and the density curve, establishing a corresponding relationship between the seismic reflection interface in the time domain and the geological stratification in the depth domain, and obtaining the time-depth relationship of the preselected well.

[0014] Optionally, the three-dimensional seismic data body and the well logging data are processed to obtain stratigraphic interpretation data, including: analyzing the well logging data of the preselected well to obtain the velocity change characteristics of the target layer; determining T geological layers of the target layer by analyzing the velocity change characteristics of the target layer and based on the velocity difference of the velocity interface; calibrating the T geological layers of the target layer by making synthetic seismic records, and determining the characteristics of the seismic reflection sync axes corresponding to each geological layer; and completing the tracking and interpretation of the seismic reflection sync axes corresponding to each geological layer based on the three-dimensional seismic data body to obtain stratigraphic interpretation data of the T geological layers of the target layer.

[0015] Optionally, the three-dimensional seismic data volume is subjected to frequency band inversion through a preset frequency band inversion strategy to obtain velocity bodies of different frequency bands, including: inverting the three-dimensional seismic data volume through the Kriging interpolation method to obtain a low-frequency band velocity body; inverting the three-dimensional seismic data volume through the spectral simulation inversion method to obtain a mid-frequency band velocity body; and inverting the three-dimensional seismic data volume through the waveform indication inversion method to obtain a high-frequency band velocity body.

[0016] Optionally, before performing frequency band inversion on the three-dimensional seismic data body through a preset frequency band inversion strategy, the method for constructing a velocity model also includes: determining T interpretation layers of the target layer from bottom to top; between each two adjacent interpretation layers, according to the contact relationship between the upper and lower layers, and according to a preset subdivision ratio C, interpolating between each two adjacent interpretation layers to obtain C interpolated sublayers; and constructing a layer grid model of the target layer from the T interpretation layers and the C interpolated sublayers between each two adjacent interpretation layers.

[0017] Optionally, the three-dimensional seismic data body is inverted by the Kriging interpolation method to obtain a low-frequency velocity body, including: obtaining a velocity curve in the logging data of the preselected well, and performing low-pass filtering on the velocity curve to retain the low-frequency velocity curve within a first cutoff frequency value; based on the obtained time-depth relationship of the preselected well, the filtered low-frequency velocity curve is converted from the depth domain to the time domain to obtain a low-frequency velocity curve in the time domain; and based on the stratigraphic grid model of the target layer, a lateral constraint condition of the low-frequency velocity curve is established, so as to perform logging interpolation on the three-dimensional seismic data body under the lateral constraint condition by the Kriging interpolation method to obtain the low-frequency velocity body.

[0018] Optionally, the stratum grid model based on the target layer establishes lateral constraints of the low-frequency velocity curve, and performs well logging interpolation under the lateral constraints on the three-dimensional seismic data body through the Kriging interpolation method to obtain the low-frequency band velocity body, including: for each layer in the stratum grid model, based on the observation value of each well of the preselected wells, obtaining the estimated value at each CDP of the layer, wherein the low-frequency velocity value read at the each well is used as the observation value of the each well, and the low-frequency velocity value of the layer is obtained through the estimated value at each CDP of the layer; and the low-frequency band velocity body is obtained through the low-frequency velocity value of each layer, and each layer includes the T interpretation layers and the C interpolated sublayers between each two adjacent interpretation layers.

[0019] Optionally, the three-dimensional seismic data volume is inverted by the spectral simulation inversion method to obtain a mid-frequency velocity volume, including: obtaining the velocity curve in the logging data of the preselected well and constructing a corresponding wave impedance curve; based on the obtained time-depth relationship of the preselected well in the target layer, the wave impedance curve of the preselected well is converted from the depth domain to the time domain to obtain a corresponding time domain wave impedance curve; and the time domain wave impedance curve is superimposed on the three-dimensional seismic data volume, and inverted by the spectral simulation inversion method to obtain the mid-frequency velocity volume.

[0020] Optionally, the time domain wave impedance curve is superimposed on the three-dimensional seismic data volume, and inverted by a spectral simulation inversion method to obtain the mid-frequency velocity volume, including: for each preselected well in the target layer, a spectral analysis is performed on the seismic trace beside the well and the time domain wave impedance curve of the preselected well, a frequency matching operator is calculated, the frequency matching operator is applied to the three-dimensional seismic data volume to obtain a spectral simulation inversion result, and according to the corresponding relationship between the reflection coefficient and the wave impedance, a wave impedance inversion result is obtained from the inversion result of the spectral simulation; the mid-frequency velocity volume is obtained by the ratio of the wave impedance inversion result of the target layer to the first density mean.

[0021] Optionally, the three-dimensional seismic data body is inverted by the waveform indication inversion method to obtain a high-frequency velocity body, including: for the wellbore seismic trace of each well in n pre-selected wells, under the stratigraphic grid model of the target layer, the seismic waveform to be predicted is compared with the waveform of the wellbore seismic trace of each well, and K waveform simulation samples of the wellbore seismic traces are selected according to two factors: waveform similarity and spatial distance, and an initial impedance model is constructed through the corresponding logging data; based on the characteristics of the initial impedance model and according to a preset high-pass frequency parameter value and a preset high-cutoff frequency parameter value, high-frequency band information is continuously filtered out, and comparisons are made between the waveform simulation samples of each wellbore seismic trace to retain a preset deterministic frequency band component; and based on the waveform simulation samples and spatial structure characteristics of the K wellbore seismic traces, an estimated value of an unknown point is determined to obtain a wave impedance inversion result of the waveform indication inversion, and the high-frequency velocity body is obtained by the ratio of the wave impedance inversion result of the target layer to the second density mean.

[0022] Optionally, fusing the velocity bodies of the different frequency bands to obtain the velocity model of the target layer includes: obtaining the velocity model of the target layer by weighted sum of the low-frequency band velocity body, the mid-frequency band velocity body and the high-frequency band velocity body.

[0023] An embodiment of the present invention also provides a control device for constructing a velocity model, the control device comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, the processor executing the computer program to implement the above-mentioned method for constructing a velocity model.

[0024] An embodiment of the present invention further provides a machine-readable storage medium having instructions stored thereon, wherein the instructions enable a machine to execute the above-mentioned method for constructing a velocity model.

[0025] Through the above technical solution, the embodiment of the present invention provides a method for constructing a velocity model for more accurate comprehensive seismic geological research related to special geological bodies, such as time-depth conversion, thickness calculation, reservoir prediction, and well location implementation. Based on the analysis of the geological characteristics of the special geological body, a detailed time-depth relationship and horizon interpretation data are obtained through fine seismic geological horizon calibration. Based on this, the well logging data and three-dimensional seismic data volume are fully utilized. Through a preset frequency band inversion strategy and data fusion technology, the three-dimensional seismic data volume is subjected to frequency band inversion to obtain velocity bodies of different frequency bands. The frequency band velocities of low-frequency, medium-frequency, and high-frequency bands are used to construct a velocity model that can accurately characterize the velocity variation characteristics within and at the boundaries of the special geological body. This provides reliable velocity data for comprehensive geological research such as time-depth conversion, fine structural implementation, and well location deployment, thereby facilitating the expansion and efficient advancement of exploration in related fields of special geological bodies.

[0026] Other features and advantages of the embodiments of the present invention will be described in detail in the subsequent detailed description. BRIEF DESCRIPTION OF THE DRAWINGS

[0027] The accompanying drawings are used to provide a further understanding of the embodiments of the present invention and constitute a part of the specification. Together with the following detailed description, they are used to explain the embodiments of the present invention, but do not constitute a limitation of the embodiments of the present invention. In the accompanying drawings:

[0028] FIG1 is a schematic flow chart of a method for constructing a velocity model according to an embodiment of the present invention;

[0029] FIG2 is a schematic diagram illustrating velocity variation characteristics of an example geological body.

[0030] FIG3 is a schematic cross-sectional view of example migrated and stacked seismic data;

[0031] FIG4 is a schematic diagram of a velocity profile of an example low-frequency velocity body of the target layer;

[0032] FIG5 is a schematic diagram of a velocity profile of a mid-frequency velocity body of an exemplary target layer;

[0033] FIG6 is a schematic diagram of a velocity profile of an example high-frequency velocity body of a target layer;

[0034] FIG7 is a schematic diagram of a velocity profile of an example velocity model; and

[0035] FIG8 is a schematic diagram of a velocity profile illustrating an example conventional velocity model. DETAILED DESCRIPTION

[0036] The following describes the specific implementation of the embodiment of the present invention in detail with reference to the accompanying drawings. It should be understood that the specific implementation described herein is only used to illustrate and explain the embodiment of the present invention and is not used to limit the embodiment of the present invention.

[0037] Existing velocity models are suitable for geological conditions with relatively stable lateral distribution of strata. However, for geological conditions with relatively complex geological structures, such as those with the development of special geological bodies such as salt domes, striped reef structures, and carbonate rock anomalies, existing velocity modeling cannot accurately characterize the velocity variation characteristics within and at the boundaries of special geological bodies, affecting the accuracy of structural characteristics and restricting exploration in areas related to special geological bodies. Therefore, for exploration areas with the development of special geological bodies, a velocity model that can accurately characterize the velocity variations within and at the boundaries of special geological bodies is needed.

[0038] FIG1 is a flow chart of a method for constructing a velocity model according to an embodiment of the present invention. Referring to FIG1 , the method for constructing a velocity model may include the following steps:

[0039] Step S110: Acquire a 3D seismic data volume of a preset area and logging data of preselected wells.

[0040] The three-dimensional seismic data volume is offset stacked seismic data after amplitude preservation processing.

[0041] Preferably, before step S110, the method for constructing a velocity model further includes: determining the range of the target layer; for wells drilled through the target layer, screening wells with acoustic wave curves, density curves and wellbore curves as the preselected wells based on corresponding logging data.

[0042] For example, wells that penetrate the target layer are selected in a preset area. Among these wells, n wells that have acoustic wave curves, density curves, and wellbore curves are preferably used as the preselected wells.

[0043] Preferably, after step S110, the method for constructing a velocity model further comprises: converting the acoustic wave curve of the preselected well into a velocity curve.

[0044] Preferably, the acoustic wave curves corresponding to the n preselected wells are converted into velocity curves by the following formula: Vel = 1000000 / DT (1)

[0045] Wherein, DT represents the sound wave curve, the unit is microsecond / meter, and Vel represents the velocity curve, the unit is meter / second.

[0046] As described above, a three-dimensional seismic data volume of a preset area and logging data of a preselected well are obtained.

[0047] Step S120: Processing the three-dimensional seismic data volume and the well logging data to obtain the time-depth relationship and horizon interpretation data of the pre-selected wells of the target layer.

[0048] Preferably, before obtaining the time-depth relationship and stratigraphic interpretation data of the preselected wells of the target layer, the method for constructing a velocity model also includes: for the preselected wells of the target layer, analyzing the histogram of the single well through the logging data to determine the longitudinal variation characteristics of the velocity of the target layer; comparing and analyzing the velocity curves of connected wells through the logging data to determine the lateral variation characteristics of the velocity of the target layer; and determining the velocity variation characteristics and geological body type of the target layer through the longitudinal variation characteristics of the velocity of the target layer and the lateral variation characteristics of the velocity of the target layer.

[0049] By way of example, the velocity variation characteristics of the target layer and the types of special geological bodies can be determined through comparative analysis of single wells and well connections within the target layer. Histogram analysis of single wells can determine the lithology and vertical velocity variation characteristics of the target layer. Comparative analysis of well connections of velocity curves of the target layer can analyze the lateral velocity variation characteristics of the target layer and determine the types of special geological bodies within the target layer (e.g., those associated with salt domes, zonal reef formations, and carbonate anomalies). These velocity variation characteristics can then be obtained, providing a foundation for establishing a velocity model.

[0050] Preferably, step S120 may include: 1. analyzing the logging data of the preselected well to obtain the time-depth relationship of the preselected well in the target layer, including: obtaining the acoustic wave curve and density curve in the logging data of the preselected well; and making a synthetic seismic record through the acoustic wave curve and the density curve, establishing a corresponding relationship between the seismic reflection interface in the time domain and the geological stratification in the depth domain, and obtaining the time-depth relationship of the preselected well.

[0051] Preferably, step S120 may include: 2. Analyzing the three-dimensional seismic data body and the well logging data to obtain stratigraphic interpretation data, including: analyzing the well logging data of the preselected well to obtain the velocity change characteristics of the target layer; determining the T geological layers of the target layer by analyzing the velocity change characteristics of the target layer and based on the velocity difference of the velocity interface; calibrating the T geological layers of the target layer by making synthetic seismic records, and determining the characteristics of the seismic reflection phase axis corresponding to each geological layer; and completing the tracking and interpretation of the seismic reflection phase axis corresponding to each geological layer based on the three-dimensional seismic data body to obtain the stratigraphic interpretation data of the T geological layers of the target layer.

[0052] 1. The fine time-depth relationship of the pre-selected wells in the target layer can be achieved through the method of synthetic seismic records. For example, synthetic seismic records are made using the acoustic wave curves and density curves of the pre-selected wells in the target layer, and the corresponding curves in the depth domain are converted to the time domain. That is, a one-to-one correspondence is established between the seismic reflection interface in the time domain and the geological stratification in the depth domain, so as to obtain a more detailed time-depth relationship for each pre-selected well.

[0053] 2. Based on the analysis results of the velocity change characteristics of the target layer determined above, and based on the velocity difference of the velocity interface, T geological layers in the target layer with obvious velocity interface changes are determined; through the synthetic seismic records produced in 1), the T geological layers of the target layer segment are calibrated, the characteristics of the seismic reflection phase axes corresponding to each geological layer are determined, and the seismic reflection phase axes corresponding to the T geological layers are tracked and interpreted based on the three-dimensional seismic data volume to obtain the stratigraphic interpretation data of the T geological layers of the target layer.

[0054] Step S130: Based on the obtained time-depth relationship of the preselected wells of the target layer and the horizon interpretation data, the 3D seismic data volume is subjected to frequency band inversion by a preset frequency band inversion strategy to obtain velocity volumes of different frequency bands.

[0055] Preferably, before step S130, the method for constructing a velocity model further includes: determining T interpretation horizons of the target layer from bottom to top; between each two adjacent interpretation horizons, according to the contact relationship between the upper and lower layers, and in accordance with a preset subdivision ratio C, interpolating between each two adjacent interpretation horizons to obtain C interpolated sublayers; and constructing a stratigraphic grid model of the target layer from the T interpretation horizons and the C interpolated sublayers between each two adjacent interpretation horizons.

[0056] For example, establishing a framework model for a target layer using T interpretation horizons involves selecting T interpretation horizons (e.g., T = 5) of the target layer segment from bottom to top. Between each pair of adjacent interpretation horizons, based on the contact relationship between the upper and lower layers and at a preset subdivision ratio C (e.g., C = 25), C interpolated sublayers are created between each pair of adjacent interpretation horizons, resulting in a total of C*(T-1) (sub)horizons. The T interpretation horizons and the interpolated sublayers between each pair of adjacent interpretation horizons constitute the layer framework model (Layer_Model) of the target layer, which can serve as a lateral constraint.

[0057] Preferably, step S130 may include: 1. Inverting the three-dimensional seismic data volume by the Kriging interpolation method to obtain a low-frequency velocity volume; 2. Inverting the three-dimensional seismic data volume by the spectral simulation inversion method to obtain a mid-frequency velocity volume; 3. Inverting the three-dimensional seismic data volume by the waveform indication inversion method to obtain a high-frequency velocity volume.

[0058] 1. Inverting the three-dimensional seismic data volume using the Kriging interpolation method to obtain a low-frequency velocity volume includes: obtaining a velocity curve from the well logging data of the preselected well, performing low-pass filtering on the velocity curve to retain the low-frequency velocity curve within a first cutoff frequency value; converting the filtered low-frequency velocity curve from the depth domain to the time domain based on the obtained time-depth relationship of the preselected well to obtain a low-frequency velocity curve in the time domain; and establishing lateral constraints for the low-frequency velocity curve based on a stratigraphic grid model of the target layer, and performing well logging interpolation on the three-dimensional seismic data volume under the lateral constraints using the Kriging interpolation method to obtain the low-frequency velocity volume.

[0059] To illustrate, the velocity curves of n preselected wells in a preset area (as described above, the velocity curves converted from the acoustic wave curves) are low-pass filtered. That is, the curve data with frequency components greater than a first cutoff frequency value FL (for example, FL = 15 Hz) are filtered out, and only the low-frequency velocity curves within FL are retained. Based on the time-depth relationship obtained above for each preselected well, the filtered low-frequency velocity curves are converted from the depth domain to the time domain to obtain low-frequency velocity curves in the time domain. Lateral constraints are established using the seismic interpretation horizon of the target layer, that is, the layer grid model Layer_Model of the target layer is used as the lateral constraint. Based on the filtered low-frequency velocity curves and the layer grid model Layer_Model, the ordinary Kriging interpolation method can be used to obtain the low-frequency (layer) velocity volume of the target layer (represented by Vel_Low).

[0060] Preferably, the stratum grid model based on the target layer establishes lateral constraints of the low-frequency velocity curve, so as to perform well logging interpolation under the lateral constraints on the three-dimensional seismic data body through the Kriging interpolation method to obtain the low-frequency band velocity body, including: for each layer in the stratum grid model, based on the observation value of each well of the preselected wells, obtaining the estimated value at each CDP of the layer, wherein the low-frequency velocity value read at each well is used as the observation value of each well, and the low-frequency velocity value of the layer is obtained through the observation value at each CDP of the layer and the estimated value at each CDP of each layer; and the low-frequency band velocity body is obtained through the low-frequency velocity value of each layer, and each layer includes the T interpretation layers and the C interpolated sublayers between each two adjacent interpretation layers.

[0061] Continuing with the above example, for each layer LayerM in the layer grid model (including the T interpretation layers and the C interpolated sublayers between each two adjacent interpretation layers), the low-frequency velocity value at each CDP along the layer LayerM in the preset area is regarded as a regional variable Z(x). For n preselected wells (for example, n=25), a low-frequency velocity value can be read at each preselected well along the layer, and the low-frequency velocity values ​​read from the n preselected wells can be regarded as n observation values, that is, Z(x1), Z(x2), Z(x3), ..., Z(x n ), the estimated value of the low-frequency velocity at a certain CDP along the layer LayerM in the preset area is recorded as Z(x0), which can be obtained by the weighted sum of the observation values ​​of n known sampling points, that is, it can be expressed by the following formula:

[0062] Among them, α i (i = 1, 2, …, n) are the weight coefficients to be determined. According to the principles of kriging interpolation, under the premise of ensuring that the estimator is unbiased and the estimated variance is minimized, the n weight coefficient values ​​can be obtained by solving the following system of equations.

[0063] Among them, γ(x i ,x j ) is the sampling point x i and x j The linear variation function value between them. From equation group (3), we can get the weighting coefficient α i By substituting the values ​​of (i=1,2,…,n) into equation (2), we can obtain the estimated value of the low-frequency velocity at a certain CDP, which is recorded as Z(x0).

[0064] In this way, the estimated value at each CDP along each layer LayerM in the preset area can be calculated, and the low-frequency velocity value along each layer LayerM in the layer grid model can be obtained in turn, thereby obtaining the low-frequency band velocity body Vel_Low of the target layer.

[0065] 2. Inverting the three-dimensional seismic data volume using a spectral simulation inversion method to obtain a mid-frequency velocity volume includes: obtaining a velocity curve from the logging data of the preselected well and constructing a corresponding wave impedance curve; based on the obtained time-depth relationship of the preselected well in the target layer, converting the wave impedance curve of the preselected well from the depth domain to the time domain to obtain a corresponding time-domain wave impedance curve; and superimposing the time-domain wave impedance curve on the three-dimensional seismic data volume and inverting it using a spectral simulation inversion method to obtain the mid-frequency velocity volume.

[0066] The inverted interval velocity volume within the mid-frequency band of the target layer can be obtained through spectral simulation inversion. Spectral simulation inversion is a frequency-domain well logging-constrained impedance inversion technique. Spectral analysis is performed on the well logging data and the 3D seismic data volume. A matching operator is designed in the frequency domain to match the seismic spectrum with the well impedance spectrum. The frequency matching operator is then applied to the 3D seismic data volume and converted into a time-domain operator for inversion, resulting in the inverted data volume.

[0067] As an example, based on the velocity curves of n preselected wells within a preset area, an impedance curve for each preselected well is constructed and used as input data for the spectral simulation inversion method. Preferably, to facilitate conversion of the inverted impedance data volume into a velocity data volume, when constructing the impedance curve, the density is preferably taken as a constant, i.e., the average density logging value of the target layer, Den_average. The velocity curve values ​​of the preselected wells are multiplied by the constant density value, Den_average, to obtain the impedance curve for the preselected well. Using the time-depth relationship already obtained for each preselected well, the impedance curve for the preselected well is converted from the depth domain to the time domain, obtaining a time-domain impedance curve, IMP_time. By inverting the time-domain impedance curves of the n preselected wells and the 3D seismic data volume (i.e., amplitude-preserving migration stack data), a mid-frequency velocity volume (represented by Vel_Medium) of the target layer can be obtained through spectral simulation.

[0068] Preferably, the time domain wave impedance curve is superimposed on the three-dimensional seismic data volume, and inverted by a spectral simulation inversion method to obtain the mid-frequency velocity volume, including: for each preselected well of the target layer, a spectrum analysis is performed on the seismic trace beside the well and the time domain wave impedance curve of the preselected well, a frequency matching operator is calculated, the frequency matching operator is applied to the three-dimensional seismic data volume to obtain a spectrum simulation inversion result, and according to the corresponding relationship between the reflection coefficient and the wave impedance, a wave impedance inversion result is obtained from the inversion result of the spectrum simulation; and the mid-frequency velocity volume is obtained by the ratio of the wave impedance inversion result of the target layer to the first density mean.

[0069] Following the above example, n wellside seismic traces are extracted correspondingly at n pre-selected wells, and each seismic trace is recorded as s i (t), (i=1,2,…,n), the seismic wavelet is denoted as b(t), and the reflection coefficient is denoted as r i (t)(i=1,2,…,n), the convolution model of the seismic trace can be expressed as follows: i (t) = r i (t)*b(t) (4)

[0070] Transforming it into the frequency domain, we get the following formula:

[0071] Among them, S i (ω) is the amplitude spectrum of the seismic trace, Ф i (ω) is the phase spectrum of the seismic trace, R i (ω) is the amplitude spectrum of the reflection coefficient, Ψ i (ω) is the phase spectrum of the reflection coefficient, B i (ω) is the amplitude spectrum of the wavelet, φ i (ω) is the phase spectrum of the wavelet, then according to formula (5), we can get the following formula: S i (ω)=R i (ω)B i (ω) (6)

[0072] Taking the logarithm of formula (6), we can get the following formula: lnS i (ω)=lnR i (ω)+lnB i (ω) (7)

[0073] Assume that the wave impedance at each preselected well is z i (t), (i=1,2,…,n), then the corresponding wave impedance z i (t) is transformed into the frequency domain and the following formula can be obtained:

[0074] Among them, Z i (ω) is the wave impedance amplitude spectrum at the preselected well, θ i (ω) is the wave impedance phase spectrum of the preselected well. According to the general relationship between wave impedance and reflection coefficient, the reflection coefficient amplitude spectrum R can be obtained by using the wave impedance amplitude spectrum i (ω), the amplitude spectrum of the reflection coefficient R i (ω) and the amplitude spectrum S of the seismic trace i Substituting (ω) into formula (7), we can get the frequency matching operator B i (ω).

[0075] Applying the matching operator obtained in the frequency domain to the 3D seismic data volume, we can get the inversion result of spectral simulation. i (ω)=1 / B i (ω), according to formula (6), we can get the following formula:

[0076] Perform inverse Fourier transform of equation (9) to the time domain to obtain the calculation formula of the reflection coefficient, which is expressed by the following equation: i (t) = a i (t)*s i (t) (10)

[0077] Among them, a i(t) is A i The time domain function of (ω) can be obtained by frequency matching operator B i (ω) is obtained.

[0078] For a seismic trace at a CDP in the preset area, through equation (10), the time domain operator a i (t) and seismic trace s i The convolution of (t) yields the inversion result of the reflection coefficient. Using the universal relationship between reflection coefficient and impedance, the reflection coefficient is then converted into an impedance value, yielding the impedance inversion result for the seismic trace. The impedance obtained here is the relative impedance. To more intuitively analyze the inversion effect of the spectrum simulation impedance and facilitate inversion quality control using wells, the mid-frequency impedance inversion result is added to the low-frequency background to obtain the final absolute impedance for the spectrum simulation.

[0079] The above calculation is performed on each seismic trace in the preset area to obtain the spectral simulated wave impedance inversion volume of the target layer, and then the wave impedance inversion volume is divided by the density average Den_average to obtain the mid-frequency velocity volume Vel_Medium of the target layer.

[0080] 3. Inverting the three-dimensional seismic data volume using a waveform indication inversion method to obtain a high-frequency velocity volume includes: for the wellbore seismic trace of each of the n pre-selected wells, comparing the seismic waveform to be predicted with the waveform of the wellbore seismic trace of each well under the stratigraphic grid model of the target layer, selecting K (K≤n) waveform simulation samples of the wellbore seismic traces by analyzing two factors: waveform similarity and spatial distance, and constructing an initial impedance model using the corresponding well logging data; based on the characteristics of the initial impedance model and according to preset high-pass frequency parameter values ​​and preset high-cutoff frequency parameter values, continuously filtering out high-frequency band information, comparing the waveform simulation samples of each wellbore seismic trace, and retaining preset deterministic frequency band components; and determining an estimated value of an unknown point based on the waveform simulation samples of the K wellbore seismic traces and spatial structure characteristics to obtain a wave impedance inversion result of the waveform indication inversion, and obtaining the high-frequency velocity volume by the ratio of the wave impedance inversion result of the target layer to the second density mean.

[0081] Waveform indication inversion is a new high-precision inversion method developed on the basis of traditional geostatistical inversion. Its theoretical basis is that the amplitude, phase, frequency and other information contained in the seismic waveform represent the characteristics of the reservoir. Although the effective seismic frequency band is narrow and the high-frequency components cannot be directly obtained, the lateral changes in the seismic waveform reflect the changing characteristics of the reservoir. Therefore, the lateral changes in the seismic waveform can be fully utilized to reflect the phase change characteristics of the reservoir space, and then the high-frequency structural characteristics of the reservoir lithologic combination can be analyzed to achieve high-precision prediction of well-seismic combination. The basic idea of ​​waveform indication inversion is: based on the two basic factors of reference spatial distribution distance and seismic waveform similarity, all wells are sorted according to the correlation, and the well with the highest correlation with the predicted point is selected as the initial model to perform unbiased optimal estimation of the high-frequency components, and ensure that the final inverted seismic waveform is consistent with the original seismic waveform. The process of obtaining the high-frequency band velocity body of the target layer by the waveform indication inversion method in the embodiment of the present invention is as follows:

[0082] The number of effective samples is K (i.e., waveform simulation samples of K near-well seismic traces selected from the seismic waveforms of n near-well seismic traces mentioned above, for example, K = 5). The waveform indication inversion requires the following three processes in sequence: 1) sample optimization, 2) sample curve (i.e., waveform simulation sample) analysis, and 3) spatial valuation.

[0083] 1) For sample optimization, the seismic waveform to be predicted is compared with the well bypass waveforms of n pre-selected wells under the layer framework model Layer_Model. Based on the dual variables of waveform similarity and spatial distance, the K well bypasses with the most similar waveforms are selected as simulation samples, and their logging data are used to construct the initial impedance model.

[0084] 2) For sample curve analysis, the impedance characteristics of the sample wells are statistically analyzed. Based on the two parameter values ​​of high-pass frequency R1 and high-cutoff frequency R2, the high-frequency band information is continuously filtered out using the frequency division method and compared between samples to retain the deterministic frequency band components.

[0085] 3) For spatial valuation, a linear, unbiased, and optimal estimation of unknown points is performed based on the original data and spatial structure characteristics of the sample wells.

[0086] Let M(x0) be the estimated value of the unknown point, M(x i )(i=1,2,…,K) is the value of the waveform’s preferred sample point, λ i is the weight of the i-th sample point to the unknown point, K is the number of valid samples, and the estimated value of the unknown point is calculated by the following formula:

[0087] Equation (11) shows the wave impedance volume obtained by waveform inversion. Dividing it by the density mean Den_average yields the interval velocity volume. However, based on the principle of waveform inversion, low- and medium-frequency information from earthquakes is applied during the inversion process. Therefore, the inversion result contains certain low- and medium-frequency components. To remove these components and retain only the high-frequency components from the waveform inversion, the waveform inversion can be high-pass filtered to obtain a data volume reflecting only the high-frequency components, i.e., the high-frequency velocity volume obtained from the waveform inversion (represented by Vel_High).

[0088] Preferably, based on the three-dimensional seismic data body, the layer grid model and the time domain impedance curve, the three key inversion parameters of smoothing radius, effective sample number and cutoff frequency are tested to optimize the inversion result. The smoothing radius mainly affects the lateral connectivity of the inversion body. The larger the smoothing radius, the better the connectivity of the inversion body. The effective sample number mainly characterizes the degree of influence of the spatial variation of the seismic waveform on the reservoir. When the effective sample number increases to a certain extent, the characterization degree of the seismic waveform on the reservoir no longer changes, and the sample number at this time is the effective sample number. The cutoff frequency controls the frequency band of the inversion, including high-pass and high-cutoff frequencies in the low-frequency band, and low-cut, low-pass, high-pass and high-cutoff frequencies in the high-frequency band. Among them, the value of the cutoff frequency can be selected according to the amplitude spectrum of the three-dimensional seismic data body, and its value principle is to effectively compensate for the low-frequency and high-frequency parts of the seismic data. Through repeated parameter testing, the values ​​of the three parameters are constantly fine-tuned to find the parameters that can minimize the waveform indication inversion error of the target layer and optimize the inversion result. For example, the key parameters are determined as follows: a smoothing coefficient of 2, a valid sample number of 5, high-pass and high-cut frequencies of 5 Hz and 10 Hz in the low-frequency band, and low-cut, low-pass, high-pass, and high-cut frequencies of 40 Hz, 50 Hz, 90 Hz, and 100 Hz in the high-frequency band. Based on the inversion parameters determined above, waveform-indicating inversion is performed to obtain a wave impedance data volume for waveform-indicating inversion, which is then divided by the density average Den_average (for example, Den_average = 2.5) to obtain a layer velocity volume for waveform-indicating inversion. However, since waveform-indicating inversion includes certain low-frequency and medium-frequency information during the calculation process, in order to obtain the inverted high-frequency velocity value, the inversion result needs to be high-pass filtered, that is, the frequency components below FH (FH can be determined based on the amplitude spectrum of the earthquake, for example, FH = 60 Hz) are filtered out, and only the high-frequency components are retained, thereby obtaining the high-frequency velocity volume Vel_High.

[0089] Step S140: fusing the velocity bodies of the different frequency bands to obtain a velocity model of the target layer.

[0090] Preferably, step S140 may include: obtaining the velocity model of the target layer by performing a weighted sum of the low-frequency velocity volume, the mid-frequency velocity volume, and the high-frequency velocity volume.

[0091] As described in step S130, the velocity model of the target layer can be obtained by fusing the inverted velocity volumes of the low-frequency, medium-frequency, and high-frequency segments. That is, the (layer) velocity model Vel_Model is calculated by using the low-frequency segment velocity volume Vel_Low, the medium-frequency segment velocity volume Vel_Medium, and the high-frequency segment velocity volume Vel_High of the target layer obtained through segmented inversion.

[0092] Preferably, during the mid-frequency inversion of the target layer (using spectral simulation inversion to invert the 3D seismic data volume to obtain a mid-frequency velocity volume), a certain amount of low-frequency information is added to more intuitively analyze the vertical and horizontal variation characteristics of the reservoir and perform effective inversion quality control. Therefore, the low-frequency information in the mid-frequency velocity volume Vel_Medium can be removed. That is, the mid-frequency velocity volume Vel_Medium is high-pass filtered (the parameters of the high-pass filter can be determined based on the spectrum of the 3D seismic data volume, for example, the high-pass frequency is 15 Hz). Through high-pass filtering, a (layer) velocity volume Vel_Mediumfilter that retains only the mid-frequency information is obtained.

[0093] Next, the low-frequency velocity body Vel_Low, the medium-frequency velocity body Vel_Medium, and the high-frequency velocity body Vel_High of the target layer are weighted according to preset weights to obtain a fused velocity model. The preset weights can be determined by analyzing the fit between the inverted interval velocity data and the wellbore. For example, they can be values ​​between 1 and 2, with higher weights assigned to higher fits and lower weights assigned to lower fits. For example, by analyzing the fit between the inversion results and the wellbore, the weights W_Low, W_Medium, and W_High for the low-frequency velocity body, medium-frequency velocity body, and high-frequency velocity body are determined. For any sampling point i in the target layer, the layer velocity values ​​of the corresponding low-frequency velocity body Vel_Low(i), medium-frequency velocity body Vel_Mediumfilter(i) and high-frequency velocity body Vel_High(i) are read in sequence. The layer velocity value Vel_Model(i) after data fusion of the point is obtained by the following formula: Vel_Model(i)=Vel_Low(i)×W_Low+Vel_Mediumfilter(i)×W_Medium+Vel_High(i)×W_High (12)

[0094] Based on this, the data fusion of each sampling point in the target layer can be realized to obtain the target layer velocity model Vel_Model.

[0095] Preferably, the embodiment of the present invention can also analyze the reliability of the velocity model obtained in steps S110-S140 through 1) velocity model comparison and / or 2) depth error comparison at the well point after time-depth conversion.

[0096] 1) The (layer) velocity body Vel_tradition obtained by conventional model-based inversion is compared with the (layer) velocity model Vel_Model of the embodiment of the present invention: the (layer) velocity body Vel_tradition cannot accurately characterize the vertical and horizontal variation characteristics of the formation velocity in the preset area, and the (layer) velocity body Vel_tradition cannot identify the velocity changes inside and at the boundaries of special geological bodies; while the (layer) velocity model Vel_Model of the embodiment of the present invention can not only reflect the velocity characteristics between layers and between different lithofacies in the vertical direction, and accurately reflect the geological laws of the area, but also can well characterize the development zone boundaries of special geological bodies, and has a high degree of consistency with the velocity on the drilling well, thereby verifying the reliability of the velocity model Vel_Model.

[0097] 2) Time-depth conversion of the target layer was performed using the (layer) velocity body Vel_tradition and the (layer) velocity model Vel_Model of the embodiment of the present invention respectively; the depth error at each well point after the conversion between the (layer) velocity body Vel_tradition and the (layer) velocity model Vel_Model of the embodiment of the present invention was calculated respectively; a comparative analysis of the depth errors at the well points obtained after the conversion between the two different velocity models showed that the depth error value at the well point obtained after the time-depth conversion of the (layer) velocity model Vel_Model of the embodiment of the present invention was significantly lower than the depth error value after the conversion using the (layer) velocity body Vel_tradition, verifying the reliability of the (layer) velocity model Vel_Model of the embodiment of the present invention.

[0098] Accordingly, embodiments of the present invention provide a method for constructing a velocity model for more accurate comprehensive seismic geological research related to special geological bodies, including time-depth conversion, thickness calculation, reservoir prediction, and well location implementation. Based on an analysis of the geological characteristics of the special geological body, a detailed drilling time-depth relationship and horizon interpretation data are obtained through precise seismic geological horizon calibration. Based on this, well logging data and three-dimensional seismic data are fully utilized. A frequency-band inversion strategy and data fusion technology are used to perform frequency-band inversion on the three-dimensional seismic data to obtain velocity bodies in different frequency bands. Frequency-band velocity bodies in low-frequency, medium-frequency, and high-frequency bands are used to construct a velocity model that can accurately characterize the velocity variation characteristics within and at the boundaries of the special geological body. This provides reliable velocity data for comprehensive geological research, such as time-depth conversion, structural refinement, and well location deployment, thereby facilitating the expansion and efficient advancement of exploration in related fields of special geological bodies.

[0099] Taking a layer segment with a special geological body developed in a certain exploration area of ​​a certain basin as an example, an implementation process of an embodiment of the present invention is described in detail.

[0100] As shown in Figure 3, amplitude-preserving offset stacked seismic data is acquired. Well logging data for the target layer in the area is also obtained. All wells penetrating the target layer are selected, and based on the corresponding well logging data, n wells with simultaneous acoustic, density, and caliper curves are selected as preselected wells. For example, 28 wells penetrating the target layer in the current area are selected. The logging data of these 28 wells is analyzed, and 25 wells with simultaneous acoustic, density, and caliper curves are selected, resulting in n = 25 preselected wells.

[0101] Then, the acoustic wave curve of each preselected well is converted into a velocity curve by using formula (1). That is, the acoustic wave curve, density curve, velocity curve and caliper curve of the 25 preselected wells are obtained, thereby completing the acquisition of the logging data of the preselected wells.

[0102] By comparing and analyzing single wells and well connections in the target layer, the velocity variation characteristics and the type of special geological bodies in the target layer are obtained. By analyzing the single-well histograms of preselected wells, the lithologic characteristics and vertical velocity variation characteristics of the target layer are obtained. In this example, the target layer in this area is mainly carbonate sedimentary strata, with lithologies mainly consisting of interbedded mudstone, limestone, and marl-limestone. The mudstone velocity is low, the limestone velocity is high, and the vertical velocity of the target layer shows interbedded high-speed and low-speed zones. By comparing the velocity curves of the target layer in a well connection, the lateral velocity variation characteristics of the target layer are analyzed, the type of special geological bodies within the target layer and their velocity characteristics are determined, and the geological foundation for the establishment of the velocity model is provided.

[0103] Among them, through well-to-well comparison of some preselected wells, as shown in Figure 2, we can obtain the thickness and velocity variation characteristics of the three target layers (P1, KT1, and MKT) from top to bottom. The KT1 and MKT segments are relatively stable in the lateral distribution of the formation, with no significant changes in velocity and thickness. However, the thickness and velocity of the top target layer, P1, vary significantly. The velocity of P1 segment varies significantly both vertically and horizontally, with low velocity in the lower part and high velocity in the middle and upper parts. The velocity in the lower part of P1 segment is relatively stable in the lateral distribution, with an average velocity of approximately 3200 m / s. However, the velocity in the middle and upper part of P1 segment varies greatly in the lateral distribution, ranging from 3300 m / s to 5600 m / s. The low velocity zone is around 3300 m / s, while the high velocity zone can reach 5600 m / s. The formation thickness of P1 segment in the high velocity zone increases significantly. After analysis, it can be determined that the velocity in the upper and middle parts of the P1 section of the target layer varies greatly laterally and is highly heterogeneous. There are obvious special geological bodies in some areas, and the type of special geological bodies is limestone high-speed anomaly bodies. The velocity changes in the lower P1 section, KT1 and MKT sections of the target layer are relatively stable.

[0104] By calibrating synthetic seismic records, the detailed time-depth relationship of the preselected wells is determined. In this example, synthetic records are generated for 25 preselected wells, and detailed well-seismic calibration is performed on the target layer to obtain the detailed time-depth relationship of the 25 preselected wells.

[0105] Next, based on the analysis results of the velocity variation characteristics of the target layer determined above, T seismic interpretation horizons of the target layer are determined. In this example, based on the analysis results of the velocity variation characteristics of the target layer determined above, it is determined that there are mainly five geological layers with obvious velocity variations in the target layer, namely T=5, which are the top of P1, the top of mudstone, the top of KT1, the top of MKT and the top of KT2. Based on the synthetic seismic records of 25 pre-selected wells, the reflection characteristics of each geological layer on the seismic profile are calibrated, and the seismic reflection horizons corresponding to the five geological layers are tracked and interpreted within the three-dimensional seismic data volume. In order to conduct seismic geological analysis, the naming of each seismic interpretation horizon is consistent with its corresponding geological layer. As shown in Figure 3, the horizon interpretation data of the five seismic reflection horizons corresponding to the target layer and the geological layer are determined through tracking and interpretation.

[0106] Based on well logging data, time-depth relationships, and target horizon interpretation data, a low-frequency velocity volume is generated using Kriging interpolation. Velocity curves from preselected wells in the area are low-pass filtered. Based on spectral analysis of the 3D seismic data volume, the first cutoff frequency FL for the low-frequency band is determined. This filter removes frequency components above the first cutoff frequency FL, retaining only velocity values ​​within the cutoff frequency. This results in a filtered low-frequency velocity curve. This filtered low-frequency velocity curve is then converted from the depth domain to the time domain based on the time-depth relationship of each preselected well.

[0107] Next, a layer grid model (Layer_Model) of the target layer is established based on the layer interpretation data as a lateral constraint for the velocity model. In this example, T interpretation layers of the target layer are selected from bottom to top. Between each two adjacent interpretation layers, C interpolated sublayers are interpolated according to a certain subdivision ratio C (for example, C = 25) based on the contact relationship between the upper and lower layers. This results in C*(T-1) (sub)layers. The layer grid model (Layer_Model) of the target layer is composed of T interpretation layers and the interpolated sublayers between each two adjacent interpretation layers. This layer grid model can be used as a lateral constraint for the low-frequency velocity body.

[0108] Based on this, the low-frequency velocity curve in the filtered time domain and the layer grid model Layer_Model are used to obtain the low-frequency velocity body of the target layer segment through the ordinary Kriging interpolation method. In this example, for each layer LayerM in the layer grid model (including the T interpretation layers and the C interpolation sublayers between each two adjacent interpretation layers), the low-frequency velocity value at each CDP along the layer in the region is regarded as the regional variable Z(x). For n pre-selected wells in the region, a low-frequency velocity value can be read at each well along the layer. The low-frequency velocity values ​​read from the n pre-selected wells can be regarded as n observations. Through the above formula (2) and equation group (3), the estimated value of the low-frequency velocity at a certain CDP can be obtained and recorded as Z(x0), that is, the estimated value of a certain layer LayerM at a certain CDP along the layer grid model Layer_Model. According to the above calculation process, the estimated value at each CDP along the layer LayerM in the area is calculated to obtain the low-frequency velocity value along the layer LayerM; for each layer in the layer grid model Layer_Model, according to the calculation process of the low-frequency velocity value along the layer LayerM, the low-frequency velocity value along each layer in the layer grid model Layer_Model can be obtained in turn, thereby obtaining the low-frequency segment velocity body Vel_Low of the target layer.

[0109] In this example, velocity logs (denoted as Vp) from 25 preselected wells are low-pass filtered. Spectral analysis of the 3D seismic data volume for the target layer reveals a frequency band of 15Hz-45Hz, with a dominant frequency of 30Hz. This means that the 3D seismic data volume lacks frequency components below 15Hz. Therefore, low-frequency components above the wells can be used to supplement these frequencies. A first cutoff frequency value, FL = 15Hz, is determined, extracting frequency components below 15Hz from the wells. Low-pass filtering of the velocity curves yields low-frequency velocity curves within 15Hz. Based on the time-depth relationships of the 25 preselected wells, these low-frequency velocity curves in the depth domain are converted to the time domain. Next, based on the T seismic interpretation horizons (T = 5) in the target layer of the region, a subdivision ratio, C = 25, is selected. 25 sublayers are interpolated between each two adjacent interpretation horizons, working from bottom to top. This results in a layer grid model, Layer_Model, consisting of five interpretation horizons and C interpolated sublayers between each two adjacent interpretation horizons. Based on the time domain low-frequency velocity curve and the layer grid model Layer_Model, for a certain layer in the layer grid model Layer_Model, the low-frequency velocity values ​​of 25 pre-selected wells at the layer are read in sequence. Then, the Kriging interpolation method is used to calculate the low-frequency velocity value at each CDP along the layer through the above formula (2) and equation group (3) to obtain the velocity value of the layer. According to the same calculation process, the velocity value of each layer in the layer grid model Layer_Model is calculated, and thus the low-frequency interval velocity body Vel_Low of the target layer is obtained. As shown in Figure 4, the example profile reflects the overall change trend of the target layer velocity and reflects the background velocity characteristics of the high-speed anomaly body of the target layer limestone.

[0110] The inversion layer velocity body within the mid-frequency band of the target layer can be obtained by the spectral simulation inversion method. According to the velocity curve of the preselected well, the impedance curve of the preselected well is constructed and used as the input data of the spectral simulation inversion. In order to facilitate the conversion of the impedance data body obtained by inversion into the velocity data body, when constructing the impedance curve, the density value can be taken as a constant, that is, the average value Den_average of the density logging value of the target layer segment (in this example, Den_average can be set to 2.5 according to the statistical results of the density logging value of the target layer segment). The velocity curve value of the preselected well is multiplied by the constant density value Den_average to obtain the impedance curve of the preselected well, and according to the obtained time-depth relationship of each preselected well, the impedance curve of the preselected well is converted from the depth domain to the time domain to obtain the time domain impedance curve IMP_time.

[0111] Next, the time domain impedance curves of n pre-selected wells and the 3D seismic data volume (i.e., amplitude-preserving migration stacking data) are used for spectral simulation inversion to obtain the mid-frequency velocity volume Vel_Medium of the target layer. By spectral analysis of the well logging data and the 3D seismic data volume, a frequency matching operator is obtained that achieves the best match between the seismic spectrum and the well logging wave impedance spectrum. At the n pre-selected wells, n wellside seismic traces are extracted, each of which is denoted as s i (t), (i=1,2,…,n), the seismic wavelet is denoted as b(t), and the reflection coefficient is denoted as r i (t)(i=1,2,…,n), then according to the above equations (4) and (5), the amplitude spectrum S of the earthquake can be obtained i (ω), the amplitude spectrum of the reflection coefficient R i (ω) and the amplitude spectrum of the wavelet B i That is, the amplitude spectrum S of the earthquake revealed in the above formula (6) i (ω) is equal to the amplitude spectrum of the reflection coefficient R i (ω) and the amplitude spectrum of the wavelet B i Taking the logarithm of equation (6), we can get the amplitude spectrum S of the earthquake which is more convenient to calculate. i (ω), the amplitude spectrum of the reflection coefficient R i (ω) and the amplitude spectrum of the wavelet B i For equation (7), since the amplitude spectrum of the earthquake is known, we only need to find the amplitude spectrum R of the reflection coefficient. i (ω), we can get the amplitude spectrum B of the wavelet i (ω), where B i (ω) is the frequency matching operator for spectrum simulation inversion. Since the wave impedance of the preselected well is known, the reflection coefficient amplitude spectrum R can be obtained based on the corresponding relationship between the general wave impedance and the reflection coefficient and the amplitude spectrum of the wave impedance. i (ω), and the amplitude spectrum of the reflection coefficient R i (ω) and the amplitude spectrum S of the seismic trace i Substituting (ω) into formula (7), we can get the frequency matching operator B i (ω).

[0112] Apply the frequency matching operator to the 3D seismic data volume to obtain the spectrum simulation inversion data volume. i (ω), the amplitude spectrum of the reflection coefficient R i (ω) and the amplitude spectrum of the wavelet B i (ω), the reciprocal of the frequency matching operator is taken as A i (ω), that is, A i (ω)=1 / B i(ω), we can get Equation (9), and perform inverse Fourier transform of Equation (9) to the time domain to obtain the calculation formula of the reflection coefficient, namely Equation (10). According to Equation (10), for a seismic trace at a CDP in the area, through the time domain operator a i (t) and seismic trace s i The convolution of (t) yields the inversion result of the reflection coefficient for the seismic trace. Using the universal correspondence between reflection coefficient and impedance, the reflection coefficient is then converted to an impedance value, yielding the impedance inversion result for the seismic trace. The resulting impedance is relative. To more intuitively analyze the inversion effect of spectral simulation impedance and facilitate inversion quality control using wells, the mid-frequency relative impedance inversion result can be added to a low-frequency background. This means adding the relative impedance inversion value obtained from spectral simulation to the low-frequency background value from the well logging. This means assigning a low-frequency impedance value to each layer based on the well logging analysis results to yield the final absolute impedance value from spectral simulation.

[0113] The above calculation is performed on each seismic trace in the area to obtain the spectral simulated wave impedance inversion volume of the target layer in the area. Then, the wave impedance inversion volume is divided by the density average Den_average to obtain the mid-frequency velocity volume Vel_Medium of the target layer.

[0114] In this example, based on the density curves of the 25 pre-selected wells, the density average of the target layer Den_average = 2.5 is calculated. The velocity curve values ​​of the 25 pre-selected wells are multiplied by the constant density value Den_average to obtain the impedance curve of the pre-selected wells. Then, based on the obtained time-depth relationship of each pre-selected well, the impedance curve of the pre-selected wells is converted from the depth domain to the time domain. Then, based on the universal correspondence between reflection coefficient and wave impedance, the logging wave impedance amplitude spectrum of the 25 pre-selected wells is converted into a reflection coefficient amplitude spectrum, and the well side channel amplitude spectrum and logging wave impedance amplitude spectrum R of the 25 selected wells are extracted. i (ω), calculate the amplitude spectrum S of the well bypass channel of 25 pre-selected wells i (ω). The frequency matching operator B for spectrum simulation inversion is obtained according to formula (7): i (ω). Then, take the inverse of the frequency matching operator and transform it to the time domain through the inverse Fourier transform to obtain the operator a in the time domain i (t), using formula (10), through the time domain operator a i (t) and each seismic trace s in the 3D seismic data volume i(t) is convolved to obtain the inversion result of the reflection coefficient. Using the universal correspondence between reflection coefficient and impedance, the reflection coefficient is then converted to an impedance value, resulting in the relative impedance inversion result. Based on the statistical results of the well logging data, a low-frequency background value is set for each layer of the target layer. For example, the low-frequency impedance background values ​​from top to bottom are 10050, 8250, 11250, and 6875, respectively. This yields the final absolute impedance volume for spectral simulation. The impedance inversion volume is then divided by the density mean of 2.5 to obtain the mid-frequency velocity volume, Vel_Medium, within the mid-frequency band of the target layer. As shown in Figure 5, the example profile reflects the velocity variation characteristics of the target layer in the same frequency band as the seismic data.

[0115] The high-frequency velocity body of the target layer is obtained using the waveform indication inversion method. According to the basic principle of the waveform indication inversion method, the key parameter test of the waveform indication inversion is carried out based on the obtained three-dimensional amplitude-preserving offset stacking data, the obtained layer grid model Layer_Model, and the time domain logging impedance curve IMP_time. The key inversion parameters include smoothing radius, effective sample number and cutoff frequency. Among them, the smoothing radius mainly affects the lateral connectivity of the inversion body. The larger the smoothing radius, the better the connectivity of the inversion body. The value range is, for example, between 2-3. The effective sample number mainly represents the degree of influence of the spatial variation of the seismic waveform on the reservoir. As the effective sample number increases, the seismic waveform's representation of the reservoir becomes better. However, when the effective sample number increases to a certain extent, the seismic waveform's representation of the reservoir no longer changes. The number of samples at this time is the effective sample number. The cutoff frequency controls the frequency band of the inversion, including high-pass and high-cutoff frequencies in the low-frequency band, and low-cutoff, low-pass, high-pass, and high-cutoff frequencies in the high-frequency band. The cutoff frequency value should be selected based on the amplitude spectrum of the three-dimensional seismic data volume. The principle of its value is to effectively compensate for the low-frequency and high-frequency parts of the seismic data. Under this principle, the cutoff frequency value is fine-tuned to find the parameters that minimize the inversion error. Through repeated parameter testing, the key parameters that can minimize the waveform indication inversion error of the target layer and achieve the best inversion results are found. Then, using the key inversion parameters determined in the parameter test, according to the above-mentioned waveform indication inversion method, the optimal wave impedance estimation is performed on the data of each sampling point of the target layer, thereby obtaining the wave impedance volume of the waveform indication inversion. Dividing it by the density average Den_average, the layer velocity volume of the waveform indication inversion can be obtained. However, based on the principle of the waveform indication inversion method, since the low-frequency and medium-frequency information of the three-dimensional seismic data body is applied in its inversion process, the obtained waveform indication inversion result contains certain low-frequency and medium-frequency components. In order to remove the low-frequency and medium-frequency components and retain only the high-frequency components of the waveform indication inversion, the waveform indication inversion can be high-pass filtered, that is, the frequency components with a frequency value below FH (FH can be determined based on the amplitude spectrum of the three-dimensional seismic data body) are filtered out, and only the high-frequency part is retained, thereby obtaining a data body that only reflects the high-frequency components, that is, the high-frequency band velocity body Vel_High obtained by the waveform indication inversion.

[0116] In this example, inversion parameters were tested based on the time-domain impedance curves IMP_time from 25 preselected wells, a 3D seismic data volume, and a layer grid model Layer_Model created from five interpreted horizons. Based on the principle of minimizing inversion error, the key inversion parameters were determined as follows: a smoothing coefficient of 2, a valid sample number of 5, high-pass and high-cut frequencies of 5Hz and 10Hz for the low-frequency band, and low-cut, low-pass, high-pass, and high-cut frequencies of 40Hz, 50Hz, 90Hz, and 100Hz for the high-frequency band. Next, waveform inversion was performed using these inversion parameters to obtain a waveform-indicator inversion impedance data volume. This data volume was then divided by the density mean of 2.5 to obtain a velocity volume. This velocity volume was high-pass filtered to remove frequency components below 60Hz (FH = 60Hz) and retain only the high-frequency components, resulting in the high-frequency velocity volume Vel_High. As shown in Figure 6, this example profile reflects the high-frequency velocity variations of the target layer and is able to fully identify the vertical and lateral velocity variations of some thin layers.

[0117] The low-frequency velocity body Vel_Low, medium-frequency velocity body Vel_Medium and high-frequency velocity body Vel_High of the target layer are obtained by frequency band inversion, and the (layer) velocity model Vel_Model is calculated. The (layer) velocity model Vel_Model of the target layer is established by using data fusion technology.

[0118] During the mid-frequency inversion of the target layer (using spectral simulation inversion to invert the 3D seismic data volume and obtain a mid-frequency velocity volume), a certain amount of low-frequency information is added to more intuitively analyze the vertical and lateral variations of the reservoir and perform effective inversion quality control. Therefore, the low-frequency information in the mid-frequency velocity volume Vel_Medium can be removed. In other words, the mid-frequency velocity volume Vel_Medium is high-pass filtered (the parameters of the high-pass filter can be determined based on the spectrum of the 3D seismic data volume, for example, the high-pass frequency is 15 Hz). Through high-pass filtering, a (layer) velocity volume Vel_Mediumfilter is obtained that retains only the mid-frequency information.

[0119] Then, according to the preset weight values, the low-frequency velocity body Vel_Low, the medium-frequency velocity body Vel_Medium, and the high-frequency velocity body Vel_High of the target layer are weighted to obtain the fused velocity model. The preset weight values ​​can be determined by analyzing the consistency between the inverted layer velocity data and the well. For example, they can be values ​​between 1 and 2, with a high weight for a high degree of consistency and a low weight for a low degree of consistency. For example, by analyzing the consistency between the inversion results and the well, the weight values ​​W_Low for the low-frequency velocity body, W_Medium for the medium-frequency velocity body, and W_High for the high-frequency velocity body are determined. For any sampling point i in the target layer, the corresponding layer velocity values ​​of the low-frequency velocity body Vel_Low(i), the medium-frequency velocity body Vel_Mediumfilter(i), and the high-frequency velocity body Vel_High(i) are read in sequence. According to the above formula (12), the layer velocity value Vel_Model(i) after the data fusion of the point can be obtained.

[0120] Based on this, the data fusion of each sampling point in the target layer can be realized to obtain the target layer velocity model Vel_Model.

[0121] In this example, the obtained mid-frequency velocity body Vel_Medium is high-pass filtered, that is, the frequency components below 15 Hz are filtered out, and only the frequency components above 15 Hz in the Vel_Medium layer velocity body are retained. The high-pass filtered mid-frequency velocity body Vel_Mediumfilter is obtained by high-pass filtering. By analyzing the inversion results and the degree of fit between the wells, the weight values ​​of the low-frequency W_Low, the mid-frequency W_Medium, and the high-frequency W_High are determined (for example, all are 1). Then, using formula (12), the low-frequency velocity body Vel_Low, the mid-frequency velocity body Vel_Mediumfilter, and the high-frequency velocity body Vel_High are fused to obtain the velocity model Vel_Model of the target layer. As shown in Figure 7, the example profile can not only clearly reflect the longitudinal and lateral variation characteristics of the target layer velocity, but also can accurately describe the velocity variation characteristics of the interior and boundary of the limestone high-speed anomaly body, and can also accurately reveal the velocity changes of some thin layers.

[0122] The reliability of the velocity model is analyzed by comparing the velocity model Vel_Model and the depth error at the well point after time-depth conversion. The (layer) velocity model Vel_tradition obtained by conventional methods is compared with the (layer) velocity model Vel_Model of the embodiment of the present invention. As shown in the comparison of Figures 7 and 8, the (layer) velocity model Vel_tradition cannot identify the velocity changes of special geological bodies (in this example, the limestone high-speed anomaly body), and the vertical and horizontal changes in the velocity of the target layer are not finely characterized; while the (layer) velocity model Vel_Model of the embodiment of the present invention can reflect the fine changes in velocity in the vertical and horizontal directions, reflect the geological laws of the study area, and can also well characterize the internal and boundary velocity change characteristics of the special geological body (in this example, the limestone high-speed anomaly body), and has a high degree of consistency with the velocity on the drilling, thereby verifying the reliability of the velocity model Vel_Model.

[0123] Next, based on the obtained (layer) velocity model Vel_tradition and the (layer) velocity model Vel_Model of the embodiment of the present invention, the target layer is converted into a time-depth. In this example, the KT1 top surface interpretation layer of the target layer is selected for time-depth conversion, and the depth error at each well point after the conversion between the layer velocity model Vel_tradition and the (layer) velocity model Vel_Model of the embodiment of the present invention is calculated, as shown in Table 1.

[0124] Table 1 - Depth error comparison of each well point at the target layer after time-depth conversion

[0125] Table 1 shows the errors between the drilling depths above sea level and the time-to-depth conversion depths for 14 wells in the KT1 geological layer. A comparative analysis of the time-to-depth conversion errors at each well point using two different velocity models shows that the depth error at the well point obtained after time-to-depth conversion using the (layer) velocity model Vel_Model according to the present embodiment of the present invention is 0.71%, significantly lower than the 2.86% depth error obtained using the (layer) velocity model Vel_tradition. This further validates the reliability of the velocity model established using frequency-band inversion according to the present embodiment of the present invention.

[0126] Based on this, the method for constructing a velocity model in an embodiment of the present invention effectively achieves the precise establishment of a velocity model for formations with special geological bodies, improving the velocity model's representation of the vertical and lateral variations in reservoir velocity. The velocity model not only reflects geological laws but also effectively characterizes the velocity variations within and at the boundaries of special geological bodies (in this example, limestone high-velocity anomalies), providing reliable velocity data for comprehensive geological studies such as time-depth conversion, detailed structural analysis, and well placement. The embodiments of the present invention are easy to implement and highly reliable, and have positive implications for the expansion and efficient advancement of exploration related to special geological bodies.

[0127] An embodiment of the present invention also provides a control device for constructing a velocity model, the control device comprising: a memory, a processor, and a computer program stored in the memory and executable on the processor, the processor executing the computer program to implement the above-mentioned method for constructing a velocity model.

[0128] An embodiment of the present invention further provides a machine-readable storage medium having instructions stored thereon, wherein the instructions enable a machine to execute the above-mentioned method for constructing a velocity model.

[0129] Those skilled in the art will appreciate that the embodiments of the present application can be provided as methods, systems, or computer program products. Therefore, the present application can adopt the form of a complete hardware embodiment, a complete software embodiment, or an embodiment in combination with software and hardware. Moreover, the present application can adopt the form of a computer program product implemented on one or more computer-usable storage media (including but not limited to magnetic disk storage, CD-ROM, optical storage, etc.) that contain computer-usable program code.

[0130] The present application is described with reference to the flow chart and / or block diagram of the method, device (system), and computer program product according to the embodiment of the present application. It should be understood that each flow process and / or box in the flow chart and / or block diagram and the combination of the flow process and / or box in the flow chart and / or block diagram can be realized by computer program instructions. These computer program instructions can be provided to a processor of a general-purpose computer, a special-purpose computer, an embedded processing machine or other programmable data processing device to produce a machine, so that the instructions executed by the processor of the computer or other programmable data processing device produce a device for realizing the function specified in one flow chart flow or multiple flows and / or one box or multiple boxes of the block diagram.

[0131] These computer program instructions may also be stored in a computer-readable memory that can direct a computer or other programmable data processing device to operate in a specific manner, so that the instructions stored in the computer-readable memory produce a product including an instruction device that implements the functions specified in one or more processes in the flowchart and / or one or more boxes in the block diagram.

[0132] These computer program instructions can also be loaded onto a computer or other programmable data processing device so that a series of operating steps are executed on the computer or other programmable device to produce a computer-implemented process, so that the instructions executed on the computer or other programmable device provide steps for implementing the functions specified in one or more processes in the flowchart and / or one or more boxes in the block diagram.

[0133] In a typical configuration, a computing device includes one or more processors (CPUs), input / output interfaces, network interfaces, and memory.

[0134] The memory may include non-permanent memory in a computer-readable medium, random access memory (RAM) and / or non-volatile memory in the form of read-only memory (ROM) or flash RAM. The memory is an example of a computer-readable medium.

[0135] Computer-readable media includes permanent and non-permanent, removable and non-removable media that can be implemented by any method or technology to store information. The information can be computer-readable instructions, data structures, program modules or other data. Examples of computer storage media include, but are not limited to, phase change memory (PRAM), static random access memory (SRAM), dynamic random access memory (DRAM), other types of random access memory (RAM), read-only memory (ROM), electrically erasable programmable read-only memory (EEPROM), flash memory or other memory technology, compact disc read-only memory (CD-ROM), digital versatile disc (DVD) or other optical storage, magnetic cassettes, magnetic tape, magnetic disk storage or other magnetic storage devices or any other non-transmission media that can be used to store information that can be accessed by a computing device. As defined herein, computer-readable media does not include transitory computer-readable media (transitory media), such as modulated data signals and carrier waves.

[0136] It should also be noted that the terms "comprises," "includes," or any other variations thereof are intended to encompass non-exclusive inclusion, such that a process, method, commodity, or apparatus that includes a series of elements includes not only those elements but also other elements not explicitly listed, or includes elements inherent to such process, method, commodity, or apparatus. In the absence of further limitations, an element defined by the phrase "comprises a ..." does not exclude the presence of other identical elements in the process, method, commodity, or apparatus that includes the element.

[0137] The above are merely embodiments of the present application and are not intended to limit the present application. For those skilled in the art, the present application may have various changes and variations. Any modifications, equivalent replacements, improvements, etc. made within the spirit and principles of the present application should all be included within the scope of the claims of the present application.

Claims

1. A method for constructing a velocity model, characterized in that: The method for constructing a velocity model comprises: Acquire 3D seismic data volume of a preset area and logging data of pre-selected wells; Processing the three-dimensional seismic data volume and the well logging data to obtain the time-depth relationship and layer interpretation data of the pre-selected wells of the target layer; Based on the obtained time-depth relationship of the pre-selected wells of the target layer and the stratum interpretation data, the three-dimensional seismic data volume is subjected to frequency-band inversion through a preset frequency-band inversion strategy to obtain velocity volumes in different frequency bands; and The velocity bodies of the different frequency bands are fused to obtain a velocity model of the target layer.

2. The method for constructing a velocity model according to claim 1, characterized in that: Before acquiring the three-dimensional seismic data volume of the preset area and the logging data of the pre-selected wells, the method for constructing the velocity model further includes: Determining a range of the target layer; Acquiring well logging data of the wells in the preset area; and For the wells drilled through the target layer, wells having acoustic wave curves, density curves and well diameter curves are screened as the pre-selected wells according to the corresponding well logging data.

3. The method for constructing a velocity model according to claim 2, characterized in that: After acquiring the three-dimensional seismic data volume of the preset area and the logging data of the pre-selected wells, the method for constructing the velocity model further includes: The acoustic wave curve of the preselected well is converted into a velocity curve.

4. The method for constructing a velocity model according to claim 1, characterized in that: Before obtaining the time-depth relationship and the stratigraphic interpretation data of the preselected well of the target layer, the method for constructing the velocity model further includes: For the pre-selected wells of the target layer, the histogram of the single well is analyzed by the well logging data to determine the vertical variation characteristics of the lithology and velocity of the target layer; Comparing and analyzing the velocity curves of the connected wells through the logging data, to determine the lateral variation characteristics of the velocity of the target layer; and The velocity variation characteristics and the geological body type of the target layer are determined through the longitudinal velocity variation characteristics of the target layer and the lateral velocity variation characteristics of the target layer.

5. The method for constructing a velocity model according to claim 1, characterized in that: The logging data of the pre-selected wells are processed to obtain the time-depth relationship of the pre-selected wells of the target layer, including: Acquiring an acoustic wave curve and a density curve from the logging data of the preselected well; and A synthetic seismic record is produced by using the acoustic wave curve and the density curve, a corresponding relationship is established between the seismic reflection interface in the time domain and the geological stratification in the depth domain, and the time-depth relationship of the preselected well is obtained.

6. The method for constructing a velocity model according to claim 5, characterized in that: The processing of the three-dimensional seismic data volume and the well logging data to obtain horizon interpretation data includes: Analyzing the logging data of the preselected well to obtain the velocity variation characteristics of the target layer; By analyzing the velocity variation characteristics of the target layer and determining T geological layers of the target layer according to the velocity difference of the velocity interface; By making synthetic seismic records, the T geological layers of the target layer are calibrated, and the characteristics of the seismic reflection event axis corresponding to each geological layer are determined; and Based on the three-dimensional seismic data volume, the seismic reflection event axes corresponding to each geological layer are tracked and interpreted to obtain the stratigraphic interpretation data of the T geological layers of the target layer.

7. The method for constructing a velocity model according to claim 1, characterized in that: The method of performing frequency-band inversion on the three-dimensional seismic data volume by using a preset frequency-band inversion strategy to obtain velocity volumes in different frequency bands includes: Inverting the three-dimensional seismic data volume by Kriging interpolation method to obtain a low-frequency velocity volume; Inverting the three-dimensional seismic data volume by spectral simulation inversion method to obtain a mid-frequency velocity volume; and The three-dimensional seismic data volume is inverted by a waveform indication inversion method to obtain a high-frequency band velocity volume.

8. The method for constructing a velocity model according to claim 7, characterized in that: Before performing frequency-band inversion on the three-dimensional seismic data volume using a preset frequency-band inversion strategy, the method for constructing a velocity model further includes: Determine T interpretation horizons of the target layer from bottom to top; Between each two adjacent interpretation horizons, according to the contact relationship between the upper and lower layers and the preset subdivision ratio C, C interpolation sublayers are obtained by interpolation between each two adjacent interpretation horizons; and A stratigraphic framework model of the target layer is constructed by using the T interpreted layers and the C interpolated sublayers between each two adjacent interpreted layers.

9. The method for constructing a velocity model according to claim 8, characterized in that: The three-dimensional seismic data volume is inverted by the Kriging interpolation method to obtain a low-frequency velocity volume, including: Acquire a velocity curve in the logging data of the preselected well, and perform low-pass filtering on the velocity curve to retain a low-frequency velocity curve within a first cutoff frequency value; Based on the obtained time-depth relationship of the pre-selected well, the filtered low-frequency velocity curve is converted from the depth domain to the time domain to obtain a low-frequency velocity curve in the time domain; and Based on the stratigraphic grid model of the target layer, the lateral constraint conditions of the low-frequency velocity curve are established, and the low-frequency velocity body is obtained by performing logging interpolation under the lateral constraint conditions on the three-dimensional seismic data body through the Kriging interpolation method.

10. The method for constructing a velocity model according to claim 9, characterized in that: The lateral constraint condition of the low-frequency velocity curve is established based on the stratigraphic frame model of the target layer, so as to perform well logging interpolation under the lateral constraint condition on the three-dimensional seismic data body by the Kriging interpolation method to obtain the low-frequency velocity body, including: For each horizon in the horizon grid model, based on the observed value of each well of the preselected wells, an estimated value at each CDP of the horizon is obtained, wherein the low-frequency velocity value read at each well is used as the observed value of each well. The low-frequency velocity body is obtained by the estimated value at each CDP of each layer. Each of the horizons includes the T interpretation horizons and C interpolation sub-layers between each two adjacent interpretation horizons.

11. The method for constructing a velocity model according to claim 7, characterized in that: The three-dimensional seismic data volume is inverted by the spectrum simulation inversion method to obtain a mid-frequency velocity volume, including: Obtaining a velocity curve from the logging data of the preselected well, and constructing a corresponding wave impedance curve; Based on the obtained time-depth relationship of the pre-selected well of the target layer, converting the wave impedance curve of the pre-selected well from the depth domain to the time domain to obtain a corresponding time domain wave impedance curve; and The time domain wave impedance curve is superimposed on the three-dimensional seismic data volume, and inversion is performed through a spectral simulation inversion method to obtain the mid-frequency band velocity volume.

12. The method for constructing a velocity model according to claim 11, characterized in that: The step of superimposing the time domain wave impedance curve onto the three-dimensional seismic data volume and inverting the data using a spectrum simulation inversion method to obtain the mid-frequency velocity volume comprises: For each preselected well of the target layer, a spectrum analysis is performed on the seismic traces near the well and the time domain wave impedance curve of the preselected well, a frequency matching operator is calculated, and the frequency matching operator is applied to the three-dimensional seismic data volume to obtain an inversion result of spectrum simulation. According to the corresponding relationship between the reflection coefficient and the wave impedance, the wave impedance inversion result is obtained from the inversion result of the spectrum simulation; and The mid-frequency velocity body is obtained by the ratio of the wave impedance inversion result of the target layer to the first density mean.

13. The method for constructing a velocity model according to claim 8, characterized in that: The three-dimensional seismic data volume is inverted by the waveform indication inversion method to obtain a high-frequency velocity volume, including: For the wellside seismic traces of each well in the n pre-selected wells, under the stratigraphic frame model of the target layer, the seismic waveform to be predicted is compared with the waveform of the wellside seismic traces of each well, and K waveform simulation samples of the wellside seismic traces are selected according to two factors: waveform similarity and spatial distance, and an initial impedance model is constructed through the corresponding logging data; Based on the characteristics of the initial impedance model and according to a preset high-pass frequency parameter value and a preset high-cut frequency parameter value, high-frequency band information is continuously filtered out, and waveform simulation samples of each well-side seismic trace are compared to retain a preset deterministic frequency band component; and According to the waveform simulation samples and spatial structure characteristics of the K well-side seismic traces, the estimated value of the unknown point is determined to obtain the wave impedance inversion result of the waveform indication inversion, and the high-frequency band velocity body is obtained by the ratio of the wave impedance inversion result of the target layer to the second density mean.

14. The method for constructing a velocity model according to claim 7, characterized in that: The fusing the velocity bodies of the different frequency bands to obtain the velocity model of the target layer includes: The velocity model of the target layer is obtained by weighted sum of the low-frequency velocity body, the mid-frequency velocity body and the high-frequency velocity body.

15. A control device for constructing a velocity model, characterized in that: The control device comprises: a memory, a processor, and a computer program stored in the memory and executable on the processor, wherein the processor executes the computer program to implement the method for constructing a velocity model according to any one of claims 1 to 14.

16. A machine-readable storage medium, characterized in that: The machine-readable storage medium stores instructions, which enable the machine to execute the method for constructing a velocity model according to any one of claims 1-14.

Images