目前，浅地层剖面探测技术已越来越多地应用于海洋、湖泊、河流多种地质调查中(孙顺才等，1988；张金城等，1995；金翔龙，2007).浅地层剖面仪利用声波在水下沉积物内的传播和反射特性，探测地层内沉积物的性质及地层间的声阻抗界面，得到的浅地层剖面图谱能够反映水下地形、浅部地层、沉积物构造和基岩构造等水下沉积特征，为近海、湖泊沉积演变研究、航道整治、码头选址论证等提供必要的基础数据(Park et al.，1996；Palmer et al.，1997；Lafferty et al.，2006；曹敬贺等，2014).

振幅是地震波非常重要的属性，地层介质对地震波的透射与反射特性能体现在振幅的变化上，因此地震波的振幅变化包含了地层的分层信息(袁子龙等，2006；滕吉文，2014；司芗等，2014).利用地震波振幅进行地层分析在煤层厚度预测、油气勘探等方面已有较为广泛的研究和应用(董守华等，2004；李文辉等，2012；石玉梅等，2014).在利用地震波振幅进行浅地层地震剖面自动解译方面，刘秀娟等做了尝试性的初步研究，根据地震波振幅在海底表面的反射特征，采用最大振幅法成功提取了海底地形(刘秀娟等，2009)，以及对浅层地震剖面中的风浪形变效应进行了校正(刘秀娟等，2007).但是从原始数据出发，结合地震波振幅的浅地层剖面的自动解译目前还没有成熟的方法.传统的剖面解译主要采用人工目视解译法，这种方法虽然灵活，但速度慢，需要解译者具有比较丰富的经验，解译结果受解译者的主观因素影响较大.因此，从浅地层剖面原始数据出发，通过数字图像处理技术，实现对浅地层剖面高效客观地自动解译，有必要进行探索与研究.

本文以青藏高原黑海的浅地层地震剖面数据为研究对象，根据其反射地震波的振幅特征，结合均值曲线法，利用优化的能量比法判读剖面中的地层分界线，达到计算机自动解译浅地层剖面的目的. 1 资料与数据 1.1 研究区概况

本文选用的数据是2011年在青藏高原东北部的湖泊 ——黑海所采集的浅地层地震剖面数据.黑海(35.9859°N，93.2638°E)位于青藏高原北部东昆仑山的北麓昆仑河-野牛沟断裂谷地内.平均海拔4477 m，沿谷地NE-SW向延伸，面积38.7 km²，湖中心最大水深21.9 m.地质构造上，北部为两组东西走向的大断裂，即东昆仑中断裂和南沟下断裂；南部为东西走向的昆仑河-野牛沟大断裂，再向南为西大滩断裂和东昆仑南断裂，构造控制明显(周连福，2014).湖盆南缘，为近东西走向的东昆仑山，系加里东晚期-燕山早期的花岗岩山体，平均海拔4500~5000 m，雪线高度约5000 m；西南部为三叠纪石灰岩，山前冲积扇发育，自山前一直伸入湖中，为磨圆度较差的砾石与细颗粒物质混杂堆积.北部和东部的山体相对较低，顶部无冰川发育，主要为三叠纪中下统洪水川组、闹仓坚沟组和希里可特组砂岩(曹凯等，2007). 1.2 数据采集

与在海上进行剖面采集不同，湖中风浪相对较弱，无其它船只活动，使用8人橡皮小艇，小艇马达产生的噪声以及周围环境中的噪声影响较弱，接收单元能够接收到更多的有效地震波反射信号，进行剖面解译时，易于提取到更多的真实边界点.同时，黑海处在高原地区，人类活动少，总体悬浮物浓度低，利于数据采集.

野外数据采集使用英国GeoAcoustic公司制造的GeoPulse 3.5 kHz浅地层剖面系统，理想状态下，有效探测深度可达70 m.系统由4个 3.5 kHz的换能器与拖体、发射机、接收机和GeoPro II工作站组成，数据采集时设定的主要参数值为：脉冲发射周期与记录长度为100 ms，信号采集率为60 μs，高低通滤波300 Hz~7 kHz，增益为6 dB或9 dB，数据存储格式为SEG-Y格式.探测时导航定位采用美国天宝公司的GPS SPS351. 1.3 数据处理

SEG-Y文件的数据存储格式比较复杂，一般通用图像分析处理软件如ArcGIS、ERDAS Imagine、MapInfo等不能直接将其打开读取其中的数据，而GeoPro工作站专用的NOMOS软件也只能对其进行回放与图像调节处理，以及进行简单的标注，无法对其原始数据进行进一步的操作.

SEG-Y格式是石油勘探行业地震数据最为普遍的格式.其标准文件结构中，前3600字节(包括3200字节的文本文件和400字节的二进制文件头)一般看作是一个整体，随后是若干个3200字节的扩展文本文件头记录，最后是若干条道的数据记录(苏欣和杜鑫，2011).经检测与分析，GeoPro工作站采集与记录的SEG-Y文件格式基本与标准文件相同，但亦存在不符合标准的部分.研究工作中，对SEG-Y文件的处理方式是：(1)在Visual C++平台下先将文件中的二进制数据转换为十六进制的数据，通过反复的对比验证，找出SEG-Y文件的数据存储方式.(2)再将二进制数据转换为十进制的txt文本文件.

经实验，最终设计的txt文本文件格式为：每一道数据按次序存于一行中，方便程序在读取此数据文件时能依次序进行；每一道数据前存入道编号、样点数两个固定的信息，一方面可以统计每一道在整个剖面中的位置，另一方面可以在程序读取时能够预先得知每行样点数的数量，无需再通过其它途径进行判别. 2 基于能量比法的地层分界线提取方法

在进行浅地层剖面自动解译时，首先用经过变长滑动时窗优化的能量比法计算每道地震记录的能量比值；采用均值曲线法将得到的能量比值曲线转换为一条能量比均值曲线，选取合适长度的时窗对曲线进行峰值提取，即边界点提取；得到每条地震记录上的边界点后，根据实际地层界线的连续性特征，结合相邻地震记录进行边界点的连续性检验，剔除不满足条件的孤立点，进而获得目标界线.在进行能量比值计算和能量比均值曲线峰值提取时，窗口大小的设定是关键，通过设计不同的参数化方案实验，能够确定适合剖面解译的最优窗口大小，提高自动解译的精度. 2.1 能量比法的基本原理

完整的一条地震剖面是由若干道地震记录组成，每一道地震记录是浅剖仪发射一次声波后，接收单元接收的反射信号的相应记录，反映的是剖面中的一个地理位置点下的沉积地层信息.单道数据存储为离散的采样点的地震波振幅值，样点数由仪器的发射频率与采样频率共同决定，因此采样点与传播时间、探测深度有一定的对应关系(图 1).时窗定义在单道记录的采样点上，其大小以采样点的数量为单位.能量比值的计算公式为(刘喜武等，2008)

图 1

Fig. 1

图 1 采样点、传播时间(双程)、探测深度之间的对应关系 Fig. 1 The correspondence between sample points，travel time(two-way) and investigation depth

能量比法的时窗设计方法目前主要有三种：固定时窗法、滑动时窗法和变长滑动时窗法(张伟等，2009).固定时窗法比较简单、便于计算，但时窗长度固定，无法得到边界点的精确位置；滑动时窗法增加了计算样本，提高了边界点的提取精度，但固定的时窗长度仍是制约其应用的最大因素.变长滑动时窗法(图 2)是以上两种方法的综合与改进，能够以不同的时窗长度，滑动计算整条地震记录的能量比值，克服了前两种方法中边界点识别准确率低与位置不精确的缺点，是目前最周全的方法，这也是本文采用的时窗设计方法.

图 2

Fig. 2

图 2 变长滑动时窗示意图 Fig. 2 The diagram for transform sliding time-window

根据地层对地震波的反射特征可知，当边界点位于时窗中点时，时窗左右两部分的能量比值最大，在能量比值曲线中表现为峰值.选择不同时窗长度，运用一定的统计方法，将众多能量比值曲线同时出现峰值的位置提取出来，该位置即边界点的位置.

本文采用均值曲线法提取能量比值曲线的峰值，为了更能突显每个点的能量比值特征，在计算能量比值均值之前，先计算方差，剔除与均值相差较大的若干个值，直到剩余值的方差在设定范围以内或者剔除值的个数达到设定的最大个数为止，最后再计算剩余数据的平均值，得到均值曲线.

在均值曲线法中，同样采用变长滑动时窗法来判别峰值点.因此，提取时窗窗口大小(L_win)的设定依然是需要着重考虑的问题.若L_win太大(图 3 W1窗口)，某些峰值较小的点可能被其附近峰值更大的点所掩盖(H1、H2)，降低了峰值检测的灵敏度；反之若L_win太小(W2窗口)，可能会把由噪音导致的曲线中极微小的突起检测为峰值，峰值检测过度灵敏，降低了准确率.

图 3

Fig. 3

图 3 单道地震记录的能量比值曲线、均值曲线、原始波形 Fig. 3 The curves of energy ratio，mean energy ratio and original waveform of single channel seismic record

在实际地层中，地层界线在一定范围内的分布是连续的，因此当某一道记录中有地层边界点出现时，与其相邻的若干道记录中应该有与之连续的边界点.在地震剖面数据的获取过程中，可能会受到外界因素的干扰，造成一个剖面中部分道数据与正常测量结果有差别.因此，针对每一道记录中识别出的地层边界点，需要根据地层连续性原则，将相对离散的孤立点剔除，在理想状态下，最终保留下来的应该是真实、连续的边界点，在地震剖面上表现为连续分布的地层界线.

边界点连续性检验的基本原理为：以自动提取得到的剖面为操作对象，遍历其中所有的边界点，判断一个点其在指定深度范围内、指定若干道相邻记录范围内是否存在一定数量的其它点，若满足条件则认为此点属于边界点，将其保留；若不满足条件则认为此点是孤立点，将其剔除. 3 参数化方案实验

剖面的地层信息主要集中在一次、二次反射内，三次反射以下的地层信息因受多次反射信号影响已较为混乱，很难分辨出真实的地层结构(图 4)，因此在统计人工解译边界点时，只统计二次及二次反射之上的边界点.人工解译边界点的数目是根据每道地震记录中穿过的地层界线的条数确定的，统计时记录边界点的垂直样点序号，最后得到一条剖面中的所有的边界点及相应的样点序号.

图 4

Fig. 4

图 4 黑海3sn浅地层剖面人工目视解译结果 Fig. 4 The result of MVI method for the 3sn seismic profile in the Heihai Lake

自动解译时，解译精度由多个参数共同决定(表 1).参数大小以样点数计算，例如计算能量比值时时窗的最小值L_min=2，是指变长滑动时窗的初始长度为2个采样点长度.

表 1(Table 1)

表 1 剖面解译时的主要参数 Table 1 Themainparameters for the interpretation of the seismic profile 参数 含义 Lmin 计算能量比值时时窗的最小值 Lmax 计算能量比值时时窗的最大值 Cwin 时窗长度变化步长 Swin 时窗滑动步长 Lwin 均值曲线法峰值提取窗口长度

表 1 剖面解译时的主要参数 Table 1 Themainparameters for the interpretation of the seismic profile

L_min-L_max和L_win是影响解译精度的主要参数，对解译精度的影响较其他参数更大，参数化方案实验以这两组参数为基础进行设计.根据窗口大小对提取效果的影响，按一定间隔给L_min-L_max、L_win取值，对提取效果较好的取值方案做进一步分析，基于这样的人机交互过程，通过与人工解译的剖面图像对比，选定提取效果最好的参数化方案进行剖面的自动解译.

提取效果以边界点密集度和地层界线相似度来衡量.边界点密集度是地层界线上下一定地层厚度内的边界点数目占提取的边界点总数的比例.在实际地层中，地层界线是有一定厚度的过渡地层，为了便于分辨地层结构，人工解译时将其划为一条线，自动解译得到的边界点会聚集在人工解译界线附近的一定范围内，边界点密集度越大，地层界线越清晰，边界点密集度越小，地层界线越模糊.地层界线相似度为人工解译得到的边界点数量与自动解译得到的边界点数量的比值，相似度越大，说明两者趋势符合越好，自动解译地层界线少有间断；相似度越小，两者符合越差，自动解译地界层线多有间断.

表 2(Table 2)

表 2 参数化方案(Lmin-Lmax/Lwin)及提取边界点数量 Table 2 The parameterization scheme(Lmin-Lmax/Lwin) and the number of the extracted boundary points Lwin Lmin-Lmax 2~20 4~40 6~60 8~80 10~100 2 16509 17080 16587 16360 15740 8 17094 17690 17154 17215 16916 20 17114 17837 17088 17330 17188 50 15840 16916 16053 16105 16060 70 14298 15045 14516 14599 14463 100 11415 11782 11051 11215 11206

表 2 参数化方案(Lmin-Lmax/Lwin)及提取边界点数量 Table 2 The parameterization scheme(Lmin-Lmax/Lwin) and the number of the extracted boundary points

从图 5中可以看出，边界点数量随L_min-L_max的变大而减小，当L_win取值2、8、20、50时，提取的边界点数量较多；当L_win取值100时，提取的边界点数量最少；当L_win取值70时，提取的边界点数量处于中间水平.从图 6中可以看出，边界点数量的变化趋势随L_win的变大而减小，当L_min-L_max取值2~20、6~60、8~80时，提取的边界点数量较为接近，且变化趋势一致.

图 5

Fig. 5

图 5 不同Lwin的提取效果 Fig. 5 The extracting effects with the different Lwin

图 6

Fig. 6

图 6 不同Lmin-Lmax的提取效果 Fig. 6 The extracting effects with the different Lmin-Lmax

野外工作时，环境中的噪声信号和地震反射信号会被同时接收，存储在原始数据中，部分噪声信号在后期处理中很难被剔除.在自动解译时，一些噪声信号会被错误地识别为“边界点”而被提取出来.分析发现，由噪声产生的“边界点”多具有离散性(不连续性)特点，而真实的地层界线边界点在局部范围内具有连续性，因此可以通过设计不同的参数化方案来降低对错误边界点的识别率. 3.2 L_min-L_max取值方案

表 3、4、5中，t：人工解译的边界点数量；a1：自动解译的边界点数量；a2：人工解译地层界线上下25个采样点的厚度地层内自动解译的边界点数量；a1/t：界线相似度；a2/a1：边界点密集度.

表 3(Table 3)

表 3 Lmin-Lmax取2~20时不同Lwin的提取效果 Table 3 The extracting effects with the different Lwin(Lmin-Lmax：2~20) 2~20 2 8 20 50 70 100 t 8210 8210 8210 8210 8210 8210 a1 7879 7919 7897 7399 6904 5694 a2 5923 5941 5935 5793 5584 4744 a1/t 0.7214 0.7236 0.7229 0.7056 0.6801 0.5778 a2/a1 0.7517 0.7502 0.7516 0.7829 0.8088 0.8332

表 3 Lmin-Lmax取2~20时不同Lwin的提取效果 Table 3 The extracting effects with the different Lwin(Lmin-Lmax：2~20)

表 4(Table 4)

表 4 Lmin-Lmax取6~60时不同Lwin的提取效果 Table 4 The extracting effects with different Lwin(Lmin-Lmax：6~60) 6~60 2 8 20 50 70 100 t 8210 8210 8210 8210 8210 8210 a1 8294 8308 8265 7733 7170 5615 a2 6127 6050 6036 5911 5621 4578 a1/t 0.7465 0.7369 0.7352 0.72 0.6847 0.5576 a2/a1 0.739 0.7282 0.7303 0.7644 0.784 0.8153

表 4 Lmin-Lmax取6~60时不同Lwin的提取效果 Table 4 The extracting effects with different Lwin(Lmin-Lmax：6~60)

表 5(Table 5)

表 5 Lmin-Lmax=8~80时不同Lwin的提取效果 Table 5 The extracting effects with different Lwin(Lmin-Lmax：8~80) 8~80 2 8 20 50 70 100 t 8210 8210 8210 8210 8210 8210 a1 8459 8431 8392 7721 7125 5708 a2 5995 5809 5781 5558 5264 4429 a1/t 0.7303 0.7079 0.7041 0.677 0.6412 0.5395 a2/a1 0.7088 0.689 0.6889 0.7199 0.7388 0.7759

表 5 Lmin-Lmax=8~80时不同Lwin的提取效果 Table 5 The extracting effects with different Lwin(Lmin-Lmax：8~80)

当L_min-L_max取2~20、6~60、8~80时，提取的边界点数量及a1/t、a2/a1的变化趋势相同，但当L_min-L_max取2~20时的a2/a1更大(图 7)，即当L_min-L_max取2~20时的解译剖面中边界点的密集度更大，界线更加清晰，而且在L_min-L_max取6~60、8~80的解译剖面中，局部的地层界线出现不连续现象(图 8)，L_min-L_max取2~20时的解译剖面中，地层界线的不连续现象明显改善，解译效果更好.

图 7

Fig. 7

图 7 不同Lmin-Lmax取值，a2/a1的分布 Fig. 7 The distribution of a2/a1 under different conditions of Lmin-Lmax

图 8

Fig. 8

图 8 Lwin=70时，不同Lmin-Lmax的解译效果 Fig. 8 The extracting effects with different values of Lmin-Lmax(Lwin=70)

当L_min-L_max取2~20时，随L_win的增大，a1与a2逐渐变小(图 9)，这与图 5的情况相吻合；同时，a1与a2之间的差距逐渐减小，a2/a1逐渐增大(图 10).即随着L_win的增大，边界点的密集度增大，解译剖面中的地层界线更加清晰；但a1/t的值在逐渐减小，即随着L_win的增大，自动解译地层界线与人工解译地层界线的相似度变差，自动解译的地层界线出现不连续现象(图 11).当L_win=70时，a1/t=0.6801，a2/a1=0.8088，解译剖面中可以分辨出清晰的地层界线，且其连续性较好，在L_win的所有取值中，解译效果最好.

图 9

Fig. 9

图 9 Lmin-Lmax取2~20时，a1、a2的分布 Fig. 9 The distribution of a1 and a2(Lmin-Lmax：2~20)

图 10

Fig. 10

图 10 Lmin-Lmax取2~20时，a1/t、a2/a1的分布 Fig. 10 The distribution of a1/t and a2/a1(Lmin-Lmax：2~20)

图 11

Fig. 11

图 11 Lmin-Lmax取2~20时，不同Lwin值的解译效果 Fig. 11 The extracting effects with different Lwin(Lmin-Lmax：2~20)

为了寻找更合适的参数化方案，在表 3的基础上将L_win的取值进一步细分(表 6).结合剖面解译效果(图 12)分析发现，当L_win=70、80时，a1/t>0.65，a2/a1>0.8，提取得到的边界点的密集度大，地层界线连续性好，剖面中可以较为清晰的分辨出地层界线，且与人工解译的地层界线吻合较好.

图 12

Fig. 12

图 12 Lmin-Lmax取2~20时，不同值Lwin的解译效果 Fig. 12 The extracting effects with different Lwin(Lmin-Lmax：2~20)

表 6(Table 6)

表 6 Lmin-Lmax取2~20时不同Lwin的提取效果 Table 6 The extracting effects with different Lwin(Lmin-Lmax：2~20) 2~20 50 60 70 80 90 100 t 8210 8210 8210 8210 8210 8210 a1 7399 7222 6904 6502 6111 5694 a2 5793 5716 5584 5337 5070 4744 a1/t 0.7056 0.6962 0.6801 0.6501 0.6175 0.5778 a2/a1 0.7829 0.7915 0.8088 0.8208 0.8291 0.8332

表 6 Lmin-Lmax取2~20时不同Lwin的提取效果 Table 6 The extracting effects with different Lwin(Lmin-Lmax：2~20)

综上分析，当L_min-L_max取2~20，L_win的取值范围在70~80之间时，自动解译剖面能够得到明显且连续的地层界线，解译效果最理想，是变长滑动时窗能量比法用于自动解译黑海浅地层剖面数据的最优参数化方案. 3.5 实例验证

为验证以上参数的可靠性，选取黑海另外两条地层分层较为明显的剖面(4ns、4we)进行解译.解译参数L_min-L_max取2~20，L_win取76.

图 13为2条剖面的自动提取的地层界线与人工识别的地层界线的对比图.剖面4 ns中自动提取的地层界线与人工识别的地层界线吻合较好，25 m以浅的地层中，地层界线S2、S3吻合度非常好；25 m~45 m深度范围内，S4、S5的走向一致；45 m深度以下，地震反射信号揭示的沉积物分层已不明显，难以识别明显的地层界线.剖面4we的一段连续起伏的埋藏谷地，自动解译能够提取出不整合面S2，30~45 m深度范围内S3的趋势也基本吻合.

图 13

Fig. 13

图 13 自动提取的地层界线与人工识别的地层界线对比图(a1、a2为剖面4ns；b1、b2为剖面4we) Fig. 13 The comparison of the extracted boundary lines between MVI and CAI methods(a1、a2：profile 4ns in Heihai Lake；b1、b2： profile 4we in Heihai Lake)

本文在传统的能量比法的基础上，选取变长滑动时窗计算能量比值，在边界点提取时采用均值曲线法，选取合适的窗口长度进行边界点的判别；进行剖面地层界线提取时再进行边界点的连续性检验，剔除不满足条件的离散孤立点：

(1)经过对传统能量比法的优化，提高了浅地层剖面地层界线自动提取的精度.

(2)通过设计不同的参数化方案进行实验分析，能够确定适合于研究区域的浅地层剖面数据进行自动解译的关键参数的最优取值，提高了自动解译中边界点的密集度，以及与人工解译获得的地层界线的相似度.

(3)运用这一方法能够快速处理同一地区获得的大量浅地层剖面数据，自动解译提取得到的地层界线，可以帮助建立研究区的基本地层结构框架，为进一步的人工判读和详细分析提供重要基础.

在不同地区进行野外工作时，使用的浅地层剖面探测仪可能不同，各项主要参数需按不同地区的地质条件、水深、动力等不同情况进行设定，发射功率也可能存在差异，因此所获得的剖面数据也不尽相同.本文所得到的关键参数(L_min-L_max，L_win)的最优取值，是以黑海地震剖面数据为对象，这组参数是否适合于其他地区的剖面解译还有待进一步研究.

在采集反射信号时，有效信号中不可避免的会掺入周围环境中的噪声信号(庄杰枣等，1996)，这对剖面的自动解译会产生很大影响，降低解译精度.本文尝试用黑海的参数(L_min-L_max=2~20，L_win=76)处理在南黄海辐射沙脊群北部西洋潮流通道采集的剖面数据时，由于周围环境中船舶以及水体中较多的噪声一同被水听器接收，并掺杂到有效的地震波反射信号中(图 14)，噪声影响大，结果无法提取出有效的地层分层信息(图 15).

图 14

Fig. 14

图 14 西洋剖面与黑海剖面单条波形对比图 Fig. 14 The comparison of single record oscillograms between Xiyang profiles and Heihai lake profiles

图 15

Fig. 15

图 15 西洋剖面的自动解译图 Fig. 15 he result of CAI method for Xiyang profile

因此，如果能在进行能量比值计算之前，先对原始信号进行降噪处理(李肃义等，2013；吴招才和刘天佑，2008)，去除原始剖面数据中的部分环境噪声信号，不仅会拓展该方法的应用范围，而且也将进一步提高解译精度，这也是下一步将要进行的研究工作.