深地震反射剖面構造信息識別研究

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

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

2 中國地質科學院地質研究所，北京 100037

1 引 言

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

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

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

2 原理與設計

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

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

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

3 實現與效果

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

3.1 強振幅提取

圖2 構造信息識別實驗（a）數據灰度圖像；（b）強振幅提取；（c）中值濾波；（d）連續性濾波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.

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

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

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

3.2 中值濾波

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

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

3.3 對象識別

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

3.4 連續性計算及長度濾波

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

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

圖4 對象長度計算原理圖中灰色區域為對象，虛線方框為最小外接矩陣，黑實線為近似對象長度.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道，縱向上則代表一個時間采樣點，因此該長度不具備實際物理量綱意義.通過連續性計算，將所有對象的長度信息存儲在對象關系表中，便可按照需求對其通過長度條件判斷實現連續性濾波形成線條圖.圖2d為按長度40進行連續性濾波的結果.

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

3.5 復雜區域屬性分析

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

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

其中xi，yi為對象中元素的坐標位置.

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

4 結論與認識

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

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

圖5 廬樅地區深地震反射剖面識別結果（a）長度≥20濾波結果；（b）長度≥40濾波結果；（c）長度≥60濾波結果.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 廬樅地區深地震反射原始波形剖面Fig.6 The Original Luzong deep seismic reflection profile

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

圖8 對象傾角計算原理示意圖灰色區域為對象，虛線為對象標準差橢圓，黑實線與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 復雜區域傾角分析結果Fig.9 Dip angle analysis of a complex area

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

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

（References）

［1］王海燕，高銳，盧占武等.深地震反射剖面揭露大陸巖石圈精細結構.地質學報，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］王海燕，高銳，薛愛民等.適用于造山帶深地震反射資料的動校正方法.吉林大學學報（地球科學版），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］盧占武，高銳，李秋生等.橫過青藏高原羌塘地體中央隆起區的深反射地震試驗剖面.地球物理學報，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］酆少英，高銳，龍長興等.銀川地塹地殼擠壓應力場：深地震反射剖面.地球物理學報，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］于常青，趙殿棟，楊文采.塔里木盆地結晶基底的反射地震調查.地球物理學報，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］徐明才，高景華.深部地震資料的處理和解釋方法.物探化探計算技術，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］王海燕，高銳，馬永生等.若爾蓋盆地—西秦嶺造山帶結合部位深反射資料的靜校正和去噪技術.地球物理學進展，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］陳滬生.下揚子地區HQ－13線的綜合地球物理調查及其地質意義.石油與天然氣地質，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分布處理及其深部地質特征.石油地球物理勘探，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］高美娟，田景文，高興友等.利用邊緣檢測法檢測地震反射同相軸.大慶石油學院學報，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］李紅星，劉財，陶春輝.圖像邊緣檢測方法在地震剖面同相軸自動檢測中的應用研究.地球物理學進展，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］陳學華，賀振華，黃德濟.地震資料的高階偽希爾伯特變換邊緣檢測.地球物理學進展，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］郭良輝，孟小紅，薛愛民.地震剖面線條化的一種簡單算法.地球科學—中國地質大學學報，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技術及其在深地震反射剖面解釋中的應用.地球物理學進展，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］姚姚.地震波場與地震勘探.北京：地質出版社，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］王寶軍.基于標準差橢圓法SEM圖像顆粒定向研究原理與方法.巖土工程學報，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.