CN103777241A - Three-dimensional seismic data quick edge detection method based on time domain generalized Hilbert conversion - Google Patents

Abstract

The invention discloses a three-dimensional seismic data quick edge detection method based on time domain generalized Hilbert conversion. The method includes the following steps that an objective interval of three-dimensional seismic data serves as a research object, the three-dimensional seismic data are input, a three-dimensional equitime slice calculation method is adopted, and a two-dimensional equitime slice of the objective interval is calculated; two one-dimensional Ehrlich's factors of a time domain are calculated to obtain an edge detection result of the current equitime slice; an ultimate edge detection result of the current slice is obtained; equitime slice numbers are replaced, and all the equitime slices are calculated; according to the positions of the equitime slices in the objective interval, edge detection results of the equitime slices are stored to obtained an ultimate three-dimensional edge detection result of the objective interval. The method has high edge detection accuracy, is high in calculation efficiency and guarantees that all sampling points are located on the respective equitime slices; underground geological conditions closest to reality are obtained, and meanwhile, restoration of the edge detection result is facilitated.

Description

Based on the three dimensional seismic data rapid edge-detection method of time domain broad sense Hilbert conversion

Technical field

The invention belongs to oil geophysical exploration field and a kind of three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion.

Background technology

In oil gas geophysical survey field, edge detection method is mainly used in finding tomography or Fractured Zone.Wherein the prediction of Fractured Zone is significant to the exploration and development of oil gas, particularly carbonate reservoir, and crack, hole and hole is wherein main migration pathway and the reservoir space of oil gas.

In natural rock or rock stratum, ubiquity hole, crack (gap) and corrosion hole, and it comes in every shape, and scale size also differs greatly.For seismic prospecting, owing to being subject to the restriction of resolution, None-identified goes out single hole, hole, seam, the Fractured Zone can only identification scale acquiring a certain degree.In fact, single hole, hole, seam are very little to the gathering role of oil gas, and what really have exploration, exploitation value is seam hole development belt of certain scale.

Fractured Zone refers to relative containment body rock stratum, and its seam hole density obviously increases, and has the rock mass of certain expanded range.Therefore, it is upper at seismic horizontal slice (or horizon slice), has certain distribution range and bearing of trend.Therefore, can adopt various edge detection method crack identification development belt.

1, Hilbert conversion

(1-1) frequency field Hilbert conversion

For 1 dimension discrete signal x (n), do Discrete Fourier Transform (DFT), obtain X (k), k=0,1 ..., N-1, wherein corresponding negative frequency, N is the sampling number of discrete signal x (t), then makes:

$H\left(k\right)=\left\{\begin{array}{cc}X\left(k\right)& k=0\\ 2X\left(k\right)& k=\mathrm{1,2},...,\frac{N}{2}-1\\ 0& k=\frac{\mathrm{<figure-callout\; id="2"\; label="N"\; filenames="HDA0000459904830000012.png,HDA0000459904830000022.png"\; state="\{\{state\}\}">N</figure-callout>}}{2},...,N-1\end{array}\right.---\left(1\right)$

H (k) is done to contrary DFT, and (this analytic signal is plural number, and its real part equals echo signal x (n) to obtain the analytic signal h (n) of x (n); Imaginary part is the result of the HT of x (n), is used to rim detection);

(1-2) time domain Hilbert conversion:

The Hilbert conversion of time domain shows as convolution relation;

$h\left(n\right)=x\left(n\right)+i\&CenterDot;\stackrel{\&OverBar;}{x}\left(n\right)---\left(2\right)$

$\stackrel{\&OverBar;}{x}\left(n\right)=x\left(n\right)*\stackrel{\&OverBar;}{h}=x\left(n\right)*\frac{1}{\mathrm{\&pi;n}}---\left(3\right)$

In formula, h (n) is the analytic signal of 1 dimension discrete signal x (n),

for the Xi Shi factor,

it is the result of the HT of echo signal x (n);

In theory, the Xi Shi factor

computing formula be:

$\stackrel{\&OverBar;}{h}\left(n\right)=\left\{\begin{array}{cc}0& n=2m\\ \frac{2}{\mathrm{\&pi;n}}& n=2m+1\end{array}\right.---\left(4\right)$

In formula, m is natural number;

2, broad sense Hilbert conversion:

LUOYi etc. (2003) have proposed broad sense Hilbert conversion (GHT) method, and it has introduced window function and exponent number, from two aspects, conventional HT is expanded, and for 1 dimension discrete signal x (t), its GHT can be expressed as follows;

h(t)＝h _r(t)+i·h _i(t) (5)

$h_r\left(t\right)={\{2\&CenterDot;{\mathrm{\&Sigma;}}_{\mathrm{\&omega;}}{\left\{\mathrm{Re}\right[X(t,\mathrm{\&omega;})\left]\right\}}^n+{\left\{\mathrm{Re}\right[X(t,0)\left]\right\}}^n\}}^{\frac{1}{n}}---\left(6\right)$

$h_i\left(t\right)={\{2\&CenterDot;{\mathrm{\&Sigma;}}_{\mathrm{\&omega;}}{\left\{\mathrm{Im}\right[X(t,\mathrm{\&omega;})\left]\right\}}^n\}}^{\frac{1}{n}}---\left(7\right)$

In formula, h (t) is the analytic signal of x (t); h _r(t) be the real number road (equaling echo signal x (t)) of h (t); h _i(t) be the imaginary number road (for rim detection) of h (t); X (t, ω) is the windowed Fourier transform of x (t); N is exponent number, works as n=1, and when box-shaped function that window function is endless, just change in quality for conventional HT in formula (5), (6), (7);

The technical scheme of the current rim detection for seismic data:

(1) Chen Xuehua, He Zhenhua, Huang Deji.The High-Order Pseudo Hilbert Transform rim detection of seismic data, in August, 2008, Advances in Geophysics, the 23rd volume the 4th gas, 1106-1110 page;

Technical scheme:

1. the object of processing: for three-dimensional time section (being horizontal time slice);

2. core calculations formula: broad sense Hilbert conversion (GHT) computing formula in proportion territory, adopts formula (5), (6), (7) to calculate;

3. calculation process: 2 dimension sections are considered as to 2 dimension groups; Again by 2 dimension groups according to laterally or being longitudinally decomposed into several 1 dimension groups; Adopt 1 dimension frequency field GHT to calculate respectively the GHT of each 1 dimension group, obtain analytic signal, and the imaginary part of extracting is for rim detection; The complete 1 all dimension group of cycle calculations, obtains the edge detection results array that 2 final dimensions are cut into slices.

(2) Xie Jing, He Zhenhua, Huang Deji, Xiong Xiaojun.Based on the high precision rim detection of windowing Hilbert conversion, in April, 2007, oil gas geophysics, the 5th the 1st phase of volume, 27-29 page

Technical scheme:

1. the object of processing: for three-dimensional time section (being horizontal time slice), 2 dimension groups;

2. core calculations formula: adopt the Hilbert transformation calculations formula of time domain, adopt formula (3), (4) to calculate;

3. calculation process: conventional time domain Hilbert transform method, in order to ensure computational accuracy, the Xi Shi factor need to be obtained very longly, but this can increase computing time and reduce resolving power in actual applications.If but that the factor obtains is too short, can produce this effect of jeep.And often get the shorter Xi Shi factor in order to improve counting yield in actual production, therefore, for overcoming this effect of jeep due to the shorter generation of the factor, adopt added-time window method.Thank and wait (2007) add Hamming window processing to the Xi Shi factor quietly.

Calculation process 1: determine that Xi Shi factor optimum length is 65, and the Xi Shi factor is added to Hamming window (2 multiply each other);

Calculation process 2: 2 dimension sections are considered as to 2 dimension groups; Again by 2 dimension groups according to laterally or being longitudinally divided into several 1 dimension groups; Adopt 1 dimension time domain HT(formula (2) and (3)) calculate respectively each 1 dimension group HT(for rim detection); The complete 1 all dimension group of cycle calculations, obtains the edge detection results array that 2 final dimensions are cut into slices.

(3) Xie Jing.Based on the earthquake rim detection research of time domain windowing Hilbert conversion, 2007, Chengdu University of Technology's Master's thesis

Technical scheme:

1. the object of processing: for three-dimensional time section (being horizontal time slice), 2 dimension groups;

2. core calculations formula: adopt the Hilbert transformation calculations formula of time domain, adopt formula (3), (4) to calculate;

3. calculation process: determine that Xi Shi factor optimum length is 65, and to the Xi Shi factor add Hamming window (2 multiply each other; 2 dimension sections are considered as to 2 dimension groups; Again by 2 dimension groups respectively according to horizontal and vertical several 1 dimension groups that is divided into; Adopt 1 dimension time domain HT(formula (2) and (3)) calculate respectively each 1 dimension group HT(for rim detection); The complete 1 all dimension group of cycle calculations, obtains the result of calculation of vertical and horizontal; Get again the root mean square of vertical and horizontal result of calculation, obtain the edge detection results array of 2 final dimension sections.

4. analyze: proposed the computing method of 2 directions, but still adopted conventional HT to calculate, computational accuracy still has much room for improvement.

(4) Xiong Xiaojun, He Zhenhua, Zhao Mingjin etc.A Crack Detection new method based on GHT, in August, 2009, geophysical prospecting for oil, the 44th the 4th phase of volume, 442-444 page

Technical scheme:

1. the object of processing: for three-dimensional time section (being horizontal time slice), 2 dimension groups;

2. core calculations formula: (a) the broad sense Hilbert in proportion territory conversion (GHT) computing formula, adopts formula (5), (6), (7) to calculate; (b) adopt 2 repair and maintenance limit denoising methods (LUOYi etc., 2002), carry out pre-service calculating.

3. calculation process: 2 dimension sections are considered as to 2 dimension groups, carry out 2 repair and maintenance limit denoisings; Again by 2 dimension groups according to laterally or being longitudinally divided into several 1 dimension groups; Adopt 1 dimension frequency field GHT to calculate respectively the GHT of each 1 dimension group, obtain analytic signal, and the imaginary part of extracting is for rim detection; The complete 1 all dimension group of cycle calculations, obtains the edge detection results array that 2 final dimensions are cut into slices.

The shortcoming of above-mentioned 4 kinds of technical schemes:

(1) 4 kind of scheme is all for three-dimensional time section (i.e. 2 dimension horizontal time slices), there is no the processing scheme for the actual formation (actual formation is non-level often) of three dimensional seismic data;

(2) 3 kinds of schemes (the 1st, 2,4 kinds) are all that 2 dimension sections are decomposed into several 1 dimension groups, only detect (horizontal or longitudinal) for 1 direction of 2 dimension sections, and the precision of its rim detection is lower;

The GHT method in (3) the 1st kinds of methods---the scheme proportion territory of Chen Xuehua etc. (2008) is calculated, and counting yield is relatively low, and suppresses noise ability poor (not having consideration in calculation process to suppress noise);

(4) the 2nd kinds of methods---thank to the scheme employing windowing HT method calculating of waiting (2007) quietly, counting yield is high, but suppresses noise ability (not having consideration in calculation process to suppress noise);

The scheme of (5) the 3rd kinds of methods---Xie Jing (2007) has been considered to suppress noise and the processing to 2 of two dimension slicing directions in calculation process, but still adopts conventional HT to calculate, and computational accuracy is low;

(6) the 4th kinds of methods---the scheme of Xiong Xiaojun etc. (2009) has been considered inhibition noise in calculation process, but still the GHT method in proportion territory is calculated, and counting yield is relatively low.

The existing edge detection method that solves 3D seismic data that lacks; The counting yield of rim detection is low, and computational accuracy is low.

Summary of the invention

The object of the embodiment of the present invention is to provide a kind of three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion, is intended to solve the existing edge detection method that solves 3D seismic data that lacks; The counting yield of rim detection is low, the problem that computational accuracy is low.

The embodiment of the present invention is achieved in that a kind of three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion, should the three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion comprise the following steps:

Take the objective interval of three dimensional seismic data as research object, input 3D seismic data, the computing method of section while adopting three-dimensional grade, section while calculating the two dimension etc. of objective interval;

Computing time 21 of territory dimensions the Xi Shi factors, section while grade for the 1st, carries out 2 repair and maintenance limit denoisings, obtains the array after denoising;

Choose length and be 65 the Xi Shi factor, the situation that order is 1, calculating and choosing length is 33 the Xi Shi factor, the situation that order is 2, until the edge detection results of section while obtaining current grade;

The result of calculation of comprehensive order 1 and order 2, obtains the final edge detection results of current slice;

Slice number while changing grade, section while calculating all grades; The position of cut into slices in objective interval when waiting, the result of the rim detection of cutting into slices while depositing each grade, obtains the three-dimensional edges testing result of final objective interval.

Further, should the three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion specifically comprise the following steps:

Step 1, take the objective interval of three dimensional seismic data as research object, input 3D seismic data, the interlayer position, top of target interval and interlayer position, the end, target interval is D, and the interlayer position, top of objective interval is Hor_1, and interlayer position, the end is Hor_2;

Step 2, the computing method of section while adopting three-dimensional grade, section while calculating the two dimension etc. of objective interval: 1. when interlayer position, the top of objective interval Hor_1 is parallel with interlayer position, the end Hor_2, the sum (N1) of cutting into slices when two dimension etc.:

N1＝L/dt （8）

In formula, L represents the time gap of Hor_2 and Hor_1, and dt represents the time sampling interval of 3D seismic data, when the top of objective interval interlayer position Hor_1 and interlayer position, end Hor_2 are when not parallel, and the number (N1) of section when two dimension waits:

N1＝Max_L/dt （9）

In formula, Max_L represents the maximum time distance of Hor_2 and Hor_1, dt represents the time sampling interval of 3D seismic data, Deng totally 15 of time sections, wherein the section of the 2nd, 3,4 deciles is in Inline or the non-sampling completely of Crossline direction, unified section mended to 0 processing, mend 0 net point of processing, do not participate in follow-up calculating, guarantee to wait the complete of time section, section is consistent in the length of vertical and horizontal when all grade;

Step 3, computing time 21 of territory dimensions the Xi Shi factors

with

its computing formula is identical,

$h_{01}^{-}\left(n\right)=h_{02}^{-}\left(n\right)\left\{\begin{array}{cc}0& n=2\mathrm{<figure-callout\; id="2"\; label="m"\; filenames="HDA0000459904830000012.png,HDA0000459904830000022.png"\; state="\{\{state\}\}">m</figure-callout>}\\ \frac{2}{\mathrm{\&pi;n}}& n=2m+1\end{array}\right.---\left(10\right)$

In formula, m is natural number, the Xi Shi factor

total length be set as 65, the Xi Shi factor

total length be set as 33; And the Xi Shi factor to time domain

with

apply Gaussian window w (n), the length of window array is 65;

${\stackrel{\&OverBar;}{h}}_1\left(n\right)={\stackrel{\&OverBar;}{h}}_{01}\left(n\right)\&times;w\left(n\right),n=1,...,65---\left(11\right)$

${\stackrel{\&OverBar;}{h}}_2\left(n\right)={\stackrel{\&OverBar;}{h}}_{02}\left(n\right)\&times;w\left(n\right),n=1,...,33---\left(12\right)$

Obtain 2 Xi Shi factors after windowing process

with

length is that 65 the Xi Shi factor is equivalent to the situation that the order in GHT is 1, i.e. conventional HT; Length is that 33 the Xi Shi factor is equivalent to the situation that the order in GHT is 2;

Step 4, section while grade for the 1st, be designated as array Slice (nx, ny), the value of nx is from 1 to Nx, Nx is the InLine sum of three dimensional seismic data, the value of ny is from 1 to Ny, and the CrossLine sum that Ny is three dimensional seismic data, carries out 2 repair and maintenance limit denoisings, obtain the array Slice_N (nx, ny) after denoising;

Step 5, chooses length and is 65 the Xi Shi factor

order is 1 situation, calculates; Concrete grammar is:

The first step, for 2 dimension group Slice_N (nx, ny) (nx=1 ..., Nx; Ny=1 ..., Ny), carry out cycle calculations in Inline direction:

Extract No. Inline=1 1 dimension group, i.e. Slice_N (1, ny), ny=1 ..., Ny; Making this 1 dimension group is Slice_X (ny), ny=1 ..., Ny;

Calculate Slice_X (ny) with convolution,

$\mathrm{Slice}\_X\_I\left(\mathrm{ny}\right)=\mathrm{Slice}\_X\left(\mathrm{ny}\right)*{\stackrel{-}{h}}_{01}\left(n\right)---\left(13\right)$

Slice_X_I in formula (ny) is the Hilbert conversion of Slice_X (ny), and the length of 1 dimension group Slice_X_I (ny) is Ny(1 dimension group Slice_X_I (ny) and 1 dimension group

the result array length of convolution is (Ny+65-1), each 32 points before and after result array is given up, Ny point in the middle of only choosing); And be array Slice_X_E (1, ny) ny=1 by 1 for array Slice_X_I (ny) writes 2 ..., Ny;

Change No. Inline, repeat the first step and second step, until all No. Inline calculated, result array Slice_X_E (nx, the ny) nx=1 of the rim detection of cutting into slices while obtaining current grade ..., Nx; Ny=1 ..., Ny;

Second step, for 2 dimension group Slice_N (nx, ny), carries out cycle calculations in Crossline direction:

Extract No. Crossline=1 1 dimension group, i.e. Slice_N (nx, 1), nx=1 ..., Nx; Making this 1 dimension group is Slice_Y (nx), nx=1 ..., Nx;

Calculate Slice_Y (nx) with convolution:

$\mathrm{Slice}\_Y\_I\left(\mathrm{ny}\right)=\mathrm{Slice}\_Y\left(\mathrm{ny}\right)*{\stackrel{-}{h}}_{01}\left(n\right)---\left(14\right)$

Slice_Y_I in formula (nx) is the Hilbert conversion of Slice_Y (nX), and the length of 1 dimension group Slice_Y_I (nx) is Nx(1 dimension group Slice_Y_I (nx) and 1 dimension group

the result array length of convolution is (Nx+65-1), by each 32 points before and after giving up, and Nx point in the middle of only choosing); And be array Slice_Y_E (nx, 1) nx=1 by 1 for array Slice_Y_I (ny) writes 2 ..., Nx;

Change No. Crossline, repeating step (6-1)～step (6-3), until all No. Crossline calculated, thus result array Slice_Y_E (nx, the ny) nx=1 of the rim detection of section while obtaining current grade ..., Nx; Ny=1 ..., Ny;

The 3rd step, the edge detection results of section while calculating current grade:

$\mathrm{Slice}\_E1(\mathrm{nx},\mathrm{ny})=\sqrt{\mathrm{Slice}\_X\_E^2(\mathrm{nx},\mathrm{ny})+\mathrm{Slice}\_Y\_E^2(\mathrm{nx},\mathrm{ny})},\mathrm{nx}=1,...,\mathrm{Nx};\mathrm{ny}=1,...,\mathrm{Ny}---\left(15\right)$

Step 6, chooses length and is 33 the Xi Shi factor

(situation that order is 2), repeats the first step to the four steps, until the edge detection results Slice_E2 (nx, ny) of section while obtaining current grade;

Step 7, the result of calculation of comprehensive order 1 and order 2, obtains the final edge detection results of current slice;

$\mathrm{Slice}\_E1(\mathrm{nx},\mathrm{ny})=\sqrt{\mathrm{Slice}\_{\mathrm{<figure-callout\; id="1"\; label="E"\; filenames="HDA0000459904830000012.png"\; state="\{\{state\}\}">E</figure-callout>}1}^2(\mathrm{nx},\mathrm{ny})+\mathrm{Slice}\_{\mathrm{<figure-callout\; id="2"\; label="E"\; filenames="HDA0000459904830000012.png,HDA0000459904830000022.png"\; state="\{\{state\}\}">E</figure-callout>}2}^2(\mathrm{nx},\mathrm{ny})},\mathrm{nx}=1,...,\mathrm{Nx};\mathrm{ny}=1,...,\mathrm{Ny}---\left(16\right)$

Step 8, changes slice number while grade, repeating step four～step 7, and while grade until all, section has been calculated;

Step 9, the position of cut into slices in objective interval when waiting, the result of the rim detection of cutting into slices while depositing each grade, obtains the three-dimensional edges testing result of final objective interval.

Further, in step 5, choose length and be 65 the Xi Shi factor

order is 1 situation, and the concrete grammar calculating is:

The first step, for 2 dimension group Slice_N (nx, ny) (nx=1 ..., Nx; Ny=1 ..., Ny), carry out cycle calculations in Inline direction:

Second step, for 2 dimension group Slice_N (nx, ny), carries out cycle calculations in Crossline direction:

The 3rd step, the edge detection results of section while calculating current grade:

$\mathrm{Slice}\_E1(\mathrm{nx},\mathrm{ny})=\sqrt{\mathrm{Slice}\_X\_E^2(\mathrm{nx},\mathrm{ny})+\mathrm{Slice}\_Y\_E^2(\mathrm{nx},\mathrm{ny})},\mathrm{nx}=1,...,\mathrm{Nx};\mathrm{ny}=1,...,\mathrm{Ny}---\left(15\right).$

Further, in the first step, extract No. Inline=1 1 dimension group, i.e. Slice_N (1, ny), ny=1 ..., Ny; Making this 1 dimension group is Slice_X (ny), ny=1 ..., Ny; Concrete grammar is:

Step 1, calculate Slice_X (ny) with

convolution,

$\mathrm{Slice}\_X\_I\left(\mathrm{ny}\right)=\mathrm{Slice}\_X\left(\mathrm{ny}\right)*{\stackrel{-}{h}}_{01}\left(n\right)---\left(13\right)$

Slice_X_I in formula (ny) is the Hilbert conversion of Slice_X (ny), and the length of 1 dimension group Slice_X_I (ny) is Ny(1 dimension group Slice_X_I (ny) and 1 dimension group

the result array length of convolution is (Ny+65-1), each 32 points before and after result array is given up, Ny point in the middle of only choosing); And be array Slice_X_E (1, ny) ny=1 by 1 for array Slice_X_I (ny) writes 2 ..., Ny;

Step 2, changes No. Inline, has calculated for all No. Inline, result array Slice_X_E (nx, the ny) nx=1 of the rim detection of section while obtaining current grade ..., Nx; Ny=1 ..., Ny.

Further, in second step, for 2 dimension group Slice_N (nx, ny), carry out cycle calculations concrete grammar in Crossline direction and be:

Step 1, extracts No. Crossline=1 1 dimension group, i.e. Slice_N (nx, 1), and nx=1 ..., Nx; Making this 1 dimension group is Slice_Y (nx), nx=1 ..., Nx;

Calculate Slice_Y (nx) with convolution:

$\mathrm{Slice}\_Y\_I\left(\mathrm{nx}\right)=\mathrm{Slice}\_Y\left(\mathrm{nx}\right)*h_{01}^{-}\left(n\right)---\left(14\right)$

Slice_Y_I in formula (nx) is the Hilbert conversion of Slice_Y (nX), and the length of 1 dimension group Slice_Y_I (nx) is Nx(1 dimension group Slice_Y_I (nx) and 1 dimension group

the result array length of convolution is (Nx+65-1), by each 32 points before and after giving up, and Nx point in the middle of only choosing); And be array Slice_Y_E (nx, 1) nx=1 by 1 for array Slice_Y_I (ny) writes 2 ..., Nx;

Step 2, changes repeating step (6-1)～step (6-3) No. Crossline; until all No. Crossline calculated, thus result array Slice_Y_E (nx, the ny) nx=1 of the rim detection of section while obtaining current grade; ..., Nx; Ny=1 ..., Ny.

Three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion provided by the invention, adopt 1 dimension time domain broad sense HT to process 3D seismic data, comprise windowing process, comprise again order processing, in calculation process, comprise squelch processing, and 2 directions of section are processed simultaneously during from 2 dimensions etc., have the precision of higher rim detection; Core calculations formula of the present invention is the HT of time domain, and counting yield is high; Adopt slicing treatment 3D seismic data while grade, and when adjacent etc., the interval of section equals 1 time sampling interval, guarantee to cut into slices when all sampling points are all positioned at grade separately upper, not only approached most real subsurface geology situation, and be conducive to the reduction of edge detection results; Adopt the Xi Shi factor of 2 different lengths to process, be equivalent to the situation of 2 orders in frequency field, and combine the result of calculation of 2 orders, computational accuracy is high.The present invention is directed to three dimensional seismic data and process, the three-dimensional edges that obtains target interval detects data volume, for the identification of Fractured Zone or seam hole development belt, effectively instructs the reservoir prediction of oil gas geophysical survey; The contact bed position, top of given 3-d seismic data set and objective interval and bottom boundary layer position, obtain three-dimensional rim detection data volume; Carry out the optimal combination of calculation process, consider squelch, counting yield and computational accuracy simultaneously.

Accompanying drawing explanation

Fig. 1 is the three dimensional seismic data rapid edge-detection method flow diagram based on time domain broad sense Hilbert conversion that the embodiment of the present invention provides;

Fig. 2 is the schematic diagram of the three-dimensional objective interval that provides of the embodiment of the present invention;

The calculating schematic diagram of section when Fig. 3 is waiting of the parallel formation that provides of the embodiment of the present invention;

The calculating schematic diagram of section when Fig. 4 is waiting of the non-parallel stratum that provides of the embodiment of the present invention.

Embodiment

In order to make object of the present invention, technical scheme and advantage clearer, below in conjunction with embodiment, the present invention is further elaborated.Should be appreciated that specific embodiment described herein, only in order to explain the present invention, is not intended to limit the present invention.

Below in conjunction with drawings and the specific embodiments, application principle of the present invention is further described.

As shown in Figure 1, the three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion of the embodiment of the present invention comprises the following steps:

S101: take the objective interval of three dimensional seismic data as research object, input 3D seismic data, the computing method of section while adopting three-dimensional grade, section while calculating the two dimension etc. of objective interval;

S102: computing time 21 of territory dimensions the Xi Shi factors, section while grade for the 1st, carries out 2 repair and maintenance limit denoisings, obtains the array after denoising;

S103: choose length and be 65 the Xi Shi factor (situation that order is 1), calculating and choosing length is 33 the Xi Shi factor (situation that order is 2), until the edge detection results of section while obtaining current grade;

S104: the result of calculation of comprehensive order 1 and order 2, obtains the final edge detection results of current slice;

S105: slice number while changing grade, section while calculating all grades; The position of cut into slices in objective interval when waiting, the result of the rim detection of cutting into slices while depositing each grade, obtains the three-dimensional edges testing result of final objective interval.

Concrete steps of the present invention comprise:

Step 1, take the objective interval of three dimensional seismic data as research object, (X-axis is Inline survey line to input 3D seismic data, Y-axis is CrossLine survey line, Z axis is the time), the interlayer position, top of target interval and interlayer position, the end, schematic diagram is shown in that (the target interval in Fig. 2 is D to Fig. 2, the interlayer position, top of objective interval is Hor_1, and interlayer position, the end is Hor_2);

Step 2, the computing method (ZengHongliu, 2010) of section while adopting three-dimensional grade, section while calculating the two dimension etc. of objective interval: 1. when interlayer position, the top of objective interval Hor_1 is parallel with interlayer position, the end Hor_2, the sum (N1) of cutting into slices when two dimension etc.:

N1＝L/dt （8）

In formula, L represents the time gap (unit: ms) of Hor_2 and Hor_1, dt represents the time sampling interval (unit: ms) of 3D seismic data, see Fig. 3 (1 time-sampling point in the grid representative longitudinally in Fig. 3, transversely 1 No. Inline or No. Crossline) Deng the calculating schematic diagram of time section; When the top of objective interval interlayer position Hor_1 and interlayer position, end Hor_2 are when not parallel, the number (N1) of section when two dimension waits:

N1＝Max_L/dt （9）

In formula, Max_L represents the maximum time distance (unit: ms) of Hor_2 and Hor_1, dt represents the time sampling interval (unit: ms) of 3D seismic data, see Fig. 4 (1 time-sampling point in the grid representative longitudinally in Fig. 4, transversely 1 No. Inline or No. Crossline Deng the calculating schematic diagram of time section, in Fig. 4, the maximum time of top interlayer position and interlayer position, the end is apart from being 16 times employing intervals, while reaching adjacent grade by 4 deciles, the time interval of section is all 1 time to adopt interval), in Fig. 4, wait totally 15 of time sections, wherein the 2nd, 3, the section of 4 deciles is in Inline or the non-sampling (only have on part grid node and have value) completely of Crossline direction, in actual computation, unification is mended 0 processing by these sections and (is mended 0 net point of processing, do not participate in follow-up calculating), while guaranteeing to wait complete (when all grade, section is consistent in the length of vertical and horizontal) of section,

Step 3, computing time 21 of territory dimensions the Xi Shi factors

with its computing formula is identical,

$h_{01}^{-}\left(n\right)=h_{02}^{-}\left(n\right)\left\{\begin{array}{cc}0& n=2\mathrm{<figure-callout\; id="2"\; label="m"\; filenames="HDA0000459904830000012.png,HDA0000459904830000022.png"\; state="\{\{state\}\}">m</figure-callout>}\\ \frac{2}{\mathrm{\&pi;n}}& n=2m+1\end{array}\right.---\left(10\right)$

In formula, m is natural number, the Xi Shi factor

total length be set as 65, the Xi Shi factor

total length be set as 33; And the Xi Shi factor to time domain

with

apply Gaussian window w (n) (length of window array is 65);

${\stackrel{\&OverBar;}{h}}_1\left(n\right)={\stackrel{\&OverBar;}{h}}_{01}\left(n\right)\&times;w\left(n\right),n=1,...,65---\left(11\right)$

${\stackrel{\&OverBar;}{h}}_2\left(n\right)={\stackrel{\&OverBar;}{h}}_{02}\left(n\right)\&times;w\left(n\right),n=1,...,33---\left(12\right)$

Obtain 2 Xi Shi factors after windowing process

with

(the Xi Shi factor that length is 65 is equivalent to the situation that the order in GHT is 1, i.e. conventional HT; Length is that 33 the Xi Shi factor is equivalent to the situation that the order in GHT is 2)

Step 4, while grade for the 1st, section (is designated as array Slice (nx, ny), the value of nx is from 1 to Nx, the InLine sum that Nx is three dimensional seismic data, the value of ny is from 1 to Ny, Ny is the CrossLine sum of three dimensional seismic data), carry out 2 repair and maintenance limit denoisings (LUOYi etc., 2002), obtain the array Slice_N (nx, ny) after denoising;

Step 5, chooses length and is 65 the Xi Shi factor

(situation that order is 1), calculates; Concrete grammar is:

The first step, for 2 dimension group Slice_N (nx, ny) (nx=1 ..., Nx; Ny=1 ..., Ny), carry out cycle calculations in Inline direction:

Extract No. Inline=1 1 dimension group, i.e. Slice_N (1, ny), ny=1 ..., Ny; Making this 1 dimension group is Slice_X (ny), ny=1 ..., Ny;

Calculate Slice_X (ny) with

convolution,

$\mathrm{Slice}\_X\_I\left(\mathrm{ny}\right)=\mathrm{Slice}\_X\left(\mathrm{ny}\right)*{\stackrel{-}{h}}_{01}\left(n\right)---\left(13\right)$

Slice_X_I in formula (ny) is the Hilbert conversion of Slice_X (ny), and the length of 1 dimension group Slice_X_I (ny) is Ny(1 dimension group Slice_X_I (ny) and 1 dimension group

the result array length of convolution is (Ny+65-1), each 32 points before and after result array is given up, Ny point in the middle of only choosing); And be array Slice_X_E (1, ny) ny=1 by 1 for array Slice_X_I (ny) writes 2 ..., Ny;

Change No. Inline, repeat the first step and second step, until all No. Inline calculated, result array Slice_X_E (nx, the ny) nx=1 of the rim detection of cutting into slices while obtaining current grade ..., Nx; Ny=1 ..., Ny;

Second step, for 2 dimension group Slice_N (nx, ny), carries out cycle calculations in Crossline direction:

Extract No. Crossline=1 1 dimension group, i.e. Slice_N (nx, 1), nx=1 ..., Nx; Making this 1 dimension group is Slice_Y (nx), nx=1 ..., Nx;

Calculate Slice_Y (nx) with

convolution:

$\mathrm{Slice}\_Y\_I\left(\mathrm{nx}\right)=\mathrm{Slice}\_Y\left(\mathrm{nx}\right)*h_{01}^{-}\left(n\right)---\left(14\right)$

Slice_Y_I in formula (nx) is the Hilbert conversion of Slice_Y (nX), and the length of 1 dimension group Slice_Y_I (nx) is Nx(1 dimension group Slice_Y_I (nx) and 1 dimension group

the result array length of convolution is (Nx+65-1), by each 32 points before and after giving up, and Nx point in the middle of only choosing); And be array Slice_Y_E (nx, 1) nx=1 by 1 for array Slice_Y_I (ny) writes 2 ..., Nx;

Change No. Crossline; repeating step (6-1)～step (6-3); content (seeing the reference position of annotations and comments) under step 2 until all No. Crossline calculated; thereby the result array Slice_Y_E (nx of the rim detection of section while obtaining current grade; ny) nx=1; ..., Nx; Ny=1 ..., Ny;

The 3rd step, the edge detection results of section while calculating current grade:

$\mathrm{Slice}\_E1(\mathrm{nx},\mathrm{ny})=\sqrt{\mathrm{Slice}\_X\_E^2(\mathrm{nx},\mathrm{ny})+\mathrm{Slice}\_Y\_E^2(\mathrm{nx},\mathrm{ny})},\mathrm{nx}=1,...,\mathrm{Nx};\mathrm{ny}=1,...,\mathrm{Ny}---\left(15\right)$

Step 6, chooses length and is 33 the Xi Shi factor

(situation that order is 2), repeats the first step to the four steps, until the edge detection results Slice_E2 (nx, ny) of section while obtaining current grade;

Step 7, the result of calculation of comprehensive order 1 and order 2, obtains the final edge detection results of current slice;

$\mathrm{Slice}\_E1(\mathrm{nx},\mathrm{ny})=\sqrt{\mathrm{Slice}\_{\mathrm{<figure-callout\; id="1"\; label="E"\; filenames="HDA0000459904830000012.png"\; state="\{\{state\}\}">E</figure-callout>}1}^2(\mathrm{nx},\mathrm{ny})+\mathrm{Slice}\_{\mathrm{<figure-callout\; id="2"\; label="E"\; filenames="HDA0000459904830000012.png,HDA0000459904830000022.png"\; state="\{\{state\}\}">E</figure-callout>}2}^2(\mathrm{nx},\mathrm{ny})},\mathrm{nx}=1,...,\mathrm{Nx};\mathrm{ny}=1,...,\mathrm{Ny};---\left(16\right)$

Step 8, changes slice number while grade, repeating step four～step 7, and while grade until all, section has been calculated;

Step 9, the position of cut into slices in objective interval when waiting, the result of the rim detection of cutting into slices while depositing each grade, obtains the three-dimensional edges testing result of final objective interval.

The present invention adopts 1 dimension time domain broad sense HT to process 3D seismic data, and GHT is similar with frequency field, has both comprised windowing process, has comprised again order processing; In addition, in calculation process, comprised squelch processing, and 2 directions of section are processed simultaneously during from 2 dimensions etc., have the precision of higher rim detection; In addition, core calculations formula of the present invention is the HT of time domain, and counting yield is high.Adopt slicing treatment 3D seismic data while grade, and the interval of cutting into slices when adjacent grade equals 1 time sampling interval.The present invention has guaranteed to cut into slices when all sampling points are all positioned at grade separately upper, not only approaches most real subsurface geology situation, and is conducive to the reduction (reconstituting 3-D data volume) of edge detection results; Adopt the Xi Shi factor of 2 different lengths to process, be equivalent to the situation of 2 orders in frequency field, and combine the result of calculation of 2 orders, computational accuracy is high.

The foregoing is only preferred embodiment of the present invention, not in order to limit the present invention, all any modifications of doing within the spirit and principles in the present invention, be equal to and replace and improvement etc., within all should being included in protection scope of the present invention.

Claims (5)

1. the three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion, is characterized in that, should the three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion comprise the following steps:

Take the objective interval of three dimensional seismic data as research object, input 3D seismic data, the computing method of section while adopting three-dimensional grade, section while calculating the two dimension etc. of objective interval;

Computing time 21 of territory dimensions the Xi Shi factors, length is respectively 65 and 33, section while grade for the 1st, carries out 2 repair and maintenance limit denoisings, obtains 2 dimension groups after denoising;

Choose length and be 65 the Xi Shi factor, the situation that order is 1, carries out Hilbert transformation calculations; Choose length and be 33 the Xi Shi factor, the situation that order is 2, carries out Hilbert transformation calculations, until the edge detection results of section while obtaining current grade;

The result of calculation of comprehensive order 1 and order 2, obtains the final edge detection results of current slice;

Slice number while changing grade, section while calculating all grades; The position of cut into slices in objective interval when waiting, the result of the rim detection of cutting into slices while depositing each grade, obtains the three-dimensional edges testing result of final objective interval.

2. the three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion as claimed in claim 1, it is characterized in that, should the three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion specifically comprise the following steps:

Step 1, take the objective interval of three dimensional seismic data as research object, input 3D seismic data, the interlayer position, top of target interval and interlayer position, the end, target interval is D, and the interlayer position, top of objective interval is Hor_1, and interlayer position, the end is Hor_2;

Step 2, the computing method of section while adopting three-dimensional grade, section while calculating the two dimension etc. of objective interval: 1. when interlayer position, the top of objective interval Hor_1 is parallel with interlayer position, the end Hor_2, the sum (N1) of cutting into slices when two dimension etc.:

N1＝L/dt

In formula, L represents the time gap of Hor_2 and Hor_1, and dt represents the time sampling interval of 3D seismic data, when the top of objective interval interlayer position Hor_1 and interlayer position, end Hor_2 are when not parallel, and the number (N1) of section when two dimension waits:

N1＝Max_L/dt

In formula, Max_L represents the maximum time distance of Hor_2 and Hor_1, dt represents the time sampling interval of 3D seismic data, Deng totally 15 of time sections, wherein the section of the 2nd, 3,4 deciles is in Inline or the non-sampling completely of Crossline direction, unified section mended to 0 processing, mend 0 net point of processing, do not participate in follow-up calculating, guarantee to wait the complete of time section, section is consistent in the length of vertical and horizontal when all grade;

Step 3, computing time 21 of territory dimensions the Xi Shi factors

with

its computing formula is identical,

$h_{01}^{-}\left(n\right)=h_{02}^{-}\left(n\right)\left\{\begin{array}{cc}0& n=2m\\ \frac{2}{\mathrm{\&pi;n}}& n=2m+1\end{array}\right.$

In formula, m is natural number, the Xi Shi factor

total length be set as 65, the Xi Shi factor total length be set as 33; And the Xi Shi factor to time domain

with

apply Gaussian window w (n), the length of window array is 65;

${\stackrel{\&OverBar;}{h}}_1\left(n\right)={\stackrel{\&OverBar;}{h}}_{01}\left(n\right)\&times;w\left(n\right),n=1,...,65$

${\stackrel{\&OverBar;}{h}}_2\left(n\right)={\stackrel{\&OverBar;}{h}}_{02}\left(n\right)\&times;w\left(n\right),n=1,...,33$

Obtain 2 Xi Shi factors after windowing process

with

length is that 65 the Xi Shi factor is equivalent to the situation that the order in GHT is 1, i.e. conventional HT; Length is that 33 the Xi Shi factor is equivalent to the situation that the order in GHT is 2;

Step 4, section while grade for the 1st, be designated as array Slice (nx, ny), the value of nx is from 1 to Nx, Nx is the InLine sum of three dimensional seismic data, the value of ny is from 1 to Ny, and the CrossLine sum that Ny is three dimensional seismic data, carries out 2 repair and maintenance limit denoisings, obtain the array Slice_N (nx, ny) after denoising;

Step 5, chooses length and is 65 the Xi Shi factor

order is 1 situation, calculates; Concrete grammar is:

The first step, for 2 dimension group Slice_N (nx, ny) (nx=1 ..., Nx; Ny=1 ..., Ny), carry out cycle calculations in Inline direction:

Extract No. Inline=1 1 dimension group, i.e. Slice_N (1, ny), ny=1 ..., Ny; Making this 1 dimension group is Slice_X (ny), ny=1 ..., Ny;

Calculate Slice_X (ny) with

convolution,

$\mathrm{Slice}\_X\_I\left(\mathrm{ny}\right)=\mathrm{Slice}\_X\left(\mathrm{ny}\right)*{\stackrel{-}{h}}_{01}\left(n\right)$

Slice_X_I in formula (ny) is the Hilbert conversion of Slice_X (ny), and the length of 1 dimension group Slice_X_I (ny) is Ny(1 dimension group Slice_X_I (ny) and 1 dimension group

the result array length of convolution is (Ny+65-1), each 32 points before and after result array is given up, Ny point in the middle of only choosing); And be array Slice_X_E (1, ny) ny=1 by 1 for array Slice_X_I (ny) writes 2 ..., Ny;

Change No. Inline, repeat the first step and second step, until all No. Inline calculated, result array Slice_X_E (nx, the ny) nx=1 of the rim detection of cutting into slices while obtaining current grade ..., Nx; Ny=1 ..., Ny;

Second step, for 2 dimension group Slice_N (nx, ny), carries out cycle calculations in Crossline direction:

Extract No. Crossline=1 1 dimension group, i.e. Slice_N (nx, 1), nx=1 ..., Nx; Making this 1 dimension group is Slice_Y (nx), nx=1 ..., Nx;

Calculate Slice_Y (nx) with

convolution:

$\mathrm{Slice}\_Y\_I\left(\mathrm{ny}\right)=\mathrm{Slice}\_Y\left(\mathrm{ny}\right)*h_{01}^{-}\left(n\right)$

Slice_Y_I in formula (nx) is the Hilbert conversion of Slice_Y (nX), and the length of 1 dimension group Slice_Y_I (nx) is Nx(1 dimension group Slice_Y_I (nx) and 1 dimension group

the result array length of convolution is (Nx+65-1), by each 32 points before and after giving up, and Nx point in the middle of only choosing); And be array Slice_Y_E (nx, 1) nx=1 by 1 for array Slice_Y_I (ny) writes 2 ..., Nx;

Change No. Crossline, repeating step (6-1)～step (6-3), until all No. Crossline calculated, thus result array Slice_Y_E (nx, the ny) nx=1 of the rim detection of section while obtaining current grade ..., Nx; Ny=1 ..., Ny;

The 3rd step, the edge detection results of section while calculating current grade:

$\mathrm{Slice}\_E1(\mathrm{nx},\mathrm{ny})=\sqrt{\mathrm{Slice}\_X\_E^2(\mathrm{nx},\mathrm{ny})+\mathrm{Slice}\_Y\_E^2(\mathrm{nx},\mathrm{ny})},\mathrm{nx}=1,...,\mathrm{Nx};\mathrm{ny}=1,...,\mathrm{Ny}$

Step 6, chooses length and is 33 the Xi Shi factor order is 2 situation, repeats the first step to the four steps, until the edge detection results Slice_E2 (nx, ny) of section while obtaining current grade;

Step 7, the result of calculation of comprehensive order 1 and order 2, obtains the final edge detection results of current slice:

$\mathrm{Slice}\_E1(\mathrm{nx},\mathrm{ny})=\sqrt{\mathrm{Slice}\_{E1}^2(\mathrm{nx},\mathrm{ny})+\mathrm{Slice}\_{E2}^2(\mathrm{nx},\mathrm{ny})},\mathrm{nx}=1,...,\mathrm{Nx};\mathrm{ny}=1,...,\mathrm{Ny};$

Step 8, changes slice number while grade, repeating step four～step 7, and while grade until all, section has been calculated;

Step 9, the position of cut into slices in objective interval when waiting, the result of the rim detection of cutting into slices while depositing each grade, obtains the three-dimensional edges testing result of final objective interval.

3. the three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion as claimed in claim 1, is characterized in that, in step 5, chooses length and be 65 the Xi Shi factor order is 1 situation, and the concrete grammar calculating is:

The first step, for 2 dimension group Slice_N (nx, ny) (nx=1 ..., Nx; Ny=1 ..., Ny), carry out cycle calculations in Inline direction;

Second step, for 2 dimension group Slice_N (nx, ny), carries out cycle calculations in Crossline direction;

The 3rd step, the edge detection results of section while calculating current grade:

$\mathrm{Slice}\_E1(\mathrm{nx},\mathrm{ny})=\sqrt{\mathrm{Slice}\_X\_E^2(\mathrm{nx},\mathrm{ny})+\mathrm{Slice}\_Y\_E^2(\mathrm{nx},\mathrm{ny})},\mathrm{nx}=1,...,\mathrm{Nx};\mathrm{ny}=1,...,\mathrm{Ny}.$

4. the three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion as claimed in claim 3, is characterized in that, in the first step, extract No. Inline=1 1 dimension group, i.e. Slice_N (1, ny), ny=1 ..., Ny; Making this 1 dimension group is Slice_X (ny), ny=1 ..., Ny; Concrete grammar is:

Step 1, calculate Slice_X (ny) with

convolution,

$\mathrm{Slice}\_X\_I\left(\mathrm{ny}\right)=\mathrm{Slice}\_X\left(\mathrm{ny}\right)*{\stackrel{-}{h}}_{01}\left(n\right)---\left(13\right)$

Slice_X_I in formula (ny) is the Hilbert conversion of Slice_X (ny), and the length of 1 dimension group Slice_X_I (ny) is Ny(1 dimension group Slice_X_I (ny) and 1 dimension group

the result array length of convolution is (Ny+65-1), each 32 points before and after result array is given up, Ny point in the middle of only choosing); And be array Slice_X_E (1, ny) ny=1 by 1 for array Slice_X_I (ny) writes 2 ..., Ny;

Step 2, changes No. Inline, has calculated for all No. Inline, result array Slice_X_E (nx, the ny) nx=1 of the rim detection of section while obtaining current grade ..., Nx; Ny=1 ..., Ny.

5. the three dimensional seismic data rapid edge-detection method based on time domain broad sense Hilbert conversion as claimed in claim 3, it is characterized in that, in second step, for 2 dimension group Slice_N (nx, ny), carrying out cycle calculations concrete grammar in Crossline direction is:

Step 1, extracts No. Crossline=1 1 dimension group, i.e. Slice_N (nx, 1), and nx=1 ..., Nx; Making this 1 dimension group is Slice_Y (nx), nx=1 ..., Nx;

Calculate Slice_Y (nx) with

convolution:

$\mathrm{Slice}\_X\_I\left(\mathrm{ny}\right)=\mathrm{Slice}\_X\left(\mathrm{ny}\right)*{\stackrel{-}{h}}_{01}\left(n\right)$

Slice_Y_I in formula (nx) is the Hilbert conversion of Slice_Y (nX), and the length of 1 dimension group Slice_Y_I (nx) is Nx(1 dimension group Slice_Y_I (nx) and 1 dimension group

the result array length of convolution is (Nx+65-1), by each 32 points before and after giving up, and Nx point in the middle of only choosing); And be array Slice_Y_E (nx, 1) nx=1 by 1 for array Slice_Y_I (ny) writes 2 ..., Nx;

Step 2, changes repeating step (6-1)～step (6-3) No. Crossline; until all No. Crossline calculated, thus result array Slice_Y_E (nx, the ny) nx=1 of the rim detection of section while obtaining current grade; ..., Nx; Ny=1 ..., Ny.

Images