1 引言

地震相分析是进行储层预测和储层特征描述的重要基础工作之一(Nivlet，2007；John et al.，2008)，通过分析地震剖面的反射结构、连续性、振幅、频率、相位和速度等地震反射参数来进行定性地质解释(Payton，1980)．现在一般直接利用地震反演、地震属性及其聚类和分类结果进行定量地质解释，通过井震标定建立地震相与地质体之间的联系，进而对地质体的构造特征、沉积相带、岩石成分、储层物性和流体特征进行精细描述．地震相分析已成为储层预测和油藏描述的有效手段(Coléou et al.，2003；Chopra and Marfurt，2005；de Matos et al.，2007；Bagheri and Riahi，2013)．地震属性是地震相分析中常用的方法，目前已有上百种属性(Brown，1996；Chen and Sidney，1997)，其中一些参数对油气识别具有较好的指示意义(林昌荣和王尚旭，2011)．基于地震属性的地震相分析也面临着一定的困难：地震属性种类繁多，地震属性与地质特征之间并不存在简单的对应关系，有些地震属性有明确的物理意义，有些则是冗余甚至无效．因此，从众多地震属性中优选出对储层特征敏感的地震属性是地震相分析重要的环节，实际中常采用多属性交汇和降维方法选择与储层相关的信息(Dumay and Fournier，1988；Kalkomey，1997；de Matos et al.，2007；John et al.，2008; Roweis and Saul，2000；龚灏等，2008；Liu et al.，2010)．多属性分析通过属性优选构建地震相分析的高维属性，采用模式识别、图像处理、神经网络和统计学等领域中的技术进行地震相划分(Saggaf et al.，2003；杨勇等，2005；唐耀华等，2009；Bagheri and Riahai，2013；田玉昆等，2013；张等，2015；袁成等，2016)，可以提高地震相分析效率．

地震数据中包含丰富的地质信息，从波形上常常不能直接识别．提取隐藏其中的地质特征，需要对地震数据进行分解，以提取各种地震属性．例如，时频分析将数据分解到不同的频段以获取各频段中不同的特征．地震属性分析采用多种属性来反映隐藏在地震数据中的地质特征，但这些属性依据不同的理论知识，采用不同方法计算，彼此之间缺少紧密的内在联系，属性选择过程中人为因素较多，易导致预测结果的多解性(Coléou et al.，2003；印兴耀和周静毅，2005)．在弹性介质中传播的声波和地震波遵从相同的微分方程，二者本质上是同一种波．实际上在地震激发时(例如炸药激发)是有声音的，由于地震波在地层的传播过程中高频能量损失，使得地震记录主体能量位于人的听觉频段范围之外或边缘，从而忽略了地震波的“听觉特征”．语音识别技术经过多年的发展，目前已经从实验室走向实际应用(Cano et al.，2005；韩纪庆等，2013)，已形成了一套基于语音高维特征参数的语音识别方法，在孤立词、说话人和连续语音识别等方面取得了较好的效果(Reynolds and Rose，1995；Pierre-Yves，2003；Muda et al.，2010)．人类语音含有丰富的信息，语音识别技术能对这些信息进行识别，说明语音特征参数对语音信号具有很强的信息表征能力．因而借鉴语音识别中的特征参数提取方法，用以对地震数据进行分解，将数据中相互叠加的信息展开到不同维度，找到隐藏其中的地质特征，将有助于地震资料的处理和解释．曹俊兴等(2011)、Tian和Cao(2011)采用语音识别中的倒谱域参数对过井剖面地震数据进行分析，发现含气井存在一阶倒谱域参数低值、二阶参数高值的现象．地震数据中包含的地质信息较多，采用高维特征参数对地震数据进行分解将更为有效．线性预测倒谱系数(Linear Prediction Cepstrum Coefficient，LPCC)是目前语音识别中广泛应用的一类特征参数，文中将借鉴LPCC系数的提取方法，对目的层段地震数据进行分解，将数据中的信息展开到多维空间，然后采用聚类方法进行地震相分析，挖掘地震数据中隐藏的地质特征．

2 地震线性预测倒谱系数

语音由声音和语言构成，声带激励的声音经过声道时，人的神经系统控制声道的形状和运动对声音进行调制，最终形成语音，因而声道的响应构成语言．基于声道模型，从语音信号中获得声带激励和声道响应信息对语音识别有着非常重要的意义(韩纪庆等，2013)．震源激发的子波经过地层时，受到地层吸收衰减、干涉、反射、透射等物理过程的改造，最终经检波器接收形成地震记录．地震记录和语音的形成过程存在一定的相似性，即震源激发可类比声带激励，地层对地震波的响应可类比声道对声音的响应，语音识别关注与声带激励对应的声音和与声道响应对应的语言，地震勘探则更多地关注与震源关联的地震子波和与地层响应对应的反射系数．因此根据地震记录将子波和反射系数信息进行分离对地震资料处理也有着重要的意义．在褶积模型假设下，地震记录可视为地震子波和地层反射系数的褶积．反褶积算法一般有两大类，一类为基于模型参数的线性预测分析，另一类为基于同态处理的倒谱域分析．对于语音信号，在线性预测分析计算信号频谱的基础上进行倒谱域分析，构成LPCC系数，作为语音识别的特征参数(屈丹等，2010；韩纪庆等，2013)．

根据地震记录褶积模型，地震记录 s(t) 可视为地震子波 w(t) 和反射系数 r(t) 的褶积：

(1)

计算(1)式两边信号的频谱，则时间域的褶积信号在频率域中变为相乘信号：

(2)

式中 S(ω)、W(ω)和R(ω) 分别表示地震记录、子波和反射系数的频谱．对(2)式两端取对数，则频率域中的相乘转变为相加：

(3)

对(3)式两端进行傅里叶逆变换，就得到地震记录的倒谱域参数：

(4)

式中FT^－1表示傅里叶逆变换．由此可见，地震记录倒谱域参数可视为是地震记录对数频谱的“频谱”，从频率域经傅里叶逆变换后为时间域，因而倒谱域参数是时间域信息．在倒谱域中，地震记录倒谱域参数为子波和反射系数倒谱域参数之和．

下面采用褶积模型来分析地震数据倒谱域参数是如何实现子波和反射系数信息分离的．假定地震子波为主频为30 Hz的Ricker子波(图 1a)，反射系数为随机反射系数序列(图 1b)，将地震子波和反射系数进行褶积运算可得到合成地震记录(图 1c)，地震子波、反射系数和地震记录的频谱与倒谱域参数见图 1．对比子波频谱(图 1d)、反射系数频谱(图 1e)和地震记录频谱(图 1f)，可以发现从地震记录频谱中直接分离子波和反射系数信息是困难的．转换至倒谱域后，子波倒谱域参数主要位于低时域段且幅值较大(图 1g)，反射系数倒谱域参数分布在整个时域段但幅值相对较小(图 1h)，地震记录倒谱域参数(图 1i)可视为子波倒谱域参数与反射系数倒谱域参数之和．因此，可以从一定位置将倒谱域分为低时域段和高时域段，低时域段主要为地震子波信息，高时域段反射系数信息占主体地位，从而很大程度上地实现了子波和反射系数信息的分离．同态反褶积就是基于这一构想，将低时域段倒谱域参数逆变换至时间域获得子波，然后对地震记录进行反褶积处理(凌云等，1995；蒋立等，2012)．在褶积模型假设下，地层特征由反射系数来表征，如果按目的层段截取储层地震数据并计算其倒谱域参数，那么这些参数本身也可以作为储层的特征参数．

倒谱域分析是在地震数据频谱的基础上进行的．数据频谱分析通常采用傅里叶变换，该方法在计算中隐含着窗函数之外未观测到的数据为零，只有当待分析信号为平稳信号且具有一定观测长度时，才能保证频率分辨率、谱估计稳定性和足够多的频点．按目的层段截取的地震数据为非平稳信号，且数据长度一般较短，采用傅里叶变换不能为倒谱域分析提供足够的频点和频谱稳定性，因而采用基于模型的谱估计方法更为合适．线性预测分析是模型谱估计中常用的一种方法，将反射系数视为一个全极点AR模型H(z)，地震记录s(n)则是子波w(n)激励该模型产生的输出，模型的传递函数为(韩纪庆等，2013)：

(5)

式中G为增益常数，P为模型阶数，a_i为模型系数．用模型系数可构建一个P阶预测系统，用地震记录的前P个数据来预测当前数据的值 s(n)：

(6)

预测均方误差为：

(7)

在求解出模型系数 a_i 后可以计算地震记录的频谱：

(8)

式中T为地震记录采样周期，E_P为模型阶数为P时的最小预测误差能量，频率f可在满足Nyquist采样定量范围内(0≤f≤1/(2T))任意取值，从而为倒谱域分析提供足够的频点．线性预测分析从观测数据出发，依据最小预测均方误差准则，估计出模型系数a_i，这些系数称为线性预测系数(LPC系数)．计算LPC系数a_i和最小预测误差能量E_P常用的方法有基于自相关的Levinson-Durbin递推算法、基于协方差的Choleskey分解法和直接由观测数据递推的Burg算法(王世一，2013)．

由于地震记录信号的极点未知，因此模型的阶数P不是先验已知的．模型阶数过低时，对地震记录的频谱估计过于平滑，不能提供足够的分辨率；而模型阶数过高时，频谱估计结果会产生震荡，出现谱分裂现象．为了消除谱分裂，提高谱估计的分辨率，需要选择合适的模型阶数，这里采用自回归传递函数准则(CAT)来选择模型阶数(Parzen，1976)：

(9)

式中N为地震记录长度，E_k为模型阶数为k时的最小预测误差能量，可取使(9)式为最小值时的阶数P作为模型最优阶数．在计算出地震记录频谱S(f)后可采用傅里叶逆变换计算LPCC系数，实际中一般采用如下递推关系直接予以计算(屈丹等，2010; 韩纪庆等，2013)：

(10)

式中n称为LPCC系数的阶数，P是谱估计模型的阶数．LPCC系数提供了一种将地震数据分解成语音特征参数的方法，可以较好地反映地震数据中隐藏的地质特征，为后续的地震相划分提供较好的特征参数．而传统的地震属性分析多达上百种属性，优选的多属性组合因人而异．

3 K均值聚类

在聚类分析中，K均值聚类算法(K-means)是无监督分类中应用范围最为广泛的一种基本方法(Selim and Ismail，1984)，能够使聚类域中所有样本到聚类中心距离的平方和最小．其实现原理是通过迭代，逐次计算每个样品到K个聚类中心的距离，将样品归入与之距离最小的类别中，更新类别的中心，直到所有样品不再改变类别或满足迭代终止条件．实现流程为：

(1) 设定聚类中心个数K；

(2) 从所有样品中随机选择K个样品作为初始聚类中心；

(3) 依次计算每个样品到聚类中心的距离；

(4) 将样品归入与之距离最小的类别中；

(5) 更新类别的聚类中心；

(6) 如果满足迭代终止条件，停止迭代，输出聚类结果；否则重复(3)~(5)直至满足迭代终止条件．

聚类中心通常定义为各类别中所有样品特征参数的平均值，样品与聚类中心距离多采用欧式距离．K均值算法聚类的结果易受到聚类中心个数和初始聚类中心选择的影响，如果样品在特征参数空间能形成相对孤立的分布，K均值聚类能取得较好的结果(Brown and Huntley，1992)．

4 模型正演分析

为了分析地震语音参数在地震相分析中的有效性，我们设计了一个四层地震模型，考虑三个方面的因素，包括反射系数、地震子波主频和地震信号信噪比的变化．为简单起见，模型厚度用时间表示，模型的参数见图 2a，第一、第三和第四层为横向均匀地层，速度分别为4000 m·s^-1、4000 m·s^-1和5000 m·s^-1．第二层横向上分为三种速度差异较大的介质，纵向上速度存在小幅度随机变化，三种介质的速度分别为2900+50×rand(-1，1)m·s^-1、3300+50×rand(-1，1)m·s^-1和3700+50×rand(-1，1)m·s^-1，其中rand(-1，1)表示位于-1至1之间的随机数，相应模型的反射系数见图 2b，模型剖面在横向上可分为3种反射系数序列组合．地震子波采用Ricker子波，主频从20 Hz开始，以0.25 Hz为间隔上升至50 Hz，对每一种反射系数序列，合成得到121道具有不同主频子波的地震记录，模型包含3组不同的反射系数序列组合，最终得到363道合成记录. 按第二层时间范围截取地震数据，对所有截取的地震数据计算24阶LPCC系数并做归一化处理．不失一般性，图 2中只显示了子波主频从20~50 Hz偶数主频的结果，共48道，图 2c、图 2e和图 2g分别是不含噪声、含10%噪声和含20%噪声时的合成地震记录，图 2d、2f、2h分别是相应地震记录计算得到的归一化LPCC系数．

分别以4至24阶(间隔为2阶)LPCC系数作为地震道的特征参数，采用K均值方法对地震记录不含噪声、含10%噪声和含20%噪声三种情况进行聚类，并计算对应聚类结果的正确率，结果如图 3a所示．合成地震记录不含噪声时，在阶数从4阶增加至12阶的过程中，聚类结果正确率上升很快；在12阶以上，随着阶数的增加，正确率高于95%，且变化很小；地震记录含10%和20%噪声时，从4阶增加至12阶的过程中，正确率也上升很快，12阶以上随阶数增加，正确率分别约为90%和80%左右，且起伏变化幅度也较小．对比不含噪声、含10%噪声和20%噪声的结果，地震记录含噪声越少，聚类结果越好．总体来看，在阶数大于12阶之后，正确率随阶数增加变化较小，表明采用12阶以上的LPCC系数作为地震数据特征参数能取得较好的聚类效果．

为了分析地震子波主频对不同阶数LPCC系数聚类结果的影响，以2 Hz为频段间隔统计子波主频不同时各频段聚类结果的正确率，图 3b、图 3c和图 3d分别是地震记录不含噪声、含10%噪声和含20%噪声时的统计结果．从图中可以看出，在没有噪声的情况下，6阶及以上聚类结果中高频段正确率为100%，错误主要出现在低频段．对于含10%和20%噪声的情况，在阶数较低时正确率分布规律不是很明显，但随着阶数的增加，总体而言高频段聚类结果的正确率显著高于低频段．对比不含噪声、含10%噪声和20%噪声的结果可以看出，采用LPCC系数作为地震道特征参数时，子波主频越高，越有利于地震相划分．

子波主频影响着地震合成记录对反射系数信息的分辨率，合成记录的LPCC系数与反射系数LPCC系数愈接近，就愈能反映地层信息，更有利于地震相分析．理想情况下，如果子波为单位脉冲信号，则合成记录将真实地反映地层反射系数，此时合成记录的LPCC系数也将全部反映反射系数信息．从图 2b中截取第二层横向3种介质的反射系数序列并计算其各自的LPCC系数，采用主频15至60 Hz的Ricker子波分别与各个反射系数褶积形成地震合成记录，计算合成记录的LPCC系数，采用欧式距离和相关系数来分析子波主频不同时合成记录LPCC系数与反射系数LPCC系数之间的接近程度，结果如图 4所示．尽管某些反射系数与合成地震记录LPCC系数之间的距离随频率升高并非严格的单调下降，但是总体而言，主频越高二者之间的距离愈小，相关系数越高，二者越接近，这进一步说明子波主频愈高，LPCC系数越有利于地震相划分．

5 实际资料处理

研究区块位于塔里木盆地中部A井区，目的层段主要为良里塔格组的礁滩相碳酸盐岩储层．储层类型主要为裂缝-孔洞型，岩性主要为泥质灰岩，埋深约为5000 m，储集层非均质性强，油气藏相态及分布复杂．良里塔格组总体构造格局为南东高、北西低的北倾斜坡与鼻状构造，该斜坡被北东向走滑断裂切割，早期存在平行于坡折带呈NW—SE向断裂．图 5a为区块内一条Xline地震剖面，良里塔格组位于层位mfs3与SB3之间，其中塔中I号台缘坡折带形成于奥陶系中晚期，发育了一套纵向上多旋回、横向上多期加积的较大规模礁滩相沉积体系的主力产层(刘嘉庆等，2012；赵文智等，2013)．图 5b为区块内一条Inline地震剖面，目的层段中间低凹区域为潮汐水道．图 6为区块的古地形图，礁滩相储层主要发育在台缘带波浪作用较强的水下低凸起区，水下低凸起之间的低势区为潮汐水道发育部位，潮汐水道通往台缘带内的低能带为泻湖沉积．受波浪作用与潮汐作用共同控制，良里塔格组礁滩体与潮汐水道相间分布．据钻探井资料显示，区块内有8口井，其中B2，B21和B23为工业油气井．

由于礁滩孔隙型储层主要受沉积相带和成岩作用控制，本文主要是在沉积层序解释的基础上，结合古地貌和地震数据LPCC系数来分析储层和沉积相的宏观展布，地震数据采用时间偏移后的三维数据体，具体流程如下：

(1) 以目的层段为窗口，截取层位mfs3与SB3之间的地震波形数据，剔除地震数据中的空值和野值；

(2) 对截取的每一个地震道波形数据计算24阶LPCC系数，并做归一化处理；

(3) 分别选取地震数据的12阶和24阶LPCC系数作为该道地震数据的特征参数，采用K均值方法对所有地震道进行聚类．

按目的层段截取地震波形数据可减少相邻地层的影响，使计算得到的LPCC系数尽可能地只反映目的层段的相关信息．采用K均值聚类方法是想检测即便使用最简单的聚类算法时，LPCC系数在地震相分析中能否取得较好的结果．图 7是24阶LPCC系数的空间分布，地震数据被渐进地分解成不同的地质特征，使得我们可以更细致地观察到断裂、礁滩相带、深水扇和储层的分布特征，特别是对一些地质现象边界的描述有独特的优势．1~4阶LPCC系数更多地反映坡折带边界、潮汐水道和断层分布，也可见台缘带内礁滩体储层的分布；5~7阶系数除礁滩体分布不明显外，坡折带边界、潮汐水道和断层分布可分辨，台缘带外深水扇边界逐渐呈现；8~16阶台缘带内礁滩体的分布开始显现，台缘带礁滩体的边界依稀可见；10~16阶对台缘带、坡折带、潮汐水道和深水扇边界描绘更为清楚，断层分布也有一定呈现，与图 6所示古地形图十分吻合；从边界的分布来看，台缘带泻湖一侧似乎存在另一个边界，9、11和12阶系数中台缘带内礁滩体分布也较为清晰；17~24阶边界清晰度下降，但依然可辨，礁滩体分布仍然明显．从1至24阶的逐渐演化来看，LPCC系数之间存在差异，同时不同阶次的LPCC系数在从不同维度描述储层特征时还存在不同程度的联系．

为便于比较，这里也对目的层段提取了地震属性分析中常用的瞬时频率、瞬时振幅和瞬时相位，如图 8所示．瞬时频率对台缘带内礁滩体储层刻画较好，坡折带边界可辨(图 8a)；瞬时振幅对坡折带和潮汐水道边界刻画较为清晰，也可见泻湖内断层分布(图 8b)；瞬时相位则对台缘带边界和礁滩体有较好的反映，也可见断层分布(图 8c)，整体形态与古地貌图也较为吻合．但是与LPCC系数相比，三瞬属性对台缘带的边界刻画不如LPCC系数清晰，且LPCC系数展现的不同边界是三瞬属性所没有的；三瞬属性对波折带下的深水扇分布反映较差，而从LPCC系数中可以清楚地看出深水扇的分布范围．从前文中对倒谱域参数的论述可知LPCC系数中的低阶部分更多地反映地震子波的特征，因而1~7阶LPCC系数与瞬时频率和振幅有些类似；高阶部分更多地反映地层信息，分布形态与反映地层不连续性的瞬时相位之间存在一定相似性．

采用K均值方法分别对前12阶和全部24阶LPCC系数进行聚类，结果示于图 9，其中图 9a是12阶聚类结果，图 9b是24阶聚类结果．可以看出两种结果对坡折带、台缘带、深水扇、潮汐水道和有利油气聚集区预测的几乎一样，较为明显的差别仅在于24阶结果在泻湖靠近台缘带附近有一带状分布，推测可能是随着LPCC系数阶次增加，高阶系数中包含的冗余信息对聚类结果产生干扰所致．从钻探井资料获悉，B2，B21和B23为工业油气井，而这三口井则位于基于LPCC系数地震相分析预测的有利油气聚集区(图 9中台缘带内红色区域)，从而也说明了利用LPCC系数进行地震相分析的可行性．

6 讨论

地震数据中包含丰富的地质信息，但从波形数据上通常难以直接识别．提取隐藏在地震数据中的信息，需要对数据进行分解，将数据中相互叠加的信息展开到不同维度，找到隐藏其中的地质特征．人类语音含有丰富的信息，语音识别技术能对这些信息进行有效识别，说明语音特征参数具有很强的信息表征能力．在弹性介质中传播的声波和地震波本质上是同一种波，可以借鉴语音识别的相关技术用于地震数据分析．

在语音识别中，一般情况下不同人发同一音节的语音时，声道的形状和运动是相似的，因而可以建立语音特征参数和发声音节之间的对应关系；不同人的声纹特征有差别，也可以通过提取说话人语音样本建立声纹特征参数和说话人之间的对应关系．语音识别是通过提取已知发声音节和说话人的特征参数并建立对应的信息库，然后结合一定的语言模式对语音进行识别．因此，语音识别是对已知对象的识别，即待识别语音信息包含在信息库中，信息库中没有的语音信息将无法予以识别．地震记录可视为一个褶积模型，从理论上讲，可以采用LPCC系数的计算方法提取对应的地震语音特征参数，分离子波和反射系数信息．但是，地震勘探所面临的则是地下未知的信息，信息的约束能力不如语音识别，很难系统地建立起地质现象与地震数据及其提取的特征参数之间对应关系．因而，地震相分析是一个不断完善的过程，需要结合尽可能多的已知信息，对地震相分析过程进行约束和对结果进行评估；在地震相分析中，更需要强调选用信息表征能力强、能反映地质现象的特征参数．地震LPCC系数提供了一种新的途径，将地震数据分解成语音特征参数，可以较好地反映地震数据中隐藏的地质特征，呈现地层细微的变化，为地震相划分提供较好的特征参数．在地震相分析阶段，采用LPCC系数作为特征参数的分析过程与采用其他地震属性或参数一样，需要结合已知的地质构造、钻井、古地形、古沉积环境和沉积层序等资料为分析过程提供先验信息，并对分析结果进行评估．

在本文中，我们采用的是经叠加和偏移处理后的时间域三维数据体，并以目的层段为窗口截取地震数据，进而计算每道地震数据的LPCC系数．在应用过程中也可以采用固定长度窗口，并逐步滑动时窗计算，从而每一地震道每个时间点有n阶系数，最终可得到对应的三维体LPCC系数．地震LPCC系数是基于褶积模型对子波和反射系数信息进行分离，因而从理论上而言，对可视为褶积模型的地震数据，我们可以采用LPCC系数的计算方法提取对应的地震语音特征参数．

7 结论

地震LPCC系数是基于线性预测分析的倒谱域高维特征参数，可在倒谱域实现子波和反射系数信息的分离，特别是对地质现象边界具有较好的描述能力，使得我们可以在不同的维度更细致地观察隐藏在地震数据中的地质特征．与三瞬属性相比，LPCC系数对储层特征(尤其是边界)的描述更为细致．与常规地震属性分析不同，不同阶次LPCC系数基于相同理论基础，在反映地震数据不同特征的同时参数之间还保持内在的联系．按目的层段截取储层地震数据计算LPCC系数，提取的特征参数主要包含储层信息，减少了储层以外地层对预测结果的影响．LPCC系数是对每一道地震数据单独计算的，不受相邻道或其他道集数据的影响．理论模型分析和实际资料应用表明基于LPCC系数分解的地震特征参数具有较高的地震相划分能力，因此，在地震相分析中采用LPCC系数作为特征参数，将有助于地震数据的处理和解释．

致谢

中国石油勘探开发研究院李劲松、孙夕平和李艳东三位老师对分析工作提出了许多建议；两位匿名审稿专家提出了许多宝贵的评审意见，对文章的修改和完善有很大的帮助．在此我们一并表示衷心的感谢