基于测井频谱分析的松科二井登娄库组地层沉积速率研究

张浩东 , 邹长春 *, 彭 诚 , 杨玉卿

1)中国地质大学(北京)地下信息探测技术与仪器教育部重点实验室, 北京 100083;2)中国地质大学(北京)地球物理与信息技术学院, 北京 100083;3)中海油田服务股份有限公司油田技术事业部, 河北廊坊 065201

沉积记录是揭示白垩纪地层结构与气候变化的基础, 目前利用海相沉积记录开展研究的成果丰富, 例如碳与氧同位素作为古气候替代指标揭示中新世气候规律(Holbourn et al., 2005), 但陆相沉积记录较少。而松辽盆地科学钻探工程松科二井钻孔钻穿了整个白垩纪地层, 获取了连续、完整的白垩纪陆相沉积记录, 为本次研究提供了丰富的数据。

沉积速率的变化可作为地质活动的间接证据,是研究沉积环境演化的重要参数(郑绵平等, 2013;丁莹莹等, 2018)。沉积速率的计算普遍采用古生物、古地磁和同位素测年资料(刘广山等, 2008; 王华等,2008; 张会领和蒲晓强, 2011), 但这些资料都存在成本高、鉴定周期长、纵向不连续等缺点(张占松等,1999)。近些年旋回地层学研究的不断发展, 认识到沉积地层中普遍记录着具有天文周期信号的地层旋回, 这为我们提供一种新的思路来计算沉积速率(吴怀春等, 2011; 邓胜徽等, 2015)。测井曲线能完整地记录地层中岩性、地球物理、地球化学等参数的旋回性信号, 具有连续、纵向分辨率高、成本低的特点, 为基于天文旋回方法进行沉积速率的研究提供了重要数据(Peng et al., 2020)。前人研究表明湖泊、海洋和河流相环境中的地层可能记录了米兰科维奇旋回信息(Zhang et al., 2019)。松辽盆地登娄库组主要发育河流、湖泊与三角洲相(Li et al., 2015),登娄库组四段孢粉化石与氧同位素证据指示该时期发生了降温事件(黄清华等, 1999; Wang et al.,2013)。但是目前登娄库组地层沉积速率的研究较少。针对早白垩世登娄库组地层, 松科二井采集了丰富的测井资料, 但仅获得了少量岩心。因此, 对测井资料的分析尤为重要。

本文采用松科二井自然伽马能谱测井数据进行频谱分析, 识别地层存在的米兰科维奇旋回信号,优选对天文旋回较为敏感的测井曲线并计算出地层沉积速率。希望从沉积速率的角度寻找松辽盆地白垩纪时期发生地质活动事件的证据。

1 地质背景

松辽盆地位于中国东北辽宁、吉林与黑龙江境内(图 1), 处于蒙古—华北板块东北部边缘带上(王璞珺等, 2015)。登娄库时期由于深部断陷活动逐渐减弱以及上地幔热流衰减, 盆地处于断坳转化的构造格局(刘招君等, 1992)。松辽盆地东南隆起区登娄库组的物源由四周向中心汇聚, 从盆地边缘到中心依次发育冲积体系、三角洲体系、滨浅湖体系、深湖体系(李占东等, 2014)。

图1 井位分布及松辽盆地一级构造单元图(改自张淑霞等, 2018)Fig.1 Distribution of well location and first-order tectonic units in the Songliao Basin(modified from ZHANG et al., 2018)

松科二井位于松辽盆地徐家围子断陷带(邹长春等, 2016), 该井钻穿整个白垩纪地层, 深度达到7018 m。同时松辽盆地科学钻探工程使用先进的测井设备获取了丰富的地球物理参数(邹长春等,2018), 其中登娄库组地层的深度段为 2533—2966 m, 地层缺失登一段, 从下往上分为登二段至登四段, 登娄库组地层与营城组地层呈不整合接触,岩性主要为细砂岩、粉砂岩、泥岩(王璞珺等, 2017)。松科二井登娄库组发育辫状河相、河流相、三角洲相与湖泊相共四种沉积相, 并发育水平、槽状交错、板状交错层理, 在登娄库组二段通过同位素测年法确定岩石年龄为(102±4) Ma(刘硕等, 2019)。

2 数据与方法

2.1 测井数据

自然伽马能谱测井记录了放射性元素U、Th、K衰变的强度, 放射性元素一般容易被黏土矿物吸附。在沉积岩中一般富集U、Th、K。Th常见于含黏土、长石、云母的沉积物中, Th的化学性质稳定,Th/K通常是作为反映古气候变化的替代指标(Schnyder et al., 2006)。而U受次生还原作用与氧化还原条件的影响, 不适合用来识别天文旋回信号。故本文选取松科二井登娄库组对天文旋回较为敏感的GR、Th、K、Th/K测井数据开展沉积速率研究。

2.2 沉积速率重建方法

本文利用松科二井连续、多参数的自然伽马能谱测井数据, 针对登娄库组地层进行频谱分析, 建立地层沉积速率剖面。步骤与方法在Acycle v2.2软件上实现(Li et al., 2019)。具体流程如下:

(1)数据预处理。测井曲线受地层压力的影响,随深度增加曲线整体呈现增加或降低的趋势。去趋势化处理是频谱分析前的关键步骤, 可通过LOWESS滤波器实现, 设置窗长参数为35%总数据长度。

(2)频谱分析。若频谱图中识别的谱峰周期比与地球轨道参数的长、短偏心率、斜率和岁差之比(约20:5:2:1)近似, 则认为地层中存在米兰科维奇旋回信息。本文采用 Multitaper Method of Spectral Analysis(MTM)与Evolutive Harmonic Analysis(EHA)这两种频谱分析方法来识别天文旋回信息。MTM方法使用多个正交锥形窗的方式, 获得整个地层深度内的平均谱结构, 具有估计方差较低的特点(Thomson, 1982); EHA方法则是采取滑动窗口的方式进行傅里叶分析, 最终形成二维频率谱。该方法常用来识别旋回地层中主频率随深度的变化信息,有利于识别可能存在的沉积间断事件(Molinie and Ogg, 1990; Meyers et al., 2001)。

(3)计算沉积速率。核心思路是通过频谱分析识别出米兰科维奇旋回厚度信息, 再结合理论天文周期中的时间信息, 获取地层沉积速率。本文选用了相关系数法(COCO)与进化相关系数法(eCOCO)来计算登娄库组地层沉积速率。COCO方法可以评估整个地层的沉积速率; eCOCO方法则是以滑动窗口的方式来计算不同深度下沉积速率(Li et al., 2018)。这些方法通过使用 Pearson乘积矩相关系数来反映实际数据记录下的功率谱与给定沉积速率下的理论天文周期功率谱之间的相关性, 同时用蒙特卡罗方法来检验没有天文频率信号的零假设, 使用102 Ma下的Laskar04年目标天文参数(Laskar et al., 2004),来测试登娄库组地层的沉积速率。相关系数越高、零假设检验显著性水平越小、且使用天文周期参数数量多的位置代表最优地层沉积速率。Pearson乘积矩相关系数的公式如下:

X为实际数据的功率谱;与sX分别为实际数据的均值与标准差;Y为给定沉积速率下理论天文周期的功率谱;与sY分别为理论天文周期序列的均值与标准差;

(4)滤波。采取高斯带通滤波器提取出 405 ka长偏心率周期信号(Hinnov and Hilgen, 2012), 来显示信号在深度域的变化情况。

3 结果

3.1 测井响应与频谱特征

由于测井曲线的分辨率较高, 能够识别出分米级尺度的岩性。因此, 测井资料可为识别地层旋回性提供数据基础。松科二井登娄库组地层岩性以砂泥岩为主, 结合录井资料, 同时依据不同岩性在测井曲线上有不同的响应值, 可细分为泥岩、粉砂质泥岩、泥质粉砂岩、粉砂岩、泥质细砂岩、细砂岩、泥质粗砂岩与粗砂岩(表1)。不同岩性GR、Th、K测井曲线响应差异性明显, 同一种岩性整体上具有一致性(图2)。GR、Th、K测井响应值随着黏土矿物含量的增加而增加, 泥岩的GR、Th、K测井数据响应高值, 砂岩的GR、Th、K测井数据响应低值。这些测井曲线能较好地反映岩性变化, 可用来识别地层的沉积旋回。

图2 松科二井登娄库组自然伽马能谱测井曲线图Fig.2 Natural gamma-ray spectral logs of Denglouku Formation from CCSD SK-2

表1 登娄库组地层不同岩性的测井响应值Table 1 Log response values of different lithologies in the Denglouku Formation

对 GR、Th、K、Th/K测井数据采用 MTM与EHA方法进行频谱分析。图中频率的倒数表示为地层识别的旋回厚度, 根据天文循环比率法, 识别的旋回厚度比满足约 20:5:2:1, 与登娄库时期的主要的长偏心率周期(405 ka)、短偏心率周期(95 ka)、斜率周期(37.6 ka)、岁差周期(22.5 ka)(Laskar et al.,2004)之比相似, 表明登娄库组地层记录着米兰科维奇旋回信号(图3)。其中EHA分析结果显示上部谱峰发生明显右移, 表明地层沉积速率发生变化。在登娄库组上下层段沉积速率变化大的情况下, 可采用MTM方法在2643 m处分上下段进行频谱分析。

对原始Th数据进行去长趋势化处理后(去长趋势化Th数据)的EHA频谱图显示在约2643 m深度以上长偏心率谱峰发生明显右移。上部地层 MTM频谱图上识别出23.9 m与7.1 m的旋回厚度分别对应长偏心率(405 ka), 短偏心率(125 ka), 并且置信水平在 99%以上, 而斜率、岁差信号较弱未能明显识别; 下部地层 MTM 功率谱识别出 57.7 m、15.8 m、5.3 m、3.4 m旋回厚度, 分别对应长偏心率(405 ka)、短偏心率(95 ka)、斜率(37.6 ka)、岁差(22.5 ka) (图 3a)。

去长趋势化GR数据的EHA频谱图显示上部地层长偏心率谱峰发生明显右移。上部MTM频谱图上识别19.2 m、8.1 m的旋回厚度被解释为长偏心率(405 ka)、短偏心率周期(125 ka); 下部MTM频谱图上识别57.8 m、21.5 m、13.9 m、5.0 m、3.7 m旋回厚度被解释为长偏心率(405 ka)、短偏心率(125 ka)、短偏心率(95 ka)、斜率(37.6 ka)、岁差(22.5 ka)(图 3b)。

去长趋势化K数据的EHA频谱图表明偏心率信号能量发生变化, 在2850 m以下, 长偏心率信号能量变弱。上部 MTM 功率谱图识别出 19.1 m、7.5 m旋回厚度分别对应长偏心率(405 ka)、短偏心率(125 ka); 下部 MTM 频谱图识别出 55.7 m、22.1 m、14.0 m、5.0 m、3.7 m旋回厚度分别对应长偏心率(405 ka)、短偏心率(125 ka)、短偏心率(95 ka)、斜率(37.6 ka)、岁差(22.5 ka)(图 3c)。

去长趋势化Th/K数据的EHA频谱图显示频谱形态发生偏移, 在2850 m以下, 斜率、岁差信号能量较弱。上部MTM功率谱图识别出23.0 m、8.7 m旋回厚度被解释为长偏心率(405 ka)、短偏心率周期(125 ka); 下部MTM频谱图识别出57.7 m、15.8 m、8.3 m、5.7 m、3.4 m旋回厚度被解释为长偏心率(405 ka)、短偏心率(125 ka)、短偏心率(95 ka)、斜率(37.6 ka)、岁差(22.5 ka)(图 3d)。

图3 登娄库组地层测井数据的频谱分析结果Fig.3 Spectrum analysis results of logging data in Denglouku Formation

3.2 沉积速率特征

Th数据的COCO测试结果显示在3.8 cm/ka、5.9 cm/ka、14.5 cm/ka、15.6 cm/ka处存在多个高相关系数的峰, 同时这些谱峰位置的零假设检验显著性水平都小于 0.1%(天文信号存在的概率大于99.9%), 受7个天文周期参数约束(图4a); GR数据测试结果显示在4.0 cm/ka、5.8 cm/ka、14.6 cm/ka、17.0 cm/ka处存在高相关系数的峰, 同时满足零假设检验显著性水平都小于1%, 受7个天文周期参数约束(图 4b); K数据测试结果表明在 4.0 cm/ka、5.8 cm/ka、16.9 cm/ka处存在高相关系数的峰, 且满足零假设检验显著性水平都小于1%, 受7个天文周期参数约束(图 4c); Th/K数据测试结果显示在4.0 cm/ka、5.6 cm/ka、16.5 cm/ka处存在高相关系数的峰, 同时满足零假设检验显著性水平都小于1%, 受7个天文周期参数约束(图4d)。这些测试结果表明登娄库组地层沉积速率大致在4～17 cm/ka范围内变化。

图4 登娄库组地层测井数据的COCO分析结果Fig.4 COCO analysis results of logging data of Denglouku Formation

对登娄库组Th测井数据进行eCOCO分析, 可以反映出地层沉积速率随深度的变化情况(图5)。提取图中满足高相关系数、零假设检验显著性水平低于1%且使用7个天文周期参数的位置为最优地层沉积速率, 结果表明登娄库组地层的沉积速率从下往上整体呈现为由高到低的趋势。该方法的测试结果与COCO分析整个登娄库组地层内存在5.9 cm/ka、14.5 cm/ka、15.6 cm/ka的沉积速率基本符合。

图5 登娄库组去长趋势化Th数据的eCOCO分析结果Fig.5 eCOCO analysis results of detrended Th data in the Denglouku Formation

在登四段地层, EHA分析结果显示出长偏心率发生明显右移, 且 eCOCO分析方法计算出沉积速率为低值, 故分两段并采用不同的高斯带通滤波器来提取地层的长偏心率信号。结果显示登娄库组共约有 10个长偏心率周期(图 6), 从而估算登娄库组沉积持续时间约为 4050 ka, 同时结合登娄库组地层厚度为 433 m, 可大致计算登娄库组平均沉积速率约为 10.7 cm/ka。在 2533—2643 m 内识别4.5个长偏心率周期, 大致计算该段地层平均沉积速率为6.0 cm/ka; 在2643—2966 m内识别5.5个长偏心率周期, 大致计算该层段平均沉积速率为14.5 cm/ka。这与 COCO 方法识别地层存在5.9 cm/ka、14.5 cm/ka沉积速率峰的结果一致。

图6 登娄库组地层长偏心率信号提取结果Fig.6 Extraction results for the long eccentricity signal in Denglouku Formation

4 讨论

4.1 自然伽马能谱测井资料中的米兰科维奇周期识别

地球轨道参数周期性变化改变地球表面接受的日照量, 导致气候的周期性变化, 最终在地层记录中表现出旋回性。利用自然伽马能谱测井数据进行频谱分析能识别出米兰科维奇周期的信号, 频谱图上识别出存在类似米兰科维奇周期比(约 20:5:2:1)的旋回厚度的信息, 表明登娄库组地层沉积过程受到天文轨道驱动力的影响, 测井数据记录着米兰科维奇周期的信息。GR、Th、K、Th/K数据的MTM频谱图上识别的天文旋回厚度信息基本相近, Th数据的 EHA分析结果显示, 在登娄库组四段深度约2643 m以上谱峰整体发生右移, 但仍然满足约20:5:2:1, 代表旋回厚度变小, 反映出沉积速率变小。登娄库组地层沉积过程受天文轨道驱动力影响, 地层记录中存在地球轨道参数的信息, 这些具有天文周期的地层旋回, 为计算地层沉积速率提供了机会。

4.2 替代指标选择对沉积速率重建精度的影响

从米兰科维奇旋回的角度上, 松科二井登娄库组地层 GR、Th、K、Th/K测井数据的频谱分析结果存在异同性(图3)。相同之处, 多种数据的频谱结果显示都存在类似于米兰科维奇周期比例的沉积旋回厚度信息, 表明登娄库组地层记录着米兰科维奇旋回周期的信号。而不同之处, 第一, 不同测井数据 MTM 分析结果中谱峰形态存在差异, 其中

GR数据测试结果中的谱峰数量最多, 反映出记录的信息多, 但同时容易导致天文信号受噪音的干扰,而K数据测试结果包含的谱峰数量最少, 反映出记录的天文旋回信息较少。第二, 在 2850 m以下,K数据测试的长偏心率信号能量较弱, Th/K测试的长偏心率信号能量较强但斜率与岁差信号弱, 表明记录的天文旋回信号不明显, 而GR与Th数据在该深度下都能较为完整地识别反映米兰科维奇旋回周期信息的谱峰。Th数据测试效果能较好的突出记录长、短偏心率、斜率和岁差的谱峰, 尤其长偏心率信号较为明显。结合上述分析, 登娄库组Th数据对天文旋回信息比较敏感。

从矿物元素的角度上, 由于 GR记录了包含放射性元素U、Th、K衰变的强度, 岩石中的U很容易受到次生还原作用与有机质吸附的影响。K由于成岩作用会导致黏土矿物的转化, 导致在不同矿物中相差极大。这些因素在频谱分析中容易产生噪音的干扰。而Th的化学性质稳定, 在风化过程中被认为是不易溶于水并富集在矿物中(Parkinson, 1996)。这可能是登娄库组 Th相比于其它参数对天文旋回信息更敏感的原因。

综上所述, Th测井数据对天文旋回信息较为敏感。在整个登娄库组深度范围内, Th测井数据与GR、K、Th/K测井数据相比记录了较为完整的天文旋回信息。若选择其他参数, 沉积速率的计算结果可能出现较大误差。

4.3 构造活动对登娄库组沉积速率的影响

登娄库组二段至三段地层具有高沉积速率的特征, 最大达到 16.2 cm/ka, 登娄库组四段整体地层沉积速率约为5.9 cm/ka, 具有低沉积速率的特征(图7)。前人研究表明, 松辽地区坳陷湖盆的平均沉积速率为5～9.8 cm/ka, 断陷盆地的平均沉积速率为12～15 cm/ka(陈建文和王德发, 1996)。本文测试沉积速率的结果与其相似, 验证了沉积速率计算结果的可靠性。

图7 登娄库组沉积速率分析结果Fig.7 Sedimentation rate analysis results of Denglouku Formation

构造作用对地层沉积的影响表现为具有时间跨度大的低频旋回特点, 一般周期在百万年以上。构造活动控制着沉积体系与展布特征(Yao et al.,2021), 断层下陷为沉积物的供应创造了空间, 由于坡度高, 为沉积物的侵蚀、运移与沉积提供了便利,并在沉积中心区域聚集多个沉积物供应区域的沉积物(Hawie et al., 2017; Matenco and Haq, 2020)。构造作用可能是影响登娄库组沉积速率剧烈变化的主控因素。根据松辽盆地区域构造背景, 登娄库组以下依次为营城组与火石岭组地层, 火石岭组发育火山岩相(胡丁玉等, 2019), 且该时期断裂活动频繁。随后由于营城组后期火山运动事件, 徐家围子地区的地震剖面上显示出断层下陷(Feng et al., 2010)。同时松科二井营城组地层内发现U富集现象也被推测是受火山活动与断层运动影响(张淑霞等, 2018)。登二段至登三段时期地层沉积主要受断陷活动控制, 在登娄库组2698 m、2754 m、2850 m与2931 m处识别出沉积速率突然增大, 推测是登娄库组早期受偶发性断层下陷事件影响, 沉积可容空间变大, 沉积物运移、沉积到盆地的速率加快, 导致沉积速率的上升。在约 2820 m处沉积速率发生突变推测是受地层抬升影响导致沉积速率下降。随后到登三段末期, 盆地的断陷活动开始减弱, 向区域坳陷的构造格局转化, 导致整体沉积速率变小。

5 结论

利用松科二井登娄库组GR、Th、K、Th/K测井资料, 采用多种时间序列方法(MTM、EHA、COCO、eCOCO)进行频谱分析, 优选出对天文旋回敏感性强的测井数据, 最终建立了登娄库组地层沉积速率剖面, 获得了以下结论。

(1)登娄库组自然伽马能谱测井数据记录了米兰科维奇周期信号。多种测井数据频谱分析结果揭示了地层中记录的长偏心率、短偏心率、斜率与岁差信号。

(2)Th测井数据对天文旋回信息敏感。与GR、K、Th/K数据频谱分析结果相比, Th测井数据在整个地层中识别的天文周期较为完整。综合考虑多种测井数据的测试结果, 建立了登娄库组地层连续、相对准确的沉积速率剖面。

(3)登二段至登三段地层具有高沉积速率的特征, 登四段地层整体呈现低沉积速率特征。推测构造活动是影响登娄库组沉积速率剧烈变化的主控因素。依据区域构造背景, 登二段至登三段时期地层沉积主要受盆地断陷活动控制, 登四段时期盆地的断陷活动逐渐减弱, 地层沉积速率从下往上呈现出由高到低的趋势。

致谢: 本文工作得到了松辽盆地资源与环境深部钻探工程首席科学家王成善院士、总地质师王璞珺教授以及副总指挥冉恒谦教授的大力支持。中国石油集团测井有限公司天津分公司测井队成员参加了测井数据采集工作, 在此一并表示衷心的感谢。

Acknowledgements:

This study was supported by China Geological Survey (No.1212011301600), National Natural Science Foundation of China (No.41790455-1), and Central Public-interest Scientific Institution Basal Research Fund (No.2652019003).