CN111722283B - Stratum velocity model building method - Google Patents
Stratum velocity model building method Download PDFInfo
- Publication number
- CN111722283B CN111722283B CN202010594866.5A CN202010594866A CN111722283B CN 111722283 B CN111722283 B CN 111722283B CN 202010594866 A CN202010594866 A CN 202010594866A CN 111722283 B CN111722283 B CN 111722283B
- Authority
- CN
- China
- Prior art keywords
- velocity
- data
- layer
- well
- calculation
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 238000000034 method Methods 0.000 title claims abstract description 46
- 238000004364 calculation method Methods 0.000 claims description 136
- 230000015572 biosynthetic process Effects 0.000 claims description 20
- 230000008859 change Effects 0.000 claims description 13
- 238000013178 mathematical model Methods 0.000 claims description 10
- 238000012360 testing method Methods 0.000 claims description 10
- 238000010606 normalization Methods 0.000 claims description 9
- 238000009499 grossing Methods 0.000 claims description 7
- 230000008569 process Effects 0.000 claims description 7
- 238000010586 diagram Methods 0.000 claims description 6
- 238000013528 artificial neural network Methods 0.000 claims description 5
- 238000000605 extraction Methods 0.000 claims description 5
- 238000012417 linear regression Methods 0.000 claims description 4
- 238000005457 optimization Methods 0.000 claims description 3
- 238000013517 stratification Methods 0.000 claims description 3
- 238000005553 drilling Methods 0.000 abstract description 9
- 230000008901 benefit Effects 0.000 abstract description 6
- 210000001015 abdomen Anatomy 0.000 abstract description 3
- 238000012545 processing Methods 0.000 description 14
- 238000012937 correction Methods 0.000 description 13
- 238000005755 formation reaction Methods 0.000 description 12
- 239000007789 gas Substances 0.000 description 11
- 210000002569 neuron Anatomy 0.000 description 11
- 238000004458 analytical method Methods 0.000 description 8
- 238000013508 migration Methods 0.000 description 8
- 230000005012 migration Effects 0.000 description 8
- 238000005516 engineering process Methods 0.000 description 7
- 230000006870 function Effects 0.000 description 6
- 238000011160 research Methods 0.000 description 6
- 238000007796 conventional method Methods 0.000 description 3
- 238000003384 imaging method Methods 0.000 description 3
- 238000003062 neural network model Methods 0.000 description 3
- 238000005070 sampling Methods 0.000 description 3
- 238000012952 Resampling Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 230000009918 complex formation Effects 0.000 description 2
- 230000007547 defect Effects 0.000 description 2
- 238000011161 development Methods 0.000 description 2
- VNWKTOKETHGBQD-UHFFFAOYSA-N methane Chemical compound C VNWKTOKETHGBQD-UHFFFAOYSA-N 0.000 description 2
- 238000001228 spectrum Methods 0.000 description 2
- 238000007476 Maximum Likelihood Methods 0.000 description 1
- 238000010521 absorption reaction Methods 0.000 description 1
- 230000001427 coherent effect Effects 0.000 description 1
- 239000000571 coke Substances 0.000 description 1
- 238000010835 comparative analysis Methods 0.000 description 1
- 238000010276 construction Methods 0.000 description 1
- 238000007598 dipping method Methods 0.000 description 1
- 230000000694 effects Effects 0.000 description 1
- 230000004927 fusion Effects 0.000 description 1
- 238000012804 iterative process Methods 0.000 description 1
- 239000003345 natural gas Substances 0.000 description 1
- 239000003209 petroleum derivative Substances 0.000 description 1
- 238000000611 regression analysis Methods 0.000 description 1
- 239000011435 rock Substances 0.000 description 1
- 238000003786 synthesis reaction Methods 0.000 description 1
- 238000012549 training Methods 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/303—Analysis for determining velocity profiles or travel times
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V1/00—Seismology; Seismic or acoustic prospecting or detecting
- G01V1/28—Processing seismic data, e.g. for interpretation or for event detection
- G01V1/30—Analysis
- G01V1/301—Analysis for determining seismic cross-sections or geostructures
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/51—Migration
- G01V2210/512—Pre-stack
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/50—Corrections or adjustments related to wave propagation
- G01V2210/51—Migration
- G01V2210/514—Post-stack
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/62—Physical property of subsurface
- G01V2210/622—Velocity, density or impedance
- G01V2210/6222—Velocity; travel time
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/63—Seismic attributes, e.g. amplitude, polarity, instant phase
-
- G—PHYSICS
- G01—MEASURING; TESTING
- G01V—GEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
- G01V2210/00—Details of seismic processing or analysis
- G01V2210/60—Analysis
- G01V2210/66—Subsurface modeling
Landscapes
- Engineering & Computer Science (AREA)
- Remote Sensing (AREA)
- Physics & Mathematics (AREA)
- Life Sciences & Earth Sciences (AREA)
- Acoustics & Sound (AREA)
- Environmental & Geological Engineering (AREA)
- Geology (AREA)
- General Life Sciences & Earth Sciences (AREA)
- General Physics & Mathematics (AREA)
- Geophysics (AREA)
- Geophysics And Detection Of Objects (AREA)
Abstract
The invention discloses a stratum velocity model building method. The method can track the track of a complicated geological structure area with large surface relief, complicated ground belly structure, quite developed two wing side reverse masked fracture zones and the like and a page gas horizontal well with higher requirements on the buried depth and the structural form of a target layer. The method can realize the establishment of an accurate stratum velocity model, further reduce the drilling risk and improve the economic benefit of related oil and gas exploration.
Description
Technical Field
The invention belongs to the field of petroleum and natural gas exploration, and particularly relates to a formation velocity model establishing method.
Background
In recent years, exploration and development of marine shale gas in the Sichuan basin are greatly advanced. For example, in Changning-Weiyuan blocks, Weirong shale gas exploration areas, coke dams and flat bridge shale gas exploration blocks of the Sichuan basin, a large number of shale gas horizontal wells are fractured to obtain industrial airflow, thereby obtaining considerable economic benefit. In the exploration and development practice of the marine shale gas, a large number of exploration results show that the good storage condition is one of the key factors of high-yield enrichment of the shale gas. In a shale gas exploration area near a large fracture zone, due to the fact that the structure is complex, seismic reflection imaging is difficult, and related shale gas well laying and exploration progress is influenced.
A large number of exploration practices show that in complex geological structure areas with large surface relief, complex ground belly structure, quite developed reverse masked fracture zones on two wing sides and the like, a depth structure diagram obtained by a conventional technical process is used for guiding drilling, so that an oil and gas reservoir is easily missed, and huge economic loss is generated.
Under highly steep complex formation conditions, only prestack depth migration can be used to obtain correct imaging of the subsurface medium from the artificial seismic data. Whereas prestack depth migration can only image correctly if the layer velocity model is accurate. Where pre-stack depth migration needs to be used, conventional velocity analysis methods often do not meet the desired accuracy requirements. In this case, to find a more accurate formation velocity model, a velocity analysis method based on migration iteration is developed. It can now be said that the migration velocity analysis is velocity analysis using migration iterations, and such velocity analysis methods use the effect of the velocity field on the migration imaging to correct the formation velocity model.
Common drawbacks of the above-described correlation methods are assumptions about subsurface and geological conditions in the velocity correction relationship, such as that the subsurface velocity is constant, and that it varies only in the depth direction and small offset assumptions. Under these assumptions, the average velocity may be approximated as a root mean square velocity, which may then be converted to a layer velocity by the Dix (Dix) formula or other formula. However, there is a large difference between the average velocity and the root mean square velocity when the formation layer velocity varies laterally. This causes divergence in the iterative process of velocity and therefore affects the accuracy and stability of the method in the presence of many complex formations and lateral variations in velocity, which basically assumes that its application is limited to the case where the formation is a horizontal layer or a subsurface formation is relatively simple, which limits the application of the method to some extent. Common disadvantages of these methods are the lack of suitable lateral upsets in velocity and dipping formations.
In the construction explanation, the information of well drilling, regional geology and the like is mainly utilized to establish the speed of each layer which changes transversely along the layer, and under the condition that the well pattern density is large enough and the distribution is uniform, a reasonable and accurate layer speed model can be established by utilizing a variable speed modeling technology.
In the non-well control area of the work area, the conventional method is to estimate the speed of each layer according to the geological data of the area or the like, or to estimate the speed by using the time-depth relation of the well drilling of the work area. And establishing a layer speed model of the whole area based on the calculated layer speed control points of the well control area and the estimated layer speed control points of the non-well control area. However, the process from point to surface to three-dimensional body is a process, the realization of the interpolation, smoothing and other processes has great limitation, and local distortion of the velocity body is easy to occur under the condition of improper operation. Therefore, the number and size of the wells also affect the formation of the interval velocity model.
As can be seen from the analysis of the conventional modeling method for the prestack depth migration velocity and the modeling method for the structure-explained variable speed velocity, the two methods have respective advantages and disadvantages and technical difficulties. For the speed model establishment based on the offset iterative speed analysis, the method has the greatest advantages of good global property, convenient and quick processing of various speed bodies, and difficult representation of complex speed structures; for the layer velocity model building based on the variable speed modeling technology, the method has the greatest advantages that an accurate stratum velocity control point of a well control area can be built, a more complex velocity structure can be represented, and the method has the greatest defects that the globality is lacked, and the interpolation and the calculation of a velocity body are difficult to be restricted.
Disclosure of Invention
Aiming at the defects in the prior art, the method for establishing the stratum velocity model solves the problems of inaccurate establishment of the stratum velocity model and high drilling risk.
In order to achieve the purpose of the invention, the invention adopts the technical scheme that: a method for establishing a stratum velocity model comprises the following steps:
s1, obtaining related preferred attribute data volume and horizon data through geological data, well logging data and seismic data, and establishing a horizon velocity calculation model through the horizon data;
s2, calculating a first-layer speed data body through the layer speed calculation model and the related preferred attribute data body, establishing a calculation grid through the first-layer speed data body, extracting grid points, and calculating speed change gradient values on the grid points;
s3, calculating the interval velocity on the grid point by using the measured interval velocity in each calculation well and the velocity change gradient value on the grid point, and determining the optimal weighting factor of each calculation well;
s4, carrying out data weighting reconstruction, interpolation and smoothing on the stratum velocity and the optimal weighting factor on the grid points of each calculation well to obtain an optimized stratum velocity data plane graph, and establishing a stratum velocity model through the optimized stratum velocity data plane graph.
Further: the geological data in the step S1 comprises core logging data, geological stratification data and petrophysical test results; the logging data comprises acoustic and density logging curves and transverse wave data; the seismic data are conventional three-dimensional pre-stack gather or post-stack seismic data volume, and the post-stack seismic data are obtained by stacking and shifting the pre-stack gather data.
Further: the specific steps of step S1 are:
s11, performing inversion and attribute extraction on the conventional three-dimensional pre-stack gather or post-stack seismic data volume aiming at the layer velocity to obtain a seismic inversion data volume and an attribute data volume, extracting seismic inversion data values and attribute data values of a target layer at each well point from the seismic inversion data volume and the attribute data volume, optimizing and normalizing the attribute data values, and determining a related optimal attribute data volume;
s12, performing relevant well-seismic calibration by using the post-stack seismic data volume, the logging data and the geological data, determining relevant meaningful reflecting layer interfaces on the well, setting an interpretation grid, and then tracking and interpreting the reflecting layer interfaces to obtain the horizon data of the relevant reflecting layer interfaces.
Further: the calculation formula of the normalization processing in step S11 is:
in the above formula, XpiTo normalize the processed sample values, XpTo normalize the pre-processed sample values, Xmax＝max{Xp}，Xmin＝min{XpN and m are positive integers, m>n≥0。
Further: the layer velocity calculation model in the step S2 includes a multiple linear regression model and a BP neural network regression mathematical model, related sample wells and blind wells are set in the layer velocity calculation model, the sample wells train the related layer velocity calculation model, and parameters of the related layer velocity calculation model are determined; the optimization of the layer velocity calculation model is to calculate the layer velocity of the blind well point by using the attribute data value of the blind well and the test layer velocity calculation model, calculate the correlation coefficient between the calculation result and the actually measured layer velocity of the blind well, and select the layer velocity calculation model with the highest correlation coefficient.
Further: the calculation formula of the correlation coefficient is as follows:
in the above formula, r is a correlation coefficient, XiAnd YiThe interval velocity of the blind well point and the measured interval velocity of the blind well are respectively, i is 1,2 …, N is the total number of data,
and
and respectively ranking the average values of the interval velocity of the blind well point and the measured interval velocity of the blind well.
Further: the calculation formula of the velocity change gradient value at the grid point in step S2 is:
in the above formula, PiFor the gradient value of the velocity change at the i-th grid point, Δ ViCalculating a first interval velocity difference, Δ I, for the ith grid pointiCalculating the linear distance of the ith grid point and the calculated well point;
wherein, is Δ ViThe calculation formula of (2) is as follows:
ΔVi＝Vi-Vo
in the above formula, ViFirst layer velocity data value, V, for the ith grid pointoA first layer data value for the well is calculated.
Further: the specific steps of step S4 are: and performing weighted reconstruction calculation on second-layer velocity data on each grid point by using the layer velocity data values and the optimal weighting factors on each grid point of the related target layer of each calculation well, interpolating and smoothing the second-layer velocity data on each grid point to obtain an optimized layer velocity data plane diagram, and establishing a stratum velocity model through the optimized layer velocity data plane diagram.
The invention has the beneficial effects that: the method can track the track of a complicated geological structure area with large surface relief, complicated ground belly structure, quite developed two wing side reverse masked fracture zones and the like and a page gas horizontal well with higher requirements on the buried depth and the structural form of a target layer. The method can realize the establishment of an accurate stratum velocity model, further reduce the drilling risk and improve the economic benefit of related oil and gas exploration.
Drawings
FIG. 1 is a flow chart of the present invention.
Detailed Description
The following description of the embodiments of the present invention is provided to facilitate the understanding of the present invention by those skilled in the art, but it should be understood that the present invention is not limited to the scope of the embodiments, and it will be apparent to those skilled in the art that various changes may be made without departing from the spirit and scope of the invention as defined and defined in the appended claims, and all matters produced by the invention using the inventive concept are protected.
As shown in fig. 1, a method for establishing a formation velocity model includes the following steps:
s1, obtaining related preferred attribute data volume and horizon data through geological data, well logging data and seismic data, and establishing a horizon velocity calculation model through the horizon data;
the geological data comprises core logging data, geological stratification data and rock physical test results; the logging data comprises acoustic and density logging curves and transverse wave data; the seismic data are conventional three-dimensional pre-stack gather or post-stack seismic data volume, and the post-stack seismic data are obtained by stacking and shifting the pre-stack gather data.
Performing inversion and attribute extraction on the conventional three-dimensional pre-stack gather or post-stack seismic data body aiming at the layer velocity to obtain a seismic inversion and attribute data body; and extracting seismic inversion and attribute data values of a target layer at each well point from the seismic inversion and attribute data volumes, calculating correlation coefficients between the seismic inversion and attribute data values, determining P inversion and attribute data volumes, performing normalization processing, and taking the P inversion and attribute data volumes as an optimal attribute data volume to enter the next step. In principle, the calculation of these inversion and attribute is related to the calculation of the layer velocity, and generally, the related inversion and attribute data volumes of longitudinal wave velocity and transverse wave velocity, poisson ratio, wave impedance and the like are calculated; and optimizing the related attribute data, and determining the related optimized attribute data to enter the next step. In addition, the preferred attribute data volume is normalized so that the data is included in a certain value range after calculation. The seismic inversion and attribute extraction can be realized by relevant geophysical exploration commercial software, such as a PAL module of Landmark company, amplitude, frequency and instantaneous seismic attributes can be extracted from three-dimensional seismic post-stack data, such as jason software, prestack or poststack wave impedance inversion data, longitudinal waves/transverse waves, density data volumes and the like can be calculated, VVA software can extract data of curvature, coherent bodies, frequency dividing bodies, maximum likelihood bodies and the like, FRS software can utilize three-dimensional prestack gather data to calculate the anisotropic strength of P waves, extract the anisotropic strength data of P waves and various absorption attenuation attributes and the like. The correlation attribute may be obtained by performing weighted fusion calculation on a plurality of specified attribute data volumes, or by performing some mathematical calculation on the attribute data volumes. The computation of these inversion and attributes is related to the computation of the layer velocity, and typically an inversion and attribute data volume related to compressional and shear velocities, density, and fracture is computed. In principle, the number of calculations of the relevant preference attribute should be greater than or equal to three. In actual operation, a plurality of seismic inversion and attribute data volumes are obtained by performing pre-stack or post-stack inversion and attribute extraction on a conventional three-dimensional pre-stack gather or post-stack seismic data volume; performing correlation coefficient calculation on seismic inversion and attribute data values on data tracks of a target layer at each well point and corresponding well point upper layer velocity data values, and selecting M inversion and attribute data with higher correlation coefficients obtained in calculation; and P inversion attribute data with lower correlation coefficient are selected from the M inversion attribute data and the M attribute data to enter the next step, and the P inversion attribute data are used as the preferred attribute data. The method comprises the steps of selecting M inversion and attribute values with higher correlation coefficients, wherein the M inversion and attribute values refer to M inversion and attribute data volumes with correlation coefficients larger than 0.65, and selecting P inversion and attribute data volumes with lower phase relation numbers from the M inversion and attribute values, wherein the attribute values refer to P inversion and attribute data volumes with correlation coefficients smaller than 0.45.
In addition, extracting the seismic inversion and attribute data values of the target interval at each well point from the seismic inversion and attribute data body mainly means that the well-seismic synthesis record is used for calibrating results, the time-depth relation of the well points is determined, the depth range data value of the target interval at the well point, the seismic data interpretation results and the like are used for determining the double-pass reflection time period corresponding to the target interval in the well, and then the relevant well point coordinates and the double-pass reflection time period of the target interval are projected into the relevant attribute data body, so that the seismic inversion and attribute data values of the target interval can be obtained. Secondly, the time-depth relation can be utilized to convert the attribute curve data of the related attribute data body of the well from the time domain to the depth domain, and the resampling calculation is carried out on the attribute curve data of the depth domain, so that the sampling interval of the attribute curve data of the depth domain is consistent with the sampling interval of the interval velocity curve of the well, and the one-to-one corresponding relation between attribute-interval velocity data values can be formed. Conversely, the corresponding relationship between the attribute-interval velocity data at the well point can be formed in the time domain.
In addition, the optimized seismic inversion and the attribute data are respectively subjected to normalization processing to obtain an optimized attribute data volume, and respective function formulas of the normalization processing are obtained. Normalization is a dimensionless processing means that transforms the absolute value of a physical system value into a relative value relationship, i.e., it uses addition, subtraction, multiplication, division or a combination thereof to perform an operation. The calculation formula of the normalization process is as follows:
in the above formula, XpiTo normalize the processed sample values, XpTo normalize the pre-processed sample values, P is 1,2, …, P, Xmax＝max{Xp}，Xmin＝min{XpN and m are positive integers, m>n is more than or equal to 0. Generally, the interval value range of n to m is determined by taking appropriate speed increment in small and large directions in a targeted manner after counting various interval speeds in a well. The specific situation can be determined according to expert experience, calculation precision, statistical data of stratum velocities of related research areas and the like, and generally the velocity increment can be set to be 1000 m/s.
S2, calculating a first-layer speed data body through the layer speed calculation model and the related preferred attribute data body, establishing a calculation grid through the first-layer speed data body, extracting grid points, and calculating speed change gradient values on the grid points;
the layer velocity calculation model comprises a multiple linear regression model and a BP neural network regression mathematical model, related sample wells and blind wells are set in the layer velocity calculation model, the sample wells train the related layer velocity calculation model, and parameters of the related layer velocity calculation model are determined; the optimization of the layer velocity calculation model is to calculate the layer velocity of the blind well point by using the attribute data value of the blind well and the test layer velocity calculation model, calculate the correlation coefficient between the calculation result and the actually measured layer velocity of the blind well, and select the layer velocity calculation model with the highest correlation coefficient.
The related optimal layer velocity calculation model has various technical methods and processes, and relatively accurate results can be obtained possibly, and the related layer velocity calculation model provided by the technology of the invention has the following methods and processes:
(1) a mathematical model method. The technical method is mainly used for establishing a related mathematical model to calculate the interval velocity, and is suitable for the condition that the number of wells in a research area is relatively large. The related calculation method and the flow are as follows:
the technical method mainly comprises the steps of establishing a plurality of layer velocity calculation models, optimizing the related models, and selecting an optimal model to participate in the calculation of the subsequent steps. The method is characterized in that all wells are respectively set into a sample well and a blind well, a related interval velocity calculation model is established by utilizing each interval velocity data value of an upper interval velocity curve of the sample well and a data value on an optimal attribute data volume curve of a corresponding depth domain, such as establishing a multivariate linear regression, a BP neural network regression and related improved mathematical models thereof to calculate interval velocity values, and an algorithm program based on seismic inversion and attribute data prediction interval velocity values can be developed. The algorithm of the correlation calculation model is briefly described as follows:
(a) and (4) performing multiple regression analysis. And establishing a multiple high degree polynomial regression model by using the data on the stratum velocity curve of the sample well and the corresponding preferred attribute data on the well point. The calculation formula is as follows:
in the above formula, yiFor predicting the ith interval velocity value, x, of interval velocity curve in a blind welliiThe attribute data of the ith preferred attribute data curve corresponding to the ith interval velocity value on the well point of the well,i≤p；aij(i-0, 1, …, p; j-1, 2, …, m, m is the number of samples) is the regression coefficient.
According to the data value y on the speed curve of the measured layer on the wellinAnd the predicted value-y on the upper layer velocity curve of the well pointiThe sum of the squares of the residuals between the values is minimized, and each coefficient a is obtained by the least square methodijThe value of (c). The calculation formula of the sum of the squares of the residuals is as follows:
Q＝∑(yin-yi)2
in the above formula, yinIs the measured i-th layer velocity value y on a certain blind welliAnd Q is the sum of the squares of the residuals, which is the value on the predicted layer velocity curve corresponding to the measured ith layer velocity value.
(b) BP neural networks and related improvements. The algorithm mainly utilizes back propagation learning to establish a neural network model for layer velocity prediction, takes actually measured layer velocity curve data as learning training and testing samples, takes curve data on a relevant optimized attribute data body on a well point as a learning sample, and trains the network. If the learning sample is set to be (x)1i,x2i,…,xpi；tp) (P ═ 1,2, …, P; p is the number of samples). Randomly giving W (W)ij，θi，vi) Then, the output y of the p sample of the network is calculated according to the formulas (9) to (11)p。
In the above formula, n is the number of neurons in the input layer; m is the number of neurons in the hidden layer; wijThe connecting weight of the neuron i of the hidden layer and the neuron j of the input layer; thetaiThe threshold of neuron i in the hidden layer.
In the above formula, IiAn input to a neuron of the ith hidden layer; oiIs the god of the ith hidden layerAnd outputting the channel element.
In the above formula, viThe connection weight of the output layer neuron and the hidden layer neuron i is obtained; y ispIs the output of the p-th sample.
Defining the connection weights W from hidden layer neurons to input layer neuronsijHidden layer neuron threshold value thetaiAnd the connection weight v of the output layer neuron and the hidden layer neuroniThe composed vector is the connection weight vector W of the network.
For sample p, the output error of the network is defined as:
and defines an error function as:
along the error function epThe direction of the negative gradient as a function of W corrects for W. Let the correction value of W be DeltaW, take
In the above formula, η is the learning rate and is a number between 0 and 1.
After calculating DeltaW, adopting an iterative formula
W+ΔW→W
And (5) performing correction calculation on the original W to obtain a new connection weight vector W.
For all the studies, the above-described calculation process is performed in the order of sample arrangement, and then the value of W is fixed. And respectively carrying out forward calculation on the P samples so as to obtain the energy function value of the learning sample.
And through repeated iteration, correcting the network connection weight W to ensure that E meets a certain precision requirement.
In the above formula, E is a standard estimation error, and the smaller the value, the better the model is built, and R is a decision coefficient, and the larger the value, the better the model is built.
And respectively calculating a layer speed curve predicted by a certain blind well section by using each obtained mathematical model, and carrying out correlation coefficient calculation with the actually measured layer speed curve, and preferably selecting the mathematical model corresponding to the predicted layer speed curve with the highest correlation coefficient as the optimal layer speed calculation model. Firstly, extracting a time-depth relation of the blind well, converting curve data obtained by calculating each mathematical model on a blind well point into a depth domain from a time domain by using the time-depth relation, and performing resampling calculation according to a sampling rate of a layer speed curve to obtain predicted layer speed curve data of a processed depth domain; and respectively carrying out correlation coefficient calculation with data values on the actually measured layer speed curve of the target layer section of the blind well, and selecting a mathematical model corresponding to the predicted layer speed curve with the highest correlation coefficient value as an optimal layer speed calculation model. In addition, the difference between the sample well and the blind well is that the sample well participates in the calculation of the velocity calculation model of the relevant stratum, but the blind well does not participate, and the blind well is only used for testing and optimizing the velocity calculation model of the relevant stratum.
The correlation coefficient is calculated by the formula:
in the above formula, r is a correlation coefficient, XiAnd YiThe interval velocity of the blind well point and the measured interval velocity of the blind well are respectively, i is 1,2 …, N is the total dataThe number of the first and second groups is,
and
and respectively ranking the average values of the interval velocity of the blind well point and the measured interval velocity of the blind well.
In addition, a calculation model with the smallest error may be selected as the optimal layer velocity calculation model based on the error analysis result of the calculation model. Specifically, the optimal layer velocity calculation model may be selected and analyzed according to actual conditions and expert experience.
S3, calculating the interval velocity on the grid point by using the measured interval velocity in each calculation well and the velocity change gradient value on the grid point, and determining the optimal weighting factor of each calculation well;
a computational grid is built and gradient values at grid points are extracted for the first layer velocity data. The specific operation is to set the parameters of the calculation grid according to the actual situation, and the grid is called a setting grid. The grid parameters comprise grid intervals, grid number and the like, the size of the grid parameters can be determined according to the grid distribution and precision requirement condition of the layer speed to be predicted, and in principle, the larger the grid interval which is usually set is, the lower the prediction precision is, and some characteristic information can be lost; the smaller the grid spacing, the higher the accuracy of the prediction, and the more detailed the plotted result, but the more the computational effort is relatively increased. Generally, the size of the set mesh is determined according to actual and layer velocity prediction requirements, and the set mesh is usually regular.
Wherein gradient values are calculated for the first layer velocity data at each calculated well pattern point. The specific operation is to extract the calculated first layer speed data value on each grid point and establish the plane distance and layer speed difference value between the related calculation well and each grid point; and gradient values for the calculated well are calculated for the relevant grid points. And repeating the steps to complete the gradient calculation of each grid point on the plane of each calculation well. According to the calculation method, the speed change gradient value of the calculation well is calculated for the correction well. Wherein the gradient value at the i-th grid point for a certain well is calculated as follows:
in the above formula, PiFor the gradient value of the velocity change at the i-th grid point, Δ ViCalculating a first interval velocity difference, Δ I, for the ith grid pointiCalculating the linear distance of the ith grid point and the calculated well point;
wherein, is Δ ViThe calculation formula of (2) is as follows:
ΔVi＝Vi-Vo
in the above formula, ViFirst layer velocity data value, V, for the ith grid pointoA first layer data value for the well is calculated.
In principle, because the optimal weighting factor is being calculated, correction wells are used, and the number of correction wells is greater than the number of calculation wells. Therefore, the number of wells required for the method of the present invention must be two or more.
And calculating the related layer velocity data value on the grid point by using the measured layer velocity data in the well and the gradient value on the grid point. The specific operation is to calculate the relevant layer velocity data value on the grid point based on the actually measured layer velocity data in the calculation well and the gradient value on the grid point, and so on, and finish the calculation of the layer velocity data on each grid point of each calculation well. The correlation calculation formula is as follows:
Vo+Pi*ΔIi＝Vi
in the above formula,. DELTA.IiFor the ith grid point and the calculated straight-line distance, P, of the well pointiFor the gradient value, V, of the ith grid point of the calculation welliSecond layer velocity data value, V, for the ith grid pointoMeasured velocity data values for the well target formation are calculated for the well.
S4, carrying out data weighting reconstruction, interpolation and smoothing on the stratum velocity and the optimal weighting factor on the grid points of each calculation well to obtain an optimized stratum velocity data plane graph, and establishing a stratum velocity model through the optimized stratum velocity data plane graph.
The specific operation is that a related equation set is established for the gradient data on the correction well point by utilizing the actually measured layer velocity data value in each calculation well, and then the related coefficient is solved for the equation set, and the related coefficient is used as the optimal weighting coefficient. The calculation formula of the optimal weighting factor is as follows:
in the above formula, VjThe measured interval velocity data value of the correction well for the jth hole,
the weighting factor value for the nth calculation well of the jth correction well,
the velocity data value calculated over the nth calculation well for the jth correction well. If there are f computation wells, then there are f equations to compute to determine the computation well
The value is obtained.
And carrying out data reconstruction on the layer velocity data on the related calculation well grid points by using respective optimal weighting factors to obtain a grid point upper layer velocity reconstruction data value. The layer velocity data reconstruction calculation formula on the related grid points is as follows:
in the above formula, Y is an optimized layer velocity data value at a certain grid point after weighted data processing; ciCalculating a zone velocity data value for the ith well at the grid point; wherein k isiOptimal weighting factor for the ith calculated wellIn the case of a hybrid vehicle,
and 0 ≦ ki。
And calculating second-layer velocity data on related grid points by using the layer velocity data values and the optimal weighting factors on the grid points of the related target layer of each calculated well, and carrying out interpolation, rounding and other processing on the second-layer velocity data on the grid points to obtain an optimized layer velocity data plane related to the target layer. And the like, completing a velocity plane diagram of each stratum, and establishing a velocity model of each stratum in the research area.
In one embodiment of the invention, the acoustic logging data values in each well are subjected to correlation calculation to obtain a layer velocity curve of each well; carrying out well-seismic calibration by using the logging data of each well, determining the position of a shale section in a seismic data body, and obtaining a relevant time-depth relation table; the method comprises the steps of calculating and extracting wave impedance data volumes, longitudinal wave velocity, transverse wave velocity, gradient data, density attributes and the like by adopting conventional commercial software, namely VVA software, FRS and jason software aiming at attributes of layer velocity calculation, calculating related layer velocities by utilizing a velocity spectrum acquired manually to obtain layer velocity data volumes, calculating six seismic inversion and attribute data volumes, calculating related coefficients and related coefficients among the wells, and selecting three attributes as optimal attribute data volumes. In the example, the layer velocity data volume calculated by the longitudinal wave velocity, the wave impedance and the velocity spectrum is found to be the optimal attribute combination, the correlation between the calculated result and the measured layer velocity data value on the well is good, and the correlation between the attribute data is low. Therefore, in the subsequent calculation, the three attribute data volumes are selected as the optimal attribute data volume combination. And respectively carrying out normalization processing on the attribute data by using respective corresponding normalization function calculation formulas, and normalizing various data and the above-well layer velocity data to a (3160, 5600) value domain. In the explanation of the horizon data, the synthetic record calibration result in the well is utilized to determine the related stratum interfaces, namely 9 total horizon data, and the horizon data is subjected to interpolation, rounding and other processing after the related explanation to obtain the related horizon data. And the related fault is interpreted so as to obtain related fault data, thereby being beneficial to establishing an accurate horizon model.
Designing a related interval velocity calculation model, setting a sample well, a blind well and the like, bringing data on an interval velocity curve of the sample well and curve data of a related attribute data body into each interval velocity calculation model for calculation, selecting the interval velocity calculation model with the calculation result closest to the interval velocity condition of the blind well as an optimal interval velocity calculation model, and calculating a first interval velocity data body. In actual operation, 15 sample well meters are set according to the drilling condition in the research area, and 1 blind well is set. Firstly, selecting a layer velocity calculation model, and mainly selecting a BP neural network model as an optimal layer velocity calculation model of a research area according to the layer velocity calculation experience of the adjacent area, the suggestion of an expert and a related test result. And substituting the three optimal attribute data volumes into a BP neural network model for calculation to obtain a first-layer speed data volume.
The 16 wells in the study area were divided into 8 wells for calculation and calibration. And calculating the interval velocity gradient of the calculation well-grid point by using the calculation well, the first interval velocity data body and the designed calculation grid, and calculating the interval velocity gradient of the calculation well and each correction well respectively to obtain related gradient data. In actual operation, the first layer speed data body and the layer data are used for calculating the layer speed data average value of the target layer, and extracting the layer speed data average values of grid points and related well points, so that the related gradient value is calculated. In actual practice, the designed computational grid is 5 lines X5.
And calculating related layer velocity data values on the grid points by using the actually measured layer velocity data in each calculation well and the velocity change gradient values on the grid points, determining an optimal weighting factor by using the correction well, performing layer velocity data weighting reconstruction processing on each grid point, and performing interpolation, smoothing and other processing to obtain an optimized layer velocity data plan. In actual operation, the relevant measured interval velocity data values in eight calculation wells and the gradient data at each grid point are used for calculation, and the interval velocity data values at the grid points of the well are obtained. And repeating the above steps to complete the calculation of the layer velocity data on the grid points of the target layer section by the eight calculation wells. And secondly, calculating the optimal weighting factors of the eight calculation wells about the target layer-by-layer speed by using the actually measured layer speed data values of the eight correction wells, and determining the optimal weighting factors of the eight calculation wells by adopting a correlation calculation formula. Then, the layer velocity data values at the eight grid points and the optimal weighting factors are used for carrying out weighted reconstruction calculation at the relevant grid points to obtain the layer velocity data values at each grid point, and then interpolation, rounding and the like are carried out on the layer velocity data values to obtain an optimized layer velocity plane. And the rest is done in turn to complete the establishment of the stratum velocity model of the relevant layer section.
The stratum velocity model result obtained by the technology of the invention is compared and analyzed with the seismic section of the stratum velocity model obtained by the related conventional method by the well-passing depth domain of the stratum velocity model and the seismic section of the stratum velocity model obtained by the related conventional method by using the same subsequent processing parameter. Relevant research results show that the depth domain seismic section obtained by the technology is superior to the depth domain seismic section obtained by conventional processing, and the horizontal well track on the seismic section is well matched with relevant data through comparative analysis. And analyzing the track position of the cross section of the subsequent related horizontal well, and verifying that the coincidence rate reaches over 81.5 percent through related data of subsequent drilling, and testing the related drilling to obtain high-yield industrial airflow. Compared with the related results, the results of the invention are superior to the results of the conventional formation velocity model building technology, and the fact that the technology is effective in modeling the formation velocity is also proved.
Claims (8)
1. A method for establishing a stratum velocity model is characterized by comprising the following steps:
s1, obtaining related preferred attribute data volume and horizon data through geological data, well logging data and seismic data, and establishing a horizon velocity calculation model through the horizon data;
s2, calculating a first-layer speed data body through the layer speed calculation model and the related preferred attribute data body, establishing a calculation grid through the first-layer speed data body, extracting grid points, and calculating speed change gradient values on the grid points;
s3, calculating the interval velocity on the grid point by using the measured interval velocity in each calculation well and the velocity change gradient value on the grid point, and determining the optimal weighting factor of each calculation well;
s4, carrying out data weighting reconstruction, interpolation and smoothing on the stratum velocity and the optimal weighting factor on the grid points of each calculation well to obtain an optimized stratum velocity data plane graph, and establishing a stratum velocity model through the optimized stratum velocity data plane graph.
2. The method for establishing a formation velocity model according to claim 1, wherein the geological data in the step S1 includes core logging data, geological stratification data and petrophysical test results; the logging data comprises acoustic and density logging curves and transverse wave data; the seismic data are conventional three-dimensional pre-stack gathers or post-stack seismic data volumes, and the post-stack seismic data volumes are obtained by stacking and shifting the pre-stack gather data.
3. The method for establishing the formation velocity model according to claim 2, wherein the step S1 includes the following steps:
s11, performing inversion and attribute extraction on the conventional three-dimensional pre-stack gather or post-stack seismic data volume aiming at the layer velocity to obtain a seismic inversion data volume and an attribute data volume, extracting seismic inversion data values and attribute data values of a target layer at each well point from the seismic inversion data volume and the attribute data volume, optimizing and normalizing the attribute data values, and determining a related optimal attribute data volume;
s12, performing relevant well-seismic calibration by using the post-stack seismic data volume, the logging data and the geological data, determining relevant meaningful reflecting layer interfaces on the well, setting an interpretation grid, and then tracking and interpreting the reflecting layer interfaces to obtain the horizon data of the relevant reflecting layer interfaces.
4. The method for establishing a formation velocity model according to claim 3, wherein the calculation formula of the normalization process in the step S11 is as follows:
in the above formula, XpiTo normalize the processed sample values, XpTo normalize the pre-processed sample values, Xmax＝max{Xp}，Xmin＝min{XpN and m are positive integers, m>n≥0。
5. The method for establishing a formation velocity model according to claim 1, wherein the layer velocity calculation model in the step S2 includes a multiple linear regression model and a BP neural network regression mathematical model, related sample wells and blind wells are set in the layer velocity calculation model, the sample wells train the related layer velocity calculation model, and parameters of the related layer velocity calculation model are determined; the optimization of the layer velocity calculation model is to calculate the layer velocity of the blind well point by using the attribute data value of the blind well and the test layer velocity calculation model, calculate the correlation coefficient between the calculation result and the actually measured layer velocity of the blind well, and select the layer velocity calculation model with the highest correlation coefficient.
6. The method of establishing a formation velocity model according to claim 5, wherein the correlation coefficient is calculated by the formula:
in the above formula, r is a correlation coefficient, XiAnd YiThe interval velocity of the blind well point and the measured interval velocity of the blind well are respectively, i is 1,2 …, N is the total number of data,
and
and respectively ranking the average values of the interval velocity of the blind well point and the measured interval velocity of the blind well.
7. The method for establishing a formation velocity model according to claim 1, wherein the calculation formula of the velocity change gradient value at the grid point in the step S2 is:
in the above formula, PiFor the gradient value of the velocity change at the i-th grid point, Δ ViCalculating a first interval velocity difference, Δ I, for the ith grid pointiCalculating the linear distance of the ith grid point and the calculated well point;
wherein, is Δ ViThe calculation formula of (2) is as follows:
ΔVi＝Vi-Vo
in the above formula, ViFirst layer velocity data value, V, for the ith grid pointoA first layer data value for the well is calculated.
8. The method for establishing the formation velocity model according to claim 1, wherein the step S4 includes the following steps: and performing weighted reconstruction calculation on second-layer velocity data on each grid point by using the layer velocity data values and the optimal weighting factors on each grid point of the related target layer of each calculation well, interpolating and smoothing the second-layer velocity data on each grid point to obtain an optimized layer velocity data plane diagram, and establishing a stratum velocity model through the optimized layer velocity data plane diagram.
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010594866.5A CN111722283B (en) | 2020-06-28 | 2020-06-28 | Stratum velocity model building method |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202010594866.5A CN111722283B (en) | 2020-06-28 | 2020-06-28 | Stratum velocity model building method |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN111722283A CN111722283A (en) | 2020-09-29 |
| CN111722283B true CN111722283B (en) | 2021-05-25 |
Family
ID=72569039
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN202010594866.5A Active CN111722283B (en) | 2020-06-28 | 2020-06-28 | Stratum velocity model building method |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN111722283B (en) |
Families Citing this family (3)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN112946742B (en) * | 2021-03-17 | 2023-10-27 | 成都捷科思石油天然气技术发展有限公司 | Method for picking up accurate superposition velocity spectrum |
| CN113341461B (en) * | 2021-06-10 | 2023-09-01 | 中国石油大学(北京) | Earthquake velocity prediction method, earthquake velocity prediction device and server |
| CN114526052B (en) * | 2021-12-31 | 2023-09-19 | 中国石油天然气集团有限公司 | Risk prediction method and device for well drilling and completion engineering |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4982382A (en) * | 1989-06-14 | 1991-01-01 | Mobil Oil Corporation | Method of modeling subsurface formations |
| CN102967882A (en) * | 2012-11-16 | 2013-03-13 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | Method for building layer velocity model of stratum |
| CN104502997A (en) * | 2015-01-22 | 2015-04-08 | 中国石油化工集团公司 | Method for using fracture density curve to forecast fracture density body |
| CN105842736A (en) * | 2016-05-27 | 2016-08-10 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | Method for building stratum velocity model |
| CN106094032A (en) * | 2016-08-30 | 2016-11-09 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | A kind of method building formation velocity model |
| CN109839660A (en) * | 2018-11-08 | 2019-06-04 | 成都捷科思石油天然气技术发展有限公司 | A method of velocity depth model is established using prestack trace gather data |
Family Cites Families (1)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US10139509B2 (en) * | 2016-03-23 | 2018-11-27 | Chevron U.S.A. Inc. | Enhanced seismic imaging systems and methods employing modified velocity model contrasts |
-
2020
- 2020-06-28 CN CN202010594866.5A patent/CN111722283B/en active Active
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US4982382A (en) * | 1989-06-14 | 1991-01-01 | Mobil Oil Corporation | Method of modeling subsurface formations |
| CN102967882A (en) * | 2012-11-16 | 2013-03-13 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | Method for building layer velocity model of stratum |
| CN104502997A (en) * | 2015-01-22 | 2015-04-08 | 中国石油化工集团公司 | Method for using fracture density curve to forecast fracture density body |
| CN105842736A (en) * | 2016-05-27 | 2016-08-10 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | Method for building stratum velocity model |
| CN106094032A (en) * | 2016-08-30 | 2016-11-09 | 中国石油集团川庆钻探工程有限公司地球物理勘探公司 | A kind of method building formation velocity model |
| CN109839660A (en) * | 2018-11-08 | 2019-06-04 | 成都捷科思石油天然气技术发展有限公司 | A method of velocity depth model is established using prestack trace gather data |
Also Published As
| Publication number | Publication date |
|---|---|
| CN111722283A (en) | 2020-09-29 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN109709603B (en) | Seismic horizon identification and tracking method and system | |
| CN111722283B (en) | Stratum velocity model building method | |
| CN104502997B (en) | A kind of method of utilization fracture spacing curve prediction fracture spacing body | |
| CN108802812B (en) | Well-seismic fusion stratum lithology inversion method | |
| CN102508293B (en) | Pre-stack inversion thin layer oil/gas-bearing possibility identifying method | |
| CN105954804B (en) | Shale gas reservoir fragility earthquake prediction method and device | |
| CN104597490A (en) | Multi-wave AVO reservoir elastic parameter inversion method based on precise Zoeppritz equation | |
| CN104570101A (en) | AVO (amplitude versus offset) three-parameter inversion method based on particle swarm optimization | |
| CN111025387B (en) | Pre-stack earthquake multi-parameter inversion method for shale reservoir | |
| WO2016008105A1 (en) | Post-stack wave impedance inversion method based on cauchy distribution | |
| CN114723095A (en) | Missing well logging curve prediction method and device | |
| CN107861149B (en) | Based on the prestack P-S wave velocity ratio analogy method under drive waveform | |
| CN112100906B (en) | Data-driven large-scale density modeling method, computing device and storage medium | |
| CN111983683B (en) | Prediction method and system for lake-facies limestone reservoir under low-well condition | |
| CN112147677B (en) | Method and device for generating parameter tag data of oil and gas reservoir | |
| CN113552624B (en) | Porosity prediction method and device | |
| CN111273349B (en) | Transverse wave velocity extraction method and processing terminal for seabed shallow sediment layer | |
| CN111025388B (en) | Multi-wave combined prestack waveform inversion method | |
| CN106226814B (en) | Utilize converted shear wave seismic data inversion reservoir S-wave impedance and the method for density | |
| CN113608258A (en) | Self-consistent deep learning method for constructing high-resolution wave impedance inversion label | |
| CN113806674A (en) | Method and device for quantifying longitudinal dimension of ancient river channel, electronic equipment and storage medium | |
| Xie et al. | A novel genetic inversion workflow based on spectral decomposition and convolutional neural networks for sand prediction in Xihu Sag of East China Sea | |
| CN114075973B (en) | Stratum element logging curve reconstruction method and device | |
| CN113514904B (en) | Stratum parameter model establishing method and device | |
| CN104375171B (en) | A kind of High-resolution Seismic Inversion method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |