1 引言

薄互层是指不同岩性薄层交互出现的一种地质现象，其相邻两层的岩性不同，如砂泥岩薄互层.在砂泥岩薄互层发育地区，薄层的调谐作用以及薄互层组合方式空间分布的多变性，都能够影响薄互层地震波的振幅和相位特征，使得薄互层的地震反射特征更为复杂.当砂泥岩薄互层中夹杂砂岩油气储层时，薄互层的地震识别与描述技术具有重要的应用价值.

对于薄互层地震识别问题，Ricker (1953)最早开始研究了薄层地震反射问题.Widess (1973)研究了薄层的调谐效应，给出薄层定量化的具体定义：认为厚度小于入射子波在其介质中传播时波长的1/4为薄层；小于1/8时，复合反射波从波形上不可识别.Hilterman (1975)研究了薄层厚度对振幅的影响.Neidell (1979)等定义地震波长的四分之一为调谐厚度，认为薄层的厚度信息包含在振幅而不是波形中.Koefoed (1980)通过合成地震模型研究认为，薄层厚度和地震反射复合波的振幅之间存在着准线性关系.同时，模型中考虑了能量传播损失和层间多次波.Kallweit和Wood (1982)论述了地层厚度变化对振幅响应的影响；对于薄层楔状模型，当顶部和底部反射接近1/4波长时，振幅响应达到最大值，调谐厚度与频率有关.

对于各向同性介质，其分界面的反射、透射系数的精确解可由Zoeppritz (1919)方程计算.许多学者对Zoeppritz方程进行了不同形式的简化，发展了地震AVO (Amplitude versus offset) 技术，即通过分析反射波振幅随偏移距的变化规律来研究反射界面两侧介质的物性特征.传统的AVO分析，都是基于Zoeppritz线性近似方程提出的.Aki和Richard (1980)把Zoeppritz方程简化为入射角与纵波振幅相关的近似公式.之后，Shuey (1985)根据该简化公式做了进一步研究，将其改写由小角度项、适中角度项和广角度项三部分之和的近似公式.Bakken和Ursin (1998)，Wapenaar (1999)开展了关于薄层AVO效应的研究，对由于偏移距产生的调谐效应做出改进.然而，层状介质的AVO特征与单一界面是不同的，因为反射系数与入射波的频率相关.实际地层为非均匀介质，地震波在地下传播时具有传播效应和层间多次波的调谐干涉效应.对此，Aki和Richard (1980)提出的传播矩阵算法在理论上可以更好地描述这种层状介质模型.Kennett (1983)给出了地震波在层状各项异性介质中传播的描述.在此之后，有许多学者对层状介质进行了AVO正演响应研究.Carcione (2001)计算了富有机质的页岩层的反射系数.Liu和Schmitt (2003)研究了泊松比和层厚度对薄层AVO响应的影响.

李庆忠 (1987)研究了含油气砂岩层的频率特征及振幅特征，指出地震波的频率特性取决于地层的砂泥岩组合情况，子波的波形及地层的含油气性.孙树海和朱彦良 (1991)结合实际资料在详细分析砂泥岩薄互层反射特征基础上，提出气层的反射特征不仅有AVO正异常，还有AVO负异常及非正非负AVO异常.苏盛甫 (1988)通过对均匀介质和薄砂层中的不同碳酸盐岩的反射特征进行研究，提出调谐厚度是划分薄层和厚层的分界点，并利用薄层厚度与复合反射波的振幅为近似线性关系的特点，制作了薄层定量解释图版来计算薄层厚度和多层储集层的总厚度.汪恩华等 (2001)给出一种新的薄储层厚度的计算方法，综合考虑了振幅、主频和频宽等信息，计算出的厚度比单纯利用振幅时计算的精度高；同年，他们还利用反射系数谱理论对薄层反射进行了模拟，获得了薄层条件下不同频率成分的纵波反射系数谱数学关系式.Chen and Liu (2006)基于反射系数谱理论对薄层多波AVO进行了讨论.齐虹等 (2012)研究了砂岩薄互层的AVO正演模型.李雪英等 (2012, 2013) 基于波动理论，利用深度域相移法对不同互层数、不同单层厚度的等厚薄互层进行正演模拟，并采用广义S变换分析零偏移距地震道反射复合波的瞬时频谱.郭智奇等 (2016)基于广义传播矩阵理论计算地震波频变反射系数，完善了频变AVO算法，为含油气储层频变AVO响应的模拟和分析提供了方法.

本文针对研究地区致密砂泥岩薄互层地震识别与描述问题，首先开发了基于传播矩阵理论的高精度合成地震记录算法；而后，应用本文提出的地震波形反演算法进行了理论模型的试算分析；最后，开发了地震波形反演技术，反演标志层以下砂泥岩薄互层结构，反演结果可以精确到厚度为2 m砂岩分布，所得结果可用于有利储层预测.

2 研究区沉积概况与薄互层建模

图 1所示为研究区A井的测井曲线，包括伽马曲线、纵波速度、横波速度、纵横波速度比、密度以及孔隙度曲线，标志以下为砂泥岩薄互层.与泥岩相比，砂岩具有较高的纵、横波速度，较低的纵横波速度比、较低的密度，以及较高的孔隙度.根据测井数据分析，研究区为典型的砂泥岩二元结构.通过测井数据统计估计得出，研究区高速砂岩的纵波速度为3800 m·s^-1，横波速度为2000 m·s^-1，密度2450 kg·m^-3；低速泥岩纵波速度为3200 m·s^-1, 横波速度为1800 m·s^-1, 密度2500 kg·m^-3.本文使用以上这些数值建立模型，从而进行地震波形反演算法的测试，为实际应用提供理论依据.

3 基于传播矩阵理论的高精度合成地震记录算法 3.1 传播矩阵理论

由于干涉、调谐等现象的存在，来自薄互层的反射地震波呈现复合波型.常规单一界面模型的反射系数可由Zoeppritz方程计算，而具有层状结构模型的地震反射波特征不仅与入射角度、物性差异有关，还与入射波频率、地层厚度、薄互层结构、地层的不均匀性等因素有关 (Rokhlin and Wang, 1992; Kennet, 1983; Ursin and Stovas, 2002; Liu and Schmitt, 2003).

根据传播矩阵理论 (Carcione，2001)，对于P波入射，地层的反射、透射系数向量r=[R_PP, R_PS, T_PP, T_PS]^T由下式计算：

(1)

其中，矩阵A₁与A₂分别为与上、下层介质物性参数有关的传播矩阵；B_α(α=1, …, N) 为具有N层结构的中间薄互层的传播矩阵；i_P为P波入射向量，与入射介质物性参数有关；同时，上述矩阵和向量都是入射波频率、波慢度的函数.

传播矩阵A₁和A₂分别为

(2)

(3)

其中, i为虚数单位，ω为入射波频率，为薄互层总厚度；同时，变量β、γ、W和Z的下标P、S分别对应准纵波、准横波，1、2分别对应上、下介质，去掉下标简写各自表达式为

(4)

(5)

(6)

(7)

式中，p.v.意为取复数的主值.对于γ，符号“+”对应qP波 (即准P波)，符号“－”对应qS波 (即准S波).并且，s_x为水平波慢度：

(8)

(9)

s_z为垂直波慢度：

(10)

(11)

(12)

s_z表达式符号约定为：(+, －)：向下传播qP波，(+, +)：向下传播qS波; (－, －)：向上传播qP波，(－, +)：向上传播qS波.

传播矩阵B_α=T(0)T^-1(h_α), 其中

(13)

并且P波入射向量：

(14)

由薄互层反射系数的传播矩阵理论，可计算各个频率下的反射、透射系数向量r=[R_PP, R_PS, T_PP, T_PS]^T，即相应反射波的频变反射系数R_f.将频变反射系数与频率域的地震子波W_f相乘可得相应反射波的振幅谱U_f，即

(15)

对U_f做反傅里叶变换则可以获得时间域的反射波波形u_t

(16)

其中f表示角频率，i为虚数单位，t为时间.

3.2 高精度合成地震记录算法

基于传播矩阵理论进一步开发砂泥岩薄互层叠前、叠后地震高精度合成记录制作技术.与阻抗差界面反射的方法相比，该方法不仅考虑地层界面间的波阻抗差异，同时考虑了地层间地震波的传播及调谐干涉等动力学作用 (Carcione，2001；郭智奇，2008；郭智奇等, 2009, 2016)，图 2所示为计算得到的AVO反射波形，图中可观察到来自标志层界面及其下方的地震反射.

图 3给出了基于单阻抗差界面反射系数与子波褶积合成的记录、传播矩阵方法合成记录、以及与井旁地震道地震数据的对比.用实际地震数据与褶积模型和传播矩阵分别做互相关得到互相关系数分别为0.8609和0.8657，由此可知传播矩阵方法计算的合成记录在地震振幅、相位的细节上优于褶积模型的合成记录.

4 基于波形对比的薄互层地震反演技术 4.1 反演算法基本原理

如图 4所示为反演算法的基本流程：首先，根据测井解释资料确定砂岩、泥岩各自的速度和密度，建立模型空间砂泥岩的组合模式.如图 5所示，以关键参数 (N, X, Y) 定量描述砂岩层数、砂岩和泥岩空间位置的组合关系，设该模型为只有砂泥两种介质的韵律型薄互层，砂体厚度相同，而间距是变化的，两层砂岩的间距由X表示.其中N为单元反射体内砂岩的个数，Y为单元反射单元体内第一个砂岩层距离上界面的距离，并且设H为单元反射体的厚度，h为单元反射体内砂岩层的厚度, 反演是根据如图所示的流程，固定单元反射体的厚度H的值逐层进行反演，参数H的个数决定反演算法的时间，本文中根据研究地区的实际情况取H=10 m.根据测井数据给出的先验信息来确定h，本文以h=2 m厚度作为单个砂体单元，在厚度H=10 m的反射体内，最大砂地比N*h/H=40%, 根据以上分析可知，需要反演的参数为X、Y，以反射体单元为单位列出X、Y和N=0, 1, 2的所有可能模型，再用传播矩阵方法计算上述模型空间中所有模型全频带上的频变AVO响应；接下来，输入叠前或叠后地震数据，计算公式 (17) 定义的目标函数，求取实际地震数据与上述模型数据空间的最大相关，并且输出由目标函数确定的关键参数向量v=(N, X, Y)，并以此确定薄互层砂泥岩组合结构，公式 (17) 中M为叠前数据的道数，当M=1时，输入的为叠后数据；最后，逐层逐道反演地震数据，得到整个地震剖面反演结果.

(17)

4.2 理论模型测试

如图 5所示，根据测井数据设计了三个研究区典型模式的薄互层模型.模型中红色部分表示砂岩，厚度分别为2 m和4 m，V_P=3800 m·s^-1，V_S=2000 m·s^-1，ρ=2450 kg·m^-3；蓝色部分表示泥岩，V_P=3200 m·s^-1, V_S=1800 m·s^-1, ρ=2500 kg·m^-3.模式1中三层砂岩应的厚度分别为4 m、2 m、2 m，相邻两层砂岩的间距分别为6 m和8 m；模式2中三层砂岩应的厚度分别为4 m、2 m、4 m，相邻两层砂岩的间距分别为6 m和8 m；模式3中三层砂岩应的厚度分别为4 m、2 m、4 m，相邻两层砂岩的间距分别为10 m和4 m.正演算法采用传播矩阵理论，地震子波是主频为40 Hz的雷克子波.从正演结果可以看出，虽然三种模型所对应的合成记录中波形不同，但是无法直观准确的的分辨出薄层的位置和厚度.

输入正演得到的合成地震记录应用图 4所示的算法进行反演，为了更加直观的对比原速度模型与反演结果的吻合度，选取每个模式其中一道的反演结果与原模型进行对比.图 5中黄色曲线为反演所得速度曲线，曲线高幅度处为高速砂岩的速度，低幅度处为低速的泥岩速度，黄色曲线所对应的位置和原模型中砂岩泥岩分布基本一致.计算图 5中反演模型正演得到的反射波形与输入波形之间的互相关系数，验证反演模型对应的地震波形与实际数据的匹配程度.计算结果表明，模式1、模式2、模式3对应的互相关系数分别为0.9924、0.9713、0.7320.互相关系数越接近1说明匹配度越好，可知反演算法中正演模型的波形与输入波形基本吻合.由以上分析可得出，本文提出的地震波形反演算法在理论上能较好的反演出薄互层位置、结构及厚度等参数.

4.3 实际数据应用

图 6所示为过井A的叠后地震剖面，其中黄色曲线为井A的声波时差曲线.为方便选取标志层下方的薄互层结构进行研究，按标志层位进行了拉平显示.根据上述反演算法，由过井A井的叠后地震剖面反演砂泥岩分布模式，反演结果如图 7以纵波速度剖面显示，其中砂岩为高速层；可以观察到，计算得到的高速砂体分布与井中观测结果吻合度较好，能够反映出主要砂体的分布位置，反演结果具有较高空间分辨率 (图 7白色区域为测井曲线按实际深度显示，图中纵坐标0~170 m为从标志层算起).图 8分部分选取井A的声波时差曲线进行局部放大来与过井A反演结果对比，可以看到反演的砂岩位置与测井曲线解释中主要砂体位置具有较高的吻合程度.

图 9为过井A井的薄互层地震反演结果与地震剖面的联合显示，其中高亮显示部分为砂体分布，砂体空间位置根据图 7的反演结果进行了深度-时间转换.从图中可以观察到地震反射波形、相位与砂岩层位的对应关系，可通过砂体厚度、空间垂向组合关系、空间横向连续性、砂层规模等预测为区分“好、中、差”储层提供依据.

5 结论

本文结合测井数据给出的先验信息，提出了基于波形对比的砂泥岩薄互层反演方法.反演算法设计了一个以砂泥岩薄互层为反射单元的模型，采用传播矩阵算法作为正演方法，更精确的考虑了反射波动力学特征，在局部振幅细节上优于褶积方法，理论模型试算证明了反演的准确度较好；而后将该方法应用于实际工区薄互层的反演，结果表明，该方法对薄互层中砂泥岩组合结构、砂岩空间展布具有一定的预测能力，反演结果可为实际生产中寻找潜在砂体油气储层提供参考.