基于魔方矩陣的斷層檢測方法

李 軍 張軍華 龔明平 楊 勇 杜玉山 武 剛(中國石油大學(華東)地球科學與技術學院，山東青島 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.