CN103777241B - Three dimensional seismic data rapid edge-detection method based on time domain generalized Hilbert transformation - Google Patents

Abstract

The invention discloses a kind of three dimensional seismic data rapid edge-detection method based on time domain generalized Hilbert transformation, comprise the following steps: with the interval of interest of three dimensional seismic data for object of study, input 3D seismic data, when adopting three-dimensional etc., the computational methods of section, calculate section during the two dimension etc. of interval of interest; Calculate the Xi Shi factor of 21 dimensions of time domain, it is thus achieved that the edge detection results of section when currently waiting; Obtain the final edge detection results of current slice; Change etc. time slice number, calculate all of wait time section; According to section position in interval of interest when waiting, deposit the result of the rim detection of section during each grade, obtain the three-dimensional edges testing result of final interval of interest. The present invention has the precision of higher rim detection, and computational efficiency is high; On ensure that all of sampling point is cut into slices when being all located at respective grade, not only closest to real subsurface geology situation, and be conducive to the reduction of edge detection results.

Description

Three dimensional seismic data rapid edge-detection method based on time domain generalized Hilbert transformation

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 generalized Hilbert transformation.

Background technology

In oil gas field of geophysical exploration, 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, and particularly carbonate reservoir, crack therein, hole and hole are main migration pathway and the reservoir spaces of oil gas.

In the rock of nature 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, due to the restriction by resolution, it is impossible to identifying single hole, hole, seam, the scale that can only identify reaches a degree of Fractured Zone. It is true that the gathering role of oil gas is very little by single hole, hole, seam, really have exploration, Development volue be fracture-vug zone of certain scale.

Fractured Zone refers to relative containment body rock stratum, and its seam hole density significantly increases, and has the rock mass of certain expanded range. Therefore, it, on seismic horizontal slice (or horizon slice), has certain distribution and bearing of trend.It is therefore possible to use various edge detection method identification Fractured Zones.

1, Hilbert conversion

(1-1) frequency domain Hilbert conversion

Tie up discrete signal x (n) for 1, do Discrete Fourier Transform (DFT), obtain X (k), k=0,1 ..., N-1, whereinCorresponding 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=1,2,\mathrm{...},\frac{N}{2}-1\\ 0& k=\frac{N}{2},\mathrm{...},N-1\end{array}\right.---\left(1\right)$

H (k) is inverse DFT, and (this analytic signal is plural number, and its real part is equal to echo signal x (n) namely 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\·\stackrel{\‾}{x}\left(n\right)---\left(2\right)$

$\stackrel{\‾}{x}\left(n\right)=x\left(n\right)*\stackrel{\‾}{h}\left(n\right)=x\left(n\right)*\frac{1}{\mathrm{\π}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 it is the result of the HT of echo signal x (n);

In theory, the Xi Shi factorComputing formula be:

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

In formula, m is natural number;

(2) generalized Hilbert transformation:

LUOYi etc. (2003) propose generalized Hilbert transformation (GHT) method, which introduce window function and exponent number, from two aspects, conventional HT having been extended, tieed up discrete signal x (t) for 1, its GHT can be expressed as follows;

H (t)=h_r(t)+i��h_i(t)(5)

$h_r\left(t\right)={\{2\·{\mathrm{\Σ}}_{\mathrm{\ω}}{\{\mathrm{Re}\[X(t,\mathrm{\ω})\]\}}^n+{\{\mathrm{Re}\[X(t,0)\]\}}^n\}}^{\frac{1}{n}}---\left(6\right)$

$h_i\left(t\right)={\{2\·{\mathrm{\Σ}}_{\mathrm{\ω}}{\{\mathrm{Im}\[X(t,\mathrm{\ω})\]\}}^n\}}^{\frac{1}{n}}---\left(7\right)$

In formula, h (t) is the analytic signal of x (t); h_rT () is the real number road (equal to echo signal x (t)) of h (t); h_iT () is 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 window function is the box-shaped function of endless, formula (5), (6), (7) are just changed in quality for conventional HT;

Technical scheme currently for the rim detection of 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 processed: cut into slices (i.e. horizontal time slice) for three-dimensional time;

2. core calculations formula: adopt generalized Hilbert transformation (GHT) computing formula of frequency domain, namely adopts formula (5), (6), (7) to calculate;

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

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

Technical scheme:

1. the object processed: cut into slices (i.e. horizontal time slice) for three-dimensional time, 2 dimension groups;

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

3. calculation process: conventional time domain Hilbert alternative approach, in order to ensure computational accuracy, it is necessary to the Xi Shi factor is obtained very long, but this can increase the calculating time and reduce resolving power in actual applications. If but the factor obtains too short, then can produce this effect of jeep. And actual production often takes the shorter Xi Shi factor to improve computational efficiency, therefore, for overcoming this effect of jeep due to the shorter generation of the factor, adopt added-time window method. Thank and wait (2007) quietly the Xi Shi factor is added the process of Hamming window.

Calculation process 1: determine that Xi Shi factor optimum length is 65, and the Xi Shi factor is added Hamming window (namely 2 are multiplied);

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

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

Technical scheme:

1. the object processed: cut into slices (i.e. horizontal time slice) for three-dimensional time, 2 dimension groups;

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

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

4. analyze: proposing the computational methods in 2 directions, but still adopt the HT of routine to calculate, computational accuracy still has much room for improvement.

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

Technical scheme:

1. the object processed: cut into slices (i.e. horizontal time slice) for three-dimensional time, 2 dimension groups;

2. core calculations formula: (a) adopts generalized Hilbert transformation (GHT) computing formula of frequency domain, namely adopts formula (5), (6), (7) to calculate; B () adopts 2 repair and maintenance limit denoising methods (LUOYi etc., 2002), carry out pretreatment calculating.

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

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

(1) 4 kind of scheme is both for three-dimensional time section (i.e. 2 dimension horizontal time slice), it does not have for the processing scheme of the actual formation (actual formation is non-horizontal often) of three dimensional seismic data;

(2) 3 kinds of schemes (the 1st, 2,4 kind) are all tie up section by 2 to be decomposed into several 1 dimension groups, and namely 1 direction just for 2 dimension sections carries out detecting (laterally or longitudinally), and the precision of its rim detection is relatively low;

The scheme of (3) the 1st kinds of method Chen Xue China etc. (2008) adopts the GHT method of frequency domain to calculate, and computational efficiency is relatively low, and suppresses noise immune poor (do not have in calculation process and consider to suppress noise);

(4) the 2nd kinds of methods are thanked to the scheme waiting (2007) quietly and are adopted windowing HT method to calculate, and computational efficiency is high, but suppress noise immune poor (do not have in calculation process and consider to suppress noise);

The scheme that (5) the 3rd kinds of methods thank quiet (2007) considers the process suppressing noise and 2 directions to two dimension slicing in calculation process, but still adopts the HT of routine to be calculated, and computational accuracy is low;

The scheme of (6) the 4th kinds of method Xiong Xiao armies etc. (2009) considers suppression noise in calculation process, but still adopts the GHT method of frequency domain to calculate, and computational efficiency is relatively low.

The existing edge detection method lacking solution 3D seismic data; The computational efficiency of rim detection is low, and computational accuracy is low.

Summary of the invention

The purpose of the embodiment of the present invention is in that to provide a kind of three dimensional seismic data rapid edge-detection method based on time domain generalized Hilbert transformation, it is intended to solve the existing edge detection method lacking and solving 3D seismic data; The computational efficiency 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 generalized Hilbert transformation, should comprise the following steps based on the three dimensional seismic data rapid edge-detection method of time domain generalized Hilbert transformation:

With the interval of interest of three dimensional seismic data for object of study, inputting 3D seismic data, when adopting three-dimensional etc., the computational methods of section, calculate section during the two dimension etc. of interval of interest;

Calculate the Xi Shi factor of 21 dimensions of time domain, section when waiting for the 1st, carry out 2 repair and maintenance limit denoisings, obtain the array after denoising;

Choosing the Xi Shi factor that length is 65, order is the situation of 1, is calculated choosing the Xi Shi factor that length is 33, and order is the situation of 2, until obtaining the edge detection results of section during current grade;

Comprehensive order 1 and the result of calculation of order 2, obtain the final edge detection results of current slice;

Change etc. time slice number, calculate all of wait time section; According to section position in interval of interest when waiting, deposit the result of the rim detection of section during each grade, obtain the three-dimensional edges testing result of final interval of interest.

Further, following steps should be specifically included based on the three dimensional seismic data rapid edge-detection method of time domain generalized Hilbert transformation:

Step one, with the interval of interest of three dimensional seismic data for object of study, inputs 3D seismic data, the interlayer position, top of target interval and interlayer position, the end, and target interval is D, and the interlayer position, top of interval of interest is Hor_1, and interlayer position, the end is Hor_2;

Step 2, when adopting three-dimensional etc., the computational methods of section, calculate section during the two dimension etc. of interval of interest: when interlayer position, the top Hor_1 of interval of interest is parallel with interlayer position, end Hor_2, the sum (N1) cut into slices during 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 interlayer position, the top Hor_1 and interlayer position, end Hor_2 of interval of interest are not parallel, two dimension etc. time section number (N1):

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 time section totally 15, wherein the section of the 2nd, 3,4 deciles is non-fully sampled in Inline or Crossline direction, unified section carries out mending 0 process, mend the mesh point of 0 process, it is not involved in follow-up calculating, it is ensured that when waiting, section is complete, and namely during all grades, section is consistent in the length of vertical and horizontal;

Step 3, calculates the Xi Shi factor of 21 dimensions of time domainWithIts computing formula is identical,

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

In formula, m is natural number, the Xi Shi factorTotal length be set as 65, the Xi Shi factorsTotal length be set as 33; And the Xi Shi factor to time domainWithApplying Gaussian window w (n), the length of window array is 65;

$\stackrel{-}{h_1}\left(n\right)={\stackrel{\‾}{h}}_{01}\left(n\right)\×w\left(n\right),n=1,\mathrm{...},65---\left(11\right)$

$\stackrel{-}{h_2}\left(n\right)={\stackrel{\‾}{h}}_{02}\left(n\right)\×w\left(n\right),n=1,\mathrm{...},33---\left(12\right)$

Obtain 2 Xi Shi factors after windowing processWithLength is the Xi Shi factor situation that order is 1 that is equivalent in GHT of 65, namely conventional HT; Length is the Xi Shi factor situation that order is 2 that is equivalent in GHT of 33;

Step 4, section when waiting for the 1st, being 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, from the CrossLine sum that 1 to Ny, Ny are three dimensional seismic data, carries out 2 repair and maintenance limit denoisings, obtain the array Slice_N after denoising (nx, ny);

Step 5, chooses the Xi Shi factor that length is 65Order is the situation of 1, is calculated; Method particularly includes:

The first step, for 2 dimension group Slice_N (nx, ny) (nx=1 ..., Nx; Ny=1 ..., Ny), it is circulated calculating in Inline direction:

1 dimension group of extraction No. Inline=1, 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) withConvolution,

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

The Hilbert that in formula, Slice_X_I (ny) is Slice_X (ny) converts, and the length of 1 dimension group Slice_X_I (ny) is Ny (1 dimension group Slice_X_I (ny) and 1 dimension groupThe result array length of convolution is (Ny+65-1), each 32 points before and after result array is given up, and namely only chooses Ny point of centre); And 1 dimension group Slice_X_I (ny) is write 2 dimension group Slice_X_E (1, ny) ny=1 ..., Ny;

Change No. Inline, repeat the first step and second step, until all of Inline calculating completes, result array Slice_X_E (nx, the ny) nx=1 of the rim detection of section when obtaining currently waiting ..., Nx; Ny=1 ..., Ny;

Second step, for 2 dimension group Slice_N (nx, ny), is circulated calculating in Crossline direction:

Step (1), 1 dimension group of extraction No. Crossline=1, 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) withConvolution:

$Slice\_Y\_I\left(nx\right)=Slice\_Y\left(nx\right)*\stackrel{-}{h_{01}}\left(n\right)---\left(13\right)$

The Hilbert that in formula, Slice_Y_I (nx) is Slice_Y (nX) converts, and the length of 1 dimension group Slice_Y_I (nx) is Nx (1 dimension group Slice_Y_I (nx) and 1 dimension groupThe result array length of convolution is (Nx+65-1), by each 32 points before and after giving up, namely only chooses Nx point of centre); And 1 dimension group Slice_Y_I (ny) is write 2 dimension group Slice_Y_E (nx, 1) nx=1 ..., Nx;

Step (2), changes No. Crossline, repeats step (1), until all of Crossline calculating completes, thus result array Slice_Y_E (nx, the ny) nx=1 of the rim detection of section when obtaining currently waiting, ..., Nx; Ny=1 ..., Ny;

3rd step, calculates the edge detection results of section during current grade:

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

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

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

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

Step 8, slice number during replacing etc., repeat step 4��step 7, until all of when waiting section calculating complete;

Step 9, according to section position in interval of interest when waiting, deposits the result of the rim detection of section during each grade, obtains the three-dimensional edges testing result of final interval of interest.

Further, in step 5, choose the Xi Shi factor that length is 65Order is the situation of 1, is calculated method particularly includes:

The first step, for 2 dimension group Slice_N (nx, ny) (nx=1 ..., Nx;Ny=1 ..., Ny), it is circulated calculating in Inline direction:

Second step, for 2 dimension group Slice_N (nx, ny), is circulated calculating in Crossline direction:

3rd step, calculates the edge detection results of section during current grade:

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

Further, in the first step, 1 dimension group of extraction No. Inline=1, i.e. Slice_N (1, ny), ny=1 ..., Ny; Making this 1 dimension group is Slice_X (ny), ny=1 ..., Ny; Method particularly includes:

Step one, calculate Slice_X (ny) withConvolution,

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

The Hilbert that in formula, Slice_X_I (ny) is Slice_X (ny) converts, and the length of 1 dimension group Slice_X_I (ny) is Ny (1 dimension group Slice_X_I (ny) and 1 dimension groupThe result array length of convolution is (Ny+65-1), each 32 points before and after result array is given up, and namely only chooses Ny point of centre); And 1 dimension group Slice_X_I (ny) is write 2 dimension group Slice_X_E (1, ny) ny=1 ..., Ny;

Step 2, changes No. Inline, and all of Inline calculating completes, result array Slice_X_E (nx, the ny) nx=1 of the rim detection of section when obtaining currently waiting ..., Nx; Ny=1 ..., Ny.

Further, in second step, for 2 dimension group Slice_N (nx, ny), it is circulated calculating in Crossline direction method particularly includes:

Step (1), 1 dimension group of extraction No. Crossline=1, 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) withConvolution:

$Slice\_Y\_I\left(nx\right)=Slice\_Y\left(nx\right)*\stackrel{-}{h_{01}}\left(n\right)---\left(13\right)$

The Hilbert that in formula, Slice_Y_I (nx) is Slice_Y (nX) converts, and the length of 1 dimension group Slice_Y_I (nx) is Nx (1 dimension group Slice_Y_I (nx) and 1 dimension groupThe result array length of convolution is (Nx+65-1), by each 32 points before and after giving up, namely only chooses Nx point of centre); And 1 dimension group Slice_Y_I (ny) is write 2 dimension group Slice_Y_E (nx, 1) nx=1 ..., Nx;

Step (2), changes No. Crossline, repeats step (1), until all of Crossline calculating completes, thus result array Slice_Y_E (nx, the ny) nx=1 of the rim detection of section when obtaining currently waiting, ..., Nx; Ny=1 ..., Ny.

Three dimensional seismic data rapid edge-detection method based on time domain generalized Hilbert transformation provided by the invention, 1 dimension time domain broad sense HT is adopted to process 3D seismic data, contain windowing process, contain again order to process, calculation process contains noise suppressed process, and during from 2 dimensions etc., 2 directions of section process, and have the precision of higher rim detection simultaneously; The core calculations formula of the present invention is the HT of time domain, and computational efficiency is high; Adopt etc. time slicing treatment 3D seismic data, and time adjacent etc. the interval of section equal to 1 time sampling interval, on ensure that all of sampling point is cut into slices when being all located at respective grade, not only closest to real subsurface geology situation, and be conducive to the reduction of edge detection results; The Xi Shi factor adopting 2 different lengths processes, and is equivalent to the situation of 2 orders in frequency domain, and combines the result of calculation of 2 orders, and computational accuracy is high. The present invention is directed to three dimensional seismic data to process, it is thus achieved that the three-dimensional edges detection data volume of target interval, for the identification of Fractured Zone or fracture-vug zone, effectively instruct the reservoir prediction of oil gas geophysical exploration; The boundary layer position, top of given 3-d seismic data set and interval of interest and bottom interface layer position, it is thus achieved that three-dimensional edge detection data body; It is calculated the optimum organization of flow process, considers noise suppressed, computational efficiency 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 generalized Hilbert transformation that the embodiment of the present invention provides;

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

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

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

Detailed description of the invention

In order to make the purpose of the present invention, technical scheme and advantage clearly understand, below in conjunction with embodiment, the present invention is further elaborated. Should be appreciated that specific embodiment described herein is 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, the application principle of the present invention is further described.

As it is shown in figure 1, the three dimensional seismic data rapid edge-detection method based on time domain generalized Hilbert transformation of the embodiment of the present invention comprises the following steps:

S101: with the interval of interest of three dimensional seismic data for object of study, inputs 3D seismic data, and when adopting three-dimensional etc., the computational methods of section, calculate section during the two dimension etc. of interval of interest;

S102: calculate the Xi Shi factor of 21 dimensions of time domain, section when waiting for the 1st, carry out 2 repair and maintenance limit denoisings, obtain the array after denoising;

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

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

S105: slice number during replacing etc., calculates section during all of grade; According to section position in interval of interest when waiting, deposit the result of the rim detection of section during each grade, obtain the three-dimensional edges testing result of final interval of interest.

The concrete steps of the present invention include:

Step one, with the interval of interest of three dimensional seismic data for object of study, (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 (target interval in Fig. 2 is D to Fig. 2, the interlayer position, top of interval of interest is Hor_1, and interlayer position, the end is Hor_2);

Step 2, computational methods (the ZengHongliu of section when adopting three-dimensional etc., 2010), calculate interval of interest two dimension etc. time section: when interlayer position, the top Hor_1 of interval of interest is parallel with interlayer position, end Hor_2, two dimension etc. time cut into slices sum (N1):

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, Deng time section calculating schematic diagram see Fig. 3 (grid in Fig. 3 represents longitudinal direction on 1 time sampling point, 1 No. Inline or No. Crossline transversely); When interlayer position, the top Hor_1 and interlayer position, end Hor_2 of interval of interest are not parallel, two dimension etc. time section number (N1):

N1=Max_L/dt (9)

In formula, Max_L represents maximum time distance (unit: ms) of Hor_2 and Hor_1, dt represents the time sampling interval (unit: ms) of 3D seismic data, Deng time section calculating schematic diagram see Fig. 4 (grid in Fig. 4 represents longitudinal direction on 1 time sampling point, 1 No. Inline or No. Crossline transversely, the maximum time distance pushing up interlayer position and interlayer position, the end in Fig. 4 is 16 times adopt interval, when reaching adjacent grade by 4 deciles, the interval of section is all 1 time adopt interval), section totally 15 when waiting in Fig. 4, wherein the 2nd, 3, (only have on Partial Mesh node and have value) is non-fully sampled in the section of 4 deciles in Inline or Crossline direction, in Practical Calculation, unified these sections being carried out mends the 0 process (mesh point of benefit 0 process, it is not involved in follow-up calculating), when ensureing etc. complete (all wait for time section consistent in the length of vertical and horizontal) of section,

Step 3, calculates the Xi Shi factor of 21 dimensions of time domainWithIts computing formula is identical,

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

In formula, m is natural number, the Xi Shi factorTotal length be set as 65, the Xi Shi factorsTotal length be set as 33; And the Xi Shi factor to time domainWithApply Gaussian window w (n) (length of window array is 65);

$\stackrel{-}{h_1}\left(n\right)={\stackrel{\‾}{h}}_{01}\left(n\right)\×w\left(n\right),n=1,\mathrm{...},65---\left(11\right)$

$\stackrel{-}{h_2}\left(n\right)={\stackrel{\‾}{h}}_{02}\left(n\right)\×w\left(n\right),n=1,\mathrm{...},33---\left(12\right)$

Obtain 2 Xi Shi factors after windowing processWith(length be 65 the Xi Shi factor situation that order is 1 that is equivalent in GHT, namely conventional HT; Length is the Xi Shi factor situation that order is 2 that is equivalent in GHT of 33)

Step 4, when waiting for the 1st, section (is designated as array Slice (nx, ny), the value of nx is total from the InLine that 1 to Nx, Nx are 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 after denoising (nx, ny);

Step 5, chooses the Xi Shi factor that length is 65(order is the situation of 1), is calculated; Method particularly includes:

The first step, for 2 dimension group Slice_N (nx, ny) (nx=1 ..., Nx; Ny=1 ..., Ny), it is circulated calculating in Inline direction:

1 dimension group of extraction No. Inline=1, 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) withConvolution,

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

The Hilbert that in formula, Slice_X_I (ny) is Slice_X (ny) converts, and the length of 1 dimension group Slice_X_I (ny) is Ny (1 dimension group Slice_X_I (ny) and 1 dimension groupThe result array length of convolution is (Ny+65-1), each 32 points before and after result array is given up, and namely only chooses Ny point of centre); And 1 dimension group Slice_X_I (ny) is write 2 dimension group Slice_X_E (1, ny) ny=1 ..., Ny;

Change No. Inline, repeat the first step and second step, until all of Inline calculating completes, result array Slice_X_E (nx, the ny) nx=1 of the rim detection of section when obtaining currently waiting ..., Nx; Ny=1 ..., Ny;

Second step, for 2 dimension group Slice_N (nx, ny), is circulated calculating in Crossline direction:

Step (1), 1 dimension group of extraction No. Crossline=1, 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) withConvolution:

$Slice\_Y\_I\left(nx\right)=Slice\_Y\left(nx\right)*\stackrel{-}{h_{01}}\left(n\right)---\left(13\right)$

The Hilbert that in formula, Slice_Y_I (nx) is Slice_Y (nX) converts, and the length of 1 dimension group Slice_Y_I (nx) is Nx (1 dimension group Slice_Y_I (nx) and 1 dimension groupThe result array length of convolution is (Nx+65-1), by each 32 points before and after giving up, namely only chooses Nx point of centre); And 1 dimension group Slice_Y_I (ny) is write 2 dimension group Slice_Y_E (nx, 1) nx=1 ..., Nx;

Step (2), change No. Crossline, repeat step (1), content (original position see annotation) under step 2 is until all of Crossline calculating completes, thus the result array Slice_Y_E (nx of the rim detection of section when obtaining currently waiting, ny) nx=1 ..., Nx; Ny=1 ..., Ny;

3rd step, calculates the edge detection results of section during current grade:

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

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

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

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

Step 8, slice number during replacing etc., repeat step 4��step 7, until all of when waiting section calculating complete;

Step 9, according to section position in interval of interest when waiting, deposits the result of the rim detection of section during each grade, obtains the three-dimensional edges testing result of final interval of interest.

The present invention adopts 1 dimension time domain broad sense HT to process 3D seismic data, similar with frequency domain GHT, has both contained windowing process, contains again order and processes; Process additionally, contain noise suppressed in calculation process, and during from 2 dimensions etc., 2 directions of section process, and have the precision of higher rim detection simultaneously; Additionally, the core calculations formula of the present invention is the HT of time domain, computational efficiency is high. Adopt etc. time slicing treatment 3D seismic data, and adjacent wait time section interval be equal to 1 time sampling interval. On present invention ensures that all of sampling point is cut into slices when being all located at respective grade, not only closest to real subsurface geology situation, and be conducive to the reduction (namely reconstituting 3D data volume) of edge detection results; The Xi Shi factor adopting 2 different lengths processes, and is equivalent to the situation of 2 orders in frequency domain, and combines the result of calculation of 2 orders, and computational accuracy is high.

The foregoing is only presently preferred embodiments of the present invention, not in order to limit the present invention, all any amendment, equivalent replacement and improvement etc. made within the spirit and principles in the present invention, should be included within protection scope of the present invention.

Claims (1)

1. the three dimensional seismic data rapid edge-detection method based on time domain generalized Hilbert transformation, it is characterised in that following steps should be specifically included based on the three dimensional seismic data rapid edge-detection method of time domain generalized Hilbert transformation:

Step one, with the interval of interest of three dimensional seismic data for object of study, inputs 3D seismic data, the interlayer position, top of target interval and interlayer position, the end, and target interval is D, and the interlayer position, top of interval of interest is Hor_1, and interlayer position, the end is Hor_2;

Step 2, when adopting three-dimensional etc., the computational methods of section, calculate section during the two dimension etc. of interval of interest: when interlayer position, the top Hor_1 of interval of interest is parallel with interlayer position, end Hor_2, the sum (N1) cut into slices during 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 interlayer position, the top Hor_1 and interlayer position, end Hor_2 of interval of interest are not parallel, two dimension etc. time section number (N1):

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 time section totally 15, wherein the section of the 2nd, 3,4 deciles is non-fully sampled in Inline or Crossline direction, unified section carries out mending 0 process, mend the mesh point of 0 process, it is not involved in follow-up calculating, it is ensured that when waiting, section is complete, and namely during all grades, section is consistent in the length of vertical and horizontal;

Step 3, calculates the Xi Shi factor of 21 dimensions of time domainWithIts computing formula is identical,

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

In formula, m is natural number, the Xi Shi factorTotal length be set as 65, the Xi Shi factorsTotal length be set as 33; And the Xi Shi factor to time domainWithApplying Gaussian window w (n), the length of window array is 65;

$\stackrel{-}{h_1}\left(n\right)={\stackrel{\‾}{h}}_{01}\left(n\right)\×w\left(n\right),n=1,...,65$

$\stackrel{-}{h_2}\left(n\right)={\stackrel{\‾}{h}}_{02}\left(n\right)\×w\left(n\right),n=1,...,33$

Obtain 2 Xi Shi factors after windowing processWithLength is the Xi Shi factor situation that order is 1 that is equivalent in GHT of 65, namely conventional HT; Length is the Xi Shi factor situation that order is 2 that is equivalent in GHT of 33;

Step 4, section when waiting for the 1st, being 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, from the CrossLine sum that 1 to Ny, Ny are three dimensional seismic data, carries out 2 repair and maintenance limit denoisings, obtain the array Slice_N after denoising (nx, ny);

Step 5, chooses the Xi Shi factor that length is 65Order is the situation of 1, is calculated; Method particularly includes:

The first step, for 2 dimension group Slice_N (nx, ny) (nx=1 ..., Nx; Ny=1 ..., Ny), it is circulated calculating in Inline direction:

1 dimension group of extraction No. Inline=1, 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) withConvolution,

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

The Hilbert that in formula, Slice_X_I (ny) is Slice_X (ny) converts, and the length of 1 dimension group Slice_X_I (ny) is Ny (1 dimension group Slice_X_I (ny) and 1 dimension groupThe result array length of convolution is (Ny+65-1), each 32 points before and after result array is given up, and namely only chooses Ny point of centre); And 1 dimension group Slice_X_I (ny) is write 2 dimension group Slice_X_E (1, ny) ny=1 ..., Ny;

Change No. Inline, repeat the first step, until all of Inline calculating completes, result array Slice_X_E (nx, the ny) nx=1 of the rim detection of section when obtaining currently waiting ..., Nx; Ny=1 ..., Ny;

Second step, for 2 dimension group Slice_N (nx, ny), is circulated calculating in Crossline direction:

Step (1), 1 dimension group of extraction No. Crossline=1, 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) withConvolution:

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

The Hilbert that in formula, Slice_Y_I (nx) is Slice_Y (nX) converts, and the length of 1 dimension group Slice_Y_I (nx) is Nx (1 dimension group Slice_Y_I (nx) and 1 dimension groupThe result array length of convolution is (Nx+65-1), by each 32 points before and after giving up, namely only chooses Nx point of centre); And 1 dimension group Slice_Y_I (ny) is write 2 dimension group Slice_Y_E (nx, 1) nx=1 ..., Nx;

Step (2), changes No. Crossline, repeats step (1), until all of Crossline calculating completes, thus result array Slice_Y_E (nx, the ny) nx=1 of the rim detection of section when obtaining currently waiting, ..., Nx; Ny=1 ..., Ny;

3rd step, calculates the edge detection results of section during current grade:

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

Step 6, chooses the Xi Shi factor that length is 33Order is the situation of 2, repeats the first step to the 3rd step, until obtaining the edge detection results Slice_E2 (nx, ny) of section during current grade;

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

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

Step 8, slice number during replacing etc., repeat step 4��step 7, until all of when waiting section calculating complete;

Step 9, according to section position in interval of interest when waiting, deposits the result of the rim detection of section during each grade, obtains the three-dimensional edges testing result of final interval of interest.