基于魔方矩阵的断层检测方法

李 军 张军华 龚明平 杨 勇 杜玉山 武 刚(中国石油大学(华东)地球科学与技术学院，山东青岛 266580；中国石化胜利油田分公司勘探开发研究院，山东东营 257015)

1 引言

断层检测与识别是地震解释的重要工作。多年来，前人在断层检测与识别方面提出了相干及其改进算法[1,2]、曲率类属性[3]、边缘检测及其他如自适应波形匹配追踪算法[4]等众多方法。边缘检测是一类应用广泛的算法，包含的检测算子众多，其中最基本的算子为Roberts[5]、Sobel[6]与Prewitt[7]算子。另外，基于算子的优缺点以及地震资料的特点，有人针对性地提出或者改进一些适用于地震处理、解释的新边缘检测算子。苟量等[8]提出了小波多尺度边缘检测技术，并获得了较好的应用效果； Al-Dossary等[9,10]利用三维Sobel滤波技术检测断层、裂缝及河道等异常体； 王清振等[11]提出了基于高维小波变换的高抗噪性边缘检测技术，能更好地检测小断层、小河道等地质体的边界信息。但大多数边缘检测技术在算法中并未考虑地震数据的倾角、方向信息，不能很好地检测一些复杂区块的断层、裂缝等地质现象。

魔方矩阵又称幻方，具有相同的行数和列数，矩阵中每个元素均不同，每行、每列、每条对角线的元素之和都相等。多年来，研究人员主要将其用于信息隐藏及加密、图像处理等方面。张萌[12]研究了基于魔方矩阵变换置乱的图像加密算法，实现了图像信息的隐藏； 赵杨[13]在魔方矩阵基础上研究了基于运动矢量的视频信息隐藏算法， 李松斌等[14]提出了一种基于整数 DCT 系数调制及N维魔方矩阵的信息隐藏方法， 刘争艳等[15]提出了基于 H.264/AVC 压缩域的视频信息隐藏算法，均提高了信息隐藏质量及效率； Kuang[16]提出了联合魔方矩阵快速编码与二进制编码的混合编码图像信息隐藏方法，可以完全解码隐藏信息，且不失真； Chen等[17]利用3×3阶二维魔方矩阵构造检测算子检测图像边缘，取得了较好的检测效果。总的来说，以上算法均利用魔方矩阵构造不同方向的检测算子，进行图像处理，信息隐藏、加密、处理等目的。将上述算法用于地震勘探领域，可以检测断层、裂缝等特殊地质现象[18-20]，在给出分布范围的同时，可定量计算倾向信息，且抗噪性好于常规方法。

2 方法原理

2.1 魔方矩阵算法简介

大多数边缘检测方法的理论基础是计算输入数据每个点的振幅梯度，在一些算子中还需要设定阈值以得到理想的计算结果。振幅梯度各个方向的计算表达式为[21]

(1)

(2)

(3)

一阶导数的3×3阶差分格式写成矩阵形式就是Sobel核的3×3模板，其沿x与y方向的表达式为

(4)

可知式(4)中对角线元素之和均为0。类似地，其他大多数边缘检测算子最后均可写成矩阵形式的模板，且对角线元素之和也为0，从而为引入魔方矩阵构造检测算子进行断层检测提供了参考。

魔方矩阵相比于其他矩阵，其对角线及每行、每列的元素之和均等于魔方值

(5)

式中N为魔方矩阵的阶数。魔方矩阵的最大特征值为魔方值，其余特征值为0或呈正、负共轭成对出现[22]。

以3×3阶魔方矩阵为例，其9个元素为数字1～9，构造的最简单的3×3阶魔方矩阵为

(6)

由上述元素构成的其他3×3阶魔方矩阵是由式(6)旋转或反转得到。将式(6)的所有元素值减去元素中值5，得到新矩阵

(7)

分析式(7)可知其对角线元素之和均为0，其特征类似于典型的边缘检测算子模板矩阵，故将魔方矩阵用于边缘检测是有数学意义的。

考虑到图像一般沿边缘方向变化幅度较小，而边缘两侧变化幅度较大，故将式(7)中元素绝对值最小的对角线[-1 0 1]用于检测边缘延伸方向，元素绝对值最大的列 [4 0 -4]用于检测边缘具体形态，其他矩阵元素取值为0，便得到魔方矩阵边缘检测算子模板

(8)

其他魔方矩阵边缘检测算子模板相当于将式(8)以间隔45°顺时针旋转得到，即

…

(9)

由于实际地震数据为三维，信噪比一般不是很高，且实际断层的连续性并不好，故需将魔方矩阵算法推广到三维，并把魔方矩阵原本的检测指示方向进行互换，即将元素绝对值最大方向用于检测断层倾向，以更好地检测断层。

2.2 算法实现流程

N×N×N(N=3,4,5,…)阶魔方矩阵的具体意义是一样的，N不同导致检测效果及计算效率不同。以3×3×3阶魔方矩阵为例，具体说明基于魔方矩阵的断层检测算法的实现方法。

(1)给定魔方矩阵。使用数字0～26作为元素构成3×3×3阶三维魔方矩阵

(10)

(2)将矩阵所有元素减去矩阵元素集的中值13，得到新矩阵

(11)

(3)类似于二维魔方矩阵检测算子构成方法，将新矩阵中非零元素中绝对值最大、最小以外的其他元素赋值为0，得到断层检测算子模板

(12)

分析式(12)可知，当元素绝对值最大对应方向为断层延伸方向时，元素绝对值最小所在方向正好与其垂直，代表断层形态，这与实际断层相符合，具有实际物理意义。

(4)将式(12)以间隔45°顺时针旋转，得到其他7个方向的检测算子k2,k3,…，k8(图1)。顺时针旋转45°后的方向检测算子为

…

(13)

图1 不同方向检测算子示意图

(5)以地震数据的某点为中心选取一个3×3×3阶数据子体，利用

f(x-a,y-b,z-c)

(14)

将选取的数据子体与步骤(4)得到的8个方向的魔方矩阵检测算子分别进行褶积运算，得到8个方向的褶积运算数据体。式中：Dx为各个方向的检测算子；f(x,y,z)为选取的以求取点为中心构造的三维地震数据子体； (a,b,c)为算子中非零元素集合。

(6)将求取的各个方向褶积运算数据体中绝对值最大的值作为该方向边缘检测结果。

为了得到更为精确的识别方向，可以将3×3×3阶魔方矩阵拓展到5×5×5阶，根据上述算法实现流程得到5×5×5阶的16个方向的魔方矩阵算子，按照类似流程得到断层检测结果。由数字1～125作为元素构成的5×5×5阶魔方矩阵为

(15)

类似地，也可得到其他阶数的魔方矩阵检测算子，但魔方矩阵的阶数越大，计算效率越低，需要根据实际情况，综合考虑选定三维魔方矩阵的阶数。

3 理论模型讨论

设计一个包含多个地垒的三维地质模型，采用30Hz主频的雷克子波，通过正演模拟得到含方差为0.2的高斯噪声的地震模型数据(图2)。文中从识别效果、抗噪性两方面讨论基于魔方矩阵断层检测方法的有效性。

3.1 识别效果

利用基于魔方矩阵的检测方法对三维模型进行处理，得到各个方向的魔方矩阵属性体数据。本文选取地垒模型Xline361剖面上断层线上的O点(图3)分析识别效果。图4为O点的8个方向的魔方矩阵计算数值玫瑰图，可见利用魔方矩阵方法可有效地识别断层，k6方向的检测结果更好。

3.2 抗噪性分析

选择相干算法检测结果与基于魔方矩阵的断层检测方法进行抗噪性分析。图5为断层方向加噪检测结果。可见基于魔方矩阵的断层检测方法的断层边界识别效果更好(图5b)，即与相干算法(图5a)相比，基于魔方矩阵的断层检测方法的抗噪性更好。

图2 三维地质模型正演模拟结果

图3 Xline361横测线剖面(添加高斯噪声)

图4 O点的8个方向的魔方矩阵计算数值玫瑰图

图5 断层方向加噪检测结果(a)相干算法； (b)k6方向

4 实际资料应用

研究区块(永3 工区)位于山东东营凹陷东北部的永安镇油田南部，是胜利探区最复杂的断块型油藏区块之一，目标层断裂众多、结构复杂，大大增加了复杂断块油藏的勘探、开发难度。因此，若能较好地检测该区的断裂结构，对勘探、开发具有重要意义。图6为沿T4地层抽取的原始沿层切片，选定数据点P(图6红色点)进行断层方位信息检测，图7为P点空间位置。根据魔方矩阵算子具体数值可知，当算子方向旋转到k4方向时(图7算子上绿色点所近似的方向)，由于该方向与断层倾向几乎一致，属性值最大，识别效果最好(图8)。

图9为实际资料检测结果。由图可见：①不同方法都能较好地检测主要大断层； ②与P点所在断层倾向接近的k4方向(图7)的算子检测结果(图9c)较其他方向的断层识别结果更清晰；③k2方向检测结果(图9a中红色箭头指示断层倾向)与k2检测算子指示的方向一致，其检测结果较其他方向的断层识别结果更清晰、轮廓更精确； ④k3算子代表水平方向，其检测结果(图9b)类似于上、下延长时间提取属性检测断层的常规方法，断层识别效果较k2(图9a)、k4(图9c)方向差； ⑤与相干算法(图9d)相比，基于魔方矩阵的断层检测方法(图9a～图9c)的抗噪效果更好，即不同方向的检测算子对不同倾向的断层检测结果不一致，当检测算子方向与断层倾向一致时，检测结果更好。

图6 沿T4地层抽取的原始沿层切片

图7 P点空间位置黄色箭头为P点所在断层倾向

图8 P点各个方向检测结果

图9 实际资料检测结果(a)k2； (b)k3； (c)k4； (d)相干算法

5 结束语

断层的检测精度对断块油田剩余油的勘探、开发具有重要意义。与传统的边缘检测方法相比，基于魔方矩阵的断层检测方法利用魔方矩阵可灵活旋转、且方向指向又很明确的特性，可有效检测断层，并定量给出断层发育的倾向信息，较常规相干算法在抗噪性上有较大优势。本文以3×3×3阶的魔方矩阵进行模型测试和应用尝试，也可得到其他阶数的魔方矩阵检测算子，但魔方矩阵的阶数越大，计算效率越低，需要根据实际情况综合考虑选定三维魔方矩阵的阶数。

参考文献

[1] 张军华,王月英,赵勇.C3相干体在断层和裂缝识别中的应用.地震学报,2004,26(5):560-564.

Zhang Junhua,Wang Yueying,Zhao Yong.Application of the third generation of coherent cube in recognizing faults and fractures.Acta Seismologica Sinica,2004,26(5):560-564.

[2] 蔡涵鹏,胡光岷,贺振华等.基于非线性变时窗相干算法的不连续性检测方法.石油地球物理勘探,2016,51(2):371-375.

Cai Hanpeng,Hu Guangmin,He Zhenhua et al.Subtle discontinuity detection with nonlinear variable-time window coherence algorithm.OGP,2016,51(2):371-375.

[3] 王雷,陈海清,陈国文等.应用曲率属性预测裂缝发育带及其产状.石油地球物理勘探,2010,45(6):885- 889.

Wang Lei,Chen Haiqing,Chen Guowen et al.Application of curvature attributes in predicting fracture-developed zone and its orientation.OGP,2010,45(6):885-889.

[4] 邓志文,赵贤正,陈雨红等．自适应波形多道匹配追踪断层识别技术.石油地球物理勘探,2017,52(3):532-537,547.

Deng Zhiwen,Zhao Xianzheng,Chen Yuhong et al.Fault identification based on multichannel adaptive wave forms matching pursuit.OGP,2017,52(3):532-537,547.

[5] Roberts L G.Machine Perception of Three-dimensional Soups[D].Massachusetts Institute of Technology,1963.

[6] Sobel I.Camera models and machine perception:Technical report.DTIC Document,1970.

[7] Prewitt J M.Object enhancement and extraction∥Picture Processing and Psychopictorics，Elsevier,1970,75-149.

[8] 苟量,彭真明.小波多尺度边缘检测及其在裂缝预测中的应用.石油地球物理勘探,2005,40(3):309-313.

Gou Liang,Peng Zhenming.Multi-scale edge detection of wavelet and application in fracture prediction.OGP,2005,40(3):309-313.

[9] Al-Dossary S and Al-Garni K.Fault detection and characterization using a 3D multidirectional Sobel filter.Presented at the SPE Saudi Arabia Section Technical Symposium and Exhibition,2013.

[10] Al-Dossary S.Preconditioning seismic data for chan-nel detection.Interpretation,2014,3(1):SD1-SD4.

[11] 王清振,张金淼,姜秀娣等.基于高维小波变换的高抗噪性边缘检测技术．石油地球物理勘探,2016,51(5):889-893.

Wang Qingzhen,Zhang Jinmiao,Jiang Xiudi et al.A robust denoise edge detection method based on high-dimensional wavelet transform.OGP,2016,51(5):889-893.

[12] 张萌.基于标签的信息嵌入与提取技术研究与实现[学位论文].北京：北京邮电大学，2014，20-22.

[13] 赵杨.基于H.264的大容量视频信息隐藏算法研究[学位论文].四川成都:西南交通大学,2010,45-46.

[14] 李松斌,付江云,刘鹏等.一种基于整数DCT系数调制及N维魔方矩阵的H.264/AVC信息隐藏方法.小型微型计算机系统,2013,34(10):2293-2297.

Li Songbin,Fu Jiangyun,Liu Peng et al.A H.264/AVC information hiding algorithm based on integer DCT coefficients modulation andN-dimensional magic matrix.Journal of Chinese Computer Systems,2013,34(10):2293-2297.

[15] 刘争艳,李絮,陈蕴.基于二维映射关系的视频信息隐藏算法.计算机工程,2010,36(22):225-227.

Liu Zhengyan,Li Xu,Chen Yun.Video information hiding algorithm based on two-dimensional mapping relationship.Computer Engineering,2010,36(22):225-227.

[16] Kuang T L.Hybrid encoding method by assembling the magic-matrix scrambling method and the binary encoding method in image hiding.Optics Communications,2011,284(7):1778-1784.

[17] Chen Zhaoxue,Nie Shengdong.Two efficient edge detecting operators derived from 3×3 magic squares.The International Conference on Wavelet Analysis and Pattern Recognition,IEEE,2007,532-534.

[18] Adetokunbo P,Al-Shuhail A A,Al-Dossary S.3D seismic edge detection using magic squares and cubes.Interpretation,2016,4(3):T271-T280.

[19] Loly P,Cameron I,Trump W et al.Magic square spectra.Linear Algebra and Its Applications,2009,430(10): 2659-2680.

[20] Andrews W S.Magic Squares and Cubes.Cosimo Classics,2004.

[21] Marques O.Practical Image and Video Processing Sing Matlab.John Wiley & Sons,2011.

[22] 单润红,高峰,宋君强.魔方矩阵的特征值分析.高等数学研究,2004,7(4):45,52.