深地震反射剖面构造信息识别研究

李文辉，高 锐＊，王海燕，李洪强

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.