CN103777241B

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

Info

Publication number: CN103777241B
Application number: CN201410027618.7A
Authority: CN
Inventors: 熊晓军; 林华伟; 吕龑
Original assignee: Chengdu Univeristy of Technology
Current assignee: Chengdu Univeristy of Technology
Priority date: 2014-01-21
Filing date: 2014-01-21
Publication date: 2016-06-08
Anticipated expiration: 2034-01-21
Also published as: CN103777241A

Abstract

The invention discloses a method for fast edge detection of three-dimensional seismic data based on time-domain generalized Hilbert transform, comprising the following steps: taking the target section of the three-dimensional seismic data as the research object, inputting the three-dimensional seismic data, and adopting a calculation method of three-dimensional isochronous slices , calculate the two-dimensional isochronous slice of the target layer; calculate two 1-dimensional Hisfactors in the time domain, and obtain the edge detection result of the current isochronous slice; obtain the final edge detection result of the current slice; replace the isochronous slice number , calculate all the isochronous slices; store the edge detection result of each isochronous slice according to the position of the isochronous slice in the target layer, and obtain the final 3D edge detection result of the target layer. The invention has high edge detection accuracy and high calculation efficiency; it ensures that all sample points are located on their respective isochronous slices, which is not only closest to the real underground geological conditions, but also facilitates the restoration of edge detection results.

Description

Fast Edge Detection Method for 3D Seismic Data Based on Generalized Hilbert Transform in Time Domain

技术领域technical field

本发明属于油气地球物理勘探技术领域及一种基于时间域广义Hilbert变换的三维地震资料快速边缘检测方法。The invention belongs to the technical field of oil and gas geophysical exploration and relates to a fast edge detection method for three-dimensional seismic data based on time-domain generalized Hilbert transform.

背景技术Background technique

在油气地球物理勘探领域，边缘检测方法主要用于寻找断层或裂缝发育带。其中裂缝发育带的预测对油气的勘探与开发具有重要意义，特别是碳酸盐岩储层，其中的裂缝、孔隙和孔洞是油气的主要运移通道和储集空间。In the field of oil and gas geophysical exploration, edge detection methods are mainly used to find faults or fracture development zones. Among them, the prediction of fracture development zones is of great significance to the exploration and development of oil and gas, especially in carbonate reservoirs, where fractures, pores and vugs are the main migration channels and storage spaces for oil and gas.

自然界的岩石或岩层中普遍存在着孔隙、裂缝(隙)和溶蚀孔洞，其形态各异，尺度大小也相差悬殊。对于地震勘探而言，由于受分辨率的限制，无法识别出单个的孔、洞、缝，只能识别规模达到一定程度的裂缝发育带。事实上，单个的孔、洞、缝对油气的聚集所起的作用是微乎其微的，真正具有勘探、开发价值的是具有一定规模的缝洞发育带。Pores, cracks (cracks) and dissolved pores are ubiquitous in rocks or rock formations in nature, and their shapes and scales are also very different. For seismic exploration, due to the limitation of resolution, it is impossible to identify individual pores, holes, and fractures, and only fracture development zones with a certain scale can be identified. In fact, individual pores, caves, and fractures play little role in the accumulation of oil and gas, and the fracture-cave belts with a certain scale are really valuable for exploration and development.

裂缝发育带是指相对围体岩层而言，其缝洞密度明显增大，并有一定延伸范围的岩体。因此，其在地震水平切片(或沿层切片)上，具有一定的分布范围和延伸方向。因此，可以采用各种边缘检测方法识别裂缝发育带。The fracture development zone refers to the rock mass whose fracture-cavity density is obviously increased compared with the surrounding rock strata, and has a certain extension range. Therefore, it has a certain distribution range and extension direction on the seismic horizontal slice (or slice along the layer). Therefore, various edge detection methods can be used to identify fracture development zones.

1、Hilbert变换1. Hilbert transform

(1-1)频率域Hilbert变换(1-1) Frequency domain Hilbert transform

对于1维离散信号x(n)，做离散Fourier变换(DFT)，得到X(k),k＝0,1,...,N-1，其中对应负频率，N为离散信号x(t)的采样点数，再令：For a 1-dimensional discrete signal x(n), do a discrete Fourier transform (DFT) to get X(k), k=0,1,...,N-1, where Corresponding to the negative frequency, N is the number of sampling points of the discrete signal x(t), and then:

$H h ((k k)) = = \{\begin{matrix} X x ((k k)) & k k = = 00 \\ 22 X x ((k k)) & k k = = 11,, 22,, ... ...,, \frac{N N}{22} - - 11 \\ 00 & k k = = \frac{N N}{22},, ... ...,, N N - - 11 \end{matrix} - - - - - - ((11))$

对H(k)做逆DFT，即得到x(n)的解析信号h(n)(该解析信号为复数，其实部等于目标信号x(n)；虚部是x(n)的HT的结果，被用于边缘检测)；Perform inverse DFT on H(k), that is, get the analytical signal h(n) of x(n) (the analytical signal is a complex number, the real part is equal to the target signal x(n); the imaginary part is the result of HT of x(n) , are used for edge detection);

(1-2)时间域Hilbert变换：(1-2) Time domain Hilbert transform:

时间域的Hilbert变换表现为褶积关系；The Hilbert transform in the time domain is represented as a convolution relationship;

$h h ((n no)) = = x x ((n no)) + + i i \cdot \cdot \overset{&OverBar; &OverBar;}{x x} ((n no)) - - - - - - ((22))$

$\overset{&OverBar; &OverBar;}{x x} ((n no)) = = x x ((n no)) * * \overset{&OverBar; &OverBar;}{h h} ((n no)) = = x x ((n no)) * * \frac{11}{π π n no} - - - - - - ((33))$

式中，h(n)为1维离散信号x(n)的解析信号，为希氏因子，是目标信号x(n)的HT的结果；In the formula, h(n) is the analytical signal of 1-dimensional discrete signal x(n), is the His factor, is the result of the HT of the target signal x(n);

理论上，希氏因子的计算公式是：Theoretically, the His factor The formula for calculating is:

$\overset{&OverBar; &OverBar;}{h h} ((n no)) = = \{\begin{matrix} 00 & n no = = 22 m m \\ \frac{22}{π π n no} & n no = = 22 m m + + 11 \end{matrix} - - - - - - ((44))$

式中m为自然数；In the formula, m is a natural number;

(2)广义Hilbert变换：(2) Generalized Hilbert transform:

LUOYi等(2003)提出了广义Hilbert变换(GHT)方法，它引入了窗函数和阶数，从两个方面对常规的HT进行了扩展，对于1维离散信号x(t)，其GHT可以表述如下；LUOYi et al. (2003) proposed the generalized Hilbert transform (GHT) method, which introduced the window function and order, and extended the conventional HT from two aspects. For a 1-dimensional discrete signal x(t), its GHT can express as follows;

h(t)＝h_r(t)+i·h_i(t)(5)h(t)=h _r (t)+i h _i (t)(5)

${h h}_{r r} ((t t)) = = {{22 \cdot &Center Dot; {Σ Σ}_{ω ω} {{Re Re [[X x ((t t,, ω ω))]]}}^{n no} + + {{Re Re [[X x ((t t,, 00))]]}}^{n no}}}^{\frac{11}{n no}} - - - - - - ((66))$

${h h}_{i i} ((t t)) = = {{22 \cdot \cdot {Σ Σ}_{ω ω} {{Im Im [[X x ((t t,, ω ω))]]}}^{n no}}}^{\frac{11}{n no}} - - - - - - ((77))$

式中，h(t)是x(t)的解析信号；h_r(t)是h(t)的实数道(等于目标信号x(t))；h_i(t)是h(t)的虚数道(用于边缘检测)；X(t,ω)是x(t)的加窗Fourier变换；n是阶数，当n＝1，窗函数为无限长的箱形函数时，公式(5)、(6)、(7)就蜕变为常规的HT；In the formula, h(t) is the analytical signal of x(t); h _r (t) is the real channel of h(t) (equal to the target signal x(t)); h _i (t) is the Imaginary number track (for edge detection); X (t, ω) is the windowed Fourier transform of x (t); n is the order, when n=1, when the window function is an infinitely long box function, the formula (5 ), (6), (7) just transform into conventional HT;

目前针对地震资料的边缘检测的技术方案：The current technical solutions for edge detection of seismic data:

(1)陈学华、贺振华、黄德济。地震资料的高阶伪希尔伯特变换边缘检测，2008年8月，地球物理学进展，第23卷第4气，第1106-1110页；(1) Chen Xuehua, He Zhenhua, Huang Deji. High-order pseudo-Hilbert transform edge detection for seismic data, August 2008, Advances in Geophysics, Vol. 23, No. 4, pp. 1106-1110;

技术方案：Technical solutions:

①处理的对象：针对三维时间切片(即水平时间切片)；①Processing objects: for three-dimensional time slices (that is, horizontal time slices);

②核心计算公式：采用频率域的广义Hilbert变换(GHT)计算公式，即采用公式(5)、(6)、(7)计算；② Core calculation formula: use the generalized Hilbert transform (GHT) calculation formula in the frequency domain, that is, use formulas (5), (6), and (7) for calculation;

③计算流程：将2维切片视为2维数组；再将2维数组按照横向或纵向分解为若干个1维数组；采用1维频率域GHT分别计算每个1维数组的GHT，得到解析信号，并提取的虚部用于边缘检测；循环计算完所有的1维数组，得到最终的2维切片的边缘检测结果数组。③ Calculation process: treat the 2-dimensional slice as a 2-dimensional array; then decompose the 2-dimensional array into several 1-dimensional arrays according to the horizontal or vertical direction; use the 1-dimensional frequency domain GHT to calculate the GHT of each 1-dimensional array to obtain the analytical signal , and the extracted imaginary part is used for edge detection; after calculating all the 1-dimensional arrays in a loop, the final 2-dimensional slice edge detection result array is obtained.

(2)谢静、贺振华、黄德济、熊晓军。基于加窗Hilbert变换的高精度边缘检测，2007年4月，油气地球物理，第5卷第1期，第27-29页(2) Xie Jing, He Zhenhua, Huang Deji, Xiong Xiaojun. High-precision edge detection based on windowed Hilbert transform, April 2007, Oil and Gas Geophysics, Vol. 5, No. 1, pp. 27-29

技术方案：Technical solutions:

①处理的对象：针对三维时间切片(即水平时间切片)，2维数组；①Processed objects: for three-dimensional time slices (that is, horizontal time slices), 2-dimensional arrays;

②核心计算公式：采用时间域的Hilbert变换计算公式，即采用公式(3)、(4)计算；② Core calculation formula: use the Hilbert transformation calculation formula in the time domain, that is, use formulas (3) and (4) for calculation;

③计算流程：常规的时间域Hilbert变换方法，为了确保计算精度，需要把希氏因子取得很长，但在实际应用中这会增加计算时间和降低分辨力。但如果因子取得过短，则会产生吉普斯效应。而实际生产中为了提高计算效率往往取较短的希氏因子，因此，为克服由于因子较短产生的吉普斯效应，采用加时窗方法。谢静等(2007)对希氏因子加汉明窗的处理。③Calculation process: In the conventional time-domain Hilbert transform method, in order to ensure the calculation accuracy, the Hilbert factor needs to be obtained very long, but this will increase the calculation time and reduce the resolution in practical applications. But if the factor is taken too short, there will be a Gipps effect. However, in actual production, a shorter Hisfactory factor is often used to improve calculation efficiency. Therefore, in order to overcome the Gibbs effect caused by the shorter factor, the method of adding time window is adopted. Xie Jing et al. (2007) dealt with His factor plus Hamming window.

计算流程1：确定希氏因子最佳长度为65，并对希氏因子加汉明窗(即2者相乘)；Calculation process 1: Determine the optimal length of the Hisfactory factor to be 65, and add a Hamming window to the Historic factor (that is, multiply the two together);

计算流程2：将2维切片视为2维数组；再将2维数组按照横向或纵向分成若干个1维数组；采用1维时间域HT(公式(2)和(3))分别计算每个1维数组的HT(用于边缘检测)；循环计算完所有的1维数组，得到最终的2维切片的边缘检测结果数组。Calculation process 2: Treat the 2-dimensional slice as a 2-dimensional array; then divide the 2-dimensional array into several 1-dimensional arrays horizontally or vertically; use the 1-dimensional time domain HT (formula (2) and (3)) to calculate each HT of 1-dimensional array (for edge detection); all 1-dimensional arrays are calculated in a loop, and the final array of edge detection results of 2-dimensional slices is obtained.

(3)谢静。基于时间域加窗Hilbert变换的地震边缘检测研究，2007，成都理工大学硕士论文(3) Xie Jing. Research on Seismic Edge Detection Based on Time Domain Windowed Hilbert Transform, 2007, Master Thesis of Chengdu University of Technology

技术方案：Technical solutions:

③计算流程：确定希氏因子最佳长度为65，并对希氏因子加汉明窗(即2者相乘，；将2维切片视为2维数组；再将2维数组分别按照横向和纵向分成若干个1维数组；采用1维时间域HT(公式(2)和(3))分别计算每个1维数组的HT(用于边缘检测)；循环计算完所有的1维数组，得到纵向和横向的计算结果；再取纵向和横向计算结果的均方根，得到最终的2维切片的边缘检测结果数组。③Calculation process: Determine the optimal length of the Hisfactory factor to be 65, and add a Hamming window to the Hisfactory factor (that is, multiply the two; treat the 2-dimensional slice as a 2-dimensional array; then divide the 2-dimensional array according to the horizontal and Divide into several 1-dimensional arrays vertically; use 1-dimensional time domain HT (formulas (2) and (3)) to calculate the HT of each 1-dimensional array (for edge detection); after looping through all 1-dimensional arrays, we get The vertical and horizontal calculation results; then take the root mean square of the vertical and horizontal calculation results to obtain the final 2-dimensional slice edge detection result array.

④分析：提出了2个方向的计算方法，但是仍采用常规的HT计算，计算精度仍有待提高。④ Analysis: Calculation methods in two directions are proposed, but the conventional HT calculation is still used, and the calculation accuracy still needs to be improved.

(4)熊晓军、贺振华、赵明金等。一种基于GHT的裂缝检测新方法，2009年8月，石油地球物理勘探，第44卷第4期，第442-444页(4) Xiong Xiaojun, He Zhenhua, Zhao Mingjin, etc. A New Method of Fracture Detection Based on GHT, August 2009, Petroleum Geophysical Exploration, Volume 44, Issue 4, Pages 442-444

技术方案：Technical solutions:

②核心计算公式：(a)采用频率域的广义Hilbert变换(GHT)计算公式，即采用公式(5)、(6)、(7)计算；(b)采用2维保边去噪方法(LUOYi等，2002)，进行预处理计算。②Core calculation formula: (a) use the generalized Hilbert transform (GHT) calculation formula in the frequency domain, that is, use formulas (5), (6), and (7) for calculation; (b) use the 2-dimensional edge-preserving denoising method (LUOYi et al., 2002), for preprocessing calculations.

③计算流程：将2维切片视为2维数组，进行2维保边去噪处理；再将2维数组按照横向或纵向分成若干个1维数组；采用1维频率域GHT分别计算每个1维数组的GHT，得到解析信号，并提取的虚部用于边缘检测；循环计算完所有的1维数组，得到最终的2维切片的边缘检测结果数组。③Calculation process: treat the 2-dimensional slice as a 2-dimensional array, and perform 2-dimensional edge-preserving and denoising processing; then divide the 2-dimensional array into several 1-dimensional arrays horizontally or vertically; use the 1-dimensional frequency domain GHT to calculate each 1 The GHT of the two-dimensional array obtains the analytical signal, and the extracted imaginary part is used for edge detection; after calculating all the one-dimensional arrays, the final array of edge detection results of the two-dimensional slice is obtained.

上述4种技术方案的缺点：The shortcoming of above-mentioned 4 kinds of technical schemes:

(1)4种方案都是针对三维时间切片(即2维水平时间切片)，没有针对三维地震资料的实际地层(实际地层往往非水平)的处理方案；(1) The four schemes are all for 3D time slices (that is, 2D horizontal time slices), and there is no processing scheme for the actual strata of 3D seismic data (the actual stratum is often non-horizontal);

(2)3种方案(第1，2，4种)都是将2维切片分解为若干个1维数组，即只针对2维切片的1个方向进行检测(横向或纵向)，其边缘检测的精度较低；(2) The three schemes (1, 2, and 4) all decompose the 2-dimensional slice into several 1-dimensional arrays, that is, only detect one direction (horizontal or vertical) of the 2-dimensional slice, and the edge detection The accuracy is low;

(3)第1种方法——陈学华等(2008)的方案采用频率域的GHT方法计算，计算效率相对较低，且抑制噪声能力较差(没有计算流程中考虑抑制噪声)；(3) The first method - the scheme of Chen Xuehua et al. (2008) uses the GHT method in the frequency domain to calculate, the calculation efficiency is relatively low, and the ability to suppress noise is poor (no noise suppression is not considered in the calculation process);

(4)第2种方法——谢静等(2007)的方案采用加窗HT方法计算，计算效率高，但是抑制噪声能力差(没有计算流程中考虑抑制噪声)；(4) The second method - the scheme of Xie Jing et al. (2007) uses the windowed HT method for calculation, which has high calculation efficiency, but poor noise suppression ability (no noise suppression is not considered in the calculation process);

(5)第3种方法——谢静(2007)的方案在计算流程中考虑了抑制噪声和对二维切片的2个方向的处理，但是仍采用常规的HT进行计算，计算精度低；(5) The third method - Xie Jing's (2007) scheme considers noise suppression and the processing of two directions of two-dimensional slices in the calculation process, but still uses conventional HT for calculation, and the calculation accuracy is low;

(6)第4种方法——熊晓军等(2009)的方案在计算流程中考虑了抑制噪声，但是仍采用频率域的GHT方法计算，计算效率相对较低。(6) The fourth method - the scheme of Xiong Xiaojun et al. (2009) considers noise suppression in the calculation process, but still uses the GHT method in the frequency domain for calculation, and the calculation efficiency is relatively low.

现有缺少解决三维地震数据的边缘检测方法；边缘检测的计算效率低，计算精度低。There is a lack of edge detection methods for 3D seismic data; the calculation efficiency and accuracy of edge detection are low.

发明内容Contents of the invention

本发明实施例的目的在于提供一种基于时间域广义Hilbert变换的三维地震资料快速边缘检测方法，旨在解决现有缺少解决三维地震数据的边缘检测方法；边缘检测的计算效率低，计算精度低的问题。The purpose of the embodiments of the present invention is to provide a fast edge detection method for 3D seismic data based on the generalized Hilbert transform in the time domain, aiming at solving the lack of existing edge detection methods for 3D seismic data; the calculation efficiency of edge detection is low, and the calculation accuracy is low The problem.

本发明实施例是这样实现的，一种基于时间域广义Hilbert变换的三维地震资料快速边缘检测方法，该基于时间域广义Hilbert变换的三维地震资料快速边缘检测方法包括以下步骤：The embodiment of the present invention is achieved in this way, a method for fast edge detection of 3D seismic data based on time domain generalized Hilbert transform, the method for fast edge detection of 3D seismic data based on time domain generalized Hilbert transform comprises the following steps:

以三维地震资料的目的层段为研究对象，输入三维地震数据，采用三维等时切片的计算方法，计算目的层段的二维等时切片；Taking the target interval of 3D seismic data as the research object, input 3D seismic data, and use the calculation method of 3D isochronous slice to calculate the 2D isochronous slice of the target interval;

计算时间域的2个1维的希氏因子，针对第1个等时切片，进行2维保边去噪处理，得到去噪后的数组；Calculate two 1-dimensional Hisfactors in the time domain, and perform 2-dimensional edge-preserving denoising processing on the first isochronous slice to obtain a denoised array;

选取长度为65的希氏因子，阶次为1的情况，进行计算选取长度为33的希氏因子，阶次为2的情况，直至获得当前等时切片的边缘检测结果；Select a Hisfactory factor with a length of 65 and an order of 1, perform calculations and select a Hisfactory factor with a length of 33 and an order of 2, until the edge detection result of the current isochronous slice is obtained;

综合阶次1与阶次2的计算结果，得到当前切片的最终的边缘检测结果；Combine the calculation results of order 1 and order 2 to obtain the final edge detection result of the current slice;

更换等时切片号，计算所有的等时切片；按照等时切片在目的层段中的位置，存放每个等时切片的边缘检测的结果，得到最终的目的层段的三维边缘检测结果。Replace the isochronous slice number, calculate all isochronous slices; store the edge detection result of each isochronous slice according to the position of the isochronous slice in the target layer, and obtain the final three-dimensional edge detection result of the target layer.

进一步，该基于时间域广义Hilbert变换的三维地震资料快速边缘检测方法具体包括以下步骤：Further, the fast edge detection method for 3D seismic data based on time-domain generalized Hilbert transform specifically includes the following steps:

步骤一，以三维地震资料的目的层段为研究对象，输入三维地震数据，目标层段的顶界层位和底界层位，目标层段为D，目的层段的顶界层位为Hor_1，底界层位为Hor_2；Step 1: Take the target interval of the 3D seismic data as the research object, input the 3D seismic data, the top boundary horizon and the bottom boundary horizon of the target interval, the target interval is D, and the top boundary horizon of the target interval is Hor_1 , the bottom layer is Hor_2;

步骤二，采用三维等时切片的计算方法，计算目的层段的二维等时切片：当目的层段的顶界层位Hor_1与底界层位Hor_2平行时，二维等时切片的总数(N1)：Step 2, using the calculation method of three-dimensional isochronous slices to calculate the two-dimensional isochronous slices of the target interval: when the top boundary horizon Hor_1 of the target interval is parallel to the bottom boundary Hor_2, the total number of two-dimensional isochronous slices ( N1):

N1＝L/dt(8)N1=L/dt(8)

式中，L代表Hor_2与Hor_1的时间距离，dt代表三维地震数据的时间采样间隔，当目的层段的顶界层位Hor_1与底界层位Hor_2不平行时，二维等时切片的个数(N1)：In the formula, L represents the time distance between Hor_2 and Hor_1, dt represents the time sampling interval of 3D seismic data, when the top boundary Hor_1 of the target layer is not parallel to the bottom boundary Hor_2, the number of 2D isochronous slices (N1):

N1＝Max_L/dt(9)N1=Max_L/dt(9)

式中，Max_L代表Hor_2与Hor_1的最大时间距离，dt代表三维地震数据的时间采样间隔，等时切片共15个，其中第2、3、4次等分的切片在Inline或Crossline方向非完全采样，统一将切片进行补0处理，补0处理的网格点，不参与后续的计算，保证等时切片的完整，即所有等时切片在纵向和横向的长度一致；In the formula, Max_L represents the maximum time distance between Hor_2 and Hor_1, dt represents the time sampling interval of 3D seismic data, and there are 15 isochronous slices, among which the 2nd, 3rd, and 4th equally divided slices are incompletely sampled in the Inline or Crossline direction , uniformly fill the slices with 0s, and fill the grid points with 0s, and do not participate in subsequent calculations to ensure the integrity of the isochronous slices, that is, the lengths of all isochronous slices in the vertical and horizontal directions are the same;

步骤三，计算时间域的2个1维的希氏因子和其计算公式相同，Step 3, calculate two 1-dimensional Hiscofactors in the time domain and Its calculation formula is the same,

$\overset{- -}{{h h}_{0101}} ((n no)) = = \overset{- -}{{h h}_{0202}} ((n no)) \{\begin{matrix} 00 & n no = = 22 m m \\ \frac{22}{π π n no} & n no = = 22 m m + + 11 \end{matrix} - - - - - - ((1010))$

式中m为自然数，希氏因子的总长度设定为65，希氏因子的总长度设定为33；并对时间域的希氏因子和施加高斯窗w(n)，窗数组的长度为65；In the formula, m is a natural number, His factor The total length is set to 65, the Hisfactor The total length of is set to 33; and the Hisfactor of the time domain and Apply a Gaussian window w(n), the length of the window array is 65;

$\overset{- -}{{h h}_{11}} ((n no)) = = {\overset{&OverBar; &OverBar;}{h h}}_{0101} ((n no)) \times \times w w ((n no)),, n no = = 11,, ... ...,, 6565 - - - - - - ((1111))$

$\overset{- -}{{h h}_{22}} ((n no)) = = {\overset{&OverBar; &OverBar;}{h h}}_{0202} ((n no)) \times \times w w ((n no)),, n no = = 11,, ... ...,, 3333 - - - - - - ((1212))$

得到加窗处理后的2个希氏因子和长度为65的希氏因子相当于GHT中的阶次为1的情况，即常规的HT；长度为33的希氏因子相当于GHT中的阶次为2的情况；Get 2 Hisfactors after windowing and A Historic factor with a length of 65 is equivalent to the case where the order in GHT is 1, that is, conventional HT; a Hisfactor with a length of 33 is equivalent to the case where the order in GHT is 2;

步骤四，针对第1个等时切片，记为数组Slice(nx,ny)，nx的取值从1至Nx，Nx为三维地震资料的InLine总数，ny的取值从1至Ny，Ny为三维地震资料的CrossLine总数，进行2维保边去噪处理，得到去噪后的数组Slice_N(nx,ny)；Step 4: For the first isochronous slice, record it as the array Slice(nx,ny), the value of nx is from 1 to Nx, Nx is the total number of InLine of 3D seismic data, the value of ny is from 1 to Ny, and Ny is The total number of CrossLines of the 3D seismic data is processed by 2D edge-preserving denoising to obtain the denoised array Slice_N(nx,ny);

步骤五，选取长度为65的希氏因子阶次为1的情况，进行计算；具体方法为：Step five, select the Hisfactory factor with a length of 65 When the order is 1, calculate it; the specific method is:

第一步，针对2维数组Slice_N(nx,ny)(nx＝1,...,Nx；ny＝1,...,Ny)，在Inline方向进行循环计算：In the first step, for the 2-dimensional array Slice_N(nx,ny) (nx=1,...,Nx; ny=1,...,Ny), the loop calculation is performed in the Inline direction:

提取Inline号＝1的1维数组，即Slice_N(1,ny),ny＝1,...,Ny；令该1维数组为Slice_X(ny),ny＝1,...,Ny；Extract the 1-dimensional array of Inline number=1, namely Slice_N(1,ny), ny=1,...,Ny; let this 1-dimensional array be Slice_X(ny), ny=1,...,Ny;

计算Slice_X(ny)与的褶积，Calculate Slice_X(ny) and the convolution,

$S S l l i i c c e e__X x__I I ((n no y the y)) = = S S l l i i c c e e__X x ((n no y the y)) * * {\overset{&OverBar; &OverBar;}{h h}}_{0101} ((n no)) - - - - - - ((1212))$

式中Slice_X_I(ny)为Slice_X(ny)的Hilbert变换，1维数组Slice_X_I(ny)的长度为Ny(1维数组Slice_X_I(ny)与1维数组褶积的结果数组长度为(Ny+65-1)，对结果数组舍弃前后各32个点，即只选取中间的Ny个点)；并将1维数组Slice_X_I(ny)写入2维数组Slice_X_E(1,ny)ny＝1,...,Ny；In the formula, Slice_X_I(ny) is the Hilbert transformation of Slice_X(ny), and the length of the 1-dimensional array Slice_X_I(ny) is Ny (1-dimensional array Slice_X_I(ny) and 1-dimensional array The length of the convolution result array is (Ny+65-1), discarding 32 points before and after the result array, that is, only selecting the Ny points in the middle); and writing the 1-dimensional array Slice_X_I(ny) into the 2-dimensional array Slice_X_E (1,ny)ny=1,...,Ny;

更换Inline号，重复第一步和第二步，直到所有的Inline号计算完成，得到当前等时切片的边缘检测的结果数组Slice_X_E(nx,ny)nx＝1,...,Nx；ny＝1,...,Ny；Replace the Inline number, repeat the first and second steps until all the Inline numbers are calculated, and get the result array of edge detection of the current isochronous slice Slice_X_E(nx,ny)nx=1,...,Nx;ny= 1,...,Ny;

第二步，针对2维数组Slice_N(nx,ny)，在Crossline方向进行循环计算：In the second step, for the 2-dimensional array Slice_N(nx,ny), perform loop calculations in the Crossline direction:

步骤(1)，提取Crossline号＝1的1维数组，即Slice_N(nx,1),nx＝1,...,Nx；令该1维数组为Slice_Y(nx),nx＝1,...,Nx；Step (1), extract the 1-dimensional array with Crossline number=1, that is, Slice_N(nx,1), nx=1,...,Nx; let the 1-dimensional array be Slice_Y(nx),nx=1,... ., Nx;

计算Slice_Y(nx)与的褶积：Calculate Slice_Y(nx) and The convolution of:

$S S l l i i c c e e__Y Y__I I ((n no x x)) = = S S l l i i c c e e__Y Y ((n no x x)) * * \overset{- -}{{h h}_{0101}} ((n no)) - - - - - - ((1313))$

式中Slice_Y_I(nx)为Slice_Y(nX)的Hilbert变换，1维数组Slice_Y_I(nx)的长度为Nx(1维数组Slice_Y_I(nx)与1维数组褶积的结果数组长度为(Nx+65-1)，通过舍弃前后各32个点，即只选取中间的Nx个点)；并将1维数组Slice_Y_I(ny)写入2维数组Slice_Y_E(nx,1)nx＝1,...,Nx；In the formula, Slice_Y_I(nx) is the Hilbert transformation of Slice_Y(nX), and the length of the 1-dimensional array Slice_Y_I(nx) is Nx (1-dimensional array Slice_Y_I(nx) and 1-dimensional array The length of the convolution result array is (Nx+65-1), by discarding 32 points before and after, that is, only selecting the middle Nx points); and writing the 1-dimensional array Slice_Y_I(ny) into the 2-dimensional array Slice_Y_E(nx ,1)nx=1,...,Nx;

步骤(2)，更换Crossline号，重复步骤(1)，直到所有的Crossline号计算完成，从而得到当前等时切片的边缘检测的结果数组Slice_Y_E(nx,ny)nx＝1,...,Nx；ny＝1,...,Ny；Step (2), replace the Crossline number, repeat step (1), until all the Crossline numbers are calculated, so as to obtain the result array Slice_Y_E(nx,ny)nx=1,...,Nx of the edge detection of the current isochronous slice ;ny=1,...,Ny;

第三步，计算当前等时切片的边缘检测结果：The third step is to calculate the edge detection result of the current isochronous slice:

$S S l l i i c c e e__E E. 11 ((n no x x,, n no y the y)) = = \sqrt{S S l l i i c c e e__X x__{E E.}^{22} ((n no x x,, n no y the y)) + + S S l l i i c c e e__Y Y__{E E.}^{22} ((n no x x,, n no y the y))},, n no x x = = 11,, ... ...,, N N x x;; n no y the y = = 11,, ... ...,, N N y the y - - - - - - ((1515))$

步骤六，选取长度为33的希氏因子(阶次为2的情况)，重复第一步至第三步，直至获得当前等时切片的边缘检测结果Slice_E2(nx,ny)；Step 6, select the Hisfactory factor with a length of 33 (the case where the order is 2), repeat the first step to the third step until the edge detection result Slice_E2(nx,ny) of the current isochronous slice is obtained;

步骤七，综合阶次1与阶次2的计算结果，得到当前切片的最终的边缘检测结果；Step 7, integrating the calculation results of order 1 and order 2 to obtain the final edge detection result of the current slice;

$S S l l i i c c e e__E E. ((n no x x,, n no y the y)) = = \sqrt{S S l l i i c c e e__E E. 11^{22} ((n no x x,, n no y the y)) + + S S l l i i c c e e__E E. 22^{22} ((n no x x,, n no y the y))},, n no x x = = 11,, ... ...,, N N x x;; n no y the y = = 11,, ... ...,, N N y the y - - - - - - ((1616))$

步骤八，更换等时切片号，重复步骤四～步骤七，直到所有的等时切片计算完成；Step 8, replace the isochronous slice number, repeat steps 4 to 7 until all isochronous slice calculations are completed;

步骤九，按照等时切片在目的层段中的位置，存放每个等时切片的边缘检测的结果，得到最终的目的层段的三维边缘检测结果。Step 9: Store the edge detection result of each isochronous slice according to the position of the isochronous slice in the target layer, and obtain the final 3D edge detection result of the target layer.

进一步，在步骤五中，选取长度为65的希氏因子阶次为1的情况，进行计算的具体方法为：Further, in step five, select the Hisfactory factor with a length of 65 When the order is 1, the specific calculation method is as follows:

$S S l l i i c c e e__E E. 11 ((n no x x,, n no y the y)) = = \sqrt{S S l l i i c c e e__X x__{E E.}^{22} ((n no x x,, n no y the y)) + + S S l l i i c c e e__Y Y__{E E.}^{22} ((n no x x,, n no y the y))},, n no x x = = 11,, ... ...,, N N x x;; n no y the y = = 11,, ... ...,, N N y the y - - - - - - ((1515)) . .$

进一步，在第一步中，提取Inline号＝1的1维数组，即Slice_N(1,ny),ny＝1,...,Ny；令该1维数组为Slice_X(ny),ny＝1,...,Ny；具体方法为：Further, in the first step, extract the 1-dimensional array with Inline number=1, that is, Slice_N(1,ny),ny=1,...,Ny; let the 1-dimensional array be Slice_X(ny),ny=1 ,...,Ny; the specific method is:

步骤一，计算Slice_X(ny)与的褶积，Step 1, calculate Slice_X(ny) and the convolution,

步骤二，更换Inline号，所有的Inline号计算完成，得到当前等时切片的边缘检测的结果数组Slice_X_E(nx,ny)nx＝1,...,Nx；ny＝1,...,Ny。Step 2, replace the Inline number, all Inline numbers are calculated, and the result array Slice_X_E(nx,ny)nx=1,...,Nx; ny=1,...,Ny of the edge detection of the current isochronous slice is obtained .

进一步，在第二步中，针对2维数组Slice_N(nx,ny)，在Crossline方向进行循环计算具体方法为：Further, in the second step, for the 2-dimensional array Slice_N(nx,ny), the specific method of loop calculation in the Crossline direction is as follows:

计算Slice_Y(nx)与的褶积：Calculate Slice_Y(nx) and The convolution of:

步骤(2)，更换Crossline号，重复步骤(1)，直到所有的Crossline号计算完成，从而得到当前等时切片的边缘检测的结果数组Slice_Y_E(nx,ny)nx＝1,...,Nx；ny＝1,...,Ny。Step (2), replace the Crossline number, repeat step (1), until all the Crossline numbers are calculated, so as to obtain the result array Slice_Y_E(nx,ny)nx=1,...,Nx of the edge detection of the current isochronous slice ;ny=1,...,Ny.

本发明提供的基于时间域广义Hilbert变换的三维地震资料快速边缘检测方法，采用1维时间域广义HT处理三维地震数据，包含了加窗处理，又包含了阶次处理，在计算流程中包含了噪声抑制处理，并从2维等时切片的2个方向进行同时处理，具有较高的边缘检测的精度；本发明的核心计算公式是时间域的HT，计算效率高；采用等时切片处理三维地震数据，且相邻等时切片的间隔等于1个时间采样间隔，确保了所有的样点都位于各自的等时切片上，不仅最接近真实的地下地质情况，而且有利于边缘检测结果的还原；采用2个不同长度的希氏因子进行处理，相当于频率域中的2个阶次的情况，并综合了2个阶次的计算结果，计算精度高。本发明针对三维地震资料进行处理，获得目标层段的三维边缘检测数据体，用于裂缝发育带或缝洞发育带的识别，有效地指导油气地球物理勘探的储层预测；给定三维地震数据体和目的层段的顶界面层位和底界面层位，获得三维的边缘检测数据体；进行计算流程的优化组合，同时考虑噪声抑制、计算效率和计算精度。The method for fast edge detection of 3D seismic data based on time-domain generalized Hilbert transform provided by the present invention uses 1-dimensional time-domain generalized HT to process 3D seismic data, including windowing processing and order processing, and includes in the calculation process Noise suppression processing, and simultaneous processing from two directions of 2-dimensional isochronous slices, has higher edge detection accuracy; the core calculation formula of the present invention is HT in time domain, and the calculation efficiency is high; using isochronous slices to process three-dimensional Seismic data, and the interval between adjacent isochronous slices is equal to one time sampling interval, ensuring that all sample points are located on their respective isochronous slices, which is not only the closest to the real underground geological situation, but also conducive to the restoration of edge detection results ; Two Hisfactors of different lengths are used for processing, which is equivalent to the case of two orders in the frequency domain, and the calculation results of the two orders are integrated, and the calculation accuracy is high. The invention processes the three-dimensional seismic data to obtain the three-dimensional edge detection data body of the target interval, which is used for the identification of the fracture development zone or the fracture-vug development zone, and effectively guides the reservoir prediction of oil and gas geophysical exploration; the given three-dimensional seismic data The top interface layer and bottom interface layer of the volume and the target interval are obtained to obtain a three-dimensional edge detection data volume; the optimization combination of the calculation process is carried out, and noise suppression, calculation efficiency and calculation accuracy are considered at the same time.

附图说明Description of drawings

图1是本发明实施例提供的基于时间域广义Hilbert变换的三维地震资料快速边缘检测方法流程图；Fig. 1 is the flow chart of the fast edge detection method of 3D seismic data based on time-domain generalized Hilbert transform provided by the embodiment of the present invention;

图2是本发明实施例提供的三维目的层段的示意图；Fig. 2 is a schematic diagram of a three-dimensional target interval provided by an embodiment of the present invention;

图3是本发明实施例提供的平行地层的等时切片的计算示意图；Fig. 3 is a schematic diagram of calculation of isochronous slices of parallel formations provided by an embodiment of the present invention;

图4是本发明实施例提供的非平行地层的等时切片的计算示意图。Fig. 4 is a schematic diagram of calculation of an isochronous slice of a non-parallel formation provided by an embodiment of the present invention.

具体实施方式detailed description

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。In order to make the object, technical solution and advantages of the present invention more clear, the present invention will be further described in detail below in conjunction with the examples. It should be understood that the specific embodiments described here are only used to explain the present invention, not to limit the present invention.

下面结合附图及具体实施例对本发明的应用原理作进一步描述。The application principle of the present invention will be further described below in conjunction with the accompanying drawings and specific embodiments.

如图1所示，本发明实施例的基于时间域广义Hilbert变换的三维地震资料快速边缘检测方法包括以下步骤：As shown in Figure 1, the method for fast edge detection of 3D seismic data based on time-domain generalized Hilbert transform of the embodiment of the present invention comprises the following steps:

S101：以三维地震资料的目的层段为研究对象，输入三维地震数据，采用三维等时切片的计算方法，计算目的层段的二维等时切片；S101: Taking the target interval of the 3D seismic data as the research object, input the 3D seismic data, and adopt the calculation method of the 3D isochronous slice to calculate the 2D isochronous slice of the target interval;

S102：计算时间域的2个1维的希氏因子，针对第1个等时切片，进行2维保边去噪处理，得到去噪后的数组；S102: Calculate two 1-dimensional Hisfactors in the time domain, and perform 2-dimensional edge-preserving denoising processing on the first isochronous slice to obtain a denoised array;

S103：选取长度为65的希氏因子(阶次为1的情况)，进行计算选取长度为33的希氏因子(阶次为2的情况)，直至获得当前等时切片的边缘检测结果；S103: Select a Historia factor with a length of 65 (in the case of order 1), perform calculations and select a Hisfactor with a length of 33 (in the case of order 2), until the edge detection result of the current isochronous slice is obtained;

S104：综合阶次1与阶次2的计算结果，得到当前切片的最终的边缘检测结果，；S104: Integrate the calculation results of order 1 and order 2 to obtain the final edge detection result of the current slice;

S105：更换等时切片号，计算所有的等时切片；按照等时切片在目的层段中的位置，存放每个等时切片的边缘检测的结果，得到最终的目的层段的三维边缘检测结果。S105: Change the isochronous slice number, calculate all isochronous slices; store the edge detection result of each isochronous slice according to the position of the isochronous slice in the target layer, and obtain the final 3D edge detection result of the target layer .

本发明的具体步骤包括：Concrete steps of the present invention include:

步骤一，以三维地震资料的目的层段为研究对象，输入三维地震数据(X轴为Inline测线，Y轴为CrossLine测线，Z轴为时间)，目标层段的顶界层位和底界层位，示意图见图2(图2中的目标层段为D，目的层段的顶界层位为Hor_1，底界层位为Hor_2)；Step 1. Taking the target layer of 3D seismic data as the research object, input 3D seismic data (X-axis is the Inline survey line, Y-axis is the CrossLine survey line, Z-axis is time), the top boundary layer and the bottom layer of the target layer Boundary horizon, schematic diagram is shown in Fig. 2 (the target stratum in Fig. 2 is D, and the top boundary stratum of purpose stratum is Hor_1, and bottom boundary stratum is Hor_2);

步骤二，采用三维等时切片的计算方法(ZengHongliu，2010)，计算目的层段的二维等时切片：当目的层段的顶界层位Hor_1与底界层位Hor_2平行时，二维等时切片的总数(N1)：Step 2, using the calculation method of 3D isochronous slice (ZengHongliu, 2010), calculate the 2D isochronous slice of the target interval: when the top boundary Hor_1 of the target interval is parallel to the bottom boundary Hor_2, the 2D isochronous slice The total number of time slices (N1):

N1＝L/dt(8)N1=L/dt(8)

式中，L代表Hor_2与Hor_1的时间距离(单位：ms)，dt代表三维地震数据的时间采样间隔(单位：ms)，等时切片的计算示意图见图3(图3中的网格代表纵向上的1个时间采样点，横向上的1个Inline号或Crossline号)；当目的层段的顶界层位Hor_1与底界层位Hor_2不平行时，二维等时切片的个数(N1)：In the formula, L represents the time distance between Hor_2 and Hor_1 (unit: ms), dt represents the time sampling interval of 3D seismic data (unit: ms), and the calculation schematic diagram of isochronous slice is shown in Fig. 3 (the grid in Fig. 3 represents the longitudinal 1 time sampling point on the top, 1 Inline number or Crossline number on the horizontal direction); when the top boundary Hor_1 of the target interval is not parallel to the bottom boundary Hor_2, the number of two-dimensional isochronous slices (N1 ):

N1＝Max_L/dt(9)N1=Max_L/dt(9)

式中，Max_L代表Hor_2与Hor_1的最大时间距离(单位：ms)，dt代表三维地震数据的时间采样间隔(单位：ms)，等时切片的计算示意图见图4(图4中的网格代表纵向上的1个时间采样点，横向上的1个Inline号或Crossline号；图4中顶界层位与底界层位的最大时间距离为16个时间采用间隔，通过4次等分达到相邻等时切片的时间间隔都为1个时间采用间隔)，图4中的等时切片共15个，其中第2、3、4次等分的切片在Inline或Crossline方向非完全采样(仅有部分网格节点上有值)，在实际计算中，统一将这些切片进行补0处理(补0处理的网格点，不参与后续的计算)，保证等时切片的完整(即所有等时切片在纵向和横向的长度一致)；In the formula, Max_L represents the maximum time distance between Hor_2 and Hor_1 (unit: ms), dt represents the time sampling interval of 3D seismic data (unit: ms), and the calculation diagram of isochronous slice is shown in Fig. 4 (the grid in Fig. 4 represents One time sampling point in the vertical direction, and one Inline number or Crossline number in the horizontal direction; the maximum time distance between the top boundary horizon and the bottom boundary horizon in Figure 4 is 16 time intervals, and the corresponding The time interval between adjacent isochronous slices is 1 time interval), and there are 15 isochronous slices in Figure 4, of which the 2nd, 3rd, and 4th equally divided slices are not completely sampled in the Inline or Crossline direction (only There are values on some grid nodes), in the actual calculation, these slices are uniformly filled with 0 (grid points with 0 filled, do not participate in subsequent calculations), to ensure the integrity of isochronous slices (that is, all isochronous slices same length in longitudinal and transverse directions);

式中m为自然数，希氏因子的总长度设定为65，希氏因子的总长度设定为33；并对时间域的希氏因子和施加高斯窗w(n)(窗数组的长度为65)；In the formula, m is a natural number, His factor The total length is set to 65, the Hisfactor The total length of is set to 33; and the Hisfactor of the time domain and Apply a Gaussian window w(n) (the length of the window array is 65);

得到加窗处理后的2个希氏因子和(长度为65的希氏因子相当于GHT中的阶次为1的情况，即常规的HT；长度为33的希氏因子相当于GHT中的阶次为2的情况)Get 2 Hisfactors after windowing and (A Hisfactor with a length of 65 is equivalent to the case where the order in GHT is 1, that is, conventional HT; a Hisfactor with a length of 33 is equivalent to the case where the order in GHT is 2)

步骤四，针对第1个等时切片(记为数组Slice(nx,ny)，nx的取值从1至Nx，Nx为三维地震资料的InLine总数，ny的取值从1至Ny，Ny为三维地震资料的CrossLine总数)，进行2维保边去噪处理(LUOYi等，2002)，得到去噪后的数组Slice_N(nx,ny)；Step 4, for the first isochronous slice (recorded as the array Slice(nx,ny), the value of nx is from 1 to Nx, Nx is the total number of InLine of 3D seismic data, the value of ny is from 1 to Ny, and Ny is The total number of CrossLines of 3D seismic data), carry out 2-dimensional edge-preserving denoising processing (LUOYi et al., 2002), and obtain the array Slice_N(nx,ny) after denoising;

步骤五，选取长度为65的希氏因子(阶次为1的情况)，进行计算；具体方法为：Step five, select the Hisfactory factor with a length of 65 (the case where the order is 1), perform the calculation; the specific method is:

提取Inline号＝1的1维数组，即Slice_N(1,ny),ny＝1,...,Ny；令该1维数组为Slice_X(ny),ny＝1,...,Ny；Extract the 1-dimensional array of Inline number=1, that is, Slice_N(1,ny), ny=1,...,Ny; let the 1-dimensional array be Slice_X(ny), ny=1,...,Ny;

计算Slice_X(ny)与的褶积，Calculate Slice_X(ny) and the convolution,

计算Slice_Y(nx)与的褶积：Calculate Slice_Y(nx) and The convolution of:

步骤(2)，更换Crossline号，重复步骤(1)，步骤2下的内容(见批注的起始位置)直到所有的Crossline号计算完成，从而得到当前等时切片的边缘检测的结果数组Slice_Y_E(nx,ny)nx＝1,...,Nx；ny＝1,...,Ny；Step (2), replace the Crossline number, repeat step (1), the content under step 2 (see the starting position of the comment) until all the Crossline numbers are calculated, so as to obtain the result array Slice_Y_E of the edge detection of the current isochronous slice ( nx,ny) nx=1,...,Nx; ny=1,...,Ny;

步骤八，更换等时切片号，重复步骤四～步骤七，直到所有的等时切片计算完成；Step 8, replace the number of isochronous slices, repeat steps 4 to 7, until all isochronous slices are calculated;

本发明采用1维时间域广义HT处理三维地震数据，与频率域GHT类似，既包含了加窗处理，又包含了阶次处理；此外，在计算流程中包含了噪声抑制处理，并从2维等时切片的2个方向进行同时处理，具有较高的边缘检测的精度；此外，本发明的核心计算公式是时间域的HT，计算效率高。采用等时切片处理三维地震数据，且相邻等时切片的间隔等于1个时间采样间隔。本发明确保了所有的样点都位于各自的等时切片上，不仅最接近真实的地下地质情况，而且有利于边缘检测结果的还原(即重新构成三维数据体)；采用2个不同长度的希氏因子进行处理，相当于频率域中的2个阶次的情况，并综合了2个阶次的计算结果，计算精度高。The present invention uses generalized HT in 1D time domain to process 3D seismic data, which is similar to GHT in frequency domain, including both windowing processing and order processing; in addition, noise suppression processing is included in the calculation process, and from 2D The two directions of the isochronous slice are processed simultaneously, which has high edge detection accuracy; in addition, the core calculation formula of the present invention is HT in the time domain, and the calculation efficiency is high. Isochronous slices are used to process 3D seismic data, and the interval between adjacent isochronous slices is equal to one time sampling interval. The invention ensures that all sample points are located on their respective isochronous slices, which is not only the closest to the real underground geological conditions, but also facilitates the restoration of edge detection results (that is, reconstructing the three-dimensional data volume); It is equivalent to the case of two orders in the frequency domain, and the calculation results of the two orders are integrated, and the calculation accuracy is high.

以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。The above descriptions are only preferred embodiments of the present invention, and are not intended to limit the present invention. Any modifications, equivalent replacements and improvements made within the spirit and principles of the present invention should be included in the protection of the present invention. within range.

Claims

1. A fast edge detection method for three-dimensional seismic data based on time-domain generalized Hilbert transform, characterized in that, the fast edge detection method for three-dimensional seismic data based on time-domain generalized Hilbert transform specifically comprises the following steps:

Step 1: Take the target interval of the 3D seismic data as the research object, input the 3D seismic data, the top boundary horizon and the bottom boundary horizon of the target interval, the target interval is D, and the top boundary horizon of the target interval is Hor_1 , the bottom layer is Hor_2;

Step 2, using the calculation method of three-dimensional isochronous slices to calculate the two-dimensional isochronous slices of the target interval: when the top boundary horizon Hor_1 of the target interval is parallel to the bottom boundary Hor_2, the total number of two-dimensional isochronous slices ( N1):

N1=L/dt

In the formula, L represents the time distance between Hor_2 and Hor_1, dt represents the time sampling interval of 3D seismic data, when the top boundary Hor_1 of the target layer is not parallel to the bottom boundary Hor_2, the number of 2D isochronous slices (N1):

N1=Max_L/dt

In the formula, Max_L represents the maximum time distance between Hor_2 and Hor_1, dt represents the time sampling interval of 3D seismic data, and there are 15 isochronous slices, among which the 2nd, 3rd, and 4th equally divided slices are incompletely sampled in the Inline or Crossline direction , uniformly fill the slices with 0s, and fill the grid points with 0s, and do not participate in subsequent calculations to ensure the integrity of the isochronous slices, that is, the lengths of all isochronous slices in the vertical and horizontal directions are the same;

Step 3, calculate two 1-dimensional Hiscofactors in the time domain and Its calculation formula is the same,

\overset{- -}{{h h}_{0101}} ((n no)) = = \overset{- -}{{h h}_{0202}} ((n no)) \{\begin{matrix} 00 & n no = = 22 m m \\ \frac{22}{ω ω 22} & n no = = 22 m m + + 11 \end{matrix}

In the formula, m is a natural number, His factor The total length is set to 65, the Hisfactor The total length of is set to 33; and the Hisfactor of the time domain and Apply a Gaussian window w(n), the length of the window array is 65;

\overset{- -}{{h h}_{11}} ((n no)) = = {\overset{&OverBar; &OverBar;}{h h}}_{0101} ((n no)) \times \times w w ((n no)),, n no = = 11,, ... ...,, 6565

\overset{- -}{{h h}_{22}} ((n no)) = = {\overset{&OverBar; &OverBar;}{h h}}_{0202} ((n no)) \times \times w w ((n no)),, n no = = 11,, ... ...,, 3333

Get 2 Hisfactors after windowing and A Historic factor with a length of 65 is equivalent to the case where the order in GHT is 1, that is, conventional HT; a Hisfactor with a length of 33 is equivalent to the case where the order in GHT is 2;

Step 4: For the first isochronous slice, record it as the array Slice(nx,ny), the value of nx is from 1 to Nx, Nx is the total number of InLine of 3D seismic data, the value of ny is from 1 to Ny, and Ny is The total number of CrossLine of the 3D seismic data is processed by 2D denoising with edge preservation, and the denoised array Slice_N(nx,ny) is obtained;

Step five, select the Hisfactory factor with a length of 65 When the order is 1, calculate it; the specific method is:

In the first step, for the 2-dimensional array Slice_N(nx,ny) (nx=1,...,Nx; ny=1,...,Ny), the loop calculation is performed in the Inline direction:

Extract the 1-dimensional array of Inline number=1, that is, Slice_N(1,ny), ny=1,...,Ny; let the 1-dimensional array be Slice_X(ny), ny=1,...,Ny;

Calculate Slice_X(ny) and the convolution,

S S l l i i c c e e__X x__I I ((n no y the y)) = = S S l l i i c c e e__X x ((n no y the y)) * * {\overset{&OverBar; &OverBar;}{h h}}_{0101} ((n no))

In the formula, Slice_X_I(ny) is the Hilbert transformation of Slice_X(ny), and the length of the 1-dimensional array Slice_X_I(ny) is Ny (1-dimensional array Slice_X_I(ny) and 1-dimensional array The length of the convolution result array is (Ny+65-1), discarding 32 points before and after the result array, that is, only selecting the Ny points in the middle); and writing the 1-dimensional array Slice_X_I(ny) into the 2-dimensional array Slice_X_E (1,ny)ny=1,...,Ny;

Replace the Inline number, repeat the first step until all Inline numbers are calculated, and get the edge detection result array of the current isochronous slice Slice_X_E(nx,ny)nx=1,...,Nx;ny=1,.. .,Ny;

In the second step, for the 2-dimensional array Slice_N(nx,ny), perform loop calculations in the Crossline direction:

Step (1), extract the 1-dimensional array with Crossline number=1, that is, Slice_N(nx,1), nx=1,...,Nx; let the 1-dimensional array be Slice_Y(nx),nx=1,... ., Nx;

Calculate Slice_Y(nx) and The convolution of:

S S l l i i c c e e__Y Y__I I ((n no x x)) = = S S l l i i c c e e__Y Y ((n no x x)) * * \overset{- -}{{h h}_{0101}} ((n no))

In the formula, Slice_Y_I(nx) is the Hilbert transformation of Slice_Y(nX), and the length of the 1-dimensional array Slice_Y_I(nx) is Nx (1-dimensional array Slice_Y_I(nx) and 1-dimensional array The length of the convolution result array is (Nx+65-1), by discarding 32 points before and after, that is, only selecting the middle Nx points); and writing the 1-dimensional array Slice_Y_I(ny) into the 2-dimensional array Slice_Y_E(nx ,1)nx=1,...,Nx;

Step (2), replace the Crossline number, repeat step (1), until all the Crossline numbers are calculated, so as to obtain the result array Slice_Y_E(nx,ny)nx=1,...,Nx of the edge detection of the current isochronous slice ;ny=1,...,Ny;

The third step is to calculate the edge detection result of the current isochronous slice:

S S l l i i c c e e__E E. 11 ((n no x x,, n no y the y)) = = \sqrt{S S l l o o c c e e__X x__{E E.}^{22} ((n no x x,, n no y the y)) + + S S l l i i c c e e__Y Y__{E E.}^{22} ((n no x x,, n no y the y))},, n no x x = = 11,, ... ...,, N N x x;; n no y the y = = 11,, ... ... N N y the y

Step 6, select the Hisfactory factor with a length of 33 In the case of order 2, repeat the first step to the third step until the edge detection result Slice_E2(nx,ny) of the current isochronous slice is obtained;

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

S S l l i i c c e e__E E. ((n no x x,, n no y the y)) = = \sqrt{S S l l i i c c e e__E E. 11^{22} ((n no x x,, n no y the y)) + + S S l l i i c c e e__E E. 22^{22} ((n no x x,, n no y the y))},, n no x x = = 11,, ... ...,, N N x x;; n no y the y = = 11,, ... ...,, N N y the y

Step 8, replace the isochronous slice number, repeat steps 4 to 7 until all isochronous slice calculations are completed;

Step 9: Store the edge detection result of each isochronous slice according to the position of the isochronous slice in the target layer, and obtain the final 3D edge detection result of the target layer.