我国西部碳酸盐岩风险勘探研究中，“串珠状”溶洞体、断裂等小尺度地质体在油藏单元的形成过程中起着重要作用.因此，研究不连续及非均质地质体检测方法在油气资源储层评价中具有重要的应用价值.在石油工业界中，目前用于断层检测的方法有多种，其中地震相干技术的出现为断层自动识别提供了坚实的基础.第一代相干体算法由 Bahorich和Farmer(1995)提出，该算法基于经典的归一化互相关理论，对地震数据的品质要求较高.为适应低信噪比数据，Marfurt等(1998)提出第二代相干体算法，通过引入协方差矩阵多道相似计算提高了抗噪声能力，但均值效应却降低了地震解释的横向分辨率.Gersztenkorn和Marfurt(1999)通过计算协方差矩阵中的特征值提出了第三代相干体算法，该算法通过子空间运算使得地震解释结果既消除了噪声又提高了横向分辨率.地震相干技术的出现同时也推动了一批不连续性属性检测技术的发展，如地震属性方法(李春峰和克里斯托弗，2005；张国栋等，2010；李建雄等，2011；方海飞等，2013；陈海清等，2014； 霍丽娜等，2014；王洪求等，2014)、谱分解方法(廖新华等，1990；魏志平，2009；朱振宇等，2009；陈波等，2010；陈红等，2011；黄捍东等，2012；于豪等，2013；蔡涵鹏等，2014；王贻朋等，2014； 张晓燕等，2014)、蚂蚁追踪法(张欣，2010；王军等，2013；霍丽娜等，2014) 等.为消除断裂带上的散射和绕射影响，一些学者引用了混沌属性(王伟，2009) 用于衡量振幅值的规律性和混乱性.廖新华等(1990)提出基于谱分解的断层识别方法，后来张亚中等(2006)将谱分解技术应用于碳酸盐岩储层预测中.然而就目前国内碳酸盐岩风险勘探而言，关注点除了小断层、断裂等不连续构造之外，还有岩溶溶洞等非均质体，实际上它们与储层特性更为紧密相关.本文研究一种能够同时考虑不连续体构造与非均质体的检测方法.鉴于该两类地质特征在数学上可用L₁非光滑范数统一表示，因此探讨一种基于地震成像数据稀疏反演的不连续及非均质地质体检测方法.

小尺度不连续及非均质地质体的地震响应为弱信号，因此只有去除强振幅反射同相轴才能有效提取隐藏的信息.平面波破坏滤波器利用局部平面波模型表征地震数据，在线性同相轴预测上具有很好的应用.本文由平面波破坏滤波器预测出振幅稳定的反射同相轴，该类数据对应地震数据中的光滑构造特征.地震成像数据减去预测出的光滑数据后，剩余的数据中包含了非线性构造特征和噪声.在此基础上，通过非线性加强算子构建L₁稀疏反演模型，求解出不连续非均质信息.

平面波破坏方法(Fomel，2002) 利用局部平面波方程描述地震数据构造特征，可用于表示光滑构造，其表达式为

(1)

其中，d表示剩余地震数据， s=[s₁，s₂，…，s_N，]^T 表示地震成像数据， R(σ) 为线性平面波算子， σ 表示地震倾角斜率，公式可详细表示为

(2)

其中， I 为单位矢量， σ 为倾角斜率矢量， P_i，j(σ_i) 表示由第 i 道地震数据预测第 j 道地震数据.

方程中， s 为给定的地震成像数据，P_i，j(σ_i)为滤波器.在地震倾角斜率求解过程中，需要给定初始倾角斜率，然后由共轭梯度法逐步迭代实现.由于本文只考虑地震数据中的光滑信息，因此选择了低阶滤波器，公式为

(3)

其中， Z_t，Z_x 分别为变量 t，x的Z 变换.

考虑到断层、裂缝、溶洞等地质异常体倾角斜率无法准确估计，因此在上述预测方法中该类信息被分离出来.经过上述平面波破坏滤波后，得出的剩余地震数据 d 中包含了地震数据中的不连续及非均质地震信息.

由于不连续及非均质地震信息为非线性构造，因此要想从剩余地震数据d中分离出有效信息，需要利用非线性构造加强滤波方法(Liu et al.，2010； 刘洋等，2014)，这里用算子D(σ) 表示，它可以是高斯相似均值滤波(Liu et al.，2009) 或下-上-中滤波(Hardie and Boncelet，1993)，本文选用下-上-中滤波(Liu et al.，2010)，公式为

(4)

其中，r^*为下-上-中滤波结果，k为平滑参数，l为锐化参数，N为样点个数，s^*为样点中间值.

由于不连续及非均质信息与噪声混叠在一起，因此分离出高信噪比的地质异常体信息，需要考虑反演方法.不连续及非均质信息具有稀疏性特点，采用L₁范数模型进行求解，其极小化问题为

(5)

其中， m 为估计模型.

引入拉格朗日算子后，该方法等价于求解如下无约束问题，公式为

(6)

其中，α为可调节的正则化参数， Ω(m) 可定义(王彦飞等，2011) 为

(7)

其中 定义为分段函数，公式为

(8)

其中，$\tilde \varepsilon $为微小的常数，当 $\tilde \varepsilon $→0 时，非光滑泛函 Ω(m)近似l₁ 范数.

设目标函数 J^α(m) 的Hessian矩阵为 H(m)， 由于地震数据规模较大，因此本文采用拟牛顿法求解.拟牛顿法只需要计算一阶偏导数，计算简单且运算量小，比解析求得Hessian 矩阵更加高效稳定.

目标函数 J^α(m) 是一个非线性方程， J^α(m) 的梯度函数为

(9)

对上式求导，可得 J^α(m) 的梯度函数为

(10)

式中 K(m) 定义为

(11)

其中，的导数定义为

(12)

记 B≈H^-1(m)，函数J^α(m) 的梯度 g_k 为第 k 次迭代的梯度 g_k=g(m_k)， 方向为 s_k=(m)_k+1－(m)_k， 则模型参数更新为(Wang et al.，2011)

(13)

其中，矩阵 B_k 迭代求解公式为

(14)

ω_k 为线性试探步长，利用下述的Armijo-Goldstein准则非精确求解获得(袁亚湘和孙文瑜，1997).

综上所述，得出逼近L₁范数的拟牛顿算法如下：

(1) 给定初始值m₀、常数r及单位矩阵B₀=I，并令k=0.

(2) 若k≠0，计算m_k及B_k.

(3) 计算 g(m_k)，若‖g(m_k)‖＜ξ×max{1，‖g(m₀)‖}，ξ ＜0.1， 则停止迭代，否则执行步骤4.

(4) 计算 .

(5) 令k=k+1，更新步骤[步骤2]中的 m_k及B_k， 并转到步骤3.

需要注意的是停机准则为反演算法一个关键问题，对于反演问题，一般算法的精度不可能超过所提供的数据误差的精度(王彦飞，2007).采用常规的梯度能量值的大小作为停机准则，可能使得迭代陷入循环，导致计算误差的累积，这里选取的停机准则是：“若 ‖g(m_k)‖ ＜ξ×max{1，‖g(m₀)‖}， 则停止迭代”，其中 ξ 为控制参数，取值一般小于0.1.

为测试基于地震成像数据稀疏反演的不连续及非均质检测方法，模型设计上包含断点、溶洞散射体、粗糙反射界面等信息.该类地质体在地震波场模拟中产生绕射波，而绕射波为不连续及非均质信息载体，经偏移后聚焦.但聚焦后的绕射信息与反射波信息在振幅和相位上均在很大差异.与反射信息相比，绕射信息很弱，因此在常规成像结果中难以突显.图 1为缝洞模型的偏移结果，主要展现层位宏观信息，而小尺度溶洞和断裂信息难以体现.利用平面波预测方法可以估计出平滑、振幅稳定的反射信息，缝洞模型平面波预测结果如图 2所示.由于该算法采用低阶滤波器，因此非光滑信息包含在剩余波场中.对剩余数据进行非线性稀疏反演，结果如图 3所示.对比图 1和图 3可见，浅层4串岩溶和深层3串岩溶在图 3中能够显现，并且图 3中提取出断点和断裂等不连续信息.在计算过程中，正则化参数 α 对反演结果起着关键作用： α 越大稀疏约束越强， α 越小平滑特征越明显.本文给出的方法利用数学上的L₁稀疏范数，因此在原理上能够更好保护不连续及非均质地质信息.与真实设计模型相比，稀疏不连续及非均质提取方法在断点、溶洞和断裂表征上具有明显优势，并且能够去除偏移噪声.虽然图 3中有些地方存在残余线性同相轴影子，但这是由于该同相轴振幅存在微小变化的原因，这也反应稀疏不连续及非均质提取方法的灵敏性.去除连续地质信息后，图 3展示沉积间断点、断层、溶洞、粗糙反射界面等信息.该类地质信息在油气藏储层研究中具有重要应用价值.

图 1

Fig. 1

图 1 缝洞模型深度偏移结果 Fig. 1 Depth migration of a fracture-cavity model

图 2

Fig. 2

图 2 反射信息平面波预测结果 Fig. 2 Reflection information predicted by plane-wave method

图 3

Fig. 3

图 3 基于稀疏反演提取出的不连续及非均质地质信息 Fig. 3 Discontinuous and inhomogeneous geological information extracted by sparsity inversion-based method

缝洞模型测试结果表明，基于稀疏反演的不连续及非均质地质体检测方法，能够清楚刻画与油气运移通道有关的地质构造及溶洞储层信息.

国外碳酸盐岩油气藏普遍以相对均质的生物礁、白云岩风化壳等为有效储层，构造较为完整，较少涉及储层非均质问题.而国内碳酸盐岩油气藏，存在年代久远，极强的非均质性是多数碳酸盐岩油气藏典型特征.尤其在塔里木盆地风险勘探中，能够预测出与古风化壳有关的大型岩溶缝洞群，对寻找规模性油气藏具有重要意义.本实例研究难点为如何清楚刻画出与奥陶系有关的断裂、溶洞等不连续及非均质地质体.

三维叠前时间偏移成像结果如图 4所示，图中3.5 s位置为风化壳界面，4.0 s为“串珠状”溶洞发育位置.虽然图 4在一定程度上体现风化壳界面振幅变化情况和断面形态，同时显示出溶洞影子，但对顺层岩溶刻画程度不够.图 5为线性预测结果，该图主要保留宏观光滑同相轴，无法刻画断裂、溶洞等地质单元.应用基于稀疏反演的不连续非均质检测方法，得到的剖面如图 5所示.去除平滑背景后，图 5中小断层表现明显，风化壳破碎描述清楚，更为重要的是能够很好显示顺层岩溶.

图 4

Fig. 4

图 4 三维叠前时间偏移剖面 Fig. 4 Three-dimensional pre-stack time migration profile

图 5

Fig. 5

图 5 平面波预测剖面 Fig. 5 Profile predicted by plane-wave method

为进一步研究断裂、溶洞等地质单元平面展布，分析奥陶系地震属性特征.图 7a为沿层相干属性，相干体技术能够将地震数据转换为相似系数，以突出断层、断裂等非相关信息.相干属性图 7a清楚描述古河道形态，但难以刻画溶洞、微断层等小尺度地质体.与图 7a对比，基于稀疏反演得到的均方根振幅属性图 7b展现出溶洞平面分布情况并清晰刻画断裂特征.

图 6

Fig. 6

图 6 基于稀疏反演方法提取的不连续及非均质信息 Fig. 6 Discontinuous and inhomogeneous geological information extracted by sparsity inversion-based method

图 7

Fig. 7

图 7 奥陶系沿层地震属性 (a) 相干属性；(b) 分离出的不连续非均质数据的均方根振幅属性. Fig. 7 Seismic attributes along the Ordovician layer (a)Coherence attribute；(b)RMS(Root-Mean-Square)amplitude attribute of separated data containing discontinuous and inhomogeneous information.

塔北四连片地震数据分析证实，基于稀疏反演的不连续及非均质地质检测方法，能够反映小尺度断层、断裂等不连续性构造，刻画溶洞等非均质体，在碳酸盐岩油气藏研究中具有重要意义.

通过构建L₁稀疏反演模型，本文给出一种不连续及非均质地质体检测方法，该方法能够分离出小断层、断裂、溶洞等小尺度构造地质单元.在方法研究与应用中，得出以下结论：

(1) 本文在线性平滑信息预测上，采用平面波破坏滤波器，其他方法如F-K滤波、中值滤波等能够达到类似技术效果；

(2) 为提高不连续及非均质地质体分离信噪比，本文构建L₁稀疏模型，并采用分段函数逼近策略和拟牛顿数值方法求解，该算法适用于三维高密度地震数据处理.

(3) 方法应用要求地震数据中不连续及非均质信息成像，如果速度模型不准确，需要进一步速度分析.

(4) 如何将小尺度不连续信息与非均质信息分开，单独应用非均质特征研究“串珠状”储层，可作进一步探讨.

感谢编辑部老师的辛勤工作，感谢审稿专家百忙之中对本文的热心审阅！