CN113671566A - Method for calculating crack parameters based on depth domain seismic data - Google Patents

Method for calculating crack parameters based on depth domain seismic data Download PDF

Info

Publication number
CN113671566A
CN113671566A CN202110937680.XA CN202110937680A CN113671566A CN 113671566 A CN113671566 A CN 113671566A CN 202110937680 A CN202110937680 A CN 202110937680A CN 113671566 A CN113671566 A CN 113671566A
Authority
CN
China
Prior art keywords
data
azimuth
attribute data
data volume
frequency
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.)
Granted
Application number
CN202110937680.XA
Other languages
Chinese (zh)
Other versions
CN113671566B (en
Inventor
王占磊
蒋裕强
谷一凡
付永红
周亚东
李杪
尹兴平
王子萌
蔡光银
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Southwest Petroleum University
Original Assignee
Southwest Petroleum University
Priority date (The priority date 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 date listed.)
Filing date
Publication date
Application filed by Southwest Petroleum University filed Critical Southwest Petroleum University
Priority to CN202110937680.XA priority Critical patent/CN113671566B/en
Publication of CN113671566A publication Critical patent/CN113671566A/en
Application granted granted Critical
Publication of CN113671566B publication Critical patent/CN113671566B/en
Active legal-status Critical Current
Anticipated expiration legal-status Critical

Links

Images

Classifications

    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/30Analysis
    • G01V1/306Analysis for determining physical properties of the subsurface, e.g. impedance, porosity or attenuation profiles
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V1/00Seismology; Seismic or acoustic prospecting or detecting
    • G01V1/28Processing seismic data, e.g. for interpretation or for event detection
    • G01V1/36Effecting static or dynamic corrections on records, e.g. correcting spread; Correlating seismic signals; Eliminating effects of unwanted energy
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/50Corrections or adjustments related to wave propagation
    • G01V2210/51Migration
    • G01V2210/512Pre-stack
    • GPHYSICS
    • G01MEASURING; TESTING
    • G01VGEOPHYSICS; GRAVITATIONAL MEASUREMENTS; DETECTING MASSES OR OBJECTS; TAGS
    • G01V2210/00Details of seismic processing or analysis
    • G01V2210/60Analysis
    • G01V2210/64Geostructures, e.g. in 3D data cubes
    • G01V2210/647Gas hydrates

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 provides a method for calculating fracture parameters based on depth domain seismic data, which comprises the following steps: step 10, acquiring a related azimuth prestack depth migration data volume; step 20, after band-pass filtering and weighting processing are carried out on the prestack depth migration data volumes of all the azimuth angles, coherent volume calculation and edge signal enhancement processing are carried out, and then a first attribute data volume of each azimuth angle is obtained; step 30, calculating the first attribute data bodies of all the azimuth angles by using an ant body algorithm respectively to generate second attribute data bodies of all the azimuth angles; then, carrying out weighted fusion processing on the second attribute data volume of each azimuth angle to obtain a third attribute data volume; and step 40, predicting and evaluating the development strength and the direction of the fracture of the target interval by using the third attribute data body. The crack prediction effect obtained by the method is superior to the result obtained by the conventional ant body calculation method.

Description

Method for calculating crack parameters based on depth domain seismic data
Technical Field
The invention relates to the technical field of seismic data interpretation processing in geophysical exploration, in particular to a method for calculating fracture parameters based on depth domain seismic data.
Background
Fracture prediction plays an important role in oil and gas exploration and is also the key point of geophysical inversion research. If a crack develops in a reservoir, the communication of the reservoir is greatly influenced, and the drilling of the crack + the pore type reservoir usually obtains high-yield industrial oil and gas flow. Therefore, the exploration of the fractured reservoir is the key point of exploration in the early exploration, and an explorationist hopes to find the reservoir with the superposed fractures and pores, for example, the prediction of the fractures in the marine reef and land tight sandstone oil and gas exploration in the Sichuan basin is also in a relatively important position.
The conventional fracture prediction technology methods are various, and can be basically divided into two main categories of pre-stack inversion and post-stack inversion, and an attribute data volume for fracture prediction can be obtained through inversion or attribute calculation. In seismic exploration, the amplitude, velocity, two-way travel time and frequency classes in seismic data of different azimuths have different response characteristics to a crack, and the amplitude, the velocity, the two-way travel time and the frequency classes are also the prerequisites for crack prediction. According to research data, the characteristic that the amplitude and the two-way travel time value of the seismic P wave change along with the azimuth angle when the underground target layer passes through the crack body is found when the crack is developed. For example, Xiangyang Li, Fracture detection using azimuthally variation of P-wave movable out from an orthogonal magnetic series circuits published in Geophysics 1999. Research data show that: the reflected P wave on the top layer of the crack body does not change along with the azimuth angle during double-pass, the amplitude slightly changes along with the azimuth angle, and the anisotropy is not obviously reflected. The reflected P wave of the bottom layer of the crack body shows strong azimuth anisotropy, the reflection amplitude and the two-way travel time of the P wave are related to the measurement line and the crack azimuth, the amplitude is strongest when the measurement line is parallel to the crack, and the two-way travel time is shortest; with the increase of the included angle between the measuring line and the crack, the amplitude is gradually reduced, and the two-way travel time is gradually lengthened; when the measuring line is vertical to the crack direction, the amplitude is weakest, and the double-pass travel time is longest. Therefore, the P wave shows strong azimuth anisotropy characteristic when passing through the crack body, the influence degree of the crack on the P wave response depends on the included angle between the crack azimuth and the measuring line and the scale of a certain underground open crack system, and the phenomenon has larger influence on the travel time of the bottom layer particularly. In the depth domain, the anisotropy of the different azimuthal depth data generated by the underlayer as the P-wave passes through the fracture volume will be more pronounced due to the addition of velocity data.
In addition, most of the conventional fracture prediction technologies are calculated by using time-domain seismic data, but the time-domain seismic data are often influenced by interval velocity changes, so that phenomena of inaccurate seismic imaging, inconsistency of reflection phase axes with actual geological conditions and the like exist. Therefore, using time domain gathers and post-stack seismic data to predict fractures introduces computational errors and risks.
Disclosure of Invention
The invention aims to provide a method for calculating fracture parameters based on depth domain seismic data, so as to solve the problems of the conventional fracture prediction technology.
The invention provides a method for calculating fracture parameters based on depth domain seismic data, which comprises the following steps:
step 10, acquiring a related azimuth prestack depth migration data volume;
step 20, after band-pass filtering and weighting processing are carried out on the prestack depth migration data volumes of all the azimuth angles, coherent volume calculation and edge signal enhancement processing are carried out, and then a first attribute data volume of each azimuth angle is obtained;
step 30, calculating the first attribute data bodies of all the azimuth angles by using an ant body algorithm respectively to generate second attribute data bodies of all the azimuth angles; then, carrying out weighted fusion processing on the second attribute data volume of each azimuth angle to obtain a third attribute data volume;
and step 40, predicting and evaluating the development strength and the direction of the fracture of the target interval by using the third attribute data body.
Further, step 1 comprises the following substeps:
step 11, setting an azimuth angle range;
step 12, carrying out azimuth division on the omnidirectional angle gather data according to a set azimuth range to obtain n azimuth gather data volumes;
and step 13, processing the n azimuth gather data volumes by respectively adopting the same prestack depth migration processing technical process and parameters to obtain the related n azimuth prestack depth migration data volumes.
Further, the method for dividing the azimuth angle of the omni-directional angle gather data according to the set azimuth angle range in step 12 is as follows: based on the symmetry relationship, converting the azimuth gather data of 0-360 degrees into the azimuth gather data of 0-180 degrees; and dividing the azimuth gather data of 0-180 degrees according to a set azimuth range to obtain n azimuth gather data volumes.
Further, in step 13, before the n azimuth gather data volumes are processed by using the same prestack depth migration processing technical process and parameters, migration preprocessing including denoising and/or resolution improvement is performed on each azimuth gather data volume.
Further, the method for performing band-pass filtering and weighting on each azimuth prestack depth migration data volume in step 20 includes the following sub-steps:
(1) calculating an azimuth prestack depth migration data volume through a layer velocity-depth model, and converting the azimuth prestack depth migration data volume into a time domain seismic data volume;
(2) according to the frequency band range of the target interval in the converted time domain seismic data volume, carrying out spectrum analysis on different crack development intensity areas in the known well of the time domain seismic data, and determining tuning frequency data sets corresponding to the related different crack development intensities;
(3) dividing tuning frequency data sets with different crack development strengths according to crack prediction precision, expert experience, geological knowledge and the like, and determining different m frequency range ranges; wherein m is more than or equal to 2;
(4) transforming the time domain seismic data volume to a frequency domain by using discrete Fourier transform, and performing band-pass filtering on the time domain seismic data volume transformed to the frequency domain by using related m frequency band ranges to obtain a series of related m time domain filtering data volumes;
(5) setting relevant weighting coefficients for a series of m time domain filtering data volumes to carry out weighting fusion processing to obtain a time domain weighting data volume;
(6) and (3) converting the time domain weighted data volume into a depth domain seismic data volume after the speed-depth model calculation.
(7) And (4) performing the steps (1) to (6) on each azimuth prestack depth migration data body by analogy, completing the band-pass filtering and weighting processing of each azimuth prestack depth migration data body, and obtaining the seismic attribute data body after the band-pass filtering and weighting processing.
Further, in step (3), the frequency data set for setting the development intensity and scale of the crack is [ f1,f2,…,fn]Determining the maximum frequency value f from the frequency data setkMinimum frequency value fbAnd Δ d is a design tolerance, the calculation formula of the frequency range is as follows:
fget up=fb-△d
fStop block=fk+△d
△d=1/n(fk-fb)
In the formula (f)Get upIs the starting frequency of the frequency band range; f. ofStop blockIs the stop frequency of the frequency band range; n is the number of frequencies of the frequency data set participating in the calculation; f. ofGet upMinimum frequency value, f, greater than the entire seismic frequency bandStop blockLess than the maximum frequency value for the entire seismic band.
Further, there may be an overlap between the frequency ranges set in step (3), but the overlap should be no more than 0.25(Δ d)1+△d2);△d1、△d2Respectively, the design tolerances of the two frequency band ranges set.
Further, in the step 20, when performing coherent body calculation and edge signal enhancement processing, coherent body calculation and edge signal enhancement processing may be performed only on the target interval of each azimuth seismic attribute data volume after band-pass filtering and weighting processing by using the layer data of the top layer and the bottom layer of the target interval.
Further, in step 30, the first attribute data volume of each azimuth is calculated by using an ant body algorithm, and the specific operation of generating the second attribute data volume of each azimuth is as follows: and after the layer leveling is carried out on the first attribute data bodies of all the azimuth angles by using the layer data of the bottom layer of the target layer section, ant body calculation is carried out on the target layer section, and then the ant body calculation result is subjected to reverse layer leveling to obtain all the second attribute data bodies.
Further, in step 30, a calculation formula for performing weighted fusion processing on the second attribute data volume of each azimuth angle is as follows:
Figure BDA0003213440640000051
in the formula, P is a third attribute data volume after weighted fusion processing; miIs the ith second attribute data volume; k is a radical ofiFor the ith second attribute data volume MiThe set weighting coefficients of (a) are set,
Figure BDA0003213440640000052
n is the number of second attribute data volumes participating in the weighted fusion processing;
the determination of the weighting coefficient is determined according to the test condition of the weighted fusion result of the data values of the development intensity of the fractures in different directions in the related well, and specifically comprises the following steps:
for the determination of the weighting coefficients of the n second attribute data volumes, the determination is made according to the following calculation formula:
P1=k1M11+k2M21+k3M31+...+knMn1
P2=k1M12+k2M22+k3M32+...+knMn2
Pj=k1M1j+k2M2j+k3M3j+...+knMnj
Pn-1=k1M1n-1+k2M2n-1+k3M3n-1+...+knMnn-1
1=k1+k2+k3+...+kn
Pj=Fj/B
in the above formula, PjA third property data value, M, corresponding to the jth point in a third property data volume of a fracture of a known well1j、M2j、M3j、...、MnjRespectively corresponding to each second attribute data value of the j point in the second attribute data body; fjFracture development intensity data value, k, of the j point of a fracture of a known welliFor the ith second attribute data volume MiB is a crack growth intensity calculation coefficient of the third attribute data volume.
In summary, due to the adoption of the technical scheme, the invention has the beneficial effects that:
the invention utilizes a sub-azimuth depth domain seismic data volume to calculate fracture development parameters. Setting a related azimuth angle range for the omnidirectional angle gather, and respectively carrying out prestack depth migration processing on related azimuth angle gather data bodies to obtain related azimuth angle prestack depth migration data bodies; performing band-pass filtering and weighting processing on the prestack depth migration data volumes of all the azimuth angles, and performing coherent body calculation and edge signal enhancement processing to obtain first attribute data volumes of all the azimuth angles; and calculating each first attribute data body by using an ant body algorithm, generating each second attribute data body, and performing weighted fusion processing on the second attribute data bodies to obtain a third attribute data body. The crack prediction effect obtained by the method is superior to the result obtained by the conventional ant body calculation method. The technical method of the invention is used for predicting the development parameters of the cracks of the marine shale reservoir in a certain area, has good effect, has high matching degree with the actual drilling well data in the area, and can guide the horizontal well layout and fracturing operation in the area.
Drawings
In order to more clearly illustrate the technical solutions of the embodiments of the present invention, the drawings in the embodiments will be briefly described below, it should be understood that the following drawings only illustrate some embodiments of the present invention, and therefore should not be considered as limiting the scope, and for those skilled in the art, other related drawings can be obtained according to the drawings without inventive efforts.
FIG. 1 is a flow chart of a method of calculating fracture parameters based on depth domain seismic data in an embodiment of the invention.
Detailed Description
In order to make the objects, technical solutions and advantages of the embodiments of the present invention clearer, the technical solutions in the embodiments of the present invention will be clearly and completely described below with reference to the drawings in the embodiments of the present invention, and it is obvious that the described embodiments are some, but not all, embodiments of the present invention. The components of embodiments of the present invention generally described and illustrated in the figures herein may be arranged and designed in a wide variety of different configurations.
Thus, the following detailed description of the embodiments of the present invention, presented in the figures, is not intended to limit the scope of the invention, as claimed, but is merely representative of selected embodiments of the invention. All other embodiments, which can be derived by a person skilled in the art from the embodiments given herein without making any creative effort, shall fall within the protection scope of the present invention.
Examples
The method adopts the depth domain seismic attribute data of the divided azimuth angle to calculate the crack development parameters, thereby achieving the purpose of accurately predicting the crack development conditions of different scales of the target layer in the research area. Therefore, as shown in fig. 1, the present embodiment provides a method for calculating fracture parameters based on depth domain seismic data, including the following steps:
step 10, obtaining a related azimuth prestack depth migration data volume:
step 11, setting an azimuth angle range; specifically, the method comprises the following steps: setting the median of the azimuth angle range as a central angle, and then setting an azimuth angle according to the relation between the central angle and the crack trend in the research area; generally, the direction line of the central angle is approximately vertical to the direction of the crack or is intersected with the direction of the crack at a large angle, and the azimuth angle range is obtained by taking the central angle as the center and adding/subtracting a certain fixed angle value. The fracture direction in the research area can be determined according to well data, regional geological data, expert experience and the like.
The calculation formula of the central angle is as follows:
Figure BDA0003213440640000071
in the formula, thetaiRepresenting the central angle of the set ith azimuth gather data;
Figure BDA0003213440640000072
representing the minimum azimuth of the ith azimuth gather data;
Figure BDA0003213440640000081
representing the maximum azimuth of the ith azimuth gather data; i is more than or equal to 2.
And step 12, carrying out azimuth division on the omnidirectional angle gather data according to a set azimuth range to obtain n azimuth angle gather data volumes. Specifically, the method comprises the following steps: based on the symmetry relationship, converting the azimuth gather data of 0-360 degrees into the azimuth gather data of 0-180 degrees; dividing the azimuth gather data of 0-180 degrees according to a set azimuth range to obtain n azimuth gather data volumes; the azimuth gather data is typically divided into equal azimuth ranges or unequal azimuth ranges, or a mixture of both.
And step 13, processing the n azimuth gather data volumes by respectively adopting the same prestack depth migration processing technical process and parameters to obtain the related n azimuth prestack depth migration data volumes. The adopted pre-stack depth migration processing technical process and parameters are conventional pre-stack depth migration processing technical process and parameters in the prior art, and are not described herein again. Further, before the n azimuth gather data volumes are processed by using the same prestack depth migration processing technical process and parameters, migration preprocessing, such as denoising and/or resolution improvement, can be performed on each azimuth gather data volume to improve the effect of subsequent processing.
Step 20, after performing band-pass filtering and weighting processing on the prestack depth migration data volume of each azimuth angle, performing coherent body calculation and edge signal enhancement processing, and obtaining a first attribute data volume of each azimuth angle:
and step 21, performing band-pass filtering and weighting processing on the pre-stack depth migration data volume of each azimuth angle to obtain the seismic attribute data volume of each azimuth angle after the band-pass filtering and weighting processing. Specifically, the method comprises the following steps:
(1) and (4) converting the azimuth prestack depth migration data volume into a time domain seismic data volume after the calculation of the interval velocity-depth model.
(2) And according to the frequency band range of the target interval in the converted time domain seismic data body, carrying out spectrum analysis on different crack development intensity areas in the known well of the time domain seismic data, and determining tuning frequency data sets corresponding to the related different crack development intensities.
(3) Dividing tuning frequency data sets with different crack development strengths according to crack prediction precision, expert experience, geological knowledge and the like, and determining different m frequency range ranges; wherein m is more than or equal to 2;
if the frequency data set of a certain set crack development intensity and scale is [ f1,f2,…,fn]Determining the maximum frequency value f from the frequency data setkMinimum frequency value fbAnd Δ d is a design tolerance, the calculation formula of the frequency range is as follows:
fget up=fb-△d
fStop block=fk+△d
△d=1/n(fk-fb)
In the formula (f)Get upIs the starting frequency of the frequency band range; f. ofStop blockIs the stop frequency of the frequency band range; n is the frequency number of the frequency data set participating in calculation, and can also be determined according to expert experience and the like; f. ofGet upShould be greater than the minimum frequency value, f, of the entire seismic frequency bandStop blockShould be less than the maximum frequency value for the entire seismic band. There may be overlap between the frequency ranges set, but the overlap should be no more than 0.25(Δ d)1+△d2),△d1、△d2Design tolerances for the two frequency range to be set, respectively, if the overlap range is greater than 0.25(Δ d)1+△d2) Then, Δ d needs to be resetj 1、△dl 2The frequency band overlapping range is less than or equal to 0.25 (delta d)j 1、△dl 2)。
(4) And transforming the time domain seismic data volume to a frequency domain by using discrete Fourier transform, and performing band-pass filtering on the time domain seismic data volume transformed to the frequency domain by using the related m frequency band ranges to obtain a series of m related time domain filtering data volumes. Generally, through testing related band-pass filtering parameters, and preferably selecting related filtering modules and parameters, the time domain seismic data volume converted into a frequency domain can relatively filter or suppress some low and high frequency interference signals (non-crack and fracture response signals) after the band-pass filtering, and the response of medium and small cracks is highlighted.
(5) And setting related weighting coefficients for a series of m time domain filtering data volumes to carry out weighting fusion processing to obtain a time domain weighting data volume. Preferably, the weighted fusion process may be performed by setting the same or substantially the same weighting coefficient for each time-domain filtered data volume.
(6) And (3) converting the time domain weighted data volume into a depth domain seismic data volume after the speed-depth model calculation.
(7) And (3) performing the steps (1) to (6) on each azimuth prestack depth migration data body, completing band-pass filtering and weighting processing on each azimuth prestack depth migration data body, and obtaining the seismic attribute data body after the band-pass filtering and weighting processing for subsequent calculation.
And step 22, performing coherent body calculation and edge signal enhancement processing on crack prediction by using the seismic attribute data volume of each azimuth angle after band-pass filtering and weighting processing to obtain a first attribute data volume of each azimuth angle.
Regarding the calculation of the seismic attributes for crack prediction, the relevant seismic attribute types may be determined according to the accuracy of crack prediction, expert experience, and the like. And performing seismic attribute calculation on the seismic data volume of each azimuth angle depth domain by using the related determined seismic attribute calculation method and parameters, thereby obtaining the related seismic attribute data volume of each azimuth angle. In general, the seismic attribute is that the response to the fracture is relatively good, the fracture and the non-fracture can be distinguished in data value, and the attribute does not relate to a time-varying parameter. If the geometric properties (such as dip angle, azimuth angle and the like) in the depth domain also reflect the attitude and the morphology of the stratum under the stress state, the local geometric properties of the fracture development in the depth domain also have abnormal response. Similarly, according to the classification of Quincy Chen (1997), among the 8 categories of amplitude, frequency, phase, energy, waveform, attenuation, correlation, ratio, there are also amplitude, waveform, coherence, etc. that respond abnormally to crack development in the deep domain. While the properties of frequency, attenuation, etc. are related to time-varying parameters, the fracture and crack cannot be directly tracked in the depth domain. The invention classifies fracture sensitivity attributes which do not relate to time-varying parameters into a class, and is uniformly named as transverse fracture attributes, such as attributes of most positive/negative body curvature, waveform coherence, chaos, variance and the like, and the fracture sensitivity attributes correspond to geological structural joints, diagenetic shrinkage joints or abnormal pressure joints instead of interlayer page-seam seams. The specific preferred attribute needs to compare the distribution of the regional fault and fracture system and the core fracture imaging, and the attribute with relatively high consistency is selected as the preferred attribute. Generally, the preferred transverse fracture properties from fault distribution and core fracture imaging are primarily one or more of bulk curvature, coherence, chaos, and variance. In the invention, the crack prediction and display of the embodiment are mainly performed by adopting a coherent body calculation algorithm, but not limited to the crack prediction method. In actual operation, which crack prediction method is adopted to participate in the subsequent step calculation can be determined according to the crack prediction precision, expert experience, test results of different crack prediction technologies and the like.
Thus, step 22 specifically includes:
(1) coherent body calculation
In this embodiment, a known well block may be set, with the development strength and direction of the fracture in the known well as known parameters; and testing related coherent body calculation modules and parameters through the known parameters, and preferably, carrying out coherent body calculation on each azimuth angle seismic attribute data body subjected to band-pass filtering and weighting processing by using the optimal coherent body calculation module and parameters to obtain coherent body attributes. In principle, the calculated coherence properties are required to highlight the response of small and medium sized fractures, suppress high frequency noise or other interference caused by non-fractures, and reflect the fracture development strength and relatively match with the measured fracture data in the well.
Preferably, the coherent body calculation algorithm may adopt all coherent body calculation related calculation methods, including a similarity coefficient algorithm, a coherent algorithm based on eigen structure analysis, a covariance matrix algorithm, and the like, and specifically, which calculation method should be selected is determined according to actual conditions, fracture prediction accuracy, and expert experience. Generally, after a coherent body calculation method is selected, the coherent body processing effect is mainly affected by two parameters, namely the number of seismic traces participating in coherent body calculation and the size of a selected window. In the test of the coherent body calculation algorithm and parameters, the test aim is to accurately identify the optimal crack shape of the known crack region, so that the coherent body calculation algorithm and the window size used in the next step are determined. The window size of the coherent body calculation algorithm is determined according to the period T of the seismic reflection wave, and in general, when the coherent window is smaller than T/2, a complete peak or trough cannot be seen due to the narrow coherent window, so that the calculation is carried out according to the complete peak or troughOften reflects noise components; when the coherent window is larger than 3T/2, the coherent window is large, the visual field is wide, a plurality of in-phase axes can be seen, the computed coherent volume data usually reflects the discontinuity of wave groups, and a plurality of fine changes are balanced; when the window is close to T, the imaging effect is best. Therefore, the window size selected by the actual coherent body calculation algorithm should be determined according to the visual period of the earthquake reflected wave or expert experience. In addition, the number of seismic traces involved in coherent body calculations can also affect the effect of coherent processing. If the number of tracks in the calculation model is more, the average effect is larger, the suppression effect on random noise is larger, and the response of the small cracks is suppressed; and the number of tracks in the computational model is reduced, which improves the response to small cracks. In general, the coherent body calculation algorithm selected mainly uses C2And carrying out corresponding tests on the apparent dip angles alpha and beta in the directions of x and y of the parameters in the coherent algorithm, wherein the parameter tests are best to clearly describe the crack form in the test area, and after the tests are finished, carrying out coherent body calculation on the seismic attribute data bodies of all azimuth angles by using the parameters. C2The coherent algorithm mainly defines a rectangular or elliptical analysis window which is centered on an analysis point and contains n channels of seismic data, and takes a local coordinate axis as the center of the analysis point, so that a calculation formula of a similarity coefficient sigma (tau, alpha, beta) is as follows:
Figure BDA0003213440640000121
where the triplet variables (τ, α, β) define a local plane at time τ, α and β are the apparent dip in the x, y directions, respectively, and H represents the hilbert transform. And changing the values of alpha and beta to adjust the inclination direction of the local plane to obtain a similarity coefficient. The maximum similarity coefficient sigma is found out to make the local plane achieve the best fit with the real reflection interface, the similarity coefficient sigma at the moment is the coherence estimation value of the analysis point, and alpha and beta are the apparent dip angles of the current analysis point in the x and y directions. Generally, a smaller coherence value indicates that a crack or larger fault may be present at the location.
(2) Edge signal enhancement processing
In practice, for the preferred relevant seismic attributes, edge signal enhancement processing must be performed before automatic tracking to highlight fracture and crack features and enhance the probability of continuous tracking of fracture and crack for subsequent automatic tracking. The specific edge signal enhancement processing method generally selects a volume attribute filtering method, and selects smooth filtering parameters in the XYZ three directions according to the dip angle and the tendency of fracture and crack.
Furthermore, the horizon data of the top layer and the bottom layer of the target interval can be utilized, and coherent body calculation and edge signal enhancement processing are only carried out on the target interval of each azimuth angle seismic attribute data body after band-pass filtering and weighting processing, so that the calculation time is conveniently shortened, and the calculation efficiency is improved. The method for determining the horizon data of the top layer and the bottom layer of the target interval comprises the following steps: and performing seismic synthetic record calibration on the seismic data volume after each azimuth angle is overlapped in a depth domain by utilizing well logging data, geological stratification data and the like of known wells, determining seismic reflection characteristics of top and bottom layers of a target interval, and performing layer interpretation and tracking on each seismic data volume according to the reflection characteristics of the top and bottom layers of the target interval respectively to obtain the related layer data of the top and bottom layers of the target interval of each azimuth angle seismic attribute data volume.
Step 30, calculating the first attribute data bodies of all the azimuth angles by using an ant body algorithm respectively to generate second attribute data bodies of all the azimuth angles; then, carrying out weighted fusion processing on the second attribute data volume of each azimuth angle to obtain a third attribute data volume;
the method comprises the following specific operations of calculating the first attribute data bodies of all the azimuth angles by using an ant body algorithm and generating the second attribute data bodies of all the azimuth angles: and after the layer leveling is carried out on the first attribute data bodies of all the azimuth angles by using the layer data of the bottom layer of the target layer section, ant body calculation is carried out on the target layer section, and then the ant body calculation result is subjected to reverse layer leveling to obtain all the second attribute data bodies. The main parameters of ant tracking used in the present invention are positive or negative tracking mode and crack inclination angle selected on rose plot. In actual operation, a related test area can be established according to the development strength and direction of small and medium-sized cracks in a known well, and the related ant body attribute calculation module and parameters are tested by utilizing the development strength and direction of different cracks on the known well. In principle, the result (plane or section) traced by the ant body is required to be able to distinguish and identify the development strength and trend of the known cracks relatively clearly.
The weighted fusion processing is performed on the second attribute data volume of each azimuth angle, so that the specific operation of obtaining the third attribute data volume is as follows: selecting related testing areas and crack data on known wells, and testing and determining weighting coefficients of the second attribute data bodies, wherein the related testing mainly comprises calculating, optimizing and testing the weighting coefficients of the second attribute data bodies according to different crack development intensity data values and directions in the known wells in the testing areas, and in principle, after weighted fusion processing of the second attribute data bodies, the acquired ant body attribute data values of the third attribute data bodies in different crack directions are consistent or approximately similar with respect to the same crack development intensity; for the cracks in different directions, when the development strengths of the cracks are consistent or approximately the same, the attribute data values of the third attribute data volume after the weighted fusion processing are consistent or approximately similar. In addition, a certain number (or pieces) of small and medium-sized fractures in different directions of a well which is relatively enough known and known fracture development intensity data values thereof are required, so that the requirements for designing and determining the weighting coefficients of the second attribute data volumes are met. The calculation formula of the weighted fusion processing is as follows:
Figure BDA0003213440640000141
in the formula, P is a third attribute data volume after weighted fusion processing; miIs the ith second attribute data volume; k is a radical ofiFor the ith second attribute data volume MiThe set weighting coefficients of (a) are set,
Figure BDA0003213440640000142
and n is the number of the second attribute data volumes participating in the weighted fusion processing.
Preferably, the design and determination of the weighting coefficient of each second attribute data volume can be determined according to the test condition of the weighted fusion result of the data values of the fracture development strength of different directions in the related well. In general, if there are j second attribute data volumes whose weighting coefficients are determined, j-1 equations must be established and j-1 fracture growth intensity data values must be known to determine the associated weighting coefficients. For example, for the determination of the weighting coefficients for the n second attribute data volumes, the determination may be made according to the following calculation formula:
P1=k1M11+k2M21+k3M31+...+knMn1
P2=k1M12+k2M22+k3M32+...+knMn2
Pj=k1M1j+k2M2j+k3M3j+...+knMnj
Pn-1=k1M1n-1+k2M2n-1+k3M3n-1+...+knMnn-1
1=k1+k2+k3+...+kn
Pj=Fj/B
in the above formula, PjA third property data value, M, corresponding to the jth point in a third property data volume of a fracture of a known well1j、M2j、M3j、...、MnjRespectively corresponding to each second attribute data value of the j point in the second attribute data body; fjFracture development intensity data value, k, of the j point of a fracture of a known welliFor the ith second attribute data volume MiB is the crack of the third attribute data volumeAnd calculating a coefficient of development strength.
And step 40, predicting and evaluating the development strength and the direction of the fracture of the target interval by using the third attribute data body.
Example (c):
according to the flow of the method for calculating the fracture parameters based on the depth domain seismic data, the target interval fracture prediction is carried out on the marine shale stratum in a certain three-dimensional work area in south China, so that the related data such as the fracture development strength, the distribution position and the like are provided for the horizontal well track design of the subsequent oil and gas development.
The associated azimuthal prestack depth migration data volume is acquired in step 10. The method comprises the steps of firstly dividing a plurality of azimuth angle ranges of omni-directional angle gather data to obtain related azimuth angle gather data bodies, and respectively carrying out pre-stack depth migration processing on each azimuth angle gather data body by using related pre-stack depth migration processing technical flows and parameters to obtain each azimuth angle pre-stack depth migration data body (namely, post-stack seismic data of each depth domain). In actual operation, the fracture development strength and direction of the target intervals on the relevant known wells are counted, and the directions of the main fractures and the secondary fractures and fractures are determined to be in the north-east-south-west direction, the north-west direction and the like. And determining the azimuth angle range and the number related to the omnidirectional angle gather data according to the development directions of the related cracks and fractures. In actual operation, six azimuth angle ranges are set for the omni-directional angle (0-180 degrees) of a gather in a research area in equal parts, the ranges are respectively 0-30 degrees, 30-60 degrees, 60-90 degrees, 90-120 degrees, 120-150 degrees and 150-180 degrees, and the six azimuth angle gather data bodies are respectively subjected to related prestack depth migration processing according to the same prestack depth migration processing flow and parameters to obtain six azimuth angle prestack depth migration data bodies.
In step 20, the first attribute data volume of each azimuth is obtained after coherent volume calculation and edge signal enhancement processing related to crack prediction are performed on the seismic attribute data volume of each azimuth after band-pass filtering and weighting processing. The method comprises the steps of utilizing logging information, stratigraphic layering information and seismic reflection characteristics to conduct logging and depth domain seismic data calibration, determining the position of a target interval in each depth domain seismic section, determining the seismic reflection characteristics of the top layer and the bottom layer of the marine shale target interval, and automatically developing tracking interpretation of the two layers for a six-azimuth pre-stack depth migration data body, so that the position of the whole target interval on the pre-stack depth migration data body is obtained. In actual operation, the relevant horizon data is used to participate in the subsequent calculation steps.
In addition, the amplitude response characteristics, the tuning frequency data and the like of the seismic data volume of the related medium and small-scale fractures in the target interval in the well in the study area (the azimuth prestack depth migration data volume of the depth domain is converted into the time domain seismic data volume) are utilized to perform band-pass filtering processing on the six azimuth time domain seismic data volumes, and then weighting fusion processing is performed to obtain the related seismic data volume highlighting the response of the fractures of different scales. In actual operation, according to the response characteristics of the amplitude and frequency of the aboveground crack, setting the frequency range of band-pass filtering to be three frequency ranges such as 23-31, 34-41, 44-51 and the like, and respectively carrying out band-pass filtering processing on six time domain seismic data bodies by utilizing the three frequency ranges; setting three weighting coefficients (0.3, 0.35 and 0.35) for the three seismic data volumes of the same azimuth angle according to the test result of the relevant weighting fusion processing to carry out weighting fusion processing, so as to obtain the relevant time domain seismic data volumes; and then according to the related interval velocity-depth model, converting the time domain seismic data volume after the weighted fusion processing into a depth domain seismic data volume to participate in subsequent calculation.
And performing coherent body calculation on the target interval on the six depth domain seismic data bodies after the weighted fusion processing, and performing edge signal enhancement processing to obtain each first attribute data body. In actual operation, coherent body calculation (C) is adopted by using the horizon data of the top layer and the bottom layer of the relevant target interval2) Performing coherent body calculation on the six depth domain seismic data bodies about the sea phase target interval (the outside of the target interval does not participate in coherent body attribute calculation) by using an algorithm and the same calculation parameter to obtain six coherent attribute data bodies; then, the six coherent attribute data volumes are subjected to edge signal enhancement processing (volume attribute filtering processing), so that the data volume with the six coherent attribute data volumes is obtainedTo six first attribute data volumes.
In step 30, the first attribute data volume of each azimuth angle is calculated by using an ant body algorithm respectively to generate a second attribute data volume of each azimuth angle; and then carrying out weighted fusion processing on the second attribute data volume of each azimuth angle to obtain a third attribute data volume. In actual operation, calculating and testing the weighting coefficients of the six second attribute data volumes according to the known fracture data on the well; and setting six weighting coefficients for the six second attribute data volumes, and performing weighting fusion processing on the six second attribute data volumes to obtain a third attribute data volume. The third attribute data volume is used for carrying out comprehensive analysis of a cross-well section, a bedding-along attribute plane graph and the like, the coincidence rate with the actually measured crack development condition on the well is better, the coincidence rate can reach 83.1 percent, and the related geological requirements are met. In addition, the comparison and analysis of the conventional crack prediction results such as ant body calculation on the depth domain seismic data body of the all-directional angle, conversion of the conventional time domain ant body calculation result into the depth domain and the like are carried out, and the result is found to be superior to the conventional ant body calculation result on the whole, the goodness of fit with the crack development condition on the well is higher than that of the conventional ant body calculation result, and the conclusion is also proved by the subsequent horizontal well drilling result.
The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention, and various modifications and changes may be made by those skilled in the art. Any modification, equivalent replacement, or improvement made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims (10)

1. A method for calculating fracture parameters based on depth domain seismic data is characterized by comprising the following steps:
step 10, acquiring a related azimuth prestack depth migration data volume;
step 20, after band-pass filtering and weighting processing are carried out on the prestack depth migration data volumes of all the azimuth angles, coherent volume calculation and edge signal enhancement processing are carried out, and then a first attribute data volume of each azimuth angle is obtained;
step 30, calculating the first attribute data bodies of all the azimuth angles by using an ant body algorithm respectively to generate second attribute data bodies of all the azimuth angles; then, carrying out weighted fusion processing on the second attribute data volume of each azimuth angle to obtain a third attribute data volume;
and step 40, predicting and evaluating the development strength and the direction of the fracture of the target interval by using the third attribute data body.
2. The method of calculating fracture parameters based on depth domain seismic data as claimed in claim 1, wherein step 1 comprises the sub-steps of:
step 11, setting an azimuth angle range;
step 12, carrying out azimuth division on the omnidirectional angle gather data according to a set azimuth range to obtain n azimuth gather data volumes;
and step 13, processing the n azimuth gather data volumes by respectively adopting the same prestack depth migration processing technical process and parameters to obtain the related n azimuth prestack depth migration data volumes.
3. The method for calculating fracture parameters based on depth-domain seismic data as claimed in claim 2, wherein the method for dividing the azimuth angle of the omni-directional angle gather data according to the set azimuth angle range in step 12 comprises: based on the symmetry relationship, converting the azimuth gather data of 0-360 degrees into the azimuth gather data of 0-180 degrees; and dividing the azimuth gather data of 0-180 degrees according to a set azimuth range to obtain n azimuth gather data volumes.
4. The method of claim 2, wherein in step 13, each azimuth gather data volume is pre-processed with migration including denoising and/or resolution enhancement before the n azimuth gather data volumes are processed using the same pre-stack depth migration processing technique and parameters.
5. The method of claim 1, wherein the step 20 of band-pass filtering and weighting each azimuthal prestack depth migration data volume comprises the sub-steps of:
(1) calculating an azimuth prestack depth migration data volume through a layer velocity-depth model, and converting the azimuth prestack depth migration data volume into a time domain seismic data volume;
(2) according to the frequency band range of the target interval in the converted time domain seismic data volume, carrying out spectrum analysis on different crack development intensity areas in the known well of the time domain seismic data, and determining tuning frequency data sets corresponding to the related different crack development intensities;
(3) dividing tuning frequency data sets with different crack development strengths according to crack prediction precision, expert experience, geological knowledge and the like, and determining different m frequency range ranges; wherein m is more than or equal to 2;
(4) transforming the time domain seismic data volume to a frequency domain by using discrete Fourier transform, and performing band-pass filtering on the time domain seismic data volume transformed to the frequency domain by using related m frequency band ranges to obtain a series of related m time domain filtering data volumes;
(5) setting relevant weighting coefficients for a series of m time domain filtering data volumes to carry out weighting fusion processing to obtain a time domain weighting data volume;
(6) and (3) converting the time domain weighted data volume into a depth domain seismic data volume after the speed-depth model calculation.
(7) And (4) performing the steps (1) to (6) on each azimuth prestack depth migration data body by analogy, completing the band-pass filtering and weighting processing of each azimuth prestack depth migration data body, and obtaining the seismic attribute data body after the band-pass filtering and weighting processing.
6. The method of claim 5, wherein the frequency data set for setting the development strength and size of the fracture in step (3) is [ f [ ]1,f2,…,fn]Determining the maximum frequency value f from the frequency data setkMinimum frequency value fbAnd Δ d is a design tolerance, the calculation formula of the frequency range is as follows:
fget up=fb-△d
fStop block=fk+△d
△d=1/n(fk-fb)
In the formula (f)Get upIs the starting frequency of the frequency band range; f. ofStop blockIs the stop frequency of the frequency band range; n is the number of frequencies of the frequency data set participating in the calculation; f. ofGet upMinimum frequency value, f, greater than the entire seismic frequency bandStop blockLess than the maximum frequency value for the entire seismic band.
7. The method of claim 5 or 6, wherein there is an overlap between the frequency ranges set in step (3), but the overlap should be no greater than 0.25(Δ d)1+△d2);△d1、△d2Respectively, the design tolerances of the two frequency band ranges set.
8. The method for calculating fracture parameters based on depth-domain seismic data according to claim 1, wherein in the step 20, when coherent body calculation and edge signal enhancement processing are performed, coherent body calculation and edge signal enhancement processing are performed only on the target interval of each azimuth seismic attribute data volume after band-pass filtering and weighting processing by using the top layer and bottom layer data of the target interval.
9. The method for calculating fracture parameters based on depth-domain seismic data as claimed in claim 1, wherein the step 30 of calculating the first attribute data volume of each azimuth angle by using an ant body algorithm respectively, and the specific operation of generating the second attribute data volume of each azimuth angle is as follows: and after the layer leveling is carried out on the first attribute data bodies of all the azimuth angles by using the layer data of the bottom layer of the target layer section, ant body calculation is carried out on the target layer section, and then the ant body calculation result is subjected to reverse layer leveling to obtain all the second attribute data bodies.
10. The method of claim 1, wherein the step 30 of performing weighted fusion on the second attribute data volumes of each azimuth angle is as follows:
Figure FDA0003213440630000041
in the formula, P is a third attribute data volume after weighted fusion processing; miIs the ith second attribute data volume; k is a radical ofiFor the ith second attribute data volume MiThe set weighting coefficients of (a) are set,
Figure FDA0003213440630000042
n is the number of second attribute data volumes participating in the weighted fusion processing;
the determination of the weighting coefficient is determined according to the test condition of the weighted fusion result of the data values of the development intensity of the fractures in different directions in the related well, and specifically comprises the following steps:
for the determination of the weighting coefficients of the n second attribute data volumes, the determination is made according to the following calculation formula:
P1=k1M11+k2M21+k3M31+...+knMn1
P2=k1M12+k2M22+k3M32+...+knMn2
Pj=k1M1j+k2M2j+k3M3j+...+knMnj
Pn-1=k1M1n-1+k2M2n-1+k3M3n-1+...+knMnn-1
1=k1+k2+k3+...+kn
Pj=Fj/B
in the above formula, PjA third property data value, M, corresponding to the jth point in a third property data volume of a fracture of a known well1j、M2j、M3j、...、MnjRespectively corresponding to each second attribute data value of the j point in the second attribute data body; fjFracture development intensity data value, k, of the j point of a fracture of a known welliFor the ith second attribute data volume MiB is a crack growth intensity calculation coefficient of the third attribute data volume.
CN202110937680.XA 2021-08-16 2021-08-16 Method for calculating crack parameters based on depth domain seismic data Active CN113671566B (en)

Priority Applications (1)

Application Number Priority Date Filing Date Title
CN202110937680.XA CN113671566B (en) 2021-08-16 2021-08-16 Method for calculating crack parameters based on depth domain seismic data

Applications Claiming Priority (1)

Application Number Priority Date Filing Date Title
CN202110937680.XA CN113671566B (en) 2021-08-16 2021-08-16 Method for calculating crack parameters based on depth domain seismic data

Publications (2)

Publication Number Publication Date
CN113671566A true CN113671566A (en) 2021-11-19
CN113671566B CN113671566B (en) 2023-05-26

Family

ID=78543351

Family Applications (1)

Application Number Title Priority Date Filing Date
CN202110937680.XA Active CN113671566B (en) 2021-08-16 2021-08-16 Method for calculating crack parameters based on depth domain seismic data

Country Status (1)

Country Link
CN (1) CN113671566B (en)

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116736374A (en) * 2023-06-20 2023-09-12 中海石油(中国)有限公司深圳分公司 Multi-azimuth seismic data fusion method, device, equipment and medium

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050219949A1 (en) * 2004-03-30 2005-10-06 Taner M T Method for detecting earth formation fractures by seismic imaging of diffractors
CN104635269A (en) * 2013-11-13 2015-05-20 中国石油化工股份有限公司 Method for predicting igneous rock fractured reservoir on basis of prestack forward azimuth trace gather
CN105629303A (en) * 2015-12-28 2016-06-01 中国石油大学(北京) Prestack crack quantitative forecast method and system based on rock physics
CN111399049A (en) * 2020-04-29 2020-07-10 西南石油大学 Crack strength prediction method based on data volume dimensionality reduction and discrete coefficient calculation
CN111399056A (en) * 2020-04-29 2020-07-10 西南石油大学 Method for predicting crack strength based on divided azimuth filtering

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20050219949A1 (en) * 2004-03-30 2005-10-06 Taner M T Method for detecting earth formation fractures by seismic imaging of diffractors
CN104635269A (en) * 2013-11-13 2015-05-20 中国石油化工股份有限公司 Method for predicting igneous rock fractured reservoir on basis of prestack forward azimuth trace gather
CN105629303A (en) * 2015-12-28 2016-06-01 中国石油大学(北京) Prestack crack quantitative forecast method and system based on rock physics
CN111399049A (en) * 2020-04-29 2020-07-10 西南石油大学 Crack strength prediction method based on data volume dimensionality reduction and discrete coefficient calculation
CN111399056A (en) * 2020-04-29 2020-07-10 西南石油大学 Method for predicting crack strength based on divided azimuth filtering

Non-Patent Citations (3)

* Cited by examiner, † Cited by third party
Title
JING ZHANG: "Seismic Azimuthal Anisotropy Analysis Applied to Natural Fracture Intensity and Azimuth Prediction: Barnett Shale Example.", 2019 AAPG ANNUAL CONVENTION AND EXHIBITION *
汪忠德;薛诗桂;李红敬;裴云龙;李子锋;: "油气地球物理解释技术研究新进展", 地球物理学进展 *
简世凯: "窄方位三维叠前地震资料裂缝预测技术研究及应用", 中国优秀硕士学位论文全文数据库 *

Cited By (1)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN116736374A (en) * 2023-06-20 2023-09-12 中海石油(中国)有限公司深圳分公司 Multi-azimuth seismic data fusion method, device, equipment and medium

Also Published As

Publication number Publication date
CN113671566B (en) 2023-05-26

Similar Documents

Publication Publication Date Title
CN109425896B (en) Dolomite oil and gas reservoir distribution prediction method and device
US8209125B2 (en) Method for identifying and analyzing faults/fractures using reflected and diffracted waves
CN101329405B (en) Simple method of multi-parameter seismic inversion
CN110231652B (en) Density-based seismic facies extraction method using spatial clustering with noise
CN111506861B (en) Method for calculating crack strength of favorable region of target layer
CN111399056B (en) Method for predicting crack strength based on divided azimuth filtering
Zhong et al. Statistical analysis of background noise in seismic prospecting
CN103869362A (en) Method and equipment for obtaining body curvature
Jesus et al. Multiattribute framework analysis for the identification of carbonate mounds in the Brazilian presalt zone
Yuan et al. 6D phase-difference attributes for wide-azimuth seismic data interpretation
CN113109875B (en) Inversion method of carbonate rock reservoir under full waveform velocity field constraint
Lay The teleseismic manifestation of pP: Problems and paradoxes
CN111090117B (en) Effective reservoir prediction method and system under phase control forward constraint
CN113484907B (en) Method for predicting distribution on different types of reservoir planes
CN113671566B (en) Method for calculating crack parameters based on depth domain seismic data
Li et al. Seismic reflection characteristics of fluvial sand and shale interbedded layers
CN111665559A (en) Method and system for describing sliding fracture zone
CN115857047A (en) Comprehensive prediction method for seismic reservoir
Shi et al. Data processing workflow for reflections from cement-formation interface in ultrasonic pitch-catch measurement in cased hole
CN111505715A (en) Method for calculating crack parameters based on central incidence angle of depth domain
Zhang et al. Identifying minor faults on top of coalfield Ordovician limestone stratum using seismic attributes derived from azimuthally stacked data
CN111352154B (en) Reservoir prediction method based on wide-azimuth earthquake
CN110850504B (en) Shale density parameter prestack inversion method based on uranium curve quasi-impedance constraint
Li et al. Fault and fracture prediction of tight gas reservoir based on seismic likelihood attribute
Baharvand Ahmadi Analysis of time-lapse 3-D VSP data for seismic monitoring of CO2 flood in Weyburn Field, Saskatchewan

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