张坦,姚威,赵永强,周雨双,黄继文,范昕禹,罗宇
(中国石化 石油勘探开发研究院 无锡石油地质研究所,江苏 无锡 214126)
塔里木盆地石炭系是油气勘探的重点层系之一,多年的勘探实践表明,石炭系具有自生、自储和自盖的能力[1]。巴麦地区自石炭系卡拉沙依组油气发现以来,不断有新的发现。近年来,随着勘探程度的持续深入,人们对石炭系卡拉沙依组层序沉积、储层特征和烃源岩生烃潜力的研究愈加重视[1-3],但地层格架建立的精度不够,限制了研究区内石炭系卡拉沙依组精细划分和对比。同时,沉积体系空间配置和烃源岩有利发育相带主要受控于沉积期古地貌的发育特征,古地貌直接决定了沉积相展布的区域差异性和平面展布特征,而剥蚀厚度的精确计算则成为恢复古地貌的一项重要的基础工作,并且地层剥蚀厚度也对恢复原型盆地、重建盆地沉积-构造演化格局和油气资源定量评价起着非常重要的作用[4]。
近些年来,米兰科维奇(Milankovitch)提出地球轨道参数偏心率、斜率和岁差发生周期性变化,会引起地球表面日照量周期性变化,进而导致地球气候系统周期性变化。这种由地球轨道驱动造成的旋回性记录称为米兰科维奇旋回[5],简称米氏旋回。旋回地层学就是通过研究天文轨道影响的地层,从而有效地识别米兰科维奇旋回,建立具有相对时间概念的“浮动”天文年代标尺和高精度的地层层序格架。随着人们对米兰科维奇旋回理论和应用方法的深入研究[6],米兰科维奇旋回能够合理地解释更多地质学领域问题,使得不同地质历史时期海相、陆相地层的沉积特征能够被揭示出来[7-8],被广泛应用于提高定年精度[9]、构建高频层序地层格架[10-13]和古环境、古气候[14-15]等方面的研究,同时该方法也被应用于恢复地层剥蚀厚度的研究。郭颖等[16]针对海相碳酸盐层的剥蚀特征提出了旋回分析法,精确地恢复了玉北地区东部中、下奥陶统鹰山组剥蚀厚度;赵军等[17]基于米兰科维奇天文旋回理论有效地分析了松辽盆地X油田青山口组剥蚀厚度,解决了盆地模拟中青山口组沉积时期构造演化问题。米兰科维奇理论的深入发展为地层剥蚀厚度的精确计算提供了一种有效的思路和方法,相较于其他方法,该方法可操作性强,适用范围较广,剥蚀厚度计算的精度较高。
本次研究针对塔里木盆地巴麦地区以陆相碎屑岩沉积为主的石炭系卡拉沙依组为研究对象,运用频谱分析法和经验模态分解法(EMD)对石炭系卡拉沙依组进行米兰科维奇旋回分析,并利用连续小波变换法对主频选取的信号进行了合理性分析,最终依据固有模态(imf)分量与长偏心率(e1)之间的对应关系,重构了卡拉沙依组具有相对时间概念的“浮动”天文年代标尺。以405.00 kyr长偏心率(e1)为划分标准,建立了高精度的层序地层格架,并进一步探讨了不同地区钻井缺失旋回数量和平均旋回厚度之间的关系,精确计算了巴麦地区石炭系卡拉沙依组顶部地层剥蚀厚度,为古地貌的精确恢复奠定了坚实的基础。“浮动”天文年代标尺的建立、高精度层序地层格架的搭建及地层剥蚀厚度的精确计算,也为下一步的油气勘探提供了新的资料。
巴麦地区地处于塔里木盆地塘古巴斯坳陷的西缘和西南坳陷的东北部,包括巴楚隆起和麦盖提斜坡两个一级构造单元(图1a)。巴楚隆起主要受到色力布亚-玛扎塔格断裂与阿恰-吐木休克断裂的挟持作用,表现为典型的大型背冲断隆特征。麦盖提斜坡位于巴楚隆起的南部地区,喜马拉雅晚期造山运动使其最后定型为向南倾的单斜坡形态[18]。
图1 巴麦地区构造位置(a)及石炭系卡拉沙依组地层综合柱状图(b)Fig.1 Tectonic location (a) and stratigraphic column (b) of the Carboniferous Kalashayi Formation in the Bamai area
巴麦地区主要发育新生界和古生界,除了三叠系、侏罗系和白垩系普遍缺失外,其余地层较为发育。石炭系是继泥盆系填平补齐的砾岩—粗碎屑岩沉积之后,在塔里木盆地构造运动平静期发育的一套稳定的碳酸盐岩-碎屑岩沉积组合[19],主要包括巴楚组、卡拉沙依组和小海子组(图1b),距今358.90 ~298.90 Ma。卡拉沙依组由下而上可分为上泥岩段、砂泥岩段和含灰岩段,上泥岩段发育泥坪相,砂泥岩段发育三角洲-潟湖相,含灰岩段发育局限台地相,沉积厚度中心主要位于巴麦地区的东部地区,从东到西呈现逐渐减薄的趋势(图2),并可划分为9 ~ 11个三级层序[20],包含有完整的沉积旋回序列。其中,巴麦地区在石炭纪末期发生了一次典型的构造抬升事件——巴楚运动,整体表现为“西强东弱”的特征,造成了研究区内石炭系卡拉沙依组存在明显的剥蚀现象,形成了典型不整合界面[21]。
图2 巴麦地区石炭系卡拉沙依组近东西向地层厚度连井对比剖面(剖面位置见图1a)Fig.2 Near east-west cross-well correlation of the Carboniferous Kalashayi Formation, Bamai area (see Fig.1 for the location)
不同沉积数据序列蕴含着不同的地质信息,对于米兰科维奇旋回信号的响应也有差异。目前应用于米氏旋回周期理论研究的数据序列包括直观岩性数据、地球化学数据、地球物理数据和古生物数据。其中,以自然伽马(GR)测井数据为代表的地球物理数据能够敏锐地反映古环境、古气候等特征的变化趋势,被广泛地应用于旋回地层学分析中。
自然伽马测井数据主要用于测量沉积物中伽马射线的强度,能够敏感地反映沉积物中泥质含量,进而反映古环境和古气候的微观变化[11],高值对应富含黏土沉积物,低值对应砂岩、碳酸盐岩沉积物[22]。本次研究采用巴麦地区4口钻井的相关自然伽马测井数据,其中YB1,YB8,HT1和H3井分别位于巴楚隆起和麦盖提斜坡东部,地层残留厚度相对较大,保留有相对完整的沉积旋回序列,采样间隔为0.125 m,GR数值范围为70 ~ 130 API,并且自然伽马数据纵向分辨率高,能够敏感地反映出岩性组合的旋回性变化。
在开展米氏旋回周期分析之前,首先要对自然伽马测井数据进行数据的预处理,去除超低频和超高频信号,从而尽量地消除地层数据中的各种环境噪声[23]。对于频谱分析,主要采用PAST(Paleonological Statistics)软件来加以实现[24],选择显著水平大于95 %置信度的峰值频率进行分析,对于90 % ~ 95 %置信度之间的数据结果选择性使用。天文年代标尺的建立,首先要搭建地层深度和沉积旋回之间的响应关系。在目前缺少连续地质取样定年数据的情况下,天文年代标尺的建立通常包括连续小波变换法、高斯带通滤波法和经验模态分解(EMD)等方法。近年来,越来越多的研究更倾向于使用经验模态分解(EMD)方法[25],按照从高到低的频率分解模式从复杂信号中提取出若干个相对简单、平稳的固有模态函数,其核心思想是逐步剔除时间序列上、下包络的平均值来获得有限数量的imf分量。本次研究主要采用EMD方法,该方法的优点是能够根据数据自身的时间尺度特征来进行信号分解,无须预先设定任何基函数,相比其他方法,该方法更简单,能够保存更多的地质信息[26]。
太阳系各星球会对地球产生一个叠加的引力场,这种引力场并不是一成不变的,而是发生准周期性变化,而这种变化产生的扰动会导致地球轨道参数周期性改变,从而进一步影响沉积物的产生、运移和累积等过程[27],并记录在地层沉积的信息中。米兰科维奇旋回周期轨道参数通常包括偏心率(e)、斜率(o)和岁差(σ),其中偏心率指轨道偏离正圆的程度,变化频率较低、周期较长。斜率是地球公转轨道面和地球赤道面的夹角,变化的频率比偏心率高,周期相对较短。岁差是最高频率的周期变化,月球和太阳对地球赤道隆起部分的吸引,使地球自转轴的方向环绕与黄道面的垂直轴作缓慢的运动,从而产生岁差。
在米兰科维奇轨道周期中,偏心率的周期比较稳定,共有3个周期长度,分别为405.00 kyr(e1),125.00 kyr(e2)和95.00 kyr(e3),而斜率和岁差周期会随着地质时代的变迁而发生变化。1992年,Berger等计算了时间间隔50.00 Myr,共计11个时间节点的长斜率周期(o1)、短斜率周期(o2)、长岁差周期(σ1)和短岁差周期(σ2)。4个天文轨道周期随着地质时代的变化表现出很好的线性变化规律(图3)。
图3 米兰科维奇轨道周期斜率和岁差随时间变化趋势(据Berger, 1992修改[28])Fig.3 Periodic variation of obliquity and precession of the Milankovitch cycle with time (modified after [28])
根据前人研究成果,塔里木盆地巴麦地区石炭系卡拉沙依组的沉积时间为350.20 ~ 323.20 Ma[20],选取沉积时间中值336.70 Ma作为石炭系卡拉沙依组沉积地质时限带入各天文轨道周期的线性方程,通过计算可得到对应的长斜率周期、短斜率周期、长岁差周期和短岁差周期分别为41.39,33.28,20.34和17.12 kyr。因此,研究区内石炭系卡拉沙依组沉积时期存在一组米兰科维奇旋回轨道周期参数,偏心率周期分别为405.00,125.00和95.00 kyr,斜率周期分别为41.39 kyr和33.28 kyr,岁差周期分别为20.34 kyr和17.12 kyr。
3.2.1 测井数据的主要频率分析
从目的层自然伽马测井曲线数据中提取出沉积旋回的主要频率是旋回分析的关键[29]。通过频谱分析判断沉积地层中是否包含米兰科维奇旋回。频谱分析是一个时域的信号在频域下的表示方式,横坐标表示为频率能量,其是频率的函数,目的就是识别出信号中(准)周期性的成分。频谱曲线中,信号强度越大该频率出现的频率越高,这种频率被称为主频,主频数值越低,代表地层的沉积周期越长。常用的频谱分析方法有快速傅里叶变换(Fast Fourier Transform,FFT)、周期图法(Periodogram)和多窗谱分析法(Multi-Taper Method,MTM)等,不同频谱分析方法的算法不同,各有优缺点,但计算结果基本一致[22]。本次研究使用Past 4.0软件基于多窗谱分析法进行统计分析,选择Spectral analysis模块中的Multi-taper程序对4口井预处理后的自然伽马数据进行主频分析,规定置信度大于95 %的频率作为主频信号(图4),需要注意的是,本次研究中岁差对应的主频明显小于置信度,未能被有效识别出来,其原因可能在于自然伽马测井曲线取样精度,也有可能在于相对短的沉积周期内,米兰科维奇轨道周期并不是沉积作用的主要控制因素。
图4 巴麦地区石炭系卡拉沙依组GR曲线多窗口法频谱分析图谱Fig.4 Multi-taper spectral analysis of GR curves of the Carboniferous Kalashayi Formation, Bamai areaa.YB1井;b.YB8井;c.HT1井;d.H3井
3.2.2 包含米兰科维奇旋回周期主频的获取
为了进一步确定哪些主频信号受米兰科维奇旋回周期的影响,需要建立单井高频旋回与米兰科维奇轨道周期之间的联系。根据前人大量的研究成果[6],如果能在地层记录中找到可以与地球轨道三要素(偏心率、斜率、岁差)之间有相同比率关系的沉积频率或周期,就认为该段地层受米兰科维奇旋回的影响。
通过前文对该地区米兰科维奇轨道旋回周期的确定,研究区内石炭系卡拉沙依组沉积时期的米兰科维奇旋回轨道周期参数为405.00,125.00,95.00,41.39.00,33.28,20.34和17.12 kyr。首先计算出卡拉沙依组沉积时的7个天文轨道周期的频率比率关系(表1),随后计算每口井所有主频信号的比率关系,选取比率关系与天文轨道周期最接近的主频信号作为对应的米兰科维奇旋回周期。
计算结果显示,4口井YB1,YB8,HT1和H3识别出的主频信号都能与米兰科维奇旋回中的长偏心率(e1)、中偏心率(e2)、短偏心率(e3)、长斜率(o1)和短斜率(o2)建立起联系(图3),轨道理论周期比值与多窗口法频谱分析图提取出来的频率比值基本接近(表2),如YB1井,从自然伽马数据多窗口法频谱分析图中可提取出一组频率为0.0195,0.0625,0.0820,0.1914和0.2344 m-1,通过频率和厚度的对应关系可得[28],5个频率对应的旋回厚度分别为6.410,2.000,1.524,0.653和0.533 m,旋回厚度的比值为12.03∶3.75∶2.86∶1.23∶1。从轨道周期比值表中可以发现,长偏心率、中偏心率、短偏心率、长斜率和短斜率对应的理论周期的比值分别为12.17∶3.76∶2.85∶1.24∶1,与频谱分析提取出来的频率比值基本相近。计算结果显示总体比率误差在0 ~ 6.40 %,可接受的平均误差为3.20 %,对应关系较好,表明研究区在石炭系卡拉沙依组沉积时受天文轨道长周期的控制,保存有完整的米兰科维奇旋回,和前人的研究结果一致[20]。
3.2.3 米兰科维奇周期识别结果的合理性分析
连续小波变换法能够将一维时间域测井信号转变为二维时频域,通过调整尺度因子的大小来确定不同尺度天文轨道周期对应的沉积旋回变化趋势,从而获得测井数据直观体现不出来的更多地质信息[30]。如图5所示,每口井中米兰科维奇循环对应的频段覆盖了更多的高能信号,表明米兰科维奇周期对应的主频具有代表性。
图5 巴麦地区石炭系卡拉沙依组GR资料连续小波变换分析Fig.5 Wavelet transform analysis of GR logging data of the Carboniferous Kalashayi Formation in the Bamai areaa.YB1井;b.YB8井;c.HT1井;d.H3井
为了分析单井米兰科维奇旋回周期的合理性,需要验证旋回对应的沉积速率是否能与巴麦地区石炭系卡拉沙依组整体的沉积速率变化范围相互吻合。黄太柱等对塔里木盆地西部海西中期巴楚运动的性质、动力学背景和时空演化进行了探讨[21],认为巴楚运动在研究区内影响范围极大,整体呈现“西强东弱”的特征,东部HT2井和H3井附近遭受剥蚀作用最弱或没有经历剥蚀作用。前人确定的卡拉沙依组沉积时间大约为27.00 Ma,HT1井的地层沉积厚度为392.88 m,H3井地层沉积厚度为397.00 m,据此可以推断巴麦地区HT1井和H3井石炭系卡拉沙依组的整体平均沉积速率分别为1.455 cm/kyr和1.470 cm/kyr。根据多窗口法频谱分析可知,HT1井长偏心率、中偏心率、短偏心率、长斜率和短斜率对应的旋回厚度分别为5.814,1.830,1.333,0.598和0.478 m。H3井长偏心率、中偏心率、短偏心率、长斜率和短斜率对应的旋回厚度分别为5.814,1.827,1.362,0.60和0.474 m。结合米兰科维奇旋回周期,可计算出米兰科维奇轨道周期控制下的地层平均沉积速率分别为1.437 cm/kyr和1.443 cm/kyr,与地层实际沉积速率基本一致,表明主频信息的提取较为准确,能够代表米兰科维奇轨道周期对研究区石炭系卡拉沙依组沉积事件的响应。
如果能够证明地层记录中的沉积旋回是天文轨道驱动所致,便可将沉积旋回或古气候的替代性指标调谐到理论天文轨道参数的目标曲线上,将深度信号校准到时间信号上,就可建立高分辨率的天文年代标尺。本次研究采用经验模态分解(EMD)方法,对研究区内YB1,YB8,HT1和H3井的自然伽马测井曲线进行经验模式分解,将复杂的信号分解为有限个本征模函数(Intrinsic Mode Function, IMF),所分解出来的各imf分量包含了原信号的不同时间尺度的局部特征信号(图6)。
图6 巴麦地区YB8井石炭系卡拉沙依组GR曲线的EMD分析结果Fig.6 EMD analysis results of GR logs of the Carboniferous Kalashayi Formation in Well YB8, Bamai area
利用经验模态分解(EMD)方法得到了研究区内4口钻井自然伽马测井曲线的各7个固有模态分量(imf1—imf7)和一个残余相(res),并根据不同固有模态分量imf曲线中波峰的位置确定包含的旋回个数(表3)。将不同imf分量中统计的旋回个数与识别出的不同级别天文轨道周期对应的平均旋回个数进行对比,发现YB1,YB8,HT1和H3井的长偏心率旋回个数与EMD方法中固有模态分量imf3统计的旋回个数匹配关系较好,因此将固有模态分量imf3的曲线作为确定天文年代标尺的依据,并以此建立了研究区内4口井的“浮动”天文时间尺度和高精度的地层层序格架(周期以405.00 kyr长偏心率周期为最小时间尺度)。需要说明的是,由于在统计固有模态分量imf3旋回个数时,在地层的顶部和底部会出现不足一个旋回的情况,在计算旋回个数时不统计在内。
表3 巴麦地区石炭系卡拉沙依组米兰科维奇旋回个数和imf分量数Table 3 Milankovitch cycle number and the imf component number, Carboniferous Kalashayi Formation, Bamai area
通过固有模态分量imf3曲线确定的天文年代标尺是以405 kyr长偏心率周期为最小时间尺度,即每个沉积旋回经历的地质事件均为405 kyr(图7)。计算结果显示,4口井中,YB8井残留地层的旋回数量最少,每个旋回的厚度为6.410 m,约有35个天文旋回数;YB1井残留的地层厚度对应的单个旋回的厚度为6.410 m,约有45个天文旋回数;HT1井残留的地层厚度对应的单个旋回的厚度为5.814 m,约有66个天文旋回数;H3井残留地层的旋回数量和HT1井一致,每个旋回的厚度为5.814 m,约有67个天文旋回数(表3)。另外,每个旋回的沉积时间对应长偏心率(e1)的时间为405 kyr,可计算出四口井石炭系卡拉沙依组沉积时间为27.13 ~ 13.36 Ma,这也与相关学者通过古生物地层对比得出的约27.00 Ma地层沉积时间基本接近[20],一定程度上验证了本文研究结果的合理性。
图7 巴麦地区YB8井石炭系卡拉沙依组 “浮动”天文年龄标尺Fig.7 Floating astronomical time scale of the Carboniferous Kalashayi Formation in Well YB8, Bamai area
根据前人对巴麦地区石炭系卡拉沙依组展布特征及巴楚运动的时空演化规律研究认为[21,31],石炭系卡拉沙依组沉积厚度中心主体位于巴麦地区的东部地区,并且巴楚运动的构造特征整体表现为西强东弱,认为巴麦地区东部石炭系卡拉沙依组经历了相对较弱的构造隆升作用,该地区地层顶部剥蚀量相对较小或没有经历剥蚀过程,因此这一地区的石炭系卡拉沙依组保留了较为完整的沉积旋回。通过对巴麦地区20口钻井进行经验模态分解(EMD)(表4),研究区内不同钻井基于EMD分解的固有模态分量imf3和e1天文旋回数结果显示,该地区存在地层缺失现象,东部地区残留旋回数量明显大于西部地区,东部地区隆升幅度相对较小,对地层的剥蚀作用相对较弱,也与前人的研究基本一致。
表4 巴麦地区石炭系卡拉沙依组天文周期和剥蚀厚度Table 4 Astronomical periods and denudation thickness of the Carboniferous Kalashayi Formation, Bamai area
结果显示,巴楚隆起东部的H3井和HT1井保留了较为完整的沉积旋回,长偏心率(e1)对应的固有模态分量imf3分别为67个旋回和66个旋回,而麦盖提斜坡东部的井保留的沉积旋回则相对不完整,如YB1井和YB8井对应的固有模态分量imf3分别为48个旋回和35个旋回(图8)。在同一稳定沉积区域内,地层记录的405.00 kyr长偏心率(e1)周期个数应相同,因此笔者认为,巴楚隆起的东部H3井和HT1井附近存在剥蚀最弱的区域,或者没有经历剥蚀作用,地层保存相对完整,可选择H3井和HT1井作为基准井,将记录68个长偏心率(e1)周期的地层定为完整的沉积序列(相对剥蚀厚度最小),而低于这个旋回数量的地层就被定义遭受到剥蚀作用,地层保存相对不完整,因而可通过缺失旋回数量和缺失地层井的平均旋回厚度计算地层剥蚀厚度,其基本原理为:
图8 巴麦地区4口井石炭系卡拉沙依组基于EMD分解的imf3分量数和e1天文旋回数Fig.8 Component of imf3 and number of e1 cycles of the Carboniferous Kalashayi Formation based on EMD analysis of four wells in the Bamai area
式中:h为地层剥蚀厚度,m;N为标准井的旋回数,个,这里为68个;Na为缺失地层井的残留旋回数,个;Ha为缺失地层井的平均旋回厚度,m。
根据上述原理对研究区内20口钻井中卡拉沙依组顶部地层剥蚀厚度进行了精确计算(图9),结果显示:巴麦地区石炭系卡拉沙依组顶部剥蚀厚度在0 ~390 m,剥蚀最强烈的地方位于K2井以北地区,代表了巴楚运动作用最强烈的区域,而HT1井和H3井附近剥蚀作用最弱,剥蚀特征整体表现为西厚东薄。
古地貌的形态控制着沉积相类型和砂体展布范围,对于寻找优质储层起着十分重要的作用[32-37]。剥蚀厚度大的地方,古地貌相对较高;剥蚀厚度小的地方,古地貌相对较低。巴麦地区西部K2井以北地区,剥蚀作用最强烈,代表了古地貌最高的区域,该地区长期处于裸露风化的状态,为沉积体提供物源供给;往东方向剥蚀强度逐渐降低,最薄弱地区位于东部H3井和HT1井附近,代表了古地貌最低的区域,该地区泥质含量较高,砂岩含量较低,有利于烃源岩优质相带的发育,而不利于储集体的形成。中部BT5井附近斜坡地带介于古地貌最高和最低区域之间,该地貌形态具备发育较好的水动力条件,易形成良好的储集体,为油气富集提供了良好的储集空间,是下一步有利的勘探区域。
1) 自然伽马测井曲线频谱分析和连续小波变换分析结果表明,巴麦地区石炭系卡拉沙依组沉积时受天文轨道周期的控制,保存有完整的米兰科维奇旋回,可检测出405.00 kyr长偏心率(e1)、125.00 kyr中偏心率(e2)、95.00 kyr短偏心率(e3)、41.39 kyr长斜率(o1)和33.28 kyr短斜率(o2)5个天文周期,米氏旋回的自然伽马测井曲线的识别,为该地区建立天文年代标尺和高精度等时旋回地层的划分和对比提供了依据。
2) 基于经验模态分解(EMD)方法计算得出的固有模态分量imf3与405.00 kyr长偏心率(e1)控制下的地层旋回个数基本一致,以固有模态imf3分量曲线作为确定天文年代标尺的依据,建立了巴麦地区石炭系卡拉沙依组具有相对时间概念的“浮动”天文年代标尺和高精度的地层层序格架。
3) 基于经验模态分解(EMD)方法计算结果,选择HT1井和H3井作为研究区的基准井,将记录68个405.00 kyr长偏心率(e1)周期的地层定义为完整的沉积序列(相对剥蚀厚度最小),建立了不同地区钻井缺失旋回数量和平均旋回厚度之间的关系,并在此基础上进行了石炭系卡拉沙依组剥蚀厚度的精确计算。
4) 研究区内石炭系卡拉沙依组顶部剥蚀厚度在0 ~ 390 m,剥蚀最强烈的地方位于K2井以北地区,古地貌形态最高,而HT1井和H3井附近剥蚀作用最弱,古地貌形态最低,整体表现为西高东低的特征,中部BT5井附近斜坡区域,易形成良好的储集体,是下一步有利的勘探区域。