深地震反射剖面構(gòu)造信息識別研究

李文輝，高 銳＊，王海燕，李洪強

1 國土資源部深部探測與地球動力學(xué)重點實驗室（建），北京 100037

2 中國地質(zhì)科學(xué)院地質(zhì)研究所，北京 100037

1 引 言

自20世紀(jì)70年代美國COCORP計劃將反射地震勘探方法的原理和技術(shù)用于地殼結(jié)構(gòu)探測以來，各國科學(xué)家以此為基礎(chǔ)相繼實施了一系列深部探測計劃，并取得了豐碩成果［1］.目前，深地震反射方法已成為探測地殼上地幔精細結(jié)構(gòu)的最有效手段之一［2－7］.與石油地震勘探相比，深地震反射剖面探測深度達20～50s，且經(jīng)常要跨越造山帶、盆山結(jié)合帶等復(fù)雜地質(zhì)條件區(qū)域，因此具有藥量大、排列長、頻率低、頻帶窄、深部信號弱、速度橫向變化大等特點.受長距離傳播能量吸收衰減和深部不同傾角的復(fù)雜地質(zhì)體等因素的影響，深地震反射剖面中下地殼的地震波組經(jīng)常表現(xiàn)為能量弱、不連續(xù)、帶狀或交織狀，給地震資料的解釋帶來了困難［8－10］.

為了清楚的反映深部地質(zhì)構(gòu)造格架，前人將地震波組以線條的形式表示，稱為線劃圖.線劃剖面視覺上一目了然，易于解釋［11］，然而人工制作線劃圖費時費力，而且受主觀因素影響很大［9，12］.隨著信息技術(shù)的進步，人們試圖利用計算機繪制線劃圖，稱為自動線條圖技術(shù).實現(xiàn)自動線條圖的方法可歸納為基于數(shù)字圖像處理和基于模式識別兩大類.其中基于數(shù)字圖像處理生成線條圖的基本原理是將地震剖面的采樣點轉(zhuǎn)換為灰度圖像像元，反射同相軸就相當(dāng)于灰度圖像的邊緣，進而應(yīng)用各種邊緣檢測算子實現(xiàn)同相軸的識別［13－15］.基于模式識別的自動線條圖技術(shù)則通過建立描述波形的模式基元、構(gòu)建描述模式基元之間關(guān)系的目標(biāo)函數(shù)、迭代計算形成三元組、連接三元組等步驟實現(xiàn)剖面線條化［8，16－22］.相比而言，數(shù)字圖像處理方法僅考慮地震波的振幅信息，易于實現(xiàn)，但其對復(fù)雜情況的識別效果一般.模式識別方法直接對波形特征進行分析，而且可以通過進一步分析識別過程中產(chǎn)生的屬性信息來輔助解釋［8］.該方法在加拿大Lithoprobe計劃中得到了廣泛應(yīng)用［23－25］，但其缺點是經(jīng)驗參數(shù)多，實現(xiàn)過程復(fù)雜.

本文在研究和總結(jié)上述兩種方法的基礎(chǔ)上，借鑒模式識別方法的基本思想，并汲取圖像處理技術(shù)相關(guān)算法，提出一種新的深地震反射剖面構(gòu)造格架識別方法.該方法通過數(shù)據(jù)預(yù)處理、振幅提取、對象識別、連續(xù)性計算和連續(xù)性濾波實現(xiàn)深地震反射剖面線條化，同時還可通過對象傾角計算對復(fù)雜區(qū)域進行屬性分析.文章最后討論了線條圖在深地震反射剖面解釋中的優(yōu)點和不足，并展望了進一步研究的方向.

2 原理與設(shè)計

深地震反射資料能夠獲取地殼尺度的精細結(jié)構(gòu)，其包含了巖性界面、構(gòu)造運動（碰撞、剪切、走滑、推覆等）、巖漿巖分布及流體等豐富的地質(zhì)信息，已成為地質(zhì)學(xué)家研究地殼內(nèi)部結(jié)構(gòu)時首選的約束依據(jù)［12，26］.由于這些信息在反射剖面上大都體現(xiàn)為同相軸，故本文所指構(gòu)造信息的識別主要指同相軸的識別.

來自同一反射界面的同相軸由于波阻抗差使其振幅大于干擾波，而具有強振幅性.另外，同一界面反射波到達相鄰接收點的射線路經(jīng)相近，其相位、頻率、到達時間也是相近的，因此在剖面上表現(xiàn)為連續(xù)平滑的曲線，稱為同相性［27］.根據(jù)以上兩個判別標(biāo)志，本文采用的研究思路是：（a）根據(jù)振幅分布特征，按照一定閾值，提取強振幅信息；（b）對提取的信息進行對象識別，生成對象關(guān)系表；（c）計算對象連續(xù)性，并通過連續(xù)性濾波獲得剖面線條圖；（d）針對復(fù)雜區(qū)域進行定量傾角分析輔助解釋.識別流程圖如圖1所示.

圖1 深地震反射剖面構(gòu)造信息識別處理流程Fig.1 Flow chart of structure information recognition from deep seismic reflection profile

3 實現(xiàn)與效果

本文選取實際剖面部分資料作為實驗數(shù)據(jù).為了便于后續(xù)處理，首先將原始疊加剖面轉(zhuǎn)為二維數(shù)字矩陣，并按照Bondar［28］的方法將每個采樣點的振幅值按線性關(guān)系將振幅范圍和灰度級建立對應(yīng)關(guān)系，形成地震數(shù)據(jù)灰度圖像（圖2a）.

3.1 強振幅提取

圖2 構(gòu)造信息識別實驗（a）數(shù)據(jù)灰度圖像；（b）強振幅提取；（c）中值濾波；（d）連續(xù)性濾波Fig.2 Recognition experiment（test data was picked out from a part of seismic profile）（a）Seismic gray level image；（b）Amplitude analysis；（c）Median filtering；（d）Continuity filtering.

地震反射剖面中有效反射波組的能量較干擾波更強，且疊加（或偏移）剖面的采樣率遠高于有效波的尼奎斯特采樣限制，因此根據(jù)振幅幅值分布特征按照合適閾值提取強振幅，不會影響剖面的基本反射特征.本文選擇保留大于一定閾值的強振幅波峰信息，其中閾值按照公式（1）計算

式中T為閾值，max，min分別對應(yīng)剖面中最大、最小振幅值，α為振幅提取因子，其取值范圍在0～1之間.經(jīng)多次測試，α取值在0.55～0.65之間較為合理.圖2b為根據(jù)圖2a的振幅分布特征（圖3）對其按α＝0.6提取強振幅的結(jié)果，結(jié)果顯示提取后的圖像保留了原數(shù)據(jù)主要反射特征.另外為了便于對象識別，在振幅提取同時將灰度圖二值化，即對振幅小于臨界閾值的樣點賦0值，大于臨界閾值部分賦1值.

圖3 剖面振幅分布Fig.3 Amplitude distribution of the profile

3.2 中值濾波

強振幅提取及二值化會使剖面圖像產(chǎn)生大量椒鹽噪音（圖2b）.中值濾波是一種非線性空間濾波器，對于濾除椒鹽噪音具有較好效果，而且在去噪時不會對連續(xù)性好的有效信息造成傷害［29］.二維中值濾波可用公式（2）（3）（4）表達

式中f（x，y）代表原始圖像x，y處的灰度值，W 為鄰域算子，它是由l和k組成的矩陣角碼集合.其原理是提取原始數(shù)據(jù)f（x，y）周圍某一鄰域W 內(nèi)的所有n個元素c1－cn，記為集合A.對A包含的n個元素排序，并取其中值g（x，y）代替原始值f（x，y）.中值濾波鄰域算子大小和形狀的選擇非常關(guān)鍵，其與剖面的采樣率、子波頻率有關(guān).為保證濾波時不影響原始有效信息形態(tài)，圖2c為采用10×5矩形鄰域算子對圖2b進行中值濾波的結(jié)果.

3.3 對象識別

定義二值圖像中連通的區(qū)域為一個獨立“對象”，為了后續(xù)連續(xù)性計算、濾波及屬性分析，需對圖像中的對象進行識別（對象矢量化）.識別過程按照經(jīng)典8方向flood－fill（漫水填充）算法通過對圖像進行逐元素連通性搜索來實現(xiàn)，識別同時對對象進行“染色”（編號）.對于該算法的具體計算原理，在此不再贅述.另外，我們引入模式識別方法中關(guān)系表的概念［8］，即建立以編號為關(guān)鍵字的關(guān)系數(shù)據(jù)庫來存儲每個對象的相關(guān)信息.本文采用的對象關(guān)系表由對象編號、元素個數(shù)、元素的位置、對象長度、對象傾角等組成.

3.4 連續(xù)性計算及長度濾波

連續(xù)性是評價反射同相軸的重要指標(biāo)，連續(xù)性好的同相軸往往刻畫了剖面的主要構(gòu)造格架.本文通過計算對象長度來衡量同相軸的連續(xù)性.由于對象往往是狹長條狀不規(guī)則圖形，因此可用對象最小外接矩形的對角線長度近似對象的長度（圖4）.最小外接矩陣可通過對原坐標(biāo)系按一個小的角度增量旋轉(zhuǎn)遍歷獲得.公式（5）（6）表示了最小外接矩陣的計算方法

式中Lm，Wm代表矩陣的長和寬，βm為最小面積外接矩陣對應(yīng)的坐標(biāo)旋轉(zhuǎn)角度，x，y分別為對象所含元素位置坐標(biāo)（矩陣角碼）組成的向量，max和min分別代表求取最大、最小值的函數(shù).公式（7）中OL即為對象長度.

圖4 對象長度計算原理圖中灰色區(qū)域為對象，虛線方框為最小外接矩陣，黑實線為近似對象長度.Fig.4 Principle of object length calculationGrey area is the object，dash rectangle is the smallest rectangle containing the region，solid black line represents the object length.

值得注意的是，以上計算均以像元為單位，在時間剖面圖像中一個像元在橫向上代表一個CDP道，縱向上則代表一個時間采樣點，因此該長度不具備實際物理量綱意義.通過連續(xù)性計算，將所有對象的長度信息存儲在對象關(guān)系表中，便可按照需求對其通過長度條件判斷實現(xiàn)連續(xù)性濾波形成線條圖.圖2d為按長度40進行連續(xù)性濾波的結(jié)果.

圖5（a，b，c）分別為對華南廬樅地區(qū)一條長92km的深地震反射疊加剖面按以上流程進行識別，并以長度≥20、≥40、≥60進行連續(xù)性濾波的結(jié)果（濾波長度選擇愈小則剖面細節(jié)保留愈多）.與原始波形剖面（圖6）對比表明圖5（a，b，c）能有效識別原始資料中的主要構(gòu)造格架信息，尤其對深部信噪比低、能量較低弱部分效果顯著。圖7為對該剖面資料使用經(jīng)圖像邊緣檢測方法處理得到的結(jié)果，相比而言本文所采用方法的效果較圖像處理方法有明顯改善.

3.5 復(fù)雜區(qū)域?qū)傩苑治?/h3>

深地震反射剖面中來自深層的地震波組往往并非由單一界面形成，而是眾多不均勻條帶疊加的結(jié)果，而且受不同期次大地構(gòu)造運動的影響，構(gòu)造產(chǎn)狀從淺到深可能出現(xiàn)傾向相反的復(fù)雜情況［6］.為了幫助解釋人員判斷這些復(fù)雜情況，本文在對象識別、連續(xù)性濾波的基礎(chǔ)上，對復(fù)雜區(qū)域中的對象進行傾角計算，并通過統(tǒng)計其分布規(guī)律達到定量分析的目的.對象傾角計算采用標(biāo)準(zhǔn)差橢圓法進行.該方法是一種常用的散點定向方法，其基本原理是在散點集的原始坐標(biāo)系下，假設(shè)存在某一方向，所有點到該方向的標(biāo)準(zhǔn)差距離最小，那么該方向與原始坐標(biāo)X軸的夾角即為散點集的方向［30］（圖8）.對象傾角按照公式（8）計算

式中θ代表對象傾角，n為對象所包含的元素個數(shù).Xi，Yi可由公式（9）（10）求得

其中xi，yi為對象中元素的坐標(biāo)位置.

為了測試其效果，我們從實際深反射資料截取一段剖面，識別結(jié)果顯示該區(qū)反射密集且關(guān)系復(fù)雜（圖9左圖）.通過對該區(qū)進行傾角分析，右圖顯示該區(qū)同相軸主要分布在－30°～30°之間，但總體以下傾方向為主.

4 結(jié)論與認識

由于深地震反射剖面反映的深部構(gòu)造信息無法由鉆井?dāng)?shù)據(jù)約束，目前深反射剖面的解釋仍以構(gòu)造解譯為主.線條圖技術(shù)作為一項專門針對深地震反射資料解釋的技術(shù)，具有特殊意義.我們通過實驗及實際剖面數(shù)據(jù)測試，得到以下結(jié)論與認識：

（1）本文提出的方法能夠識別深地震反射資料中的主要構(gòu)造格架信息，其效果較圖像處理方法有顯著改善.另外，由于去除了波形特征描述和迭代等復(fù)雜步驟，與模式識別法相比更高效.

圖5 廬樅地區(qū)深地震反射剖面識別結(jié)果（a）長度≥20濾波結(jié)果；（b）長度≥40濾波結(jié)果；（c）長度≥60濾波結(jié)果.Fig.5 Recognition result of the Luzong deep seismic reflection profile（a）Filtered with length≥20；（b）Filtered with length≥40；（c）Filtered with length≥60.

圖6 廬樅地區(qū)深地震反射原始波形剖面Fig.6 The Original Luzong deep seismic reflection profile

圖7 廬樅地區(qū)深地震反射剖面數(shù)字圖像處理邊緣檢測結(jié)果Fig.7 Edge detecting result of the Luzong deep seismic reflection profile

圖8 對象傾角計算原理示意圖灰色區(qū)域為對象，虛線為對象標(biāo)準(zhǔn)差橢圓，黑實線與X軸夾角為對象傾角.Fig.8 Principle of object angle calculatingGrey area is the object；dash ellipse is standard deviational ellipse；the angle between the long axis of ellipse and the Xaxis represents the object dip angle.

圖9 復(fù)雜區(qū)域傾角分析結(jié)果Fig.9 Dip angle analysis of a complex area

（2）本文在識別過程中對整條剖面使用相同的參數(shù)，而未考慮地震資料在橫向上和縱向上的不均一性，后續(xù)可改進為開窗口并使用不同參數(shù)處理.

（3）深部結(jié)構(gòu)特征是深地震反射剖面能夠反映的最主要的構(gòu)造信息，然而除此之外，地震波組的能量對比、頻率變化同樣包含了豐富的構(gòu)造信息.例如通過亮點分析獲取深部流體及熔融體特征［31］、利用弱反射“透明體”追蹤花崗巖基的分布等等.因此，如何有效利用除結(jié)構(gòu)特征外的其他構(gòu)造信息并進行綜合解釋，是值得進一步研究的方向.

（References）

［1］王海燕，高銳，盧占武等.深地震反射剖面揭露大陸巖石圈精細結(jié)構(gòu).地質(zhì)學(xué)報，2010，84（6）：818－838.Wang H Y，Gao R，Lu Z W，et al.Fine structure of the continental lithosphere circle revealed by deep seismic reflection profile.Acta Geolocica Sinica （in Chinese），2010，84（6）：818－838.

［2］Gao R，Li P W，Li Q S，et al.Deep process of the collision and deformation on the northern margin of the Tibetan plateau：revelation from investigation of the deep seismic profiles.Science in China （Series D），2001，44：71－78.

［3］Gao R，Lu Z W，Li Q S，et al.Geophysical probe and geodynamic study of the crust and upper Mantle in the Qinghai－Tibet plateau，China.Episodes，2005，28（4）：263－273.

［4］王海燕，高銳，薛愛民等.適用于造山帶深地震反射資料的動校正方法.吉林大學(xué)學(xué)報（地球科學(xué)版），2006，36（4）：622－626.Wang H Y，Gao R，Xue A M，et al.A NMO method adapted to deep seismic reflection in orogenic belt.Journal of Jilin University－Earth Science Edition （in Chinese），2006，36（4）：622－626.

［5］盧占武，高銳，李秋生等.橫過青藏高原羌塘地體中央隆起區(qū)的深反射地震試驗剖面.地球物理學(xué)報，2009，52（8）：2008－2014.Lu Z W，Gao R，Li Q S，et al.Testing deep seismic reflection profiles across the central uplift of the Qiangtang terrane in the Tibetan Plateau.Chinese Journal of Geophysics （in Chinese），2009，52（8）：2008－2014.

［6］酆少英，高銳，龍長興等.銀川地塹地殼擠壓應(yīng)力場：深地震反射剖面.地球物理學(xué)報，2011，54（3）：692－697.Feng S Y，Gao R，Long C X，et al.The compressive stress field of Yinchuan garben：Deep seismic reflection profile.Chinese Journal of Geophysics （in Chinese），2011，54（3）：692－697.

［7］于常青，趙殿棟，楊文采.塔里木盆地結(jié)晶基底的反射地震調(diào)查.地球物理學(xué)報，2012，55（9）：2925－2938.Yu C Q，Zhao D D，Yang W C.Seismic reflection investigations of basement in Tarim basin.Chinese Journal of Geophysics （in Chinese），2012，55（9）：2925－2938.

［8］Li Q，Vasudevan K，Cook F A.Seismic skeletonization：a new approach to interpretation of seismic reflection data.Journal of Geophysical Research （Solid Earth），1997，102（B4）：8427－8445.

［9］徐明才，高景華.深部地震資料的處理和解釋方法.物探化探計算技術(shù)，1999，21（2）：151－158.Xu M C，GAO J H.Methods of processing and interpretation about the deep seismic data.Computing Techniques for Geophysical and Geochemical Exploration （in Chinese），1999，21（2）：151－158.

［10］王海燕，高銳，馬永生等.若爾蓋盆地—西秦嶺造山帶結(jié)合部位深反射資料的靜校正和去噪技術(shù).地球物理學(xué)進展，2007，22（3）：743－749.Wang H Y，Gao R，Ma Y S，et al.Static correction and noise removal for deep seismic reflection data from the contact zone of the Zoigêand the west Qinling orogen.Progress in Geophysics，2007，22（3）：743－749.

［11］陳滬生.下?lián)P子地區(qū)HQ－13線的綜合地球物理調(diào)查及其地質(zhì)意義.石油與天然氣地質(zhì)，1988，9（3）：211－222.Chen H S.Comprehensive geophysical survey of HQ－13line in the lower Yangzi reaches and its geological significance.Oil and Gas Geology （in Chinese），1988，9（3）：211－222.

［12］陳志德，楊文采，李玲等.松遼盆地北部深反射地震χ2分布處理及其深部地質(zhì)特征.石油地球物理勘探，2003，38（6）：654－660.Chen Z D，Yang W C，Li L，et al.χ2distribution processing of deep seismic reflection in north of Songliao basin and its deep geological feature. Oil Geophysical Prospecting （in Chinese），2003，38（6）：654－660.

［13］高美娟，田景文，高興友等.利用邊緣檢測法檢測地震反射同相軸.大慶石油學(xué)院學(xué)報，2000，24（3）：8－11.Gao M J，Tian J W，Gao X Y，et al.Detection of seismic reflection event by edge detection.Journal of Daqing Petroleum Institute （in Chinese），2000，24（3）：8－11.

［14］李紅星，劉財，陶春輝.圖像邊緣檢測方法在地震剖面同相軸自動檢測中的應(yīng)用研究.地球物理學(xué)進展，2007，22（5）：1607－1610.Li H X，Liu C，Tao C H.The study of application of edge measuring technique to the detection of phase axis of the seismic section.Progress in Geophysics （in Chinese），2007，22（5）：1607－1610.

［15］陳學(xué)華，賀振華，黃德濟.地震資料的高階偽希爾伯特變換邊緣檢測.地球物理學(xué)進展，2008，23（4）：1106－1110.Chen X H，He Z H，Huang D J.Seismic data edge detection based on higher－order pseudo Hilbert transform.Progress in Geophysics（in Chinese），2008，23（4）：1106－1110.

［16］Tφdstheim D.Improved seismic discrimination using pattern recognition.Physics of the Earth and Planetary Interiors，1978，16：85－108.

［17］Lu S Y.A string－to－string correlation algorithm for image skeletonization.Proc.International Conference on Pattern Recognition，1982，6：178－190.

［18］Cheng Y C，Lu S Y.The binary consistency checking scheme and its applications to seismic horizon detection.IEEE Trans.Pattern Analysis，1989，11：439－447.

［19］Le L，Nyalnd E.Pattern analysis of seismic records.Geophysics，1990，55：20－28.

［20］郭良輝，孟小紅，薛愛民.地震剖面線條化的一種簡單算法.地球科學(xué)—中國地質(zhì)大學(xué)學(xué)報，2007，32（4）：545－548.Guo L H，Meng X H，Xue A M.A simple algorithm for seismic skeletonization.Earth Science Journal of China University of Geosciences （in Chinese），2007，32（4）：545－548.

［21］李文輝，高銳，王海燕.Skeletonization技術(shù)及其在深地震反射剖面解釋中的應(yīng)用.地球物理學(xué)進展，2010，25（4）：1161－1167.Li W H，GAO R，WANG H Y.Seismic Skeletonization and Its application in deep seismic reflection Profiles′interpretation.Progress in geophysics （in Chinese），2010，25（4），1161－1167.

［22］Lu S Y，Cheng Y C.An iterative approach to seismic skeletonization.Geophysics，1990，55：1312－1320.

［23］Eaton D，Vasudevan K.Skeletonization of aero－magnetic data.Geophysics，2004，69（2）：478－488.

［24］Vasudevan K，Eaton D，Cook F A.Adaptation of seismic skeletonization for other geosciences applications.Geophysical Journal International，2005，161：975－993.

［25］Vasudevan K.，Eaton D，Cook F A.Seismic skeletonization：a useful tool for geophysical data analysis.Cseg Recorder，2006，31（9）：36－42.

［26］Mooney W D，Brocher T M.Coincident seismic reflection／refraction studies of the continental lithosphere：global view.Review of Geophysics，1987，25（4）：723－742.

［27］姚姚.地震波場與地震勘探.北京：地質(zhì)出版社，2007：199－200.Yao Y.Seismic Wave Field and Seismic Exploration （in Chinese）.Beijing：Geological Publishing House，2007：199－200.

［28］Bondar I.Seismic horizon detection using image processing algorithms.Geophysical Prospecting，1992，40：785－800.

［29］Gonzalez R C，Woods R E.Digital Image Processing（Second Edition）.Beijing：Publishing House of Electronics Industry，2002：14－19.

［30］王寶軍.基于標(biāo)準(zhǔn)差橢圓法SEM圖像顆粒定向研究原理與方法.巖土工程學(xué)報，2009，31（7）：1082－1087.Wang B J.Theories and methods for soil grain orientation distribution in SEM by standard deviational ellipse.Chinese Journal of Geotechnical Engineering （in Chinese），2009，31（7）：1082－1087.

［31］Kelly C，Lu Z W，Klemperer S L，et al.SinoProbe seismic reflection imaging of an upper－crustal“bright spot”beneath the Karakoram Fault，West Tibet.2012，AGU 2012Fall Meeting Abstract.