全球海平面变化是层序地层学的重要理论基础和核心内容。高频海平面变化是指四、五和六级与米兰科维奇天文周期有关的海平面振荡变化(Haq et al., 1987; Vail et al., 1991)。研究表明, 短时序(20~400 ka)的米兰科维奇周期变化引起了陆地冰盖的消融和扩张, 控制了高频海平面升降旋回, 米兰科维奇旋回周期(M)反映了地球轨道参数的变化, 它包括地球轨道的长和短偏心率周期(400 ka和100 ka)、地轴斜率周期(54 ka和41 ka)和岁差周期(23 ka和19 ka)。据资料, 在冰室期由米兰科维奇旋回驱动的高频海平面升降幅度波动最高可达150m, 温室期要低, 一般在25m左右(Berger et al., 1990; Pomar, 1991; Abreu and Anderson, 1998; Pillans et al., 1998; Miller et al., 2003), 并且在地质历史时期的地层记录中保存了普遍的印迹。目前, 在古代浅海碳酸盐岩地层中不但已经识别出了四级和五级(400 ka和100 ka)的长短偏心率周期, 而且在深海泥灰岩— 泥岩韵律沉积中也发现了六级(20 ka)的岁差周期(Einsele and Seilacher, 1991; Schwarzacher, 1993; Weedon, 2003)。

层序地层学关注的是冰川或构造驱动的长周期三级海平面变化, 其周期介于0.5~5 Ma(王鸿祯和史晓颖, 1998)。一般认为三级海平面变化具有幕式旋回的特点(episodic cycles, Haq et al., 1987; Vail et al., 1991), 但也有人提出它反映了米兰科维奇长周期变化, 具有周期性波动的特点(periodic oscillation, Schwarzacher, 1993)。

这些低频长周期和高频短周期海平面波动都有其特征的振幅和频率, 它们相互叠加产生全球海平面波动。全球海平面变化和构造沉降速率、沉积供给速率三者之间相互作用并在时间和空间上形成一系列复杂的地层叠置样式。

为了直观、定量地反映全球海平面振荡变化周期、振幅变化与区域沉降速率以及沉积速率之间的关系, 下面采用数值模拟图解的方法予以说明。

按照Weedon(2003)提出的旋回地层学中时间序列分析理论, 任意一个海平面波动过程可以采用正弦波数学模型来描述, 每一个正弦波包括频率、周期和相位3个要素, 而全球海平面的变化轨迹则可以采用这些频率和振幅不同的正弦函数波的复合叠加曲线表示。图1中长周期海平面变化周期设定为1 Ma, 振荡幅度设定为± 20m, 高频海平面变化周期分别设定为100 ka和20 ka, 对应米兰科维奇旋回的偏心率和岁差周期, 振幅设定为± 5m, 则全球海平面或复合海平面随时间升降的变化曲线可通过下式计算:

复合海平面变化曲线=20 m× sin(1/1000 ka)+5 m× sin(1/100 ka)+5 m× sin(1/20 ka)

图1

Fig.1

	Figure Option ViewDownloadNew Window
	图1 高频和低频复合海平面变化曲线Fig.1 Curves of composite sea level change of low and high frequency

图1中横坐标表示时间, 以ka为单位, 纵坐标为海平面升降幅度, 正值指示海平面上升, 负值指示海平面下降。图中蓝色曲线表示频率 为1/1000、 1/100 和1/20的3 个正弦波叠加产生的复合海平面变化曲线, 红色曲线表示长周期三级海平面变化轨迹, 显示在时间长度为2 Ma的区间内, 复合叠加海平面变化化曲线上存在2个长周期、频率 为1/1000的 低频正弦波动, 升降幅度为 ± 30 m。 按照Haq等(1987)和Vail等(1991)的层序地层学理论, 复合海平面变化具有全球一致性, 是进行地层区域对比的基础; 如果海平面波动的周期不变, 则海平面振幅愈大, 海平面变化速率愈高; 而低频海平面变化曲线上的拐点, 即海平面变化速率最大值或最小值的转换点, 对应三级层序界面或层序内部单元分界线的位置。

图2

Fig.2

	Figure Option ViewDownloadNew Window
	图2 构造沉降速率与相对海平面变化速率之间的关系Fig.2 Relationship between rates of tectonic subsidence and relative sea level change

构造沉降(T)是地壳热沉降和构造载荷引起的盆地基底下沉量两者之和。当构造沉降速率增加时, 盆地中同一地点水深增加, 反之则水深减小。据资料, 被动大陆边缘盆地沉降速率低, 一般为0.01~0.025 m/ka, 而以热沉降为主的克拉通盆地到以挠曲沉降为主的前陆盆地, 沉降速率变化为0.01~0.4 m/ka(Angevine et al., 1990)。

全球海平面变化与构造沉降之和决定了相对海平面升降幅度, 进而控制了沉积物可容纳空间的大小。如果全球海平面变化速率保持不变, 区域构造沉降速率愈大, 则沉积盆地的可容纳空间增加越快, 两者之间关系如图2所示。一般来说, 在同一个沉积盆地中, 其盆地边缘的沉降速率小于盆地中心位置的沉降速率, 而且不同沉积时段构造沉降速率还可能有所差别, 因此将会产生一系列复杂的相对海平面变化模式。在三级海平面变化的时间尺度上, 由于海平面变化幅度远大于构造沉降量, 一般假定构造沉降速率为一个常数, 即沉降速率随时间呈线性变化(Read et al., 1986)。图2中分别设定沉降速率为0.025 m/ka和0.055m/ka, 全球海平面变化参数与图1相同, 从图中全球海平面变化与沉降速率叠加产生的2条相对海平面变化曲线可以看出, 随着沉降速率增加, 相对海平面升降幅度和沉积物可容纳空间明显增加。

沉积速率是指海洋或湖泊沉积物在一定时间内在单位面积上堆积的地层厚度。根据McNeill(2005)对第四纪沉积速率的研究, 碳酸盐岩台地的平均沉积速率一般为0.01~0.05 m/ka, 但变幅极大, 最大可达2 m/ka。

海平面变化速率、沉降速率和沉积速率3者之间存在以下函数关系:

相对海平面变化速率(RS)=全球海平面变化速率(ES)+沉降速率(T)

可容纳空间变化速率(A)=相对海平面变化速率(RS)

全球海平面变化速率(ES)+沉降速率(T)=水深变化速率(W)+沉积速率(S)

相对海平面升降幅度=沉积速率× 时间+水深变化速率× 时间

对于任一时间点的水深和沉积高度(地层厚度), 可通过下式计算:

水深=可容纳空间变化速率(A)× 时间— 沉积速率× 时间

沉积高度=可容纳空间变化速率(A)× 时间— 水深变化速率× 时间

当可容纳空间变化速率大于沉积速率时, 即A> S, 这时沉积厚度或沉积高度仅与沉积速率有关; 反之, 如果可容纳空间变化速率小于沉积速率时, 即A< S, 这时沉积厚度受海平面高度控制, 两者关系如图3中时间点T1、T2、T3和T4所示。

如果全球海平面变化速率和构造沉降速率恒定, 当沉积速率为0时, 可容纳空间最大, 海平面升降高度等于水深; 当沉积速率大于可容纳空间变化速率时, 所有的可容纳空间都为沉积物充填, 这时沉积高度等于海平面升降幅度; 而在沉积速率小于可容纳空间变化速率的情况下, 只有部分可容纳空间为沉积物充填, 剩余空间为水体占据, 如图3所示。

图3

Fig.3

	Figure Option ViewDownloadNew Window
	图3 沉积速率与可容纳空间变化速率之间的关系Fig.3 Relationship between sedimentation rate and rate of accommodation change

层序地层等时对比的理论核心是海平面变化具有全球同步的特点, 而低频长周期海平面变化旋回又调控着高频海平面变化的形式和升降幅度, 因此以地层记录中高频沉积旋回厚度作为一个参数, 论证高频旋回厚度与三级海平面旋回过程的关系。

要求解沉积旋回厚度与全球海平面波动周期及幅度之间的关系, 关键是要将ka级时间坐标单位(时间域数据系列)转换为m级深度域数据系列。如果引入一个变量, 称为旋回厚度偏差,

图4

Fig.4

	Figure Option ViewDownloadNew Window
	图4 旋回沉积厚度与沉积速率之间的关系Fig.4 Relationship between depositional thickness of a cycle and sedimentation rate

即每个旋回厚度与其平均厚度值之差, 标记为± C, 其值可正可负, 就可以将时间域转换为深度域, 进而评估高频旋回叠置样式与三级海平面变化之间的关系。兹以图4和图5说明计算过程和结果。首先, 假设沉降速率为常数, 取0.055 m/ka作为计算值, 复合海平面变化的周期和振幅设定与图1和图2相同。为了简化计算过程, 计算时间长度取2000 ka, 时间坐标间距设定为1 ka, 则根据共计2000个数据点计算的可容纳空间变化曲线如图4-A所示。以任一周期作为步长, 对时间单位进行等间距采样, 则可求一个高频旋回的个数。图4-B中设定采样间距为10 ka, 则在时间长度为2000 ka的区间内, 计有200个周期为10 ka的数据点。设定不同的沉积速率变量, 根据沉积高度与沉积速率和可容纳空间变化速率之间的关系, 可求得每一个数据点的沉积厚度。图4-B表示沉积速率大于、小于和等于沉降速率时沉积厚度随时间变化的曲线。

为了考察沉积物堆积速率与可容纳空间变化速率之间的关系, 按旋回个数或旋回厚度分别计算了累积旋回厚度偏差, 其值与旋回数或深度的关系如图5所示。为了对比研究累积旋回厚度偏差曲线与长周期沉积层序界面的对应关系, 图5中同时标绘了全球海平面变化曲线。

图5

Fig.5

	Figure Option ViewDownloadNew Window
	图5 累积旋回厚度偏差与全球海平面变化关系的图解Fig.5 Diagram showing relationship between cumulative cycle thickness departures and eustatic sea level changes

从图5模拟计算结果可以发现, 如果全球海平面升降幅度和区域构造沉降速率一定, 随着沉积速率的变化, 旋回厚度偏差及其随时间的变化可能出现如下几种情况:

1)当沉积速率极高, 或沉积速率远大于沉降速率时, 即S> > T时, 由于沉积堆积厚度一般小于海平面高度, 这时高频沉积旋回的累积偏差曲线与低频海平面波动曲线一致, 累积旋回偏差曲线的波动幅度等于低频海平面波动幅度。这种情况常出现在浅水碳酸盐岩台地沉积区, 在阿尔卑斯三叠纪潮坪环境发育的Latermar旋回和Lofer旋回都是这一方面的实例(Anderson et al., 1990, Goldhammer et al., 1990, Osleger et al., 1990)。相反地, 如果沉积速率极低, 或者沉降速率等于或大于沉积速率, 即T≥ S时, 如在海盆深水区, 出现低速凝缩段沉积, 则沉积节律缺失。

2)在正常的海洋环境中, 沉积物堆积速率一般大于盆地沉降速率。分别取中等沉积速率0.055 m/ka和0.085 m/ka进行模拟计算, 结果如图5所示:旋回厚度累积偏差曲线与海平面变化曲线相比较不但发生了变化, 而且旋回厚度累积偏差曲线升降的幅度也明显小于海平面变化幅度, 但值得注意的是两条曲线升降转换节点的位置相同。这个结果提供了一个重要启示, 即高频旋回厚度累积偏差曲线并不一定能直接标定低频海平面变化幅度, 但在地层序列中旋回厚度累积曲线升降拐点指示沉积层序界面的位置, 这就为应用测井曲线划分沉积旋回, 进而判别三级沉积层序界面提供了理论依据。

四川盆地西北部中三叠统雷口坡组(T₂l)厚200~1000m, 主要由一套潮坪— 潟湖相白云岩组成, 中下部夹膏盐岩, 顶部见少量的灰岩; 雷口坡组与下伏嘉陵江组(T₁j)白云岩为整合接触, 两组之间以1层层厚数米的“ 绿豆岩” 分界, 顶部与上覆的须家河组(T₃x)碎屑岩之间以不整合面分隔, 以岩性突变面作为识别标志。

根据对江油马脚坝地区雷口坡组地表剖面和钻井岩心的观察, 雷口坡组可分3段, 其中雷三段进一步细分为上亚段和下亚段。雷一段由薄层状深灰色粉细晶白云岩、黄灰色膏质白云岩组成, 夹石膏层; 雷二段岩性为深灰色粉细晶白云岩、膏质白云岩和灰白色石膏层, 局部层位见鲕粒白云岩; 雷三段下亚段以发育褐灰色针孔状白云岩、藻粒白云岩、鲕粒白云岩、云质灰岩为特征, 底部以深灰色灰岩层与雷二段分界。雷三段上亚段以深灰色、褐灰色粉细晶白云岩为主, 局部层位见砂屑和藻屑白云岩。

测井曲线是一系列深度域振动信号序列的组合, 它是时间域周期性变化的叠加响应。测井曲线具有等间距采样的特点, 而且数据序列连续、纵向分辨率高, 因此可以满足数学含义上时间序列分析的基本要求(李庆谋, 1996; 陈茂山, 1999)。在常规测井曲线中, 自然伽马曲线能最敏感地反映泥质含量变化, 有效记录地层序列中沉积旋回的特点, 因而在高频旋回地层研究中得到了广泛应用(张占松, 2000; 李霞等, 2006; 吴怀春等, 2008)。

青林1井雷口坡组深度位于3562~4056m, 地层厚度494 m。为了全面了解雷口坡组米级沉积旋回的叠置样式, 选择同一深度的自然伽马数据进行研究, 共计采集数据点3952个, 测量间距0.125m, 自然伽马值变化范围在30~65 API。 自然伽马测井是测量地层中放射性元素衰变过程中放射出的伽马射线的强度, 其值大小取决于岩石中所含放射性核素(铀、钍、钾)的种类和数量。一般伽马高值对应泥岩, 低值对应石灰岩或白云岩。地层剖面中泥质含量从低到高的变化, 相应地反映了从石灰岩或白云岩到泥灰岩或泥岩的沉积旋回。

根据青林1井雷口坡组伽马值, 通过Excel提供的函数程序, 计算沉积旋回厚度并进行旋回个数计数的步骤如下:

1)首先采用线性内插法对原始数据进行校正, 求取等间距数据序列。

2)为了消除系统误差和测井曲线的长趋势偏移, 采用最小二乘法对原始曲线进行拟合, 取其与原数据序列之差作为新数据序列。

3)采用移动平均法过滤背景噪音, 突出沉积旋回的周期波动。移动平均的周期视旋回的大小而定, 实际计算过程中采用5点移动平均值消除高频噪音。

4)由于预处理过程中采用参数的差异, 可能造成计算结果的误差, 因此必须使数据中心化或归一化。实际计算过程采用一阶差分法处理, 保证最终计算结果具有一致性。

5)最后采用逻辑判别函数, 提取归一化伽马测井数据序列正负偏差数据, 计算出沉积旋回厚度。

图5是计算结果, 根据伽马测井数据, 求出青林1井雷口坡组计有238个高频沉积旋回, 旋回厚度变化为0.5~5.5m, 旋回平均厚度为2.1 m。

应该强调, 如果测井数据操作处理方法不同, 可能造成伽马数据的计算结果有一定的差别, 但在以上处理流程中, 伽马数据仅是一个用来求取旋回厚度的过渡参数, 旋回厚度的计算结果仅涉及伽马数据的相对变化幅度而与其绝对值大小无关。

根据青林1井雷口坡组旋回个数和厚度, 取纵坐标为旋回数, 或者以井深为纵坐标, 横坐标为累积旋回厚度偏差, 以米为单位, 可编绘2条旋回厚度偏差曲线, 结果如图6所示。分析旋回组合及其叠加样式在地层剖面中的变化, 可以识别层序界面和层序结构。

图6

Fig.6

	Figure Option ViewDownloadNew Window
	图6 川西北地区青林1井应用自然伽马测井数据计算旋回数及厚度的处理流程Fig.6 Procedure of cycle number counting and thickness calculation from natural gamma ray logging data of Well Qinglin 1 in northwestern Sichuan area

从图7中可以看出, 高频旋回厚度的累积偏差曲线表现为多级次组合谱系。按照旋回厚度累积偏差曲线的波动特征, 雷口坡组中至少可以识别出2个低频变化周期, 其曲线形态表现出低缓波状起伏且近于对称的特点, 两者分界在井深为3885m的位置。根据生物地层资料和火山成因绿豆岩测年数据(胡世玲等, 1996), 初步认为这2个低频变化周期可能对应于三级海平面变化, 曲线转换点在井深3885m处对应于海平面下降最低点, 即三级层序界面的位置, 最高点对应最大海泛期即凝缩段层位。进一步研究还可以发现, 在三级低频海平面变化曲线上, 叠加有一系列尖突的锯齿状波峰, 这在海平面上升曲线上表现得最为明显。如雷口坡组上部三级层序中可以观察到8个高频旋回组合, 下部三级层序中可以观察到6个高频旋回组合, 每个旋回组合分别由13个到16个高频旋回组合而成。这些旋回组合是构成海侵— 海退沉积旋回的重要单元, 相当于一般三级层序的准层序组或高级别的四级亚旋回(王鸿祯和史晓颖, 1998)。