A method of mapping trend velocity by controlling earthquake attributes

Abstract

The present invention relates to the technical field of seismic data structural interpretation, and in particular to a method for seismic attribute controlled trend speed mapping. The method comprises: step S1, fine well seismic calibration; step S2, seismic structural interpretation; step S3, calculating well point time residuals; step S4, selecting seismic attributes; step S5, extracting seismic attribute trend surfaces; step S6, calculating residual grids; step S7, correcting T 0 Grid; Step S8, establish an empty variable velocity field; The seismic attribute control trend variable velocity mapping method provided by the present invention has high accuracy, can meet the problem of converting the actual drilling trajectory of the horizontal well from the depth domain to the time domain, and maximize the work efficiency.

Description

Technical Field

[0001] The invention relates to the technical field of seismic data structural interpretation, and in particular to a seismic attribute controlled trend speed-varying mapping method. Background Art

[0002] With the continuous deepening of oilfield development, the oilfield has entered the stage of tapping the potential of residual oil. High-angle wells and horizontal wells are the main well types for tapping the potential of residual oil. In the design of horizontal wells, accurate three-dimensional structural models are the basis and guarantee for horizontal well design. The oilfield structural model includes a layer model and a fault model. This requires that the high-precision velocity field established can simultaneously meet the requirements of high-precision time-to-depth conversion of both seismic interpretation layers and faults. In the process of on-site drilling tracking of horizontal wells, it is often necessary to convert the depth domain well trajectory data of the actual drilling of the horizontal well into time domain data on-site and in real time, and then load it into the seismic or inversion profile, and make on-site adjustments to the drilling in combination with the profile characteristics. This requires that the high-precision velocity field established can not only perform time-to-depth conversion from the time domain to the depth domain, but also perform time-to-depth conversion from the depth domain to the time domain without distinction.

[0003] At present, the most commonly used variable velocity mapping technology for oilfield development blocks with a high degree of well control is to use the wellbore velocity, seismic interpretation layer and velocity spectrum to jointly establish the velocity field. The specific method is to use the Dix formula to obtain the layer velocity from the stacked velocity of the velocity spectrum (this formula is based on the assumption that the velocity field does not change laterally). The converted layer velocity is used as the initial velocity, and the velocity information on the well is added. The stratigraphic grid of the region is established using the seismic interpretation layer. The stratigraphic grid is used as the control unit, and a three-dimensional velocity body is generated by spatial interpolation. The three-dimensional velocity body is used to perform time-to-depth conversion on the seismic interpretation time domain layer, and the converted depth layer is then error corrected (to meet the requirements of structural modeling accuracy) to obtain the depth domain layer. However, the time domain fault of the seismic interpretation is directly converted to time-to-depth to obtain the depth domain fault.

[0004] The velocity field established by the conventional variable velocity mapping method is not accurate enough, and structural error correction is required after the time-depth conversion of the layer, while the time-depth conversion is directly performed on the fault. This leads to a mismatch between the layers and faults converted to the depth domain by the conventional method, and there is an intersection or deviation problem in the spatial relationship between the two. Such a velocity field cannot meet the conversion problem of the depth domain to the time domain of the actual drilling trajectory of the horizontal well. Summary of the invention

[0005] 1. Technical issues to be resolved

[0006] The present invention provides a method for variable velocity mapping of seismic attribute controlled trend, so as to overcome the defects in the prior art that the velocity field established by the conventional variable velocity mapping method has low precision and mismatch between the layers and faults converted to the depth domain, resulting in the velocity field being unable to satisfy the conversion from the depth domain to the time domain of the actual drilling trajectory of the horizontal well.

[0007] (II) Technical solution

[0008] In order to solve the above problems, the present invention provides a method for mapping seismic attribute controlled trend speed, comprising:

[0009] Step S1: Perform fine well seismic calibration on multiple wells in the work area;

[0010] Step S2: According to the calibration results of step S1, the marked seismic phase axis is finely tracked in the whole area to complete the seismic layer and fault interpretation corresponding to the geological stratification;

[0011] Step S3: According to the altitude depth of the geological stratification, the time value of the interpretation layer is constructed and the regional time-depth relationship is fitted; using the fitted time-depth relationship, the time residual value at the well point is calculated;

[0012] Step S4: According to the quality of the seismic data in this area, a suitable attribute extraction time window is selected, and multiple seismic attributes are extracted and optimized within the time window;

[0013] Step S5: extracting trend surfaces from the seismic attributes selected in step 4. The extracted trend surfaces should have good correlation with the high and low distribution characteristics of the regional structure.

[0014] Step S6: Taking the geological layered well point time residual in step 3 as a known quantity and the trend surface extracted from the seismic attributes in step 5 as a trend, select a suitable trend-controlled interpolation algorithm to calculate the time residual grid for the entire area;

[0015] Step S7: gridding the time domain horizon obtained by the structural interpretation in step 2 to obtain an equal T0 grid, and correcting the equal T0 grid using the time residual grid calculated in step 6 to obtain a corrected equal T0 grid;

[0016] Step S8: According to the altitude depth of the well point geological stratification and the time value of the corrected equal T0 grid well point position, the regional time-depth relationship is fitted twice, and the space-varying velocity field of the whole area is established based on the fitted time-depth relationship; according to the space-varying velocity fields of multiple layers, a three-dimensional velocity field of the whole area is established to realize the time-depth conversion of layers and faults.

[0017] Preferably, the step S1 specifically includes: selecting as many wells as possible in the work area, using acoustic time difference curves and density curves, establishing artificial synthetic records, and determining the time-depth relationship of each well.

[0018] Preferably, in step S4, the principle for selecting the attribute extraction time window is: taking the time value infinitely close to time 0 as the starting time and the interpreted layer as the ending time, performing correlation analysis on the extracted seismic attributes and the time residual value at the well point obtained in step 3, and selecting a good seismic attribute with the highest correlation and matching the geological knowledge.

[0019] Preferably, in step S3,

[0020] The calculation formula of the time residual value at the well point is: Terror=FV(z)-Twell

[0021] Where: Terror is the time residual value of the geological stratification at the well point, in ms; FV(z) is the time-depth relationship formula obtained by fitting the geological stratification time value and altitude depth; z is the altitude depth of the geological stratification at the well point, in m; Twell is the time value of the corresponding layer of the geological stratification at the well point, in ms.

[0022] Preferably, in step S5, a multivariate regression analysis is performed based on the attribute data at the well point to obtain a smooth mathematical surface, which is a trend surface. For high and steep structural areas with large structural fluctuations, a high-order trend surface is selected to fit the trend; the trend surface equation can be divided into:

[0023] The linear trend surface equation is: Z = b0 + b1x + b2y

[0024] The quadratic trend surface equation is: Z = b0 + b1x + b2y + b3x 2 +b4xy+ b5y 2

[0025] The L-order trend surface equation is: Z = b0 + b1x + b2y + b3x 2 +b4xy+ b5y 2 +……+ b n-1 xy l-1 +b n y l

[0026] Where: Z is the trend surface function; b i is the coefficient, dimensionless; n is the coefficient number, dimensionless; l is the power series, dimensionless; x is the horizontal coordinate, unit is m; y is the vertical coordinate, unit is m.

[0027] Preferably, in step S8, the quadratic fitting best time-depth relationship is an exponential function, and its formula is: FV(z)=119.98*e 0.0016*z , the correlation coefficient is 0.99;

[0028] Where FV(z) is the time-depth relationship formula obtained by fitting the geological stratification time value and altitude depth; z is the altitude depth of the geological stratification at the well point, in meters.

[0029] (III) Beneficial effects

[0030] The seismic attribute controlled trend variable speed mapping method provided by the present invention fully utilizes the advantage that the number of development wells per unit area in the development zone is greater than the sampling points of the velocity spectrum, takes the well points as the basis, controls the spatial trend with seismic attributes, and establishes a variable speed mapping method that can simultaneously perform bidirectional time-depth conversion from the time domain to the depth domain and from the depth domain to the time domain on the strata and faults. The method has high accuracy, can meet the problem of converting the depth domain to the time domain of the actual drilling trajectory of the horizontal well, and maximizes work efficiency. BRIEF DESCRIPTION OF THE DRAWINGS

[0031] Figure 1 A flow chart showing the method for mapping earthquake attribute controlled trend speed changes implemented by the present invention is shown;

[0032] Figure 2 A schematic diagram showing the scope of seismic attribute extraction;

[0033] Figure 3 Well point residual contour maps are shown;

[0034] Figure 4 The regional time-depth relationship trend diagram of the quadratic fit is shown. DETAILED DESCRIPTION

[0035] The present invention is described in detail below with reference to the accompanying drawings and embodiments.

[0036] like Figure 1 As shown, the present invention provides a method for mapping earthquake attribute controlled trend speed, which specifically includes the following steps:

[0037] Step S1: Fine well seismic calibration:

[0038] Fine well-seismic calibration is performed on as many development wells as possible in the work area. The correlation coefficient between artificial synthetic records and seismic records can be improved through appropriate stretching and compression. The purpose of calibration is to determine the time-depth relationship of each well, especially the layers within the vertical range of variable speed mapping (target layer section) must have a high well-seismic correlation coefficient, which is the basis for achieving high-precision variable speed mapping. Through calibration, the seismically identifiable geological stratification interface is marked on the seismic profile. To improve accuracy, as many seismic phase axes corresponding to geological stratification as possible are marked. In this example, the calibration of 36 wells and 5 oil group-level stratifications was completed. Through calibration, the seismic reflection phase axes corresponding to 5 strata were identified on the seismic profile.

[0039] Step S2: Seismic structural interpretation:

[0040] The scope of structural interpretation is mainly developed areas with a high degree of well control. According to the calibration results of step S1, the 5 seismic reflection phase axes marked are tracked in detail in the whole area to complete the seismic layer and fault interpretation corresponding to the geological stratification. First, the well system layer is connected. By pulling multiple well connection lines, check whether each well is calibrated to the same reflection phase axis. If there are inconsistent wells, return to step S1 and re-make the synthetic record until the same geological layer of all wells is calibrated on the same seismic reflection axis. Then establish a 4×4 line interpretation skeleton, and finally make the calibrated layers of the synthetic record tracked to the same level in the three-dimensional data body to achieve the interpretation and tracking of the seismic layer.

[0041] Step S3, calculate the well point time residual:

[0042] First, according to the altitude depth of geological stratification and the time value of the interpretation layer corresponding to the geological stratification at the well point, the regional time-depth relationship formula is calculated by statistical fitting. Finally, the time residual value at the well point is calculated using the fitted time-depth relationship trend formula.

[0043] In this implementation, the time-depth relationship index correlation coefficient of the initial fitting is the highest, and its formula is: FV(z)=136.2*e 0.0015*z The initial time-depth relationship correlation coefficient is 0.93. Where: FV(z) is the time-depth relationship formula obtained by fitting the geological stratification time value and altitude depth; z is the altitude depth of the geological stratification at the well point, in meters.

[0044] The calculation formula of the residual Terror of the well point in this interpretation layer is: 136.2*e 0.0015*z -Twell, where: Terror is the time residual value of the geological stratification at the well point, in ms; z is the altitude depth of the geological stratification at the well point, in m; Twell is the time value of the corresponding layer of the geological stratification at the well point, in ms.

[0045] Step S4: Optimize earthquake attributes:

[0046] According to the quality of seismic data in this area, select a suitable attribute extraction time window, and extract multiple seismic attributes within this time window. The principle of selecting the attribute extraction time window is: the start time is as close to the time 0 as possible, and the end time is the interpretation layer. In this example, according to the effective data range of the seismic data, 150ms is selected as the time window start value, and the target layer is selected as the time window end time (about 1400ms). Figure 2As shown. Calculate the seismic attributes within this time window, give priority to energy-type seismic attributes, such as root mean square amplitude attributes, energy average attributes, etc., and filter the attributes to remove singular values ​​in the attributes. Perform a correlation analysis on the extracted seismic attributes and the time residual values ​​at the well points obtained in step 3, and select a good seismic attribute with the highest correlation and matching the geological knowledge. In this example, the root mean square amplitude attribute with the highest correlation coefficient is taken as the preferred seismic attribute.

[0047] Step S5, extracting earthquake attribute trend surface:

[0048] The trend surface extraction technology is used to extract the trend surface of the seismic attributes selected in step 4. Based on the attribute data at the well point, a multivariate regression analysis is performed to obtain a smooth mathematical surface, which is called a trend surface. For high-steep structural areas with large structural fluctuations, a high-order trend surface can be selected to better fit the trend. The trend surface equation can be divided into:

[0049] The linear trend surface equation is: Z = b0 + b1x + b2y

[0050] The quadratic trend surface equation is: Z = b0 + b1x + b2y + b3x 2 +b4xy+ b5y 2

[0051] The L-order trend surface equation is: Z = b0 + b1x + b2y + b3x 2 +b4xy+ b5y 2 +……+ b n-1 xy l-1 +b n y l

[0052] Where: Z is the trend surface function; b i is the coefficient, dimensionless; n is the coefficient number, dimensionless; l is the power series, dimensionless; x is the horizontal coordinate, unit is m; y is the vertical coordinate, unit is m.

[0053] In this implementation, the trend surface of the root mean square amplitude attribute obtained in step 4 is extracted using a cubic trend surface.

[0054] Step S6, calculate the residual grid:

[0055] Taking the time residual of the geological layer well point in step 3 as a known quantity and the trend surface extracted from the seismic attributes in step 5 as the trend, select an appropriate trend-controlled interpolation algorithm to calculate the time residual grid for the entire area. This implementation uses the least squares method to calculate the time residual grid for the entire area under the control of the tertiary trend surface calculated in step 5 and the time residual value at the well point calculated in step 3 as a known value, as shown in Figure 3 shown.

[0056] Step S7, correcting the T0 grid:

[0057] First, the time domain horizon obtained by the construction interpretation in step 2 is gridded to obtain an equal T0 grid, and then the time residual grid calculated in step 6 is added to the equal T0 grid to obtain a corrected equal T0 grid.

[0058] Step S8, establishing an empty variable velocity field:

[0059] First, according to the geological layering elevation depth of the well point and the corrected equal T0 grid time value of the well point, the regional time-depth relationship is fitted twice, as shown in Figure 4 In this example, the best time-depth relationship of quadratic fitting is an exponential function, and its formula is: FV(z)=119.98*e 0.0016*z , the correlation coefficient is 0.99.

[0060] Where FV(z) is the time-depth relationship formula obtained by fitting the geological stratification time value and altitude depth; z is the altitude depth of the geological stratification at the well point, in meters.

[0061] Based on the fitted time-depth relationship, the space-variable velocity field of the entire layer is established. Under the control of the space-variable velocity fields of multiple layers, the velocity curves calculated on the well are used to establish the three-dimensional velocity field of the entire area through interpolation to achieve the time-depth conversion of layers and faults.

[0062] The above-mentioned embodiments are only embodiments of the present invention, and the description thereof is relatively specific and detailed, but it cannot be understood as limiting the scope of the present invention. It should be pointed out that, for a person skilled in the art, without departing from the concept of the present invention, several variations, equivalent substitutions, improvements, etc. can be made, which all belong to the protection scope of the present invention. Therefore, the protection scope of the present invention patent shall be subject to the attached claims.

Claims

1. A method for mapping earthquake attribute controlled trend speed, characterized in that: include: Step S1: Perform fine well seismic calibration on multiple wells in the work area; Step S2: According to the calibration results of step S1, the marked seismic phase axis is finely tracked in the whole area to complete the seismic layer and fault interpretation corresponding to the geological stratification; Step S3: According to the altitude depth of the geological stratification, the time value of the interpretation layer is constructed and the regional time-depth relationship is fitted; using the fitted time-depth relationship, the time residual value at the well point is calculated; Step S4: According to the quality of the seismic data in this area, a suitable attribute extraction time window is selected, and multiple seismic attributes are extracted and optimized within the time window; Step S5: extracting trend surfaces from the seismic attributes selected in step 4. The extracted trend surfaces should have good correlation with the high and low distribution characteristics of the regional structure. Step S6: Taking the geological layered well point time residual in step 3 as a known quantity and the trend surface extracted from the seismic attributes in step 5 as a trend, select a suitable trend-controlled interpolation algorithm to calculate the time residual grid for the entire area; Step S7: gridding the time domain horizon obtained by the structural interpretation in step 2 to obtain an equal T0 grid, and correcting the equal T0 grid using the time residual grid calculated in step 6 to obtain a corrected equal T0 grid; Step S8: According to the altitude depth of the well point geological stratification and the time value of the corrected equal T0 grid well point position, the regional time-depth relationship is fitted twice, and the space-varying velocity field of the whole area is established based on the fitted time-depth relationship; according to the space-varying velocity fields of multiple layers, a three-dimensional velocity field of the whole area is established to realize the time-depth conversion of layers and faults.

2. The method for mapping earthquake attribute control trend by space calibration and speed change according to claim 1, characterized in that: The step S1 specifically includes: selecting as many wells as possible in the work area, using the acoustic time difference curve and density curve, establishing artificial synthetic records, and determining the time-depth relationship of each well.

3. The method for mapping earthquake attribute controlled trend by space calibration and speed change according to claim 1, characterized in that: In step S4, the principle of selecting the attribute extraction time window is: taking the time value infinitely close to time 0 as the starting time and taking the interpreted layer as the ending time, performing a correlation analysis on the extracted seismic attributes and the time residual value at the well point obtained in step 3, and selecting a good seismic attribute with the highest correlation and matching the geological knowledge.

4. The method for mapping earthquake attribute control trend by space calibration and speed change according to claim 1, characterized in that: In the step S3, The calculation formula of the time residual value at the well point is: Terror=FV(z)-Twell Where: Terror is the time residual value of the geological layer at the well point, in ms; FV(z) is the time-depth relationship formula obtained by fitting the geological layer time value and altitude depth; z is the altitude depth of the geological layer at the well point, in m; Twell is the time value of the corresponding layer of the geological stratification at the well point, in ms.

5. The method for mapping earthquake attribute control trend by space calibration and speed change according to claim 1, characterized in that: In the step S5, a multivariate regression analysis is performed based on the attribute data at the well point to obtain a smooth mathematical surface, which is a trend surface; For high and steep structural areas with large structural fluctuations, high-order trend surfaces are used to fit the trend; trend surface equations can be divided into: The linear trend surface equation is: Z = b0 + b1x + b2y The quadratic trend surface equation is: Z = b0 + b1x + b2y + b3x 2 +b4xy+ b5y 2 The L-order trend surface equation is in the form of: Z = b0 + b1x + b2y + b3x 2 +b4xy+ b5y 2 +……+ b n-1 xy l-1 +b n y l Where: Z is the trend surface function; b i is the coefficient, dimensionless; n is the coefficient number, dimensionless; l is the power series, dimensionless; x is the horizontal coordinate, unit is m; y is the vertical coordinate, unit is m.

6. The method for mapping earthquake attribute controlled trend by space calibration and variable speed according to claim 1, characterized in that: In step S8, the quadratic fitting best time-depth relationship is an exponential function, and its formula is: FV(z)=119.98*e 0.0016*z , the correlation coefficient is 0.99; Where FV(z) is the time-depth relationship formula obtained by fitting the geological stratification time value and altitude depth; z is the altitude depth of the geological stratification at the well point, in meters.

Images