基于坐标变换的三维旋转曲率属性计算方法

张晓琦, 文晓涛, 张超铭, 何易龙, 王锦涛

(成都理工大学 a.地球物理学院,b.油气藏地质及开发工程国家重点实验室，成都 610059)

0 引言

裂缝是重要的油气储集空间与运移通道，因此有效地识别裂缝已成为裂缝型油气藏研究的重点与难点。一般而言，裂缝与构造有密切的关系，在构造高部位及断层附近易发育裂缝，而在这些部位，同相轴的曲率与构造其他部位也有差异，因此曲率属性分析是一种较有效的裂缝识别方法[1-2]。Robert[3]总结了各种曲率属性的定义和计算方法，分析了曲率在三维解释中的应用，发展了一代传统的二维曲面曲率计算方法；Al-Dossary等[4]发展了第二代的体曲率计算方法，该方法直接利用地震振幅计算体曲率，验证了曲率属性对地下裂缝以及裂缝型碳酸盐岩储层有良好的反映能力；Marfurt等[5]将相干属性和曲率属性进行对比融合显示，能较好地突出地质构造特征，取得了很好的效果；陈学华[6]改进了算法，将多尺度自适应微分算子引入体曲率计算中，充分利用不同尺度的曲率来突出地质异常，因此此方法可以用来描述地震在时间和空间上的多尺度特征；刘松鸣[7]利用推导的方位滤波算子对体曲率属性进行方位加强处理，其主要特点是突出了特定方向的地质体异常信息，同时也抑制了地震背景干扰。

一般来说，常规曲率算法在描述平缓地层的几何形态上还是比较可靠的，但如果地层陡峭倾斜甚至于反转，就会有很多算法高估曲率属性，而这常常是因为传统算法对相关的方位角及倾角做出了不准确的估计，进而影响了区域曲率精度。研究者们对此也提出了诸多办法来解决这一方面的问题，如进行大倾角扫描[8]、离心窗扫描[9]，或提出为消除倾角影响采用振幅曲率加权算法[10]等各种优化改进方法。虽然上述方法在裂缝识别中可以取得一定的效果，但总体来看，常规的曲率计算在处理地震数据时均在地球直角坐标系中进行处理，计算方法复杂，且方位角倾角的计算与使用存在较大的误差，不利于裂缝断层的准确描述。为了弥补这些方法的不足，增强曲率属性的提取精度，笔者提出了基于坐标变换的三维旋转坐标变换曲率计算方法，将构造面进行三维旋转坐标变换，借助构造倾角将局部构造曲面转换为平面，间接减少角度影响。这样既简化了计算，减少了误差，又提高了曲率属性的分辨率，即增强了曲率属性的提取精度。

1 基本原理

体曲率属性[11]作为直接描述断层和预测裂缝走向及分布的一个几何属性，在地震勘探中占有重要地位。三维地震数据的采集和显示通常是在直角坐标系中，如图1所示，z轴为深度方向，θ角为真倾角，φ角为从x轴出发的方位角。

图1 倾角及方位角关系图

根据Marfurt等理论，三维空间中的特定某点处的曲率值[12-13]，是利用这一点及其相邻的采样点的视倾角的拟合方程进行计算。然后通过使用最小二乘法近似的获得二次曲面方程：

z(x,y)=ax2+by2+cxy+dx+ey+f

(1)

通过计算上述微分方程，趋势面系数计算如下：

(2)

所以倾角θ和方位角φ分别表示为：

(3)

(4)

这里是在基于地球直角坐标系旋转[14]后的新坐标系中对每个样本的二阶导数，采用复地震道分析方法进行曲率属性计算，在旋转后的新系统中，其中X1轴沿倾斜方位角方向，Z1轴与过原点的表面垂直，x0、y0、z0表示原始坐标系中坐标参数，x1、y1、z1表示坐标系更新后的坐标参数。具体操作步骤为：将Z0轴作为旋转轴，方位角φ为旋转角度，对原始的坐标系统X0-Y0-Z0沿逆时针(右侧)进行旋转，得到新的X1-Y1-Z0系统。 然后将Y1作为旋转轴，倾角θ为旋转角度对X1-Y1-Z0系统沿逆时针(右侧)进行旋转，经过两次旋转得到新的X1-Y1-Z1系统，各步骤如图2所示。按照右手坐标系约定，逆向旋转为正。原始坐标系统X0-Y0-Z0经过图2的两个步骤旋转至新的X1-Y1-Z1坐标系统，沿方位角旋转后得到的X1轴与过倾斜面的Z1轴垂直。在数学中，上述两步旋转可以用两个旋转矩阵相乘的结果来表示：

图2 坐标转换示意图

(5)

其中旋转系数分别为：

(6)

(7)

(8)

注意矩阵的右乘。

相应的坐标逆向旋转为：

(9)

(10)

(11)

(12)

其中:θy和φz分别表示沿y1轴和z0轴逆时针旋转的角度。

(13)

(14)

(15)

中间涉及的相关参数如下：

(16)

(17)

(18)

(19)

(20)

(21)

(22)

(23)

利用一阶公式并借助链式法则可用二次曲面拟合系数表示以上公式中的二阶导数。

根据Di[16]提出的理论，我们将式(13)、式(14)、式(15)代入最大正曲率kpos式(24)和最小负曲率kneg式(25)中。

(24)

(25)

常规曲率计算是基于最小二乘法近似获得二次曲面方程，并计算趋势面系数，带入常规曲率计算方程计算而出。而本次使用的三维坐标变换方法所用的处理流程是在计算完趋势面系数后，先计算出倾角，方位角，再依次以Z0为旋转轴逆时针旋转该点的方位角角度，以Y0为旋转轴逆时针旋转该点倾角角度，这样就可以计算新曲面的一阶导数，将一阶导数为零，至此即为新坐标系原点，完成三维旋转，最后计算新曲面的二阶导数，带入曲率计算方程。在相同的情况下，本文方法计算的曲率属性分辨率更高，效果更好，断层、裂缝带显示更清晰，细节更为丰富。

2 模型方法测试

笔者借助三维French模型来验证本文基于坐标变换的旋转体曲率计算方法，在提高曲率属性分辨率上的有效性(图3)。

图3 三维French模型及目的层反射切片

图4是常规体曲率和本文,基于坐标旋转后体曲率算法所分别得出的最小负曲率沿层切片与最大正曲率沿层切片。图4(a)～图4(b)表示常规体曲率算法的沿层切片计算结果，图4(c)～图4(d)表示本文基于坐标旋转的体曲率算法沿层切片模型计算结果。如图4(a)所示，经过本文方法计算后所获得的曲率结果，将河道显示地更为清晰，检测精度也更高，没出现常规算法中出现的部分河道缺失等问题。对比图4(a)和图4(b)可以发现，经过文中方法计算后，图4(d)中展示的细节更丰富，各细小构造得到了加强，可以更加突出地显示细小构造的特征，对于我们精确识别微小断层及裂缝带很有帮助。

图4 不同算法的曲率属性沿层切片

因为实际地震数据会包含许多干扰信息(噪声)，而噪声对曲率计算是有一定影响的，当噪声过大，还会出现错误结果，因而误导我们对地下真实构造的判断，所以需要检验此方法是否在现实生产中仍能保持相对稳定。关于抗噪性分析(图5、图6)，图5为模型加不同信噪比后的沿层切片图，本实验对模型添加随机噪声，添加的噪声分别为信噪比5.0、10.0、20.0，图6为当模型加不同信噪比后，用本方法处理所得到的最大正曲率切片。对比图6(a)、图6(b)、图6(c)可以发现,模型随噪声的加剧逐渐变得很难直观判断出构造分布，而经过处理后的曲率属性还可以在一定程度上识别出模型真实情况。虽然随着噪声的加大，细节开始逐渐变得模糊不可分辨，曲率效果也随之越来越差，但是走向与大致轮廓依旧清楚可现，因此总体来说，经过信噪比测试，可以肯定本方法相对稳定，具有一定抗噪性。

图5 模型加噪后沿层切片图

图6 加噪后经曲率处理的最大正曲率沿层切片图

3 实际地震资料测试

研究区位于川东北盆地，处在我国西南地区四川盆地的东北部，其北边有秦岭造山带，东部为龙门山冲断带，西部为大巴山褶皱带。研究区以位于四川A工区一定范围的实际地震数据为实验样本进行计算研究。目前已知A强形变地区分布须二段高产富集带，富集最有利区为(海陆相烃源岩+烃源断裂、陆相断裂+砂体发育区)。混源充注，输导体系与裂缝发育带控制“甜点”分布。此处M3井钻遇良好油气显示，海相断裂供烃，陆相断裂控缝，海陆断裂联合控制天然气富集带分布，井中岩石采样可见高角度裂缝(图7)。

图7 M3高角度裂缝

研究区须二段属于裂缝发育带控制着“甜点”的分布，该区域特点是“大面积含气、局部富集高产”。须家河组二段作为晚三叠世出现的第一个稳定层位，其累产气量近90×108m3，探明产量可占须家河组的70%，此段已测试出了多口高产稳定井，获得工业气流，由此可见须二段是整个地区须家河组天然气勘探生产的主要产层段。

研究区裂缝带走向呈北-东向分布，该方向裂缝形成时期较晚，与现在最大应力方向一致，对于裂缝开启是非常有利的。因为该区域须二段为大型岩性气藏，裂缝控制油气的运移与储集，所以对研究区裂缝带精确分析有助于后期的井位部署。结合以上背景原因，对须二段的部分工区地震资料进行基于坐标变换的三维体曲率属性提取，以获得更为精确的曲率结果与裂缝分布显示。图7引自中国石化勘探分公司。

图8(a)～图8(d)分别为研究区域的地层T0图、该区域的沿层振幅切片图、最大正曲率沿层切片图及最小负曲率沿层切片图，从中我们可以看出该区域断层、裂缝带宏观分布，但是很难发现一些小的断层和裂缝带，即常规方法提取的曲率属性对细节的反映不够。

图8 研究区域地层T0图、沿层振幅切片图、最大正曲率沿层切片图及最小负曲率沿层切片图

图9是使用坐标旋转体曲率方法对该区域所做的最小负曲率沿层切片和最大正曲率沿层切片。从图9中可以明显发现，最小负曲率与最大正曲率比常规方法分辨率更高，能够显示更为丰富的细节，断层、裂缝带分布情况更明显，解释结果与实际工区区域地质构造保持一致，与此地区裂缝体展布相同。而且通过测井数据可知本实验地区M3井裂缝发育，测井裂缝走向与处理所得曲率属性切片中的裂缝发育方向一致，曲率属性计算预测结果与成像测井吻合(图10)。说明了本方法所计算的曲率属性较为准确，证明了本文方法预测结果的可靠性。相关的地震剖面与本工区测井图，如图11、图12所示。图11中，原始地震剖面显示为彩色，曲率属性显示为黑白色，从图11(a)～图11(d)对比可知，传统的计算方法在剖面上的显示较为粗略，相对来说不够细致，而本方法计算的曲率属性在剖面上有更多细节体现。

图9 坐标旋转曲率属性方法沿层切片

图10 M3二段成像测井裂缝走向与裂缝倾角图

图11 曲率属性与原始地震数据叠合的地震剖面局部图

图12 M3过井剖面图

4 结论

通过模型试算和实际地震数据的分析可以得到以下结论：

1)本方法将三维坐标旋转应用到曲率属性计算中，借助构造倾角将局部构造曲面转换为平面，使目的层在新的旋转坐标系中处于水平状态，将计算中的一阶导数变为零，既符合工区实际数据的应用，又提高了计算效率，使精准度有所提高。

2)噪声对本方法也有影响，但在模型加噪后，经本方法处理所得的曲率属性切片其宏观构造形态依然清晰可见，具有一定程度的抗噪性，证明方法具有实用性与可靠性。

3)通过模型试算与实际地震数据分析可以看出，相对于传统体曲率计算方法，本方法可提取出更多细节，更能详细地描述地下构造，描述裂缝和断层等的区域分布情况。