李 軍 張軍華 龔明平 楊 勇 杜玉山 武 剛(中國石油大學(華東)地球科學與技術(shù)學院，山東青島 266580；中國石化勝利油田分公司勘探開發(fā)研究院，山東東營 257015)

1 引言

斷層檢測與識別是地震解釋的重要工作。多年來，前人在斷層檢測與識別方面提出了相干及其改進算法[1,2]、曲率類屬性[3]、邊緣檢測及其他如自適應(yīng)波形匹配追蹤算法[4]等眾多方法。邊緣檢測是一類應(yīng)用廣泛的算法，包含的檢測算子眾多，其中最基本的算子為Roberts[5]、Sobel[6]與Prewitt[7]算子。另外，基于算子的優(yōu)缺點以及地震資料的特點，有人針對性地提出或者改進一些適用于地震處理、解釋的新邊緣檢測算子。茍量等[8]提出了小波多尺度邊緣檢測技術(shù)，并獲得了較好的應(yīng)用效果； Al-Dossary等[9,10]利用三維Sobel濾波技術(shù)檢測斷層、裂縫及河道等異常體； 王清振等[11]提出了基于高維小波變換的高抗噪性邊緣檢測技術(shù)，能更好地檢測小斷層、小河道等地質(zhì)體的邊界信息。但大多數(shù)邊緣檢測技術(shù)在算法中并未考慮地震數(shù)據(jù)的傾角、方向信息，不能很好地檢測一些復(fù)雜區(qū)塊的斷層、裂縫等地質(zhì)現(xiàn)象。

魔方矩陣又稱幻方，具有相同的行數(shù)和列數(shù)，矩陣中每個元素均不同，每行、每列、每條對角線的元素之和都相等。多年來，研究人員主要將其用于信息隱藏及加密、圖像處理等方面。張萌[12]研究了基于魔方矩陣變換置亂的圖像加密算法，實現(xiàn)了圖像信息的隱藏； 趙楊[13]在魔方矩陣基礎(chǔ)上研究了基于運動矢量的視頻信息隱藏算法， 李松斌等[14]提出了一種基于整數(shù) DCT 系數(shù)調(diào)制及N維魔方矩陣的信息隱藏方法， 劉爭艷等[15]提出了基于 H.264/AVC 壓縮域的視頻信息隱藏算法，均提高了信息隱藏質(zhì)量及效率； Kuang[16]提出了聯(lián)合魔方矩陣快速編碼與二進制編碼的混合編碼圖像信息隱藏方法，可以完全解碼隱藏信息，且不失真； Chen等[17]利用3×3階二維魔方矩陣構(gòu)造檢測算子檢測圖像邊緣，取得了較好的檢測效果。總的來說，以上算法均利用魔方矩陣構(gòu)造不同方向的檢測算子，進行圖像處理，信息隱藏、加密、處理等目的。將上述算法用于地震勘探領(lǐng)域，可以檢測斷層、裂縫等特殊地質(zhì)現(xiàn)象[18-20]，在給出分布范圍的同時，可定量計算傾向信息，且抗噪性好于常規(guī)方法。

2 方法原理

2.1 魔方矩陣算法簡介

大多數(shù)邊緣檢測方法的理論基礎(chǔ)是計算輸入數(shù)據(jù)每個點的振幅梯度，在一些算子中還需要設(shè)定閾值以得到理想的計算結(jié)果。振幅梯度各個方向的計算表達式為[21]

(1)

(2)

(3)

一階導數(shù)的3×3階差分格式寫成矩陣形式就是Sobel核的3×3模板，其沿x與y方向的表達式為

(4)

可知式(4)中對角線元素之和均為0。類似地，其他大多數(shù)邊緣檢測算子最后均可寫成矩陣形式的模板，且對角線元素之和也為0，從而為引入魔方矩陣構(gòu)造檢測算子進行斷層檢測提供了參考。

魔方矩陣相比于其他矩陣，其對角線及每行、每列的元素之和均等于魔方值

(5)

式中N為魔方矩陣的階數(shù)。魔方矩陣的最大特征值為魔方值，其余特征值為0或呈正、負共軛成對出現(xiàn)[22]。

以3×3階魔方矩陣為例，其9個元素為數(shù)字1～9，構(gòu)造的最簡單的3×3階魔方矩陣為

(6)

由上述元素構(gòu)成的其他3×3階魔方矩陣是由式(6)旋轉(zhuǎn)或反轉(zhuǎn)得到。將式(6)的所有元素值減去元素中值5，得到新矩陣

(7)

分析式(7)可知其對角線元素之和均為0，其特征類似于典型的邊緣檢測算子模板矩陣，故將魔方矩陣用于邊緣檢測是有數(shù)學意義的。

考慮到圖像一般沿邊緣方向變化幅度較小，而邊緣兩側(cè)變化幅度較大，故將式(7)中元素絕對值最小的對角線[-1 0 1]用于檢測邊緣延伸方向，元素絕對值最大的列 [4 0 -4]用于檢測邊緣具體形態(tài)，其他矩陣元素取值為0，便得到魔方矩陣邊緣檢測算子模板

(8)

其他魔方矩陣邊緣檢測算子模板相當于將式(8)以間隔45°順時針旋轉(zhuǎn)得到，即

…

(9)

由于實際地震數(shù)據(jù)為三維，信噪比一般不是很高，且實際斷層的連續(xù)性并不好，故需將魔方矩陣算法推廣到三維，并把魔方矩陣原本的檢測指示方向進行互換，即將元素絕對值最大方向用于檢測斷層傾向，以更好地檢測斷層。

2.2 算法實現(xiàn)流程

N×N×N(N=3,4,5,…)階魔方矩陣的具體意義是一樣的，N不同導致檢測效果及計算效率不同。以3×3×3階魔方矩陣為例，具體說明基于魔方矩陣的斷層檢測算法的實現(xiàn)方法。

(1)給定魔方矩陣。使用數(shù)字0～26作為元素構(gòu)成3×3×3階三維魔方矩陣

(10)

(2)將矩陣所有元素減去矩陣元素集的中值13，得到新矩陣

(11)

(3)類似于二維魔方矩陣檢測算子構(gòu)成方法，將新矩陣中非零元素中絕對值最大、最小以外的其他元素賦值為0，得到斷層檢測算子模板

(12)

分析式(12)可知，當元素絕對值最大對應(yīng)方向為斷層延伸方向時，元素絕對值最小所在方向正好與其垂直，代表斷層形態(tài)，這與實際斷層相符合，具有實際物理意義。

(4)將式(12)以間隔45°順時針旋轉(zhuǎn)，得到其他7個方向的檢測算子k2,k3,…，k8(圖1)。順時針旋轉(zhuǎn)45°后的方向檢測算子為

…

(13)

圖1 不同方向檢測算子示意圖

(5)以地震數(shù)據(jù)的某點為中心選取一個3×3×3階數(shù)據(jù)子體，利用

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

(14)

將選取的數(shù)據(jù)子體與步驟(4)得到的8個方向的魔方矩陣檢測算子分別進行褶積運算，得到8個方向的褶積運算數(shù)據(jù)體。式中：Dx為各個方向的檢測算子；f(x,y,z)為選取的以求取點為中心構(gòu)造的三維地震數(shù)據(jù)子體； (a,b,c)為算子中非零元素集合。

(6)將求取的各個方向褶積運算數(shù)據(jù)體中絕對值最大的值作為該方向邊緣檢測結(jié)果。

為了得到更為精確的識別方向，可以將3×3×3階魔方矩陣拓展到5×5×5階，根據(jù)上述算法實現(xiàn)流程得到5×5×5階的16個方向的魔方矩陣算子，按照類似流程得到斷層檢測結(jié)果。由數(shù)字1～125作為元素構(gòu)成的5×5×5階魔方矩陣為

(15)

類似地，也可得到其他階數(shù)的魔方矩陣檢測算子，但魔方矩陣的階數(shù)越大，計算效率越低，需要根據(jù)實際情況，綜合考慮選定三維魔方矩陣的階數(shù)。

3 理論模型討論

設(shè)計一個包含多個地壘的三維地質(zhì)模型，采用30Hz主頻的雷克子波，通過正演模擬得到含方差為0.2的高斯噪聲的地震模型數(shù)據(jù)(圖2)。文中從識別效果、抗噪性兩方面討論基于魔方矩陣斷層檢測方法的有效性。

3.1 識別效果

利用基于魔方矩陣的檢測方法對三維模型進行處理，得到各個方向的魔方矩陣屬性體數(shù)據(jù)。本文選取地壘模型Xline361剖面上斷層線上的O點(圖3)分析識別效果。圖4為O點的8個方向的魔方矩陣計算數(shù)值玫瑰圖，可見利用魔方矩陣方法可有效地識別斷層，k6方向的檢測結(jié)果更好。

3.2 抗噪性分析

選擇相干算法檢測結(jié)果與基于魔方矩陣的斷層檢測方法進行抗噪性分析。圖5為斷層方向加噪檢測結(jié)果?？梢娀谀Х骄仃嚨臄鄬訖z測方法的斷層邊界識別效果更好(圖5b)，即與相干算法(圖5a)相比，基于魔方矩陣的斷層檢測方法的抗噪性更好。

圖2 三維地質(zhì)模型正演模擬結(jié)果

圖3 Xline361橫測線剖面(添加高斯噪聲)

圖4 O點的8個方向的魔方矩陣計算數(shù)值玫瑰圖

圖5 斷層方向加噪檢測結(jié)果(a)相干算法； (b)k6方向

4 實際資料應(yīng)用

研究區(qū)塊(永3 工區(qū))位于山東東營凹陷東北部的永安鎮(zhèn)油田南部，是勝利探區(qū)最復(fù)雜的斷塊型油藏區(qū)塊之一，目標層斷裂眾多、結(jié)構(gòu)復(fù)雜，大大增加了復(fù)雜斷塊油藏的勘探、開發(fā)難度。因此，若能較好地檢測該區(qū)的斷裂結(jié)構(gòu)，對勘探、開發(fā)具有重要意義。圖6為沿T4地層抽取的原始沿層切片，選定數(shù)據(jù)點P(圖6紅色點)進行斷層方位信息檢測，圖7為P點空間位置。根據(jù)魔方矩陣算子具體數(shù)值可知，當算子方向旋轉(zhuǎn)到k4方向時(圖7算子上綠色點所近似的方向)，由于該方向與斷層傾向幾乎一致，屬性值最大，識別效果最好(圖8)。

圖9為實際資料檢測結(jié)果。由圖可見：①不同方法都能較好地檢測主要大斷層； ②與P點所在斷層傾向接近的k4方向(圖7)的算子檢測結(jié)果(圖9c)較其他方向的斷層識別結(jié)果更清晰；③k2方向檢測結(jié)果(圖9a中紅色箭頭指示斷層傾向)與k2檢測算子指示的方向一致，其檢測結(jié)果較其他方向的斷層識別結(jié)果更清晰、輪廓更精確； ④k3算子代表水平方向，其檢測結(jié)果(圖9b)類似于上、下延長時間提取屬性檢測斷層的常規(guī)方法，斷層識別效果較k2(圖9a)、k4(圖9c)方向差； ⑤與相干算法(圖9d)相比，基于魔方矩陣的斷層檢測方法(圖9a～圖9c)的抗噪效果更好，即不同方向的檢測算子對不同傾向的斷層檢測結(jié)果不一致，當檢測算子方向與斷層傾向一致時，檢測結(jié)果更好。

圖6 沿T4地層抽取的原始沿層切片

圖7 P點空間位置黃色箭頭為P點所在斷層傾向

圖8 P點各個方向檢測結(jié)果

圖9 實際資料檢測結(jié)果(a)k2； (b)k3； (c)k4； (d)相干算法

5 結(jié)束語

斷層的檢測精度對斷塊油田剩余油的勘探、開發(fā)具有重要意義。與傳統(tǒng)的邊緣檢測方法相比，基于魔方矩陣的斷層檢測方法利用魔方矩陣可靈活旋轉(zhuǎn)、且方向指向又很明確的特性，可有效檢測斷層，并定量給出斷層發(fā)育的傾向信息，較常規(guī)相干算法在抗噪性上有較大優(yōu)勢。本文以3×3×3階的魔方矩陣進行模型測試和應(yīng)用嘗試，也可得到其他階數(shù)的魔方矩陣檢測算子，但魔方矩陣的階數(shù)越大，計算效率越低，需要根據(jù)實際情況綜合考慮選定三維魔方矩陣的階數(shù)。

參考文獻

[1] 張軍華,王月英,趙勇.C3相干體在斷層和裂縫識別中的應(yīng)用.地震學報,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] 蔡涵鵬,胡光岷,賀振華等.基于非線性變時窗相干算法的不連續(xù)性檢測方法.石油地球物理勘探,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] 王雷,陳海清,陳國文等.應(yīng)用曲率屬性預(yù)測裂縫發(fā)育帶及其產(chǎn)狀.石油地球物理勘探,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] 鄧志文,趙賢正,陳雨紅等．自適應(yīng)波形多道匹配追蹤斷層識別技術(shù).石油地球物理勘探,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] 茍量,彭真明.小波多尺度邊緣檢測及其在裂縫預(yù)測中的應(yīng)用.石油地球物理勘探,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] 王清振,張金淼,姜秀娣等.基于高維小波變換的高抗噪性邊緣檢測技術(shù)．石油地球物理勘探,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] 張萌.基于標簽的信息嵌入與提取技術(shù)研究與實現(xiàn)[學位論文].北京：北京郵電大學，2014，20-22.

[13] 趙楊.基于H.264的大容量視頻信息隱藏算法研究[學位論文].四川成都:西南交通大學,2010,45-46.

[14] 李松斌,付江云,劉鵬等.一種基于整數(shù)DCT系數(shù)調(diào)制及N維魔方矩陣的H.264/AVC信息隱藏方法.小型微型計算機系統(tǒng),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] 劉爭艷,李絮,陳蘊.基于二維映射關(guān)系的視頻信息隱藏算法.計算機工程,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] 單潤紅,高峰,宋君強.魔方矩陣的特征值分析.高等數(shù)學研究,2004,7(4):45,52.