CN104166161A - Method and device for predicating fractures based on elliptical velocity inversion of anisotropism - Google Patents

Abstract

The present invention provides a fracture prediction method and device based on anisotropic ellipse velocity inversion, by rearranging wide-azimuth or multi-azimuth seismic data into azimuth gathers according to azimuth; Velocity spectrum is obtained by warp velocity analysis; the dynamic correction velocity curve is obtained by using the velocity spectrum automatic picking algorithm, and the target layer azimuth dynamic correction velocity is screened out; the dynamic correction velocity distribution is ellipse fitted, and the fracture development direction and Crack development strength. The present invention automatically picks up the velocity spectrum and η with high precision, greatly improves the efficiency of non-hyperbolic velocity analysis, and helps to reveal the spatial distribution law of fracture development, thereby enabling more effective detection of fractured oil and gas, greatly improving Success Rate of Fractured Reservoir Prediction.

Description

A Fracture Prediction Method and Device Based on Anisotropic Ellipse Velocity Inversion

technical field

The invention belongs to the field of oil and gas exploration, in particular to a fracture prediction method and device based on anisotropic ellipse velocity inversion.

Background technique

Effective fractures are the key to obtaining industrial productivity in tight sandstone gas reservoirs. In order to effectively identify fractures, fracture prediction technology is crucial, especially in the exploration of fractured reservoirs. Based on wide-azimuth or even multi-azimuth 3D seismic data, the distribution of faults and fractures can be predicted by pre-stack inversion, extracting anisotropic velocity by azimuth, and applying ellipse fitting technology. Studies have shown that: whether fractures are developed and how to identify fracture development zones in shale reservoirs, and then predict the distribution of high-yield and enriched zones of oil and gas are the key and difficult points for successful exploration of fractured reservoirs. Therefore, it is an important requirement in exploration and production to develop technologies to identify the spatial development and distribution of reservoir fractures.

Contents of the invention

The object of the present invention is to provide a fracture prediction method and device based on anisotropic ellipse velocity inversion, aiming at improving the accuracy and reliability of oil and gas exploration, providing a reliable basis for the fracture prediction work of interpreters, thus effectively Judging the distribution of high oil and gas enrichment zones.

The present invention is achieved in this way, a fracture prediction method based on anisotropic ellipse velocity inversion, comprising the following steps:

S1. Rearrange the wide-azimuth or multi-azimuth seismic data into an azimuth gather according to the azimuth;

S2. Obtain the velocity spectrum by non-hyperbolic velocity analysis for each azimuth gather;

S3. Use the velocity spectrum automatic picking algorithm to obtain the dynamic correction velocity curve, and screen out the target layer azimuth dynamic correction velocity;

S4. Ellipse fitting is carried out from the velocity distribution of dynamic correction, and the fracture development direction and fracture development intensity are obtained by least square method inversion.

Preferably, the step S2 includes the following specific steps:

(1) Perform conventional velocity analysis on the seismic data of the short array to obtain the initial dynamic correction value;

(2) Use the Akhalifa non-hyperbolic velocity analysis formula to substitute the initial dynamic correction velocity value into the original trace seismic data to scan to determine the value of η;

(3) Substituting η into the non-double speed analysis formula to redetermine the dynamic correction speed value to iterate repeatedly;

(4) Repeat steps (2), (3) until the η value and the dynamic correction velocity value are stable, then output the velocity spectrum and η spectrum.

The present invention further provides a fracture prediction device based on anisotropic ellipse velocity inversion, including:

The azimuth gather selection module is used to rearrange the wide-azimuth or multi-azimuth seismic data into azimuth gathers according to the azimuth;

The non-hyperbolic velocity analysis module is used to obtain the velocity spectrum by non-hyperbolic velocity analysis for each azimuth gather in the azimuth gather selection module;

The automatic picking module is used to obtain the dynamic correction velocity curve by using the velocity spectrum automatic picking algorithm on the velocity spectrum obtained by the non-hyperbolic velocity analysis module, and screen out the target layer azimuth dynamic correction velocity;

The ellipse fitting module is used for ellipse fitting of the dynamic correction velocity distribution in the automatic picking module, and the fracture development direction and fracture development intensity are obtained by least square method inversion.

Preferably, the non-hyperbolic velocity analysis module includes:

The conventional velocity analysis module is used to conduct conventional velocity analysis on the seismic data of the short array to obtain the initial dynamic correction value;

The scanning module is used to use the Akhalifa non-hyperbolic velocity analysis formula to substitute the initial dynamic correction velocity value into the original channel seismic data to scan and determine the η value;

The iterative module is used for substituting η into the non-double speed analysis formula to re-determine the dynamic correction speed value and iterate repeatedly;

The output module is used to output the velocity spectrum and the η spectrum after the repeated scanning module and the iterative module process until the η value and the dynamic correction velocity value are stable.

The present invention overcomes the deficiencies of the prior art, and provides a fracture prediction method and device based on anisotropic ellipse velocity inversion, which mainly uses pre-stack wide-azimuth or multi-azimuth seismic data for processing. By using the non-hyperbolic velocity analysis method, the omnidirectional or wide azimuth dynamic correction velocity information and anisotropy parameter η of the target layer are extracted, so as to invert the fracture development intensity and development direction according to the relationship between the dynamic correction velocity and azimuth angle . This method provides a reliable prediction algorithm for crack prediction, and provides a reliable basis for the crack prediction work of interpreters.

In the present invention, the ellipse velocity inversion anisotropic fracture prediction technology is a technology that is being vigorously developed and continuously improved at home and abroad. It detects the azimuth and density of fracture development by using the variation of seismic wave velocity with azimuth when it is relayed in anisotropic media, and the prediction results are more closely related to the microscopic characteristics of the fracture development zone. Therefore, it is more conducive to revealing the spatial distribution of fracture development, so that the fractured oil and gas can be detected more effectively, and the success rate of fractured reservoir prediction can be greatly improved.

The prior art relevant to the present invention includes:

1. Conventional speed analysis:

Conventional velocity analysis is a method to estimate the velocity distribution of the underground medium based on the Dix formula assumed by the layered isotropic geophysics of the subsurface medium through pre-stack seismic data. Velocity analysis is an important key link in seismic exploration, and the accuracy of velocity analysis directly affects the follow-up seismic data time-domain processing and interpretation process.

Conventional velocity analysis methods include the following methods:

(1) t ² -x ² method

According to the Dix formula, the time-distance equation of the reflected wave in the t ² -x ² plane is linear. Therefore, the linear equation of reflected wave offset and reflected wave travel time can be fitted, and the square root of the reciprocal slope of the equation obtained is the dynamic correction speed sought in the present invention.

(2) Speed scanning method

It is to use a series of dynamic correction speeds to perform dynamic correction on the CMP gathers, and obtain a corrected image every time the dynamic correction is performed, and evaluate the leveling degree. The scanning speed that can make the event of the reflected wave the largest leveling degree is the present invention. The desired dynamic correction speed.

(3) Velocity Spectrum Scanning Method

By setting a series of scanning dynamic correction velocities and a series of zero-offset two-way travel times, using these scanning dynamic correction velocities and zero-offset two-way travel times, the CMP (CDP) gathers are dynamically corrected and superimposed, and the calculation Evaluate parameters (average energy, similarity coefficient, correlation coefficient, etc.), and place the calculated evaluation parameters in the two-dimensional plane of scanning motion correction velocity and scanning zero offset two-way travel time, and the scanning result is called velocity spectrum. The dynamic correction speed and zero-offset round-trip travel time corresponding to the maximum value of the speed spectrum are the dynamic correction speed and zero-offset round-trip travel time calculated in the present invention.

(4) CVS method

This method is to take a small section of the survey line, use a series of dynamic correction speeds to carry out dynamic correction superposition on this section, different superimposition speeds, the superimposed image is called a CVS image, and the best dynamic correction speed is obtained from a series of CVS images For the desired result of the present invention.

2. Non-hyperbolic velocity analysis:

Non-hyperbolic velocity analysis is a seismic wave velocity inversion method based on an anisotropic geological model, and a non-hyperbolic moveout formula is used on the basis of conventional velocity analysis. In the traditional velocity analysis technique, the formation is assumed to be isotropic, and when the seismic wave propagates in the horizontal layered medium, the travel time curve of the reflected wave is a hyperbolic equation. The actual strata are mostly anisotropic. If the conventional velocity analysis is used to process the seismic data, the velocity information obtained will have a huge error, so that the dynamic correction cannot level the far distance, so that reliable and satisfactory results cannot be obtained. Based on the anisotropic medium model, the non-hyperbolic velocity analysis method is established, which is not only more in line with the actual geological conditions, but also can obtain more accurate velocities, thereby obtaining higher-precision seismic data.

3. VVAZ crack prediction technology:

VVAZ (Velocity Variation with AZimuth) technology is a fracture prediction technology based on the relationship between the dynamic velocity of the HTI medium (high-angle fracture medium model) and the azimuth angle to evaluate the fracture development direction and development intensity. Fracture spatial distribution information and development are very important for oil and gas exploration. VVAZ can be used to predict the direction and strength of fracture development, so as to obtain more detailed fracture spatial distribution information and development conditions, and help determine well locations.

Description of drawings

Fig. 1 is the flow chart of the steps of the fracture prediction method based on anisotropic elliptic velocity inversion of the present invention;

Fig. 2 is the flow chart of the steps of bispectral method non-hyperbolic velocity analysis in the embodiment of the present invention;

Fig. 3 is a flow chart of the steps of the velocity spectrum automatic picking algorithm in the embodiment of the present invention;

Fig. 4 is the change figure of speed with azimuth angle in the embodiment of the present invention;

Fig. 5 is a structural schematic diagram of a fracture prediction device based on anisotropic ellipse velocity inversion of the present invention;

Fig. 6 is a schematic structural diagram of a non-hyperbolic velocity analysis module in the fracture prediction device based on anisotropic elliptic velocity inversion of the present invention.

Detailed ways

In order to make the object, technical solution and advantages of the present invention clearer, the present invention will be further described in detail below in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described here are only used to explain the present invention, not to limit the present invention.

A fracture prediction method based on anisotropic ellipse velocity inversion, as shown in Fig. 1, includes the following steps:

S1. Rearrange wide-azimuth or multi-azimuth seismic data into azimuth gathers according to azimuth

In step S1, the azimuth gather is the basis of the azimuth non-hyperbolic velocity analysis, which directly affects the subsequent processing, so it is necessary to select and arrange the azimuth gathers for the super CMP gather. There are three selection methods:

(1) Optimized gun spacing method

According to the distribution range of running distances, the maximum vertical and horizontal offsets of the survey line are consistent, and at the same time, the area outside the "inscribed circle" is cut off, so that the wide gathers become omni-directional, which is more conducive to studying the influence of azimuth on velocity analysis .

(2) borrowing the law

The coverage times in the panel and the distribution of the shot detection pairs are more reasonable by means of the method, such as the uniformization of the panel and the processing of data rules.

(3) Variable angle method

According to the coverage times of seismic data and the distribution of offset distances, the size of the sectors should be adjusted reasonably based on the principle that the coverage times of each sector are basically equal.

Different data have different effects on the same method, so when selecting azimuth angle gathers for super CMP gathers, the method needs to be selected according to the actual situation of the data.

S2. Obtain the velocity spectrum by non-hyperbolic velocity analysis for each azimuth gather

In step S2, the Akhalifah non-hyperbolic time difference analysis formula is used as the basis of anisotropic non-hyperbolic velocity, such as formula (1):

${\mathrm{<span\; class="notranslate">\; t</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; x</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}\mathrm{<span\; class="notranslate">\; 0</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; +</span>\frac{{\mathrm{<span\; class="notranslate">\; x</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}}{{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; nmo</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</span>}<span\; class="notranslate">\; -</span>\frac{\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; \&eta;</span>}}_{\mathrm{<span\; class="notranslate">\; eff</span>}}{\mathrm{<span\; class="notranslate">\; x</span>}}^{\mathrm{<span\; class="notranslate">\; 4</span>}}\mathrm{}}{{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; nmo</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; [</span>{\mathrm{<span\; class="notranslate">\; t</span>}}_{\mathrm{<span\; class="notranslate">\; p</span>}\mathrm{<span\; class="notranslate">\; 0</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}{\mathrm{<span\; class="notranslate">\; V</span>}}_{\mathrm{<span\; class="notranslate">\; nmo</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 0</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; +</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 2</span>}{\mathrm{<span\; class="notranslate">\; \&eta;</span>}}_{\mathrm{<span\; class="notranslate">\; eff</span>}}<span\; class="notranslate">\; )</span>{\mathrm{<span\; class="notranslate">\; x</span>}}^{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ]</span>}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 1</span>}<span\; class="notranslate">\; )</span>$

The hyperbolic time difference analysis formula consists of two parts: the first part is the hyperbolic time difference composed of the first two terms, and the second part is the anisotropic non-hyperbolic time difference composed of the third term. Through theoretical analysis, it can be seen that this formula has a good processing effect on seismic data with large offsets (large arrays).

The non-hyperbolic time difference formula mainly has two parameters: one is V _nmo (dynamic correction speed), and the other is η _eff (equivalent anisotropy parameter η). The influence of these two parameters on the non-hyperbolic travel time is different in the offset distribution, and the influence of the dynamic correction velocity V _nmo on the non-hyperbolic travel time is global and plays a decisive role in the entire arrangement. The equivalent anisotropy parameter η _eff will have a more obvious impact on the non-hyperbolic travel time only in the case of large offsets, so a more accurate η _eff can only be obtained by analyzing the long-distance data value. Combined with the judgment criterion of the similarity coefficient of conventional velocity analysis, an anisotropic non-hyperbolic velocity analysis is formed.

On the basis of Akhalifa (1997), the present invention has designed the bispectrum algorithm of following non-hyperbolic velocity analysis, realizes based on anisotropic non-hyperbolic velocity analysis, bispectral method non-hyperbolic velocity analysis, as shown in Figure 2, The more specific implementation steps are as follows:

(1) Perform conventional velocity analysis on the seismic data of the short array to obtain the initial dynamic correction value

In step (1), at a certain t ₀ , the short array CMP gathers are used to perform conventional velocity analysis to obtain an initial dynamic correction velocity V _nmo .

(2) Use the Akhalifa non-hyperbolic velocity analysis formula to substitute the initial dynamic correction velocity value into the original trace seismic data to determine the η value

In step (2), the initial dynamic correction velocity and t ₀ obtained in step (1) are substituted into the non-hyperbolic moveout formula of Alkhalifah, and the large array seismic data is scanned to obtain an initial value of η.

(3) Substituting η into the non-double-velocity analysis formula to re-determine the dynamic correction velocity value and iterate repeatedly

In step (3), the initial η value is substituted into the non-hyperbolic moveout formula, and the large-array seismic data is scanned to re-determine the dynamic correction velocity V _nmo .

(4) Repeat steps (2), (3) until the η value and dynamic correction velocity value are stable, output the velocity spectrum and η spectrum

In step (4), step (2) and step (3) are repeated to iterate _Vnmo and η two to three times, and stable _Vnmo and η will be obtained; next _t0 is scanned, and the preceding three steps are performed, Until all t ₀ are scanned, the velocity spectrum and η spectrum are output.

The core of bispectral non-hyperbolic velocity analysis is to use short-arranged seismic data for conventional velocity analysis to determine an initial dynamic correction velocity value, and then use long-distance seismic data to analyze the time difference formula of non-hyperbolic velocity of Akhaliafa and pass The initial dynamic correction velocity value, determine the initial η value, and then iterate repeatedly to obtain a stable dynamic correction velocity spectrum and η spectrum.

S3. Use the velocity spectrum automatic picking algorithm to obtain the dynamic correction velocity curve, and screen out the target layer azimuth dynamic correction velocity

In step S3, the automatic picking of velocity spectrum generally refers to the process of extracting the most effective velocity spectrum energy group, which is a very important content in velocity analysis. This process can be done manually or automatically by computer. Manual picking has strong flexibility, but it has some shortcomings in estimating the parameters of the velocity spectrum energy group, not only the accuracy is not high, but also the efficiency is low and the cost is high. When the seismic section is long and the information of the velocity spectrum energy group is very rich, automatic technology should be integrated into the picking process, which can not only solve the above economic losses caused by the loss of effective geological information due to only partial picking, but also greatly Greatly improve work efficiency.

The process flow of velocity spectrum automatic picking algorithm (Viterbi algorithm) is shown in Figure 3, including:

(1) Carry out mean filtering process to velocity spectrum;

(2) Optimizing the integral calculation of the velocity spectrum;

(3) Automatically pick up the velocity spectrum.

As an optimal search algorithm for finding the shortest path, the Viterbi algorithm has been used in telecommunications to decipher convolutional codes. The main principle of the algorithm is derived from the observation that if the shortest path between a point and a path passes through an intermediate point, then this path segment between a point and a sum is also the shortest path between that path segment. This algorithm consists of two steps: forward shortest path accumulation (integration) calculation and backward descending tracing. In the application of automatic velocity picking, the present invention is called "forward maximum velocity spectrum energy group accumulation calculation and backward descending tracking".

The specific principle is:

Suppose the observed sequence is:

$\mathrm{<span\; class="notranslate">\; Z</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&Center\; Dot;</span><span\; class="notranslate">\; \&Center\; Dot;</span><span\; class="notranslate">\; \&Center\; Dot;</span><span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; ,</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 2</span>}<span\; class="notranslate">\; )</span>$

Then its joint probability distribution can always be expressed as

$\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span>\underset{\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 2</span>}}{\overset{\mathrm{<span\; class="notranslate">\; T</span>}}{\mathrm{<span\; class="notranslate">\; \&Pi;</span>}}}\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 3</span>}<span\; class="notranslate">\; )</span>$

The formula shows that the conditional distribution probability of the observed sequence variable zt at time t depends on all values before zt-1. In the formula, P(zt|zt-11) is the posterior probability, and its prior probability is P(zt-1). In order to effectively record the previous sequence with the largest prior probability, the present invention introduces a variable state sequence variable of a non-observation sequence, also known as a transfer sequence variable, which carries all the information before zt-1, and this information helps to describe the following The distribution of an observation sequence zt.

If the state qt at time t is a finite number in the range of {1,...,M}, and assuming that the processing process is only from 0 to time T, and the initial state and final state are known, then the state sequence can be represented by a finite vector representation, namely:

$\mathrm{<span\; class="notranslate">\; Q</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&Center\; Dot;</span><span\; class="notranslate">\; \&CenterDot;</span><span\; class="notranslate">\; \&CenterDot;</span><span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; \&CenterDot;</span>\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&CenterDot;</span><span\; class="notranslate">\; \&Center\; Dot;</span><span\; class="notranslate">\; \&CenterDot;</span><span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 4</span>}<span\; class="notranslate">\; )</span>$

It is used to represent the probability of the state qt=k at the time t=k based on all the states from time zero to time t=k-1 as conditions, then this processing process is called a first-order Markov process. That is to say, the state qt=k at the moment t=k gives all the states up to the moment t=k-1, and these states are uniquely dependent on the state at the previous moment, such as the state qt= at the moment t=k-1 k-1, that is:

P(q _k |q ₁ ,q ₂ ,…q _k-1 )＝P(q _k |q _k-1 ) (5)

Then, the T-order Markov process is:

$\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; k</span>}}<span\; class="notranslate">\; )</span>\underset{\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; k</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; T</span>}}{\mathrm{<span\; class="notranslate">\; \&Pi;</span>}}}\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; k</span>}}^{\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 6</span>}<span\; class="notranslate">\; )</span>$

According to Bayes' rule, the aforementioned observation sequence:

$\mathrm{<span\; class="notranslate">\; Z</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; =</span><span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 2</span>}}<span\; class="notranslate">\; ,</span><span\; class="notranslate">\; \&Center\; Dot;</span><span\; class="notranslate">\; \&Center\; Dot;</span><span\; class="notranslate">\; \&CenterDot;</span><span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 7</span>}<span\; class="notranslate">\; )</span>$

In short, when it is necessary to predict the observation sequence or the next state, the state variable qt=k contains the values of all related observation sequences and state variables in the past, that is, the sum of their past values. According to the assumption of conditional independent probability distribution, the relationship between the joint probability distribution of state variables and the observation sequence can be further simplified as:

$\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; )</span>\underset{\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 2</span>}}{\overset{\mathrm{<span\; class="notranslate">\; T</span>}}{\mathrm{<span\; class="notranslate">\; \&Pi;</span>}}}\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; +</span>\mathrm{<span\; class="notranslate">\; 1</span>}}<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; )</span>\underset{\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; =</span>\mathrm{<span\; class="notranslate">\; 1</span>}}{\overset{\mathrm{<span\; class="notranslate">\; T</span>}}{\mathrm{<span\; class="notranslate">\; \&Pi;</span>}}}\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; |</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 8</span>}<span\; class="notranslate">\; )</span>$

Then, to use a given observation sequence z ^T ₁ to deduce the corresponding maximum possible state sequence q ^T ₁ , it can be obtained by the following maximization algorithm.

The Viterbi algorithm finds a relevant and effective recursive solution to find the above maximum value, first define:

V(i,t)=P(z _t |q _t =i)max{P(q _t =i|q _t-1 =j)×V(j,t-1)} (9)

Its initialization value is:

V(i,1)=P(z ₁ |q _t =i)P(q ₁ =i) (10)

Thus obtaining the final sequence:

${\mathrm{<span\; class="notranslate">\; max</span>}}_{{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}}\mathrm{<span\; class="notranslate">\; P</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; z</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; ,</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; 1</span>}}^{\mathrm{<span\; class="notranslate">\; T</span>}}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; max</span>}}_{\mathrm{<span\; class="notranslate">\; i</span>}}\mathrm{<span\; class="notranslate">\; V</span>}<span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; i</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; T</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 11</span>}<span\; class="notranslate">\; )</span>$

If the maximum increment j*(i, t) of the above recursive formula is determined, then the optimal q*T1 can be obtained through a backward recursive calculation, that is, using:

${\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; -</span>\mathrm{<span\; class="notranslate">\; 1</span>}}^{<span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; =</span>{\mathrm{<span\; class="notranslate">\; j</span>}}^{<span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; (</span>{\mathrm{<span\; class="notranslate">\; q</span>}}_{\mathrm{<span\; class="notranslate">\; t</span>}}^{<span\; class="notranslate">\; *</span>}<span\; class="notranslate">\; ,</span>\mathrm{<span\; class="notranslate">\; t</span>}<span\; class="notranslate">\; )</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 12</span>}<span\; class="notranslate">\; )</span>$

Computing backward recursively is similar: the sum is simply replaced by the maximum.

S4. Carry out ellipse fitting based on dynamic correction velocity distribution, and obtain fracture development direction and fracture development intensity through least square method inversion

In step S4, the wide-azimuth or multi-azimuth seismic data are rearranged into azimuth gathers according to the azimuth angles, and the velocity spectrum is obtained by non-hyperbolic velocity analysis for each azimuth gather, and then the dynamic correction velocity is obtained by using the velocity spectrum automatic picking algorithm The dynamic correction velocity of the azimuth of the target layer is screened out, and the ellipse fitting is carried out according to the distribution of the dynamic correction velocity. Finally, the fracture development direction and fracture development intensity are obtained by least square method inversion.

If the anisotropy in the rock is caused by a set of oriented vertical fractures, then according to the theory of seismic wave propagation, longitudinal waves have different travel velocities as they propagate parallel or perpendicular to the fractures. When propagating parallel to the fracture, it propagates at a fast wave velocity; when propagating perpendicular to the fracture, it propagates at a slow wave velocity. As with the AVAZ principle, when a P-wave passes through a fractured medium, for a fixed offset, its azimuth velocity and fracture azimuth satisfy the following relationship:

V＝V ₀ +α·cos2β (13)

In the formula, V——zimuth velocity of longitudinal wave, m/s;

V ₀ —average value of azimuth velocity, m/s;

α—modulation factor related to azimuth velocity;

radian;

φ——observation azimuth from the excitation point to the receiver point, in radians;

——crack orientation, radian;

It can be seen from equation (13) that when the observation azimuth is parallel to the fracture strike, the velocity is the largest; as the angle between the observation azimuth and the fracture strike increases, the velocity gradually decreases, and when the angle is 90 degrees, the velocity After that, the speed gradually increases with the increase of the included angle, and reaches the maximum when the included angle is 180 degrees, and the change period is 180 degrees. The relationship between velocity and azimuth angle can be represented by an ellipse as shown in Figure 4. It can be seen from Figure 4 that the velocity along the fracture direction is V ₀ +α, and the velocity perpendicular to the fracture direction is V ₀ -α.

Theoretically, as long as the equation knows the velocities of three or more azimuths, the three parameters V ₀ , α, and β of the equation can be solved to obtain the azimuth-velocity ellipse equation. For wide-azimuth or omni-directional seismic data, assuming that offsets and azimuths are evenly distributed, usually at a given CDP position, there are multiple azimuths (generally greater than 3) of seismic observation data, then solving equation (13) is becomes an overdetermined problem. If it is defined that the direction from the true north is zero degrees, and the seismic data of each observation azimuth φ _i (i=1,2,…,N) are sorted clockwise, then the velocity of the corresponding azimuth angle is

For seismic data with N (N > 3) observation azimuths, the parameter value of equation (14) can be obtained by the least square fitting method at this time, and the variable e is defined as

For V ₀ , α and Find the partial derivatives, and make them equal to zero respectively, to get the following equations

Solving equations (16) gives

According to equations (17), (18) and (19), we can get α and the exact value of the V ₀ parameter, so as to obtain the anisotropic elliptic equation. According to the major and minor axes of the ellipse, the direction of crack development can be obtained.

On this basis, the fracture development strength can be further calculated, and through the above least squares algorithm, the α and V ₀ , _the ellipse is fitted by the ellipse fitting formula ( ₁₃ ), and then the fractures and development Intensity e:

$\mathrm{<span\; class="notranslate">\; e</span>}<span\; class="notranslate">\; =</span>\frac{{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; high</span>}}<span\; class="notranslate">\; -</span>{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; low</span>}}}{{\mathrm{<span\; class="notranslate">\; v</span>}}_{\mathrm{<span\; class="notranslate">\; high</span>}}}<span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; -</span><span\; class="notranslate">\; (</span>\mathrm{<span\; class="notranslate">\; 20</span>}<span\; class="notranslate">\; )</span>$

In the present invention, the development orientation and strength of fractures are detected by utilizing the velocity and azimuth angle when the seismic wave propagates in the anisotropic medium. It proposes the application of sub-azimuth selection and sorting of super CMP gathers, and VTI medium approximation processing for each azimuth gather; based on the non-hyperbolic velocity analysis, a Gaussian window is introduced to improve the resolution and analysis accuracy of the velocity spectrum. The viterbi optimization algorithm can automatically pick up the velocity spectrum and η with high precision, which greatly improves the efficiency of non-hyperbolic velocity analysis; in addition, a group of overdetermined elliptic equations can be obtained through the velocity and corresponding azimuth angles. The present invention utilizes least squares Finally, the ellipse is fitted to determine the development direction of the fracture, and the strength development formula is used to calculate the development strength of the fracture. The present invention introduces the model concept of ellipse fitting to predict the direction and intensity of fracture development, and forms a set of systematic fracture prediction means, which is beneficial to reveal the spatial distribution law of fracture development, so that the fractured oil and gas can be detected more effectively. The success rate of fracture reservoir prediction is greatly improved.

The present invention further provides a fracture prediction device based on anisotropic ellipse velocity inversion, as shown in Figure 5, comprising:

Azimuth gather selection and arrangement module 1, used to rearrange wide-azimuth or multi-azimuth seismic data into azimuth gathers according to azimuth;

The non-hyperbolic velocity analysis module 2 is used to obtain the velocity spectrum by non-hyperbolic velocity analysis for each azimuth gather in the azimuth gather selection module 1;

The automatic picking module 3 is used to use the velocity spectrum automatic picking algorithm obtained by the non-hyperbolic velocity analysis module 2 to obtain the dynamic correction velocity curve, and screen out the target layer azimuth dynamic correction velocity;

The ellipse fitting module 4 is used to perform ellipse fitting on the distribution of the dynamic correction velocity in the automatic picking module 3, and obtain the fracture development direction and fracture development intensity through least square method inversion.

More specifically, as shown in Figure 6, the non-hyperbolic velocity analysis module 2 includes:

The conventional velocity analysis module 21 is used to perform conventional velocity analysis on the seismic data of the short array to obtain an initial dynamic correction value;

The scanning module 22 is used for utilizing the Akhalifa non-hyperbolic velocity analysis formula to substitute the initial dynamic correction velocity value into the original track seismic data to scan and determine the η value;

Iteration module 23, is used for substituting n into the non-double speed analysis formula to re-determine the dynamic correction speed value and iterate repeatedly;

The output module 24 is used to output the velocity spectrum and the n-spectrum after the repeated scanning module and the iterative module process until the n value and the dynamic correction speed value are stable.

The device in the embodiment of the present invention corresponds to the method in the above-mentioned embodiment, and the device in this embodiment is also explained with the content described in the above-mentioned method embodiment, and details are not repeated here.

Compared with the shortcomings and deficiencies of the prior art, the present invention has the following beneficial effects:

(1) The present invention detects the development orientation and strength of fractures by utilizing the velocity of seismic waves propagating in the anisotropic medium along with the azimuth angle. The application is proposed to sort super CMP gathers by azimuth, and to do VTI medium approximation for each azimuth gather.

(2) Based on the non-hyperbolic velocity analysis, the Gaussian window is introduced to improve the resolution and analysis accuracy of the velocity spectrum. At the same time, the viterbi optimization algorithm is used to automatically pick up the velocity spectrum and η with high precision, which greatly improves the non-hyperbolic Warp velocity analysis efficiency.

(3) A set of overdetermined elliptic equations can be obtained through the velocity and corresponding azimuth angles. The present invention uses least squares to solve, and finally fits the ellipse to determine the development direction of the fracture, and uses the strength development formula to calculate the development intensity of the fracture. The present invention introduces the model concept of ellipse fitting to predict the direction and intensity of fracture development, and forms a set of systematic fracture prediction means, which is beneficial to reveal the spatial distribution law of fracture development, so that the fractured oil and gas can be detected more effectively. The success rate of fracture reservoir prediction is greatly improved.

The above descriptions are only preferred embodiments of the present invention, and are not intended to limit the present invention. Any modifications, equivalent replacements and improvements made within the spirit and principles of the present invention should be included in the protection of the present invention. within range.

Claims (4)

1. the crack prediction method based on anisotropic elliptical velocity inverting, is characterized in that, comprises the following steps:

S1, wide-azimuth or multi-faceted geological data are rearranged into orientation Dao Ji according to position angle;

S2, each road, orientation collection is obtained to velocity spectrum with non-double curve velocity analysis;

S3, utilization velocity spectrum automatic Picking algorithm obtain NMO velocity curve, filter out destination layer orientation NMO velocity;

S4, distribute and carry out ellipse fitting by NMO velocity, and obtain fracture azimuth and fracture development intensity by Least-squares inversion.

2. the crack prediction method based on anisotropic elliptical velocity inverting as claimed in claim 1, is characterized in that, described step S2 comprises following concrete steps:

(1) seismic data of short arrangement is carried out to conventional speeds analysis and obtain initial normal moveout correction value;

(2) utilize Akhalifa non-double curve velocity analysis formula that initial NMO velocity value substitution Dui Yuan road seismic data is scanned and determines η value;

(3) more non-η substitution double wave speeds analyse formula being redefined to NMO velocity value iterates;

(4) repeating step (2), (3) are to η value and NMO velocity value stabilization, and output speed spectrum and η compose.

3. the FRACTURE PREDICTION device based on anisotropic elliptical velocity inverting, is characterized in that, comprising:

Road, position angle collection gather module, for being rearranged into orientation Dao Ji by wide-azimuth or multi-faceted geological data according to position angle;

Non-double curve velocity analysis module, for obtaining velocity spectrum to each road, orientation collection of road, position angle collection gather module with non-double curve velocity analysis;

Automatic Picking module, uses velocity spectrum automatic Picking algorithm to obtain NMO velocity curve for the velocity spectrum that non-double curve velocity analysis module is obtained, and filters out destination layer orientation NMO velocity;

Ellipse fitting module, for the ellipse fitting that carries out that automatic Picking module NMO velocity is distributed, obtains fracture azimuth and fracture development intensity by Least-squares inversion.

4. the FRACTURE PREDICTION device based on anisotropic elliptical velocity inverting as claimed in claim 3, is characterized in that, described non-double curve velocity analysis module comprises:

Conventional speeds analysis module, obtains initial normal moveout correction value for the seismic data of short arrangement being carried out to conventional speeds analysis;

Scan module, for utilizing Akhalifa non-double curve velocity analysis formula that initial NMO velocity value substitution Dui Yuan road seismic data is scanned and determines η value;

Iteration module, iterates for non-η substitution double wave speeds analyse formula is redefined to NMO velocity value;

Output module, for processing to η value and NMO velocity value stabilization in multiple scanning module and iteration module, output speed spectrum and η spectrum.