渤海湾盆地沾化凹陷罗69井沙三下亚段旋回地层学研究

张若琳，金思丁

(1.成都理工大学能源学院，四川成都，610059；2.油气藏地质及开发工程国家重点实验室，四川成都，610059)

在地学研究领域拓宽的过程中，很多地质现象和过程无法仅从地球系统本身得到有效解释的答案，这促使了旋回地层学即米兰科维奇理论[1]的诞生：由地球轨道驱动力造成沉积地层的旋回性[2-7]。细粒沉积岩同时具有相对不间断、分辨率高、古环境信息丰富且持续时间长的特征，是保存天文信号的理想载体，也是开展旋回地层研究的合适对象。JIN等[8-12]在细粒沉积岩中发现了米兰科维奇旋回。曾有学者认为湖泊沉积环境相对更易受构造作用、陆源输入的影响，不利于天文信号的保存，国际上细粒沉积岩的旋回地层研究多集中于海相沉积中[13-15]。自20世纪90年代以来，在多处湖相沉积中发现了典型的米兰科维奇旋回[16-18]。随非常规陆相细粒沉积岩油气勘探的开展，WU等[19]在松辽盆地晚白垩地层、鄂尔多斯盆地晚三叠延长组地层[20]、渤海湾盆地古近系地层[21-23]等湖相细粒沉积物中发现了米氏旋回。由于天文轨道周期从根本上控制着气候变化，气候变化不仅影响着湖泊的机械沉积作用，而且影响着营养物质供给、古生产力条件、氧化还原条件等[24-25]，利用旋回地层学理论不仅可以尝试以轨道周期控制沉积作用的角度进行解释，从成因上解决地层划分问题，亦能从天文角度还原湖泊沉积的古气候、古环境变化特征，有机质富集因素，对湖相细粒沉积岩的勘探开发具有重要科学意义。本次研究拟对渤海湾盆地济阳坳陷沾化凹陷罗69井的自然伽马测井序列进行旋回地层学研究，运用普通频谱(MTM)、滑动窗口频谱(LAH)和相关系数法(“COCO”)分析技术，匹配最优沉积速率，排除误差信号，提取沙三下亚段(Es3l)细粒沉积岩中的米兰科维奇旋回信息，通过滤波与调谐，建立天文年代标尺，计算不同时期的沉积速率，并将长偏心率周期(0.405 Ma)的滤波曲线作为沙三下亚段内部四级层序划分的参考曲线，定量划分细粒沉积岩高频层序地层。最后以GR曲线为主，并结合不同的能反映沉积环境、古气候的地球化学替代指标，验证天文驱动对有机质富集控制作用的模型。

1 区域地质背景

济阳坳陷位于渤海湾盆地东南部，是渤海湾盆地的1个次级构造单元，也是渤海湾盆地油气最丰富的地区之一[26]。济阳坳陷总面积约为2.65×104km2，内部发育东营、沾化、惠民、车镇4个凹陷及若干凸起，展现出“多凸多凹、凹凸相间”的构造格局(图1(a))[27-28]。始新世沙三段是1段厚度为120～700 m 的深湖相沉积，湖盆中心厚度可达1 000 m 以上，主要岩性为富有机质的细粒沉积岩，含砂岩、灰岩夹层，是古近纪的主要生油岩层段[29]。沙三段从下到上可划分为3段，分别为沙三下亚段(Es3l)、沙三中亚段(Es3m)、沙三上亚段(Es3u)(图1(b))。

图1 渤海湾盆地济阳坳陷构造简图及地层垂向充填序列[23-24,30]Fig.1 Structural sketch and stratigraphic filling histogram of Jiyang Depression,Bohai Bay Basin[23-24,30]

沾化凹陷位于济阳坳陷东北部，属于济阳坳陷的1个次级构造单元，其南部为陈家庄凸起，北部为义和庄凸起，是典型的“北断南超”的箕状盆地，面积约为2 800 km2[30]。从始新世早期开始，沾化凹陷进入断拗阶段，湖盆持续下沉，接受沙河街组沉积，其中沙三下亚段(Es3l)沉积时期处于沾化凹陷的深断陷期，盆地进入非补偿阶段，气候潮湿，处于稳定而持久的深湖，半深湖的沉积环境，发育了厚为100～500 m的暗色泥页岩夹灰岩沉积，为沾化凹陷的主力烃源层。目标井罗69 井位于沾化凹陷的罗家鼻状构造区，沙三下亚段(2 930～3 127.5 m)属1 套深灰色—灰褐色的深湖—半深湖相富有机质纹层状、层状泥页岩、灰质泥页岩沉积，累计厚度为197.5 m[31]。这口井在2 911.00～3 140.75 m井段实施了连续密闭取芯，彭丽等[28]将罗69 井沙三下亚段进一步划分为湖扩体系域(EST)、早期高位体系域(EHST)和晚期高位体系域(LHST)。

2 测试方法与数据预处理

本次研究采用的软件包为基于Matlab 平台下的Acycle2.2[32]。对比Kaleidagraph™[33]，Analyseries 2.0.8[34]，SSA-MTM[35]，基于R 软件的“Astrochron”程序包[36]以及Past3[37]等，Acycle2.2 将旋回地层学研究所要做的所有工作统一编写在同一种操作系统下的同一个程序中，无须来回切换操作系统或软件，操作方便简洁。

本研究中涉及的主要计算方法包括Multi-taper method (MTM)频谱分析、Fast Fourier transform(LAH)滑动窗口频谱分析、滤波分析、相关系数“COCO”及“eCOCO”分析等。罗69 井的GR 测井序列数值范围为20～100 Ω·m，采样间隔为0.1 m。GR曲线的高值与暗色泥页岩对应，低值与灰岩及粉砂岩相对应，肉眼可见明显的旋回特征(图2)[28]。由于GR测井曲线是由各种地质因素所引起的周期性变化的综合反映，其中也包括与地质因素无关的干扰信号，且地层记录会出现随时间或深度逐渐增加或减小的趋势。为了更好地提取GR测井序列中存在的天文旋回频段，在进行频谱分析之前，要将这种低频部分的形变影响去除(见图2中去趋势化后的Detrended GR 曲线)。去趋势化后GR 数据序列中受天文周期影响的信号将更为突出。

图2 罗69井旋回特征Fig.2 Cycle characteristic of Well Luo 69

3 旋回地层学研究

3.1 天文周期解决方案

地球轨道3 个重要参数为偏心率(eccentricity)、斜率(obliquity)及岁差(precession)，受太阳系其他行星及月球与地球间的引力影响，这3个参数在地史时期呈周期性缓慢变化。LASKAR 等[38]引入月球潮汐力的损耗等因素将米兰科维奇理论模型延伸至过去的250 Ma，提出La2004 解决方案，后采用天文历表INPOP08 进行计算[39]，得到更精确的偏心率变化曲线，提出La2010 解决方案。目前的天文周期解决方案对50 Ma以来的地球轨道参数的计算非常精确。

LIU等[22]通过对渤海湾盆地沙河街组四口科学探井的伽马测井数据进行天文旋回计算，得到了22～66 Ma 的高精度天文年代标尺，他们将渤海湾盆地济阳坳陷沙三下亚段的顶和底年龄锚定在40.20～42.47 Ma，该年龄也被ZHAO 等[23]应用在渤海湾盆地东营凹陷牛页1井的天文年代调谐的年龄控制点中，并取得了较好的对应关系。本文拟采用La2010天文周期解决方案，计算40.2～42.47 Ma之间的偏心率、斜率与岁差变化(E+O-P)的理论值，采样间隔为0.001 Ma(图3(a))。随后对E+O-P曲线进行频谱分析与滑动窗口频谱分析(图3(b)和(c))，结果显示40.20～42.47 Ma 期间的主要天文周期为0.405，0.125，0.095，0.051，0.040，0.023，0.022和0.019 Ma。其中，0.405 Ma为长偏心率(E)周期，0.125 Ma与0.095 Ma为短偏心率(e)周期；0.051 Ma与0.040 Ma 属于斜率(O)周期；0.023，0.022 和0.019 Ma 属于岁差周期。利用计算出的各轨道参数的比值，与伽马测井序列通过频谱分析之后的谱峰进行对比，若存在相等或相近的旋回比率，则可以初步判断地层中记录了米兰科维奇旋回信息[40]。

图3 目的层段天文周期解决方案频谱分析图Fig.3 Spectral analysis of target astronomical solution

3.2 深度域频谱分析

对趋势化后的伽马测井数据进行MTM(multitaper method of spectral analysis)频谱分析(图4(b))，图4(b)中，横坐标显示为频率，纵坐标为功率，功率越大表示该功率对应的旋回越明显。每个高置信度峰值(置信度大于95%)对应的厚度旋回分别为39.40，9.74，3.94，2.78，2.54，1.92 和1.70 m，显示出的频谱结构具有能够明显区分开的4 个优势频率带，分别为39.40～27.77，9.74～6.32，3.94～2.54和1.92～1.70 m(图4(b))。分别将4 个频率谱峰值对应的周期比20.44∶5.05∶2.09∶1.00 与理论周期比进行比较，发现天文轨道周期20∶5∶2∶1 的比例关系与之非常接近。

为了识别旋回地层中主频率随深度的变化，观察整体数据序列是否具有稳定性，继而对罗69井沙三下亚段伽马测井曲线进行滑动窗口频谱分析。本次滑动窗口频谱分析中的窗口步长为0.3 m，滑动窗口长度为40 m(图4(c))。谱图中峰值所对应的频率大致与MTM功率谱中高置信度峰值相对应，表明各个深度段均保存了偏心率、斜率和岁差周期的信号。然而，由图4(c)可见在2 988 m及3 040 m处附近出现了明显的峰值信号不连续现象，可能意味着某些深度段序列中包含更多噪声，高频信号被具有更高功率的噪声所抑制[41]，但更有可能的是陆源碎屑物质的输入导致沉积速率发生变化，从而使频谱的周期性发生变化[42]。

据彭丽等[28]对罗69 的层序地层学划分结果，本次研究将GR 数据序列进一步划分为3 小段(2 930～2 988，2 988～3 040和3 040～3 127.5 m)进行频谱分析(图4(d)～(f))。此处只是初步估计沉积速率发生变化的区间段，对罗69 井的沉积速率变化区间精确计算。在2 930～2 988 m 深度段，高于95%置信度的4 个频率谱峰对应的周期为44.42，10.42，4.06 和1.75 m，其比值与405∶95∶40∶22 相近；在2 988～3 040 m 段，高于95%置信度的4 个频率谱峰对应的周期为40.25，12.48，3.89和1.78 m，其比值与405∶125∶40∶19 相近；在3 040～3 127.5 m段，高于95%置信度的4个频率谱峰对应的周期为45.02，10.31，4.52和2.38 m，其比值与405∶95∶40∶19 相近。同时，图4(d)～(f)的频谱结果显示了在不同深度段，由于沉积环境的改变，即使沉积速率发生微小变化，也会导致沉积物中保存的主要天文轨道信号发生变化。

图4 罗69井沙三下亚段GR曲线频谱分析图Fig.4 Spectral analysis of Es3 l in Well Luo 69

3.3 最优沉积速率的估算

由于选用的GR测井数据是以深度为单位，尽管在95%的置信度之上检测到了功率的周期性，但并未对空间周期所代表的时间范围进行估计，因此，仅通过传统的“比值法”对天文旋回进行识别是比较粗略的。本次研究进一步采用LI 等[32]提出的“COCO(correlation coefficient)”和“eCOCO(evolutionary correlation coefficient)”分析来提高旋回识别的准确性，通过自回归模型(autoregressive-1，AR1)，利用蒙特卡洛模拟(Monte Carlo)对非天文轨道驱动的信号(噪声)进行零假设检验，增大天文轨道信号识别的精度。检验中的显著性水平(H0，significance level)表示地层记录中零假设被错误拒绝的概率[43]，例如0.1%的显著性水平表示有99.9%的频谱被检测到，被错误拒绝的零假设(非天文轨道驱动的信号)仅为0.1%。因此H0越小，分析的结果越可靠。“COCO”分析同时还可以计算研究层段的沉积速率，通过“COCO”计算得罗69井沙三下亚段的平均沉积速率为7.2×103cm/Ma 或8.2×103cm/Ma(图5(a))，两者都超过了要求的显著性水平(图5(b))，并且沉积速率为8.2×103cm/Ma，且7 个天文轨道参数分量都参与的置信水平低于0.05%，可以拒绝没有轨道信号参与的零假设(图5(b)～(c))。根据前人对研究区对应层段的沉积速率的估算，本次研究中“eCOCO”分析选取(5～50)×103cm/Ma的速率区间，滑动步长为0.7 m，滑动窗口为35 m，蒙特卡洛模拟2 000 次，得到如图6所示结果。由图6可知：数据全段与图5有良好的对应关系，且(7.2～8.2)×103cm/Ma的沉积速率在沙三下亚段整段保持稳定(图6(a))。因此，“eCOCO”分析可用于对罗69 井进行的旋回地层学分析，其在沙三下亚段保存有完整的米兰科维奇旋回。

图5 罗69井沙三下亚段GR曲线的“COCO”分析Fig.5 "COCO"analysis of detrended GR curve in Es3l of Well Luo 69

图6 罗69井沙三下亚段GR曲线的“eCOCO”分析Fig.6 "eCOCO"analysis of detrended GR curve in Es3l of Well Luo 69

3.4 滤波和调谐

基于频谱分析的结果，将代表长偏心率(E：0.0405 Ma)的沉积旋回(39.40 m)通过高斯带通滤波提取出来，其滤波频率为(0.027±0.003 3)周/m，滤波结果显示罗69 井沙三下亚段约记录了6 个长偏心率旋回(图2)；提取代表短偏心率e(0.100 Ma)的沉积旋回(9.74 m)通过高斯带通滤波，其滤波频率为(0.114±0.035) 周/m，滤波结果显示沙三下亚段约记录了24个短偏心率旋回(图2)；提取代表斜率O(0.040 Ma)的沉积旋回(3.94 m)通过高斯带通滤波，其滤波频率为(0.300±0.020) 周/m，滤波结果显示沙三下亚段约记录了59个斜率旋回(图2)。

由于50 Ma以来的天文模型已经十分完善且可靠，可尝试将滤波曲线与天文目标曲线对应。长偏心率0.405 Ma 在地质历史时期一直显示比较稳定的状态，本次研究将选择代表长偏心率的滤波曲线建立深度-时间模型，并选取长偏心率、短偏心率、斜率3 个天文旋回参数作为理论调谐曲线(图7)。由图7可见：以沙三下亚段底界面年龄42.47 Ma 作为天文年代调谐的初始锚点，沙三下亚段顶界面年龄40.20 Ma 为控制点，用来界定本次研究的时间范围并作为深度与时间转换的依据。当获得时间域数据列后，结合对应的深度域序列可获得罗69 井沙三下亚段连续变化的沉积速率，在2 963.3 m(40.552 4 Ma)，2 993.6 m(40.957 4 Ma)，3 031.6 m(41.362 4 Ma)，3 066.7 m(41.767 4 Ma)这4处出现了沉积速率的明显变化(图7)，沉积速率从下到上存在明显的减小趋势。

图7 罗69井沙三下亚段天文调谐结果Fig.7 Astronomical tuning of GR series of Es3l in Well Luo 69

3.5 调谐后时间域频谱分析

为了检验滤波和调谐结果的准确性，对调谐后得到的时间域GR数据序列再次进行频谱分析及滑动窗口频谱分析(图8)，GR曲线频谱分析图中高于95%置信度的波峰对应的频率将被用于对应轨道参数比值的计算，被识别的频率有2.5，9.4，21.8，24.9，32.3，42.6，47.8 和57.8 周/Ma，其中频率2.5，9.4，24.9，42.6和47.8 周/Ma的倒数为其所对应的天文轨道周期，分别为0.405，0.100，0.040，0.023 和0.020 Ma。以40.2～42.47 Ma 时期的各天文周期范围为参考(图3)，结果表明对罗69井进行的旋回地层学分析是可靠的。对时间域GR曲线进行滑动窗口频谱分析，发现各峰值对应的频率大致与MTM 功率谱(图8(c))中高置信度峰值相对应。

图8 罗69井时间域GR的频谱分析Fig.8 Spectral analysis of Well Luo 69 in time domain

4 天文旋回控制下的湖相细粒沉积

4.1 自然伽马测井曲线的古气候和古环境意义

自然伽马测井(nature gamma-ray)是测量岩层中放射性元素(主要是40K，232Th和238U)衰变后放射出来的伽马射线强度，相当于测量沉积地层中钾、钍、铀总含量的变化。不同种类岩石的放射性元素的含量存在明显差异，例如，泥岩中的黏土矿物和有机质对U的吸附能力较强[44]，而由钙质生物壳体构成的碳酸盐岩，具有较低的放射性元素含量，促使生物有机体发育和生长[45]。伽马曲线不仅能应用于区分基本岩性，目前已被广泛应用于反演古气候和古环境研究[46-47]。在地层中GR 值的升高在一定程度上能够表征泥质含量的增高，沉积环境水动力相对较弱，水体变深，海(湖)扩张等，从而反映相对降水多、湿润的古气候特征[48]。刘泽纯等[49]将碳酸盐岩含量、孢粉和δ18O 等与GR 曲线相结合作为古气候研究的替代指标，GR曲线和不同指标之间显示出相对一致的波动趋势及较好的可比性。

自然伽马曲线在旋回地层学中被频繁使用[50-52]。除去自然伽马对气候、环境变化相对敏感之外，GR数据序列最大的优势在于：油田的探井资料较容易获得，且采样间隔(0.125 m)具有高分辨率。本次研究以GR曲线为主进行旋回地层学分析，并通过分析GR数据与其他气候指标之间的相关性，找出影响气候变化的主要天文周期因素及其控制机理，为反演中国东部始新世气候变化提供参考。

4.2 基于米兰科维奇理论的细粒沉积物层序地层划分

米兰科维奇旋回与层序地层的定量划分有着紧密的联系。石巨业等[53]研究表明，三级层序的持续时间一般为1～10 Ma，与构造因素相关；0.405 Ma的长偏心率旋回已经被证明与四级层序存在直接成因关系，具有高幅水深变化，且彼此有成因联系的水进-水退序列；五级层序受0.100 Ma的短偏心率控制，具有低幅水深特征；六级层序受0.040 Ma的斜率周期控制[54]。

考虑到自然伽马曲线的意义，GR的极小值代表着湖平面当时的最低值，可作为四级层序界面的划分依据(图9)。综上所述，本次研究可将罗69井沙三下亚段地层划分为7个四级层序，其中Es3l1和Es3l7发育不完整。

4.3 沙三下亚段细粒沉积物中有机质富集作用的天文轨道影响

图9中使用的其他古气候替代指标数据见文献[55]，包括罗69 井沙三下亚段2 930～3 127 m 之间的沉积地层所保存的古气候古环境信息记录。以本文所建立的罗69 井的天文年代标尺为基础，将各类替代指标的深度数据对应到相对的地质时间上，得到时间域的数据序列。其中GR 曲线的0.405 Ma 滤波曲线(图9中红色实线)与TOC、黏土矿物质量分数、陆源碎屑指数、碳酸盐岩质量分数(图9蓝色实线为红色实线的镜像)均具有良好的对应关系，与化学风化强度不具有明显对应关系。

图9 时间域GR、TOC、黏土矿物质量分数、陆源碎屑指数、碳酸盐岩质量分数、化学风化强度指标数据序列与0.405 Ma理论轨道周期曲线的相位关系Fig.9 Time domained GR,TOC,clay mineral mass fraction,detrital classic index,carbonate mass fraction,chemical weathering index series and their relationships with astronomical theory period(0.405 Ma)

RACHOLD等[56]总结了轨道强迫气候过程及其对沉积环境的影响，其中长偏心率周期(0.405 Ma)对长期季风周期的调节作用有重要影响[57]，由于岁差旋回的最强和最弱振幅变化对应最大和最小长偏心率，从而形成了2 个截然不同的季风气候[58]。JIN等[8]对古生代页岩地层中有机质富集规律建立2种假设，确定有机质富集作用主要受控于长偏心率轨道周期。由于太阳系的混沌行为，天文学家还无法获得超过250 Ma的精确天文轨道周期方案，因此，对影响古生代细粒沉积岩的季风变化控制沉积过程、沉积环境的探讨仅仅停留在2种模型的假设：第一种假设认为初级生产力是细粒沉积物中有机质富集(高TOC，高GR)的主要因素。偏心率极大值对应岁差幅度变化最剧烈，从而形成了强烈的季风系统，频繁的干湿气候变化及较强的季风降水作用。剧烈的风化作用使向湖盆中输送碎屑物质含有更多的养分，有利于湖盆内生命体的发育，即偏心率极大值控制有机制的富集，对应TOC、GR 的极大值(图10(a))；第二种假设认为在长偏心率周期控制作用下，有机质的保存条件是细粒沉积物中有机质富集(高TOC，高GR)的主要因素。在偏心率极大值时，尽管大量的养分输入促进了季节性的高生产力，但在偏心率极小的干旱条件下，有机质又会被氧化分解，相当于并未真正富集。相反，在偏心率极小值时，季节的交替减少，提供了一个稳定的沉积环境，有机物不受保存条件的影响[59]，可以大量富集(图10(b))。

图10 2种有机质富集因素模型Fig.10 Two factor models of organic accumulation

本次研究获得的GR数据序列的滤波信号与理论轨道周期之间的相位关系在沙三下亚段早晚期呈不同状态(图9)：Es3l1，Es3l2和Es3l3这3个四级层序(Es3l早期)的调谐GR曲线中，0.405 Ma滤波曲线的高值与E+O-P-0.405 Ma滤波曲线的低值相对应；而Es3l4，Es3l5，Es3l6和Es3l7这4 个四级层序(Es3l 晚期)的调谐GR曲线中，0.405 Ma滤波曲线的高值与E+O-P-0.405 Ma 滤波曲线的高值相对应。即提出的模型二在罗69 井沙三下亚段早期的沉积记录更

可信，而模型一更适用于沙三下亚段晚期的沉积记录；沙三下亚段早期(Es3l1，Es3l2，Es3l3)长偏心率极小值控制作用下的季节交替减少，稳定的沉积环境降低了有机质的氧化分解效率，在沙三下亚段晚期(Es3l4，Es3l5，Es3l6，Es3l7)长偏心率极大值控制作用下，强烈的季风作用导致碎屑物质提供更充足的养分，促进了有机生命体的生长。从图9所示的化学风化指数与偏心率滤波曲线并无明显相关性也可推测，模型一中剧烈风化作用为有机质富集提供养分的假设不成立。

5 结论

1)根据GR测井曲线频谱分析结果，渤海湾盆地沾化凹陷罗69 井沙三下亚段地层中沉积旋回周期分别为0.405，0.095，0.040，0.023，0.022 和0.019 Ma。其中长偏心率周期0.405 Ma控制的地层厚度为39.4 m；短偏心率周期控制的地层厚度为9.74 m；地轴斜率周期控制的地层厚度为3.94～2.54 m；岁差周期控制的地层厚度为1.92～1.70 m。罗69井沙三下亚段最优沉积速率为8.2×103cm/Ma。

2)基于米兰科维奇旋回周期的识别结果，对去趋势化的自然伽马测井数据进行高通滤波处理，据0.405 Ma 长偏心率旋回建立“浮动”天文年代标尺，并以沙三下亚段底部磁性地层年代数据42.47 Ma 为锚点，沙三下亚段顶部磁性地层数据40.20 Ma为控制点，建立罗69井的天文年代标尺。

3)以长偏心率0.405 Ma 周期滤波曲线作为四级层序划分的参考曲线，将沙三下亚段划分为7个四级层序，实现湖相细粒沉积岩的高频地层单元划分，并精确计算沙三下亚段沉积速率的变化。

4)从天文学角度验证了天文轨道周期变化和有机质富集因素之间的联系，并建立了2个有机质富集因素的对比模式：长偏心率控制初级生产力为主的有机质富集模式；长偏心率控制氧化还原条件为主的有机质富集模式。本研究所获得的GR数据序列的滤波信号与理论轨道周期的滤波信号比对结果，证实了沙三下亚段早期GR数据的极大值与长偏心率之间的极小值具有良好的对应关系，沙三下亚段晚期GR数据的极大值与长偏心率之间的极大值具有良好的对应关系，表明沙三下亚段中有机质的富集存在保存条件和生产力条件这2个不同的主控因素。