基于魔方矩陣的斷層檢測(cè)方法

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

1 引言

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

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

2 方法原理

2.1 魔方矩陣算法簡(jiǎn)介

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

(1)

(2)

(3)

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

(4)

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

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

(5)

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

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

(6)

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

(7)

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

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

(8)

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

…

(9)

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

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

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

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

(10)

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

(11)

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

(12)

分析式(12)可知，當(dāng)元素絕對(duì)值最大對(duì)應(yīng)方向?yàn)閿鄬友由旆较驎r(shí)，元素絕對(duì)值最小所在方向正好與其垂直，代表斷層形態(tài)，這與實(shí)際斷層相符合，具有實(shí)際物理意義。

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

…

(13)

圖1 不同方向檢測(cè)算子示意圖

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

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

(14)

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

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

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

(15)

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

3 理論模型討論

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

3.1 識(shí)別效果

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

3.2 抗噪性分析

選擇相干算法檢測(cè)結(jié)果與基于魔方矩陣的斷層檢測(cè)方法進(jìn)行抗噪性分析。圖5為斷層方向加噪檢測(cè)結(jié)果?？梢?jiàn)基于魔方矩陣的斷層檢測(cè)方法的斷層邊界識(shí)別效果更好(圖5b)，即與相干算法(圖5a)相比，基于魔方矩陣的斷層檢測(cè)方法的抗噪性更好。

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

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

圖4 O點(diǎn)的8個(gè)方向的魔方矩陣計(jì)算數(shù)值玫瑰圖

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

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

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

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

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

圖7 P點(diǎn)空間位置黃色箭頭為P點(diǎn)所在斷層傾向

圖8 P點(diǎn)各個(gè)方向檢測(cè)結(jié)果

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

5 結(jié)束語(yǔ)

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

參考文獻(xiàn)

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

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

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