Method for constructing speed model, control device and storage medium

Abstract

The embodiment of the invention provides a method for constructing a speed model, a control device and a storage medium, and belongs to the technical field of geological exploration. The method comprises the following steps: acquiring a three-dimensional seismic data volume of a preset area and logging data of a pre-selected well; processing the three-dimensional seismic data volume and the logging data to obtain a time-depth relationship and horizon interpretation data of a pre-selected well of a target stratum; based on the obtained time-depth relationship and horizon interpretation data of the pre-selected well of the target stratum, through a preset multi-frequency-band inversion strategy, performing multi-frequency-band inversion on the three-dimensional seismic data volume to obtain speed volumes of different frequency bands; and fusing the velocity bodies of different frequency bands to obtain a velocity model of the target stratum. According to the method, logging data and a three-dimensional seismic data body are fully utilized, frequency-band-divided inversion is carried out on the three-dimensional seismic data body through a preset frequency-band-divided inversion strategy and a data fusion technology, speed bodies of different frequency bands are obtained, and a speed model capable of finely representing internal and boundary speed change characteristics of a special geological body is constructed.

Description

Technical Field

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

[0002] Velocity parameters are the basic physical properties of underground strata and important parameters of seismic exploration. They affect all aspects of seismic exploration and play a vital role in oil and gas exploration and development. Strata or geological bodies of different ages and lithologies usually have different layer velocities; underground reservoir pores will also show different degrees of velocity differences due to different fluids (gas, oil, water) or abnormalities and high-pressure closure of geological bodies. Therefore, an accurate velocity model that conforms to geological laws is not only the key to prestack migration imaging and time-depth conversion, but also affects many aspects of exploration and development, such as structural accuracy, well depth implementation, stratum thickness, reservoir prediction, and pressure interface prediction.

[0003] With the continuous deepening of exploration and development, the geological conditions are becoming more and more complex. We often face some localized special geological bodies accompanied by salt domes, belt reef structures, igneous intrusive bodies and carbonate anomalies. These special geological bodies are usually quite different from the surrounding rocks in terms of lithology, velocity and thickness, and have strong heterogeneity. The current velocity modeling method cannot reflect the local changes in the velocity of special geological bodies, resulting in large errors in time-depth conversion, low structural implementation accuracy, and easy formation of structural illusions, which seriously restricts the exploration and development process of geological targets related to special geological bodies. Therefore, complex geological conditions pose new challenges to velocity modeling.

[0004] At present, the methods for establishing layer velocity models mainly include: 1) using pre-stack time migration velocity or stacking velocity to calculate layer velocity through the Dix formula. This method is based on a uniform medium or horizontal layered medium model and is estimated from the seismic processing velocity spectrum. The velocity change trend and relative relationship can be obtained. This method has poor accuracy for sedimentary stratigraphic sequences with widespread spatial anisotropy. 2) Grid tomography velocity modeling method. This method uses data-driven and computer-automatic iterative solution methods to obtain more detailed velocity changes, but this method has a large amount of calculations and is difficult to control velocity errors under complex geological conditions due to the lack of human intervention. 3) Using seismic logging and acoustic logging data from drilled wells, combined with the interpretation of geological strata, layer velocity models are established. This method can obtain high-precision layer velocity data that truly reflects the lithological characteristics of the formation at the location of the control well, but it is controlled by the plane distribution of the drilled wells and the drilling density, especially in areas with few wells or no wells. The velocity prediction accuracy is poor. Summary of the invention

[0005] 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 boundary of a special geological body.

[0006] 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 of a preset area and logging data of a pre-selected well; processing the three-dimensional seismic data volume and the logging data to obtain the time-depth relationship and stratigraphic interpretation data of the pre-selected well of the target layer; based on the obtained time-depth relationship and stratigraphic interpretation data of the pre-selected 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 different frequency bands to obtain a velocity model of the target layer.

[0007] 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.

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

[0009] Optionally, before obtaining the time-depth relationship and stratigraphic interpretation data of the pre-selected wells of the target layer, the method for constructing a velocity model also includes: for the pre-selected wells of the target layer, analyzing the bar chart 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.

[0010] 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 the 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.

[0011] 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 pre-selected well to obtain velocity variation characteristics of the target layer; determining T geological layers of the target layer by analyzing the velocity variation 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 event axes corresponding to each geological layer; and completing the tracking and interpretation of the seismic reflection event 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.

[0012] Optionally, the three-dimensional seismic data volume is inverted in frequency bands by a preset frequency band inversion strategy to obtain velocity volumes in different frequency bands, including: inverting the three-dimensional seismic data volume by Kriging interpolation method to obtain a low-frequency band velocity volume; inverting the three-dimensional seismic data volume by spectral simulation inversion method to obtain a mid-frequency band velocity volume; and inverting the three-dimensional seismic data volume by waveform indication inversion method to obtain a high-frequency band velocity volume.

[0013] Optionally, before performing frequency band inversion on the three-dimensional seismic data volume 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 stratigraphic grid model of the target layer from the T interpretation layers and the C interpolated sublayers between each two adjacent interpretation layers.

[0014] 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 pre-selected well, and low-pass filtering 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 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, lateral constraints of the low-frequency velocity curve are established to obtain the low-frequency velocity body by performing logging interpolation on the three-dimensional seismic data body under the lateral constraints by the Kriging interpolation method.

[0015] Optionally, the lateral constraints of the low-frequency velocity curve are established based on the stratum grid model of the target layer, 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 stratum 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 stratum, 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 stratum is obtained through the estimated value at each CDP of the stratum; and the low-frequency band velocity body is obtained through the low-frequency velocity value of each stratum, wherein each stratum includes the T interpreted strata and the C interpolated sublayers between each two adjacent interpreted strata.

[0016] 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.

[0017] Optionally, the time domain wave impedance curve is superimposed on the three-dimensional seismic data volume, and inverted by spectral simulation inversion method to obtain the mid-frequency band velocity body, including: for each preselected well of the target layer, spectral analysis is performed on the seismic trace beside the well and the time domain wave impedance curve of the preselected well, and 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 spectral simulation, and according to the corresponding relationship between reflection coefficient and wave impedance, a wave impedance inversion result is obtained from the inversion result of the spectral simulation; and the mid-frequency band velocity body is obtained by the ratio of the wave impedance inversion result of the target layer to the first density mean.

[0018] Optionally, the three-dimensional seismic data body is inverted by a waveform indication inversion method to obtain a high-frequency velocity body, including: for the wellbore seismic traces 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 traces 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 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 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.

[0019] 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 taking the weighted sum of the low-frequency band velocity body, the mid-frequency band velocity body and the high-frequency band velocity body.

[0020] 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.

[0021] An embodiment of the present invention further provides a machine-readable storage medium, on which instructions are stored, and the instructions enable a machine to execute the above-mentioned method for constructing a velocity model.

[0022] Through the above technical scheme, the embodiment of the present invention provides a method for constructing a velocity model for more accurate comprehensive seismic geological research such as time-depth conversion, thickness calculation, reservoir prediction, and well location implementation related to special geological bodies. Based on the analysis of the geological characteristics of special geological bodies, fine time-depth relationship and layer interpretation data are obtained through fine seismic geological layer calibration. On this basis, the well logging data and three-dimensional seismic data body are fully utilized, and the three-dimensional seismic data body is subjected to frequency band inversion through the preset frequency band inversion strategy and data fusion technology to obtain velocity bodies of different frequency bands, and the frequency band velocities of low frequency, medium frequency and high frequency bands are used to construct a velocity model that can finely characterize the velocity change characteristics of the interior and boundary of the special geological body, and provide reliable velocity data for comprehensive geological research such as time-depth conversion, fine structural implementation, and well location deployment, so as to help expand and efficiently promote exploration in related fields of special geological bodies.

[0023] 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

[0024] 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 specific embodiments, they are used to explain the embodiments of the present invention, but do not constitute a limitation on the embodiments of the present invention. In the accompanying drawings:

[0025] Figure 1 It is a schematic flow chart of a method for constructing a velocity model provided by an embodiment of the present invention;

[0026] Figure 2 It is a schematic diagram of the velocity variation characteristics of an example geological body.

[0027] Figure 3 is a schematic cross-sectional view of example offset stacked seismic data;

[0028] Figure 4 It is a velocity profile diagram of a low-frequency velocity body of an example target layer;

[0029] Figure 5 It is a velocity profile diagram of a mid-frequency velocity body in an example target layer;

[0030] Figure 6 It is a velocity profile diagram of the high-frequency velocity body of the target layer;

[0031] Figure 7 is a velocity profile diagram of an example velocity model; and

[0032] Figure 8 is a schematic diagram of a velocity profile of an example conventional velocity model. DETAILED DESCRIPTION

[0033] The specific implementation of the embodiment of the present invention is described in detail below in conjunction with the accompanying drawings. It should be understood that the specific implementation described here 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.

[0034] The existing velocity model is suitable for geological conditions with relatively stable lateral distribution of strata. However, for geological conditions with relatively complex geological structures, for example, with the development of special geological bodies such as salt domes, belt reef structures and carbonate rock anomalies, the existing velocity modeling cannot accurately characterize the velocity change characteristics inside and at the boundaries of special geological bodies, which affects the accuracy of structural characteristics and restricts the exploration of related fields of special geological bodies. Therefore, for exploration areas with the development of special geological bodies, a velocity model that can accurately characterize the velocity changes inside and at the boundaries of special geological bodies is needed.

[0035] Figure 1 This is a flow chart of a method for constructing a velocity model provided by an embodiment of the present invention. Please refer to Figure 1 , the method for constructing a velocity model may include the following steps:

[0036] Step S110: Acquire a three-dimensional seismic data volume of a preset area and logging data of a pre-selected well.

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

[0038] 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 having acoustic wave curves, density curves and well diameter curves as the pre-selected wells based on corresponding logging data.

[0039] By way of example, wells that penetrate the target layer are selected in a preset area, and among these wells, n wells that have acoustic wave curves, density curves and well diameter curves are preferably used as the pre-selected wells.

[0040] 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.

[0041] Preferably, the acoustic wave curves corresponding to the n pre-selected wells are converted into velocity curves by the following formula:

[0042] Vel=1000000 / DT (1)

[0043] Among them, DT represents the sound wave curve, the unit is microseconds / meter, and Vel represents the velocity curve, the unit is meter / second.

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

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

[0046] Preferably, before obtaining the time-depth relationship and stratigraphic interpretation data of the pre-selected wells of the target layer, the method for constructing a velocity model also includes: for the pre-selected wells of the target layer, analyzing the bar chart 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.

[0047] By way of example, the velocity variation characteristics of the target layer and the type of special geological body can be obtained through comparative analysis of single wells and wells in the target layer. The lithology and vertical velocity variation characteristics of the target layer are determined through the histogram analysis of the single well; the lateral velocity variation characteristics of the target layer are analyzed through well-by-well comparative analysis of the velocity curve of the target layer, and the type of special geological bodies (for example, special geological bodies such as salt domes, belt-like reef structures, and carbonate rock anomalies) in the target layer are determined, and their velocity variation characteristics are obtained, providing a basis for the establishment of a velocity model.

[0048] 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 the 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.

[0049] Preferably, step S120 may include: 2. Analyzing the three-dimensional seismic data body and the logging data to obtain stratigraphic interpretation data, including: analyzing the 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 oscillators corresponding to each geological layer; and based on the three-dimensional seismic data body, completing the tracking and interpretation of the seismic reflection oscillators corresponding to each geological layer to obtain the stratigraphic interpretation data of the T geological layers of the target layer.

[0050] 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 through 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.

[0051] 2. According to the analysis results of the velocity change characteristics of the target layer determined above, according to the velocity difference of the velocity interface, determine T geological layers of the target layer with obvious velocity interface changes; through the synthetic seismic records produced in 1), calibrate the T geological layers of the target layer segment, determine the characteristics of the seismic reflection phase axis corresponding to each geological layer, and track and interpret the seismic reflection phase axis corresponding to the T geological layers based on the three-dimensional seismic data volume, so as to obtain the stratigraphic interpretation data of the T geological layers of the target layer.

[0052] Step S130: Based on the obtained time-depth relationship of the preselected 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.

[0053] Preferably, before step S130, 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 stratigraphic grid model of the target layer from the T interpretation layers and the C interpolated sublayers between each two adjacent interpretation layers.

[0054] To illustrate, the framework model of the target layer is established through T interpretation layers of the target layer, including: selecting T interpretation layers of the target layer segment from bottom to top (for example, T=5), and interpolating C interpolated sublayers between each two adjacent interpretation layers according to the contact relationship between the upper and lower layers and the preset subdivision ratio C (for example, C=25), so there are a total of 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, and the layer grid model Layer_Model can be used as a lateral constraint condition.

[0055] 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.

[0056] 1. The inversion of the three-dimensional seismic data volume by the Kriging interpolation method to obtain a low-frequency velocity volume comprises: obtaining a velocity curve in the well logging data of the preselected well, and performing 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 preselected well, converting the filtered low-frequency velocity curve from the depth domain to the time domain to obtain a low-frequency velocity curve in the time domain; and establishing a lateral constraint condition of the low-frequency velocity curve based on the stratigraphic framework model of the target layer, so as to perform well logging interpolation under the lateral constraint condition on the three-dimensional seismic data volume by the Kriging interpolation method to obtain the low-frequency velocity volume.

[0057] To illustrate by way of example, the velocity curves of n pre-selected 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 the first cut-off frequency value FL (for example, FL=15 Hz) are filtered out, and only the low-frequency velocity curves within FL are retained; and according to the time-depth relationship of each pre-selected well obtained as described above, 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; a lateral constraint condition is established through the seismic interpretation horizon of the target layer, that is, the horizon grid model Layer_Model of the target layer is used as the lateral constraint condition; based on the filtered low-frequency velocity curve and the horizon grid model Layer_Model, the low-frequency (layer) velocity body of the target layer (represented by Vel_Low) can be obtained through the ordinary Kriging interpolation method.

[0058] Preferably, the lateral constraints of the low-frequency velocity curve are established based on the stratum grid model of the target layer, 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 stratum 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 stratum, 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 stratum is obtained through the observation value at each CDP of the stratum and the estimated value at each CDP of each stratum; and the low-frequency band velocity body is obtained through the low-frequency velocity value of each stratum, wherein each stratum includes the T interpreted strata and the C interpolated sublayers between each two adjacent interpreted strata.

[0059] 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 pre-selected wells (for example, n=25), a low-frequency velocity value can be read at each pre-selected well along the layer. Then, the low-frequency velocity values ​​read from the n pre-selected wells can be regarded as n observations, that is, Z(x1), Z(x2), Z(x3), ..., Z(x4). n ), the estimated value of the low-frequency velocity at a 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:

[0060]

[0061] Among them, α i(i=1,2,…,n) is the weight coefficient to be determined. According to the principle of Kriging interpolation, under the premise of ensuring that the estimator is unbiased and the estimated variance is minimal, n weight coefficient values ​​can be obtained by solving the following set of equations.

[0062]

[0063]

[0064] Among them, γ(x i ,x j ) is the sampling point x i and x j The linear variation function value between. From equation group (3), we can get the weight coefficient α i By substituting the value 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).

[0065] 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.

[0066] 2. 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 pre-selected well and constructing a corresponding wave impedance curve; based on the obtained time-depth relationship of the pre-selected well in the target layer, the wave impedance curve of the pre-selected 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.

[0067] The inversion layer velocity body within the mid-frequency band of the target layer can be obtained by spectral simulation inversion method. Spectral simulation inversion is a frequency domain logging constrained wave impedance inversion technology. The well logging data and the three-dimensional seismic data body data are subjected to spectral analysis respectively, and a matching operator is designed in the frequency domain to match the seismic spectrum with the wave impedance spectrum of the well. Then, the frequency matching operator is applied to the three-dimensional seismic data body, and the time domain operator is inversely calculated to participate in the inversion operation to obtain the inversion data body.

[0068] To illustrate by way of example, based on the velocity curves of n pre-selected wells in a preset area, an impedance curve of the pre-selected well is constructed and used as input data for the spectrum simulation inversion method. Preferably, in order to facilitate the conversion of the impedance data volume obtained by inversion into a velocity data volume, when constructing the impedance curve, it is preferred that the density value is taken as a constant, that is, the average value Den_average of the density logging value of the target layer, and the velocity curve value of the pre-selected well is multiplied by the constant density value Den_average to obtain the impedance curve of the pre-selected well, and the impedance curve of the pre-selected well is converted from the depth domain to the time domain using the time-depth relationship of each pre-selected well that has been obtained to obtain the time domain impedance curve IMP_time. The time domain impedance curves of the n pre-selected wells and the three-dimensional seismic data volume (i.e., amplitude-preserving offset stacking data) can be inverted by spectrum simulation to obtain the mid-frequency band velocity volume of the target layer (represented by Vel_Medium).

[0069] Preferably, the time domain wave impedance curve is superimposed on the three-dimensional seismic data volume, and inverted by spectral simulation inversion method to obtain the mid-frequency band velocity volume, including: for each pre-selected 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 pre-selected well, and 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 spectral simulation, and according to the corresponding relationship between reflection coefficient and wave impedance, a wave impedance inversion result is obtained from the inversion result of the spectral simulation; and the mid-frequency band velocity volume is obtained by the ratio of the wave impedance inversion result of the target layer to the first density mean.

[0070] Following the above example, n wellside seismic traces are extracted corresponding to 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:

[0071] s i (t) = r i (t)*b(t) (4)

[0072] Transform it to the frequency domain and get the following formula:

[0073]

[0074] 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, is the phase spectrum of the wavelet, then according to formula (5), we can get the following formula:

[0075] S i (ω)=R i (ω)B i (ω) (6)

[0076] Taking the logarithm of formula (6), we can get the following formula:

[0077] oeLh i (ω)=lnR i (ω)+lnB i (ω) (7)

[0078] 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:

[0079]

[0080] 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 equation (7), we can get the frequency matching operator B i (ω).

[0081] 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:

[0082]

[0083] 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:

[0084] r i (t) = a i (t)*s i (t) (10)

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

[0086] 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) can obtain the inversion result of the reflection coefficient, and then use the general correspondence between the reflection coefficient and the wave impedance to convert the reflection coefficient into the wave impedance value to obtain the wave impedance inversion result of the seismic trace. The wave impedance obtained at this time is the relative wave impedance. In order to more intuitively analyze the inversion effect of the spectrum simulation wave impedance and facilitate the inversion quality control of the well, the obtained mid-frequency wave impedance inversion result is added with the low-frequency background to obtain the final spectrum simulation absolute wave impedance.

[0087] 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 medium-frequency velocity volume Vel_Medium of the target layer.

[0088] 3. 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, by analyzing the two factors of waveform similarity and spatial distance, K (K≤n) waveform simulation samples of the wellside seismic traces are selected, and an initial impedance model is constructed through 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 the waveform simulation samples of each wellside 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 wellside seismic traces, an estimated value of an unknown point is determined to obtain a wave impedance inversion result of waveform indication inversion, and the high-frequency velocity volume is obtained by the ratio of the wave impedance inversion result of the target layer to the second density mean.

[0089] 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 lithology 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:

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

[0091] 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 waveform similarity and spatial distance dual variables, K well bypasses with the most similar waveforms are selected as simulation samples, and their logging data are used to build the initial impedance model.

[0092] 2) For sample curve analysis, the impedance characteristics of the sample wells are statistically analyzed. According to 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.

[0093] 3) For spatial valuation, the original data and spatial structure characteristics of the sample wells are used to perform linear unbiased and optimal estimation of the unknown points.

[0094] 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:

[0095]

[0096] The wave impedance body obtained by waveform indication inversion can be obtained by formula (11), and it is divided by the density mean Den_average to obtain the layer velocity body. However, based on the principle of waveform indication inversion, the low-frequency and medium-frequency information of earthquakes are applied in its inversion process, so the inversion result contains certain low-frequency and medium-frequency components. In order to remove the low-frequency and medium-frequency components and only retain the high-frequency components of waveform indication inversion, preferably, the waveform indication inversion can be high-pass filtered to obtain a data body that only reflects the high-frequency components, that is, the high-frequency band velocity body obtained by waveform indication inversion (represented by Vel_High).

[0097] Preferably, based on the three-dimensional seismic data body, the stratum 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 range 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 tests, 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: the smoothing coefficient is 2, the number of valid samples is 5, the high-pass and high-cut frequencies of the low-frequency band are 5Hz and 10Hz respectively, and the low-cut, low-pass, high-pass and high-cut frequencies of the high-frequency band are 40Hz, 50Hz, 90Hz and 100Hz respectively. Based on the above determined inversion parameters, waveform indication inversion is performed to obtain the wave impedance data body of waveform indication inversion, and then it is divided by the density average Den_average (for example, Den_average = 2.5) to obtain the layer velocity body of waveform indication inversion. However, since the waveform indication inversion contains certain low-frequency and medium-frequency information in the calculation process, in order to obtain the inverted high-frequency band velocity value, it is necessary to perform high-pass filtering on the inversion result, that is, the frequency components with a frequency value below FH (FH can be determined according to the amplitude spectrum of the earthquake, for example, FH = 60Hz) are filtered out, and only the high-frequency part is retained, so that the high-frequency band velocity body Vel_High can be obtained.

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

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

[0100] As described in step S130, the velocity model of the target layer can be obtained by fusing the inverted velocity bodies of the low frequency, medium frequency and high frequency bands, that is, the (layer) velocity model Vel_Model is calculated by the low frequency band velocity body Vel_Low, the medium frequency band velocity body Vel_Medium and the high frequency band velocity body Vel_High of the target layer obtained by segmented inversion.

[0101] Preferably, in the process of inversion of the mid-frequency band of the target layer (inverting the three-dimensional seismic data body through the spectrum simulation inversion method to obtain the mid-frequency band velocity body), in order to more intuitively analyze the vertical and horizontal variation characteristics of the reservoir and perform effective inversion quality control, a certain amount of low-frequency information is added, so the low-frequency information in the mid-frequency band velocity body Vel_Medium can be removed. That is, the mid-frequency band velocity body Vel_Medium is high-pass filtered (the parameters of the high-pass filter can be determined based on the spectrum of the three-dimensional seismic data body, for example, the high-pass frequency is 15Hz), and through high-pass filtering, a (layer) velocity body Vel_Mediumfilter that only retains the mid-frequency band information is obtained.

[0102] Next, 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, it can be a value between 1-2, with a high weight for a large consistency and a low weight for a small consistency. For example, by analyzing the consistency between the inversion result and the well, the weight value W_Low of the low-frequency velocity body, the weight value W_Medium of the medium-frequency velocity body and the weight value W_High of the high-frequency velocity body are determined. For any sampling point i of the target layer segment, 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 turn, and the layer velocity value Vel_Model(i) after the data fusion of the point is obtained by the following formula:

[0103] Vel_Model(i)=Vel_Low(i)×W_Low+Vel_Mediumfilter(i)×W_Medium

[0104] +Vel_High(i)×W_High (12)

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

[0106] 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-to-depth conversion.

[0107] 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 longitudinal and lateral 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 region, 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, thereby verifying the reliability of the velocity model Vel_Model.

[0108] 2) The target layer is converted into time-depth 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 is calculated respectively; the comparative analysis of the depth errors at the well points obtained after the conversion of the two different velocity models shows 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 is significantly lower than the depth error value after the conversion of the (layer) velocity body Vel_tradition, which verifies the reliability of the (layer) velocity model Vel_Model of the embodiment of the present invention.

[0109] Accordingly, the embodiment of the present invention provides a method for constructing a velocity model for more accurate comprehensive seismic geological research such as time-depth conversion, thickness calculation, reservoir prediction, and well location implementation related to special geological bodies. Based on the analysis of the geological characteristics of special geological bodies, fine drilling time-depth relationship and layer interpretation data are obtained through fine seismic geological layer calibration. On this basis, the logging data and three-dimensional seismic data body are fully utilized, and the three-dimensional seismic data body is subjected to frequency band inversion through a preset frequency band inversion strategy and data fusion technology to obtain velocity bodies of different frequency bands. The frequency band velocity bodies of low frequency, medium frequency, and high frequency bands are used to construct a velocity model that can finely characterize the velocity change characteristics of the interior and boundary of the special geological body, and provide reliable velocity data for comprehensive geological research such as time-depth conversion, fine structural implementation, and well location deployment, so as to help expand and efficiently promote exploration in related fields of special geological bodies.

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

[0111] like Figure 3 As shown, the offset stacked seismic data after amplitude preservation processing is obtained; the logging data of the target layer in the area is obtained. All wells that penetrate the target layer are selected, and n wells that have acoustic wave curves, density curves and well diameter curves are screened according to the corresponding logging data as pre-selected wells. For example, 28 wells that penetrate the target layer in the current area are selected, and the logging data of these 28 wells that penetrate the target layer are analyzed, and 25 wells that have acoustic wave curves, density curves and well diameter curves are screened, that is, the pre-selected wells n = 25.

[0112] Next, the acoustic wave curve of each pre-selected well is converted into a velocity curve by using formula (1), that is, the acoustic wave curves, density curves, velocity curves and caliper curves of the 25 pre-selected wells are obtained, thereby completing the acquisition of the logging data of the pre-selected wells.

[0113] According to the comparative analysis of single wells and connected wells in the target layer, the velocity change characteristics of the target layer and the type of special geological body are obtained. By analyzing the single well column chart of the pre-selected wells, the lithological characteristics and longitudinal velocity change characteristics of the target layer are obtained. In this example, the target layer in this area is mainly carbonate sedimentary strata, and the lithology is mainly mudstone, limestone and marl-limestone interlayers. The mudstone velocity is low and the limestone velocity is high. The vertical velocity of the target layer shows the characteristics of interlayers of high-speed areas and low-speed areas. By comparing the velocity curves of the target layer with connected wells, the lateral variation characteristics of the velocity of the target layer are analyzed, and the types of special geological bodies in the target layer and their velocity characteristics are determined, providing a geological basis for the establishment of the velocity model.

[0114] Among them, through the comparison of well connections of some pre-selected wells, such as Figure 2As shown in the figure, the thickness and velocity variation characteristics of the three layers of P1, KT1 and MKT from top to bottom in the target layer can be obtained. Among them, the strata of KT1 and MKT are relatively stable in the lateral distribution, and the velocity and thickness changes are not very obvious, while the thickness and velocity differences of the P1 section at the top of the target layer are relatively obvious. The velocity of the P1 section has obvious differences both vertically and horizontally. The P1 section has the characteristics of low velocity in the lower part and high velocity in the middle and upper parts. The velocity of the lower part of the P1 section changes horizontally relatively steadily, with an average velocity of about 3200m / s, while the velocity of the middle and upper part of the P1 section changes horizontally greatly, with the velocity value varying between 3300m / s-5600m / s, the low-speed area is about 3300m / s, and the high-speed area can reach 5600m / s, and the stratum thickness of the P1 section in the high-speed area 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 has strong heterogeneity. There are obvious special geological bodies in some parts, 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.

[0115] By calibrating synthetic seismic records, the detailed time-depth relationship of the pre-selected wells is determined. In this example, synthetic records are made for 25 pre-selected wells, and the target layer is finely calibrated by well seismic, so that the detailed time-depth relationship of the 25 pre-selected wells can be obtained.

[0116] Next, according to the analysis results of the velocity change characteristics of the target layer determined above, the T seismic interpretation horizons of the target layer are determined. In this example, according to the analysis results of the velocity change characteristics of the target layer determined above, it is determined that there are 5 geological strata with obvious velocity changes in the target layer, namely T=5, which are P1 top, mudstone top, KT1 top, MKT top and KT2 top. According to the synthetic seismic records of 25 pre-selected wells, the reflection characteristics of each geological stratum on the seismic profile are calibrated, and the seismic reflection horizons corresponding to the 5 geological strata 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 stratum, such as Figure 3 As shown, through tracking and interpretation, the horizon interpretation data of 5 seismic reflection horizons corresponding to the target layer and the geological stratification are determined.

[0117] According to the logging data, time-depth relationship and the horizon interpretation data of the target layer, the low-frequency velocity body is obtained through the Kriging interpolation method. The velocity curves of the pre-selected wells in the area are low-pass filtered, and the first cutoff frequency value FL of the low-frequency band is determined according to the spectrum analysis of the 3D seismic data body, that is, the frequency components with frequency components greater than the first cutoff frequency value FL are filtered out, and only the velocity values ​​within the cutoff frequency value are retained to obtain the filtered low-frequency velocity curve, and according to the time-depth relationship of each pre-selected well, the filtered low-frequency velocity curve is converted from the depth domain to the time domain.

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

[0119] 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 interpolated 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 layer velocity body Vel_Low of the target layer.

[0120] In this example, the velocity logging curves (denoted as Vp) of 25 pre-selected wells are low-pass filtered. Through the spectrum analysis results of the three-dimensional seismic data body of the target layer segment, it can be obtained that the frequency band range of the target layer is 15Hz-45Hz, and the main frequency is 30Hz, that is, the three-dimensional seismic data body lacks frequency components less than 15Hz. Therefore, the low-frequency components on the well can be used to supplement. The first cutoff frequency value FL=15Hz is determined, that is, the frequency components less than 15Hz are extracted from the well. By low-pass filtering the velocity curve, a low-frequency velocity curve within 15Hz is obtained, and according to the time-depth relationship of the 25 pre-selected wells obtained above, the low-frequency velocity curve in the depth domain is converted to the time domain. Then, according to the T seismic interpretation horizons (T=5) of the target layer in the area, the subdivision ratio C=25 is selected, and 25 sub-layers are interpolated between each two adjacent interpretation horizons from bottom to top, and a layer grid model Layer_Model consisting of 5 interpretation horizons and C interpolated sub-layers between each two adjacent interpretation horizons is obtained. Based on the low-frequency velocity curve in the time domain 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 turn, and then the Kriging interpolation method is used to calculate the low-frequency velocity value at each CDP along the layer in turn 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 layer velocity body Vel_Low of the target layer is obtained. Figure 4 As shown, the example profile reflects the overall change trend of the velocity of the target layer and embodies the background velocity characteristics of the high-speed anomaly body of the limestone in the target layer.

[0121] 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 pre-selected well, the impedance curve of the pre-selected well is constructed accordingly, and it is used as the input data of the spectral simulation inversion. In order to facilitate the conversion of the impedance data body obtained by inversion to the velocity data body, when constructing the impedance curve, the density 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), and the velocity curve value of the pre-selected well is multiplied by the constant density value Den_average to obtain the impedance curve of the pre-selected well, and according to the obtained time-depth relationship of each pre-selected well, the impedance curve of the pre-selected well is converted from the depth domain to the time domain to obtain the time domain impedance curve IMP_time.

[0122] Next, the time domain impedance curves of n pre-selected wells and the three-dimensional seismic data volume (i.e., amplitude-preserving migration stacking data) are used for spectrum simulation inversion to obtain the mid-frequency velocity volume Vel_Medium of the target layer. By spectrum analysis of the well logging data and the three-dimensional seismic data volume, the frequency matching operator that makes the seismic spectrum and the well logging wave impedance spectrum achieve the best match is obtained. At the n pre-selected wells, n wellside seismic traces are correspondingly extracted, 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), 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 earthquake amplitude spectrum S revealed by the above equation (6) is i (ω) is equal to the amplitude spectrum R of the reflection coefficient 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 easier 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 obtain 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 correspondence 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 equation (7), we can obtain the frequency matching operator B i (ω).

[0123] 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 inverse 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, that is, 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) can obtain the inversion result of the reflection coefficient of the seismic trace, and then the reflection coefficient is converted to the wave impedance value through the general correspondence between the reflection coefficient and the wave impedance to obtain the wave impedance inversion result of the seismic trace. The wave impedance obtained is the relative wave impedance. In order to more intuitively analyze the inversion effect of the spectrum simulation wave impedance and facilitate the inversion quality control by wells, the relative wave impedance inversion result of the mid-frequency band can be added with the low-frequency background, that is, the relative impedance inversion value obtained by the spectrum simulation inversion is added with the low-frequency background value on the well logging, that is, according to the well logging analysis results, a low-frequency impedance value is assigned to each layer to obtain the final absolute wave impedance value of the spectrum simulation.

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

[0125] In this example, based on the density curves of the 25 pre-selected wells, the density average Den_average of the target layer is calculated as 2.5, and the velocity curve values ​​of the 25 pre-selected wells are multiplied by the constant density value Den_average, then the impedance curve of the pre-selected wells can be obtained, and the impedance curve of the pre-selected wells can be converted from the depth domain to the time domain through the obtained time-depth relationship of each pre-selected well. Then, according to the universal correspondence between reflection coefficient and wave impedance, the logging wave impedance amplitude spectrum of the 25 pre-selected wells is converted into the 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 side channel of 25 pre-selected wells i (ω). According to formula (7), the frequency matching operator B of spectrum simulation inversion is obtained: i (ω). Next, take the inverse of the frequency matching operator and transform it to the time domain through 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 iThe inversion result of the reflection coefficient can be obtained by convolution of (t), and then the corresponding relationship between the reflection coefficient and the wave impedance is used to convert the reflection coefficient to the wave impedance value to obtain the relative wave impedance inversion result. According to the statistical results of the 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, to obtain the absolute wave impedance body of the final spectrum simulation, and then the wave impedance inversion body is divided by the density mean 2.5 to obtain the mid-frequency band velocity body Vel_Medium within the mid-frequency band of the target layer. Figure 5 As shown, the example profile reflects the velocity variation characteristics of the target layer in the same frequency band as the earthquake.

[0126] The high-frequency velocity body of the target layer is obtained by 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 number of effective samples increases, the seismic waveform represents the reservoir better. However, when the number of effective samples increases to a certain extent, the degree of representation of the reservoir by the seismic waveform 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 according to the amplitude spectrum of the three-dimensional seismic data body. 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 tests, 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 each sampling point data of the target layer, thereby obtaining the wave impedance body of the waveform indication inversion, and dividing it by the density average Den_average, the layer velocity body 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 result of the waveform indication inversion contains certain low-frequency and medium-frequency components. In order to remove the low-frequency and medium-frequency components and only retain 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, so as to obtain 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.

[0127] In this example, the inversion parameters are tested based on the time domain impedance curves IMP_time of 25 pre-selected wells, the three-dimensional seismic data body, and the layer grid model Layer_Model created by 5 interpreted layers. According to the principle of minimum inversion error, the key parameters of inversion are determined as follows: the smoothing coefficient is 2, the number of effective samples is 5, the high-pass and high-cut frequencies of the low-frequency band are 5Hz and 10Hz respectively, and the low-cut, low-pass, high-pass and high-cut frequencies of the high-frequency band are 40Hz, 50Hz, 90Hz and 100Hz respectively. Then, the waveform indication inversion is performed using the above inversion parameters to obtain the wave impedance data body of the waveform indication inversion, which is then divided by the density mean of 2.5 to obtain the velocity body. The velocity body is high-pass filtered, that is, the frequency components below 60Hz (FH=60Hz) are filtered out, and only the high-frequency part is retained to obtain the high-frequency band velocity body Vel_High, as shown in Figure 6 As shown, the example profile reflects the high-frequency variation characteristics of the velocity of the target layer and can fully identify the longitudinal and lateral velocity variation characteristics of some thin layers.

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

[0129] Among them, in the process of inversion of the mid-frequency band of the target layer (through the spectrum simulation inversion method, the three-dimensional seismic data body is inverted to obtain the mid-frequency band velocity body), in order to more intuitively analyze the vertical and horizontal variation characteristics of the reservoir and perform effective inversion quality control, a certain amount of low-frequency information is added, so the low-frequency information in the mid-frequency band velocity body Vel_Medium can be removed. That is, the mid-frequency band velocity body Vel_Medium is high-pass filtered (the parameters of the high-pass filter can be determined based on the spectrum of the three-dimensional seismic data body, for example, the high-pass frequency is 15Hz), and through high-pass filtering, a (layer) velocity body Vel_Mediumfilter that only retains the mid-frequency band information is obtained.

[0130] 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. Among them, the preset weight values ​​can be determined by analyzing the consistency between the inverted layer velocity data and the well, for example, it can be a value between 1-2, the weight of the large consistency is high, and the weight of the small consistency is low. For example, by analyzing the consistency between the inversion result and the well, the weight value W_Low of the low-frequency velocity body, the weight value W_Medium of the medium-frequency velocity body and the weight value W_High of the high-frequency velocity body are determined. For any sampling point i of the target layer segment, 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 turn. Through the above formula (12), the layer velocity value Vel_Model(i) after the data of the point is fused can be obtained.

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

[0132] In this example, the obtained mid-frequency velocity body Vel_Medium is high-pass filtered, that is, the frequency components below 15Hz are filtered out, and only the frequency components above 15Hz in the Vel_Medium layer velocity body are retained. The mid-frequency velocity body Vel_Mediumfilter after high-pass filtering 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 As shown in the figure, the example profile can not only clearly reflect the longitudinal and lateral variation characteristics of the target layer velocity, but also can finely describe the velocity variation characteristics inside and at the boundary of the limestone high-speed anomaly body, and can also finely reveal the velocity variation of some thin layers.

[0133] 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. Figure 7 and Figure 8As shown in the comparison, 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 lateral changes in the velocity of the target layer are not precisely represented; 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 lateral directions, reflect the geological laws of the study area, and also can 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.

[0134] Next, according to 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 errors 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 are calculated respectively, as shown in Table 1.

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

[0136]

[0137] Table 1 shows the error values ​​of the drilling depth above sea level and the time-depth conversion depth corresponding to the KT1 geological layer of 14 wells in the area. Through the comparative analysis of the time-depth conversion errors of the two different velocity models at each well point, 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 is 0.71%, which is significantly lower than the depth error value of 2.86% after the (layer) velocity model Vel_tradition conversion. This further verifies the reliability of the velocity model established by the embodiment of the present invention based on the frequency-divided band inversion.

[0138] Accordingly, the method for constructing a velocity model in the embodiment of the present invention effectively realizes the fine establishment of a velocity model of a special geological body development stratum, improves the velocity model's representation of the vertical and horizontal variation characteristics of the reservoir velocity, and the velocity model can not only reflect geological laws, but also well characterize the velocity variation characteristics inside and at the boundary of a special geological body (in this example, a limestone high-speed anomaly body), providing reliable velocity data for comprehensive geological studies such as time-depth conversion, fine structural implementation, and well location deployment. The embodiment of the present invention is easy to implement and has high credibility, and has positive significance for the expansion and efficient advancement of exploration in related fields of special geological bodies.

[0139] 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.

[0140] An embodiment of the present invention further provides a machine-readable storage medium, on which instructions are stored, and the instructions enable a machine to execute the above-mentioned method for constructing a velocity model.

[0141] Those skilled in the art will appreciate that the embodiments of the present application may be provided as methods, systems, or computer program products. Therefore, the present application may 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 may adopt the form of a computer program product implemented in one or more computer-usable storage media (including but not limited to disk storage, CD-ROM, optical storage, etc.) that include computer-usable program code.

[0142] The present application is described with reference to the flowcharts and / or block diagrams of the methods, devices (systems), and computer program products according to the embodiments of the present application. It should be understood that each process and / or box in the flowchart and / or block diagram, as well as the combination of the processes and / or boxes in the flowchart and / or block diagram, can be implemented 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 processor, or other programmable data processing device to generate a machine, so that the instructions executed by the processor of the computer or other programmable data processing device generate instructions for implementing the processes in the flowchart and / or block diagram. Figure 1 A process or multiple processes and / or boxes Figure 1 A device that provides the functions specified in a block or multiple blocks.

[0143] These computer program instructions may also be stored in a computer-readable memory capable of directing 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 an article of manufacture comprising an instruction device, which implements the process Figure 1 A process or multiple processes and / or boxes Figure 1 A function specified in one or more boxes.

[0144] 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, thereby providing instructions for implementing the process. Figure 1 A process or multiple processes and / or boxes Figure 1The steps for the functions specified in one or more boxes.

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

[0146] 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.

[0147] Computer readable media include permanent and non-permanent, removable and non-removable media that can be implemented by any method or technology to store information. 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 disk read-only memory (CD-ROM), digital versatile disk (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 temporary computer readable media (transitory media), such as modulated data signals and carrier waves.

[0148] It should also be noted that the terms "include", "comprises" or any other variations thereof are intended to cover non-exclusive inclusion, so that a process, method, commodity or device including a series of elements includes not only those elements, but also other elements not explicitly listed, or also includes elements inherent to such process, method, commodity or device. In the absence of more restrictions, the elements defined by the sentence "comprises a ..." do not exclude the existence of other identical elements in the process, method, commodity or device including the elements.

[0149] The above are only 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 modification, equivalent replacement, improvement, etc. made within the spirit and principle of the present application should 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