研究地质历史时期的海平面变化对沉积环境恢复、砂体成因剖析和有利储层刻画等都具有重要意义。绝对海平面升降一度被认为是海相沉积环境和储集砂体发育的唯一决定因素。海平面周期性的变化造成了沉积环境的周期性变化，进而产生了具有周期性韵律的地层记录(Mitchum et al., 1977；Vail and Mitchum, 1977)。海平面变化研究成为石油地质和层序地层学工作者们的热点研究内容，世界上第一条地质历史时期的全球海平面变化曲线(也叫Haq曲线)便是在这一大背景下产生(Haq et al., 1987；Hunt and Tucker, 1992；Embry，1993)。随着研究的深入与认识的进步，人们发现除去绝对的海平面变化之外，构造升降甚至环境因素在特定的时期特定的区域也可以成为地层发育的控制因素。因此海平面的定义也被进一步改进——相对海平面被提出，用以代替绝对海平面和构造沉降的综合效果(可近似看做沉积可容纳空间，其变化可近似看做基准面变化，参见Schumm, 1993；Shanley and McCabe, 1994；Cross and Lessenger, 1996；Catuneanu, 2002, 2006；Li et al., 2015)。后文中的海平面即指的是相对海平面。

海平面变化研究方法有很多，主要有海岸上超点法、同位素地球化学法、古地磁法和古生物法以及各种方法的综合等(林孝先和侯中健，2014)。关于珠江口盆地海平面变化，前人早已有研究。秦国权(1996)通过海岸上超点法和浮游有孔虫、超微生物事件以及化石带的分析建立了珠江口盆地渐新世至今的海平面变化曲线。王春修(1996)利用浮游有孔虫和底栖有孔虫分析，结合与同位素曲线正向偏移曲线和Haq曲线的对比，讨论了珠江口盆地珠海组至万山组二级周期尺度的古水深变化和海平面升降过程。田世峰(2012)利用古地磁、古生物和放射性同位素方法，结合天文年代标尺，分析了惠州凹陷中新世的海平面和沉积速率变化。

然而，这些方法各有利弊。其中最经典的当属海岸上超点法，又称魏乐(Wheeler)图法。这种方法即Haq et al.(1987)采用的核心方法，其简单易操作，但过于依赖地震资料，因而精度普遍较低，而且对于资料的位置和范围也有要求，一般含有海岸上超的且延伸大而连续的地震资料比较适用(Zhang et al., 2015a)。古生物、古地磁和地球化学方法常常是基于采样的样本而进行实验室观察或测定的方法，具有可信度高的特点，但二者所依赖的采样数据较稀疏离散，因而在较大尺度内是比较可靠的办法，但无法针对局部进行高精度的研究。

近年来，大数据分析获得前所未有的重视(张旗和周永章，2017)。人们期望，在认识事物不全面的情况下，可以利用已有经验帮助做出大致合理的判断、决策，以后如有客观的新信息、新数据更新最初关于某个事物的信念后，就会得到一个新的、改进了的信念。贝叶斯-拉普拉斯原理因此获得新的活力(周永章等，2017a)。

在上述背景下，本文以珠江口盆地惠州凹陷南部为例，利用Fischer图解法和INPEFA曲线法研究K系列地层沉积时期的相对海平面变化，并应用于大数据贝叶斯-拉普拉斯(Bayes-Laplace)分析方法，建立新信息、新数据不断更新下的开源研究策略，获得相对海平面在时间域(geological time domain)的高精度刻画，为小尺度古环境研究和精细油藏描述提供支撑。

1 地质背景

惠州凹陷位于珠江口盆地中北部地区, 面积约10000km²，是珠江口盆地重要的油气富集区之一。其北靠北部断阶带，南接中央隆起带的东沙隆起，东西分别与惠陆低凸起和西江凹陷相邻，呈北东走向，是晚新生代古珠江三角洲-滨岸沉积体系形成、发育及演化的主要场所(国朋飞，2015)。H21油区位于惠州凹陷的南部，惠西半地堑东南角，其三面环洼，水深约116m，具有长期上隆的构造背景(图 1)。

凹陷新近系珠江组是该区油气勘探的主要目的层。珠江组厚度基本上分布在750m和350m之间。从珠江组开始，南海扩张进入冷却阶段，全盆地进入稳定的热沉降阶段。K系列地层发育于珠江组内部，位于T50地震反射截面之上，分布面积主要在惠州凹陷南部。K系列沉积时期，随着南海运动的结束，盆地由断陷向拗陷过渡，惠州凹陷整体进入平稳沉降时期(高慧，2007；吴娟, 2010, 2013)(图 2)。受海侵影响，K系列地层发育三角洲沉积和三角洲前缘改造形成的滨外沙坝、滨外沙席沉积(Zhang et al., 2017)。

在H21油区，K系列的三级层序地层界面主要有SB17.5、S7、SB18、MFS18.5，共可以划分为1.5个三级层序(图 3)。K系列包含4个四级层序，每个层序可以分为一个海侵体系域(TST)，和一个高位体系域(HST)。这些体系域自下而上分别被命名为TST1.2、HST1.2、TST1.3、HST1.3、TST2.1、HST2.1、TST2.2、HST2.2和HST3.1(丁琳等，2014；杜家元等，2014；Zhang et al., 2017)(图 4)。

2 Fischer图解与INPEFA 2.1 Fischer图解法

Fischer图解法最早是由Fischer(1964)提出，主要是通过计算高频旋回(如米氏旋回)的累计厚度偏差(偏差为相对旋回厚度平均值)来估算地层沉积时期的可容纳空间以及水进水退变化。原始的Fischer图解主要用于碳酸盐岩野外剖面，随着后来学者的补充，Fischer图解也被证明可以应用到非碳酸盐岩剖面，同时资料也不限于野外露头，利用测井曲线也可以进行研究(Day, 1997；Zhang et al., 2015b)。深度域(depth domain)的Fischer图解可以表示为公式(1)：

(1)

式中，CMDT为累计厚度偏差(单位为m)，CMDT₀和CMDT_N；TOP₀为第零个旋回的顶深，就是研究层位底深(单位为m)；TOP_i为第i个旋回的顶深(单位为m)；T为地层的总厚度(单位为m)；N为旋回的个数(一般要>50个)。

式(1)中的N可用能代表旋回变化的测井曲线通过滤波得出，也可用INPEFA曲线算出(见后文2.2)。因为计算的是地质历史时期的旋回厚度变化，上式中TOP、T等参数都要用现今地层进行压实校正得出。压实校正一般有两种假设，第一种假设是压实前后岩石骨架体积不变，另一种假设是压实前后密度不变。其中用得最多的是前者，即公式(2)(修改自林孝先和侯中健，2014)：

(2)

式中，T为沉积物原始厚度，即沉积物总量S(单位为m)；z为深度(单位为m)；z₁为现金底层的底深(单位为m)；z₂为现今地层的顶深(单位为m)；ϕ(z)为孔隙度函数(单位为%)。

2.2 INPEFA曲线法(最大熵谱分析法)

地层的泥质含量是反映沉积时水深的重要依据，因此可以通过泥质含量分析地层的沉积旋回，也可用来间接推算水体深度的的变化。但是测井曲线本身并不能直接反映水体的深度大小，必须通过一定的变换才能有效反映水体的深度变化。合成预测误差滤波分析(INPEFA)便是一种常用的变换，GR曲线经过变换后得到的INPEFA曲线能够有效反映水体深度的变化(Zhou，1999；翟玉兰，2009；薛欢欢等，2015)。本文有个重要的假设：泥质含量的高低与水体深度的大小呈正相关(现实问题中可能更复杂，比如潮坪沉积中泥质含量与水体深度很可能呈反相关，这里的假设只考虑理想情况)。

要进行INPEFA分析，首先要对测井曲线进行滤波处理(去除奇异值，进行光滑校正等)(图 5a)。然后运用最大熵谱分析(Maximum Entropy Spectral Analysis)，采用外推自相关函数的办法由部分已知点推算下一个点的最可能的值(李景哲，2013)(图 5b)。所谓“下一个点”本身具有实际观测值，这个值与预测出来的值不见得处处一致，因而也就会有误差。误差大于0，代表观测值大于预测值，也就是真实的沉积地层泥质含量比本来的变化趋势更大了，说明水体加深了；反之，误差小于0，则代表观测值小于预测值，真实的沉积地层泥质含量比本来的变化趋势小了，说明水体深度变小了。一系列的误差点组成的曲线称为预测误差滤波曲线(PEFA曲线)(图 5c)。但是PEFA曲线只是代表变化的趋势(相当于加速度)，要得到变化后的结果(相当于路程)，还需要对PEFA曲线积分。PEFA曲线积分后得到的曲线便叫做INPEFA曲线(图 5d)。

水深变化是相对海平面变化和沉积作用共同作用的结果，三者关系如公式(3)：

(3)

写成积分的形式如公式(4)(深度域)：

(4)

式中，A为相对海平面或可容纳空间(单位为m)；为深度域沉积充填速率(单位为1)；为深度域水深变化速率(单位为1)。

沉积充填速率与沉积物的总供给量密切相关。沉积物的总量可以通过现今沉积物的厚度经过压实校正得出(见上文2.1)。另一个重要参数便是水深变化速率。由前面的分析可知，INPEFA曲线可以衡量水深的变化。但INPEFA曲线本身只是定性反映水深变化，如果能将这种变化定量化，便可以计算出该参数。惠州凹陷地处陆架地区，水深不超过200m，而最小水深可能为0(考虑到该地区有过SB18沉积间断，水深很可能被消耗到0)。因此可以利用INPEFA数据构建出水深方程，如(5)：

(5)

式中，w是一个经验系数，指示最大水深。

3 海平面变化曲线

通过2.1中的Fischer方法，带入孔隙度函数(参考中海油深圳分公司内部经验公式，ϕ=40e-0.00038z)还有现今地层的顶底深度，便可算出沉积物原始沉积厚度(H21井所在位置的K系列原始沉积厚度约652.67m)。通过上述方法，本文建立了H21井K系列地层的Fischer曲线(图 6)。该曲线反映的海平面总体是先下降后上升，在后半段总体的上升过程中也有小幅度的下降过程。该曲线纵坐标为累计厚度偏差，横坐标为旋回序号，为了更好的分析也可以在深度域展示(图 7)。

而与此相对比，通过2.2中的INPEFA曲线方法也可以得到K系列沉积时期海平面变化(图 7蓝色线)。该曲线总体上是上升的，但上升的速率在不同的位置不尽相同。总体而言，四级层序海侵体系域(TST)的上升速率要稍高于与其相邻的同级别高位体系域(HST)的速率。如最下部的HST3.1的上升速率就比紧邻其上的TST2.2要高很多。

4 讨论与分析

通过海平面变化的综合分析，一共得到了两种不同的曲线(图 7中黑色线和蓝色线)。一种是通过Fischer图解得到的旋回厚度累计偏差曲线，另一种是通过最大熵谱分析得到的相对海平面变化曲线。

通过INPEFA曲线得出的相对海平面(或可容纳空间)在深度域的变化主要有以下规律：在海侵体系域时期，相对海平面快速上升；而在高位体系域时期，相对海平面缓慢上升。相对海平面总体处于上升状态，上升的总距离近似为地层压实校正后的总厚度。这与地质常识相吻合：从2770m到2406m(时间上18.5~17.5Ma)，可容纳空间变化等于构造沉降和全球海平面变化之和，二者共同创造了650多米的可容纳空间(可容纳空间在本次研究区主要是构造沉降创造，这段时期全球海平面变化的升降幅度不到50m(Haq et al., 1987)，可以忽略不计)。而且这个结果与前人研究惠州凹陷沉降历史的数据也接近(赵中贤等，2010)。经验系数w值的选择与曲线的形态直接相关。w值越大，模拟出来的曲线形态就越类似于INPEFA；w值越小，模拟出来的曲线形态则越接近直线。在实际操作中，要充分研究目标井所在的古地理环境，合理设定w取值。如果能有古生物和元素地球化学资料支撑，w的取值将更准确。需要指出的是，由于测井资料是地层记录的反映，对于没有地层记录的沉积间断或剥蚀现象，测井资料的各种变换也无法很好的表征(比如下降体系域所对应的时期，局部地区可能就没有地层记录)。

后者则是假定高频旋回的时间相当的前提下计算出来的各个旋回相对于全部旋回厚度平均值的累积偏移。这个假定条件通常在沉积作用主要受控于天文旋回的地区更适用。本次研究区由于构造是最重要的控制因素，构造因素与天文因素的混合可能会降低该假设的适用性。在本地区该曲线不是可容纳空间(或相对海平面)的近似。因为从结果上看，Fischer曲线以零开始，中间大部分是负值，最后又以零结束，总的效果是回到了原点。如果这条曲线代表可容纳空间的变化，那么长时间的负值可容纳空间条件下绝不可能有几百米厚的地层记录。从原理上说，累积厚度偏移实际上为沉积充填速率的升高或降低(Fischer图解假定条件为各旋回周期相当，那么沉积充填速率就是旋回厚度唯一的区分因素)。研究点H21井所在的位置在K系列沉积时期一直处于内陆架环境，属于水下沉积环境，沉积充填的短板一直是沉积物供给。高位体系域时期，沉积物可以源源不断供给到研究点所在位置，而海侵体系域时期，因为水深的增大，沉积物被大量限制在河口附近，难以送到研究点位置。因此，与其说这条曲线反映可容纳空间变化，不如说这条曲线反映沉积物能够供给到研究点的速率和剩余可容纳空间(水深)变化。另外，通过图 7可以看出一些大致规律：在三级高位体系域时期，曲线向上不断变小；而三级海侵体系域时期，曲线向上不断变大。图中的两个红线段是累积偏移逐渐增加的两个主要时间段，而这两个时间段恰巧对应K系列的两个主要的储集层HST1.3和HST2.2。这也说明了陆源碎屑供给能力与富砂程度的响应关系。

综上，通过INPEFA曲线计算的相对海平面变化更接近实际情况。不管是三级、四级体系域的形成，还是构造沉降历史都能很好地吻合(没有地层记录的沉积间断或者剥蚀现象除外)。Fischer图解在本次研究中并没有直接反映相对海平面(或可容纳空间)变化，而是反映水深(剩余可容纳空间)变化。原因可能是本次研究区一直处于内陆架环境，水体深度有较大变化，沉积物厚度变化并不完全代表可容纳空间变化，但是其仍然在预测砂体的富集位置方面很有价值。当然，如果在水深很浅且变化可以忽略的盆地，沉积充填的短板则是可容纳空间变化速率，可容纳空间上升多高，地层就沉积多厚，累积的旋回厚度偏移自然也就能够反映可容纳空间变化。

5 挑战和展望

随着大数据科学的兴起，多源异构数据的综合分析越来越成为地球科学研究和地质工程应用的热点(郭华东，2014；李超岭等，2015；吴冲龙和刘刚，2015；严光生等，2015; 张旗和周永章，2017; 赵鹏大，2015；周永章等，2017a)。

上述所得海平面曲线是一个针对局部的高分辨率曲线，其主要依赖的资料是测井资料，而其与基于样本测试的传统方法(古生物、古地磁和地球化学等方法)相比，所依赖的资料各有侧重，前者在深度域连续且具有较高精度，后者则可以提供部分离散点时间数据。

随着工作程度的深入，离散点时间数据会逐步增加。设想离散点时间数据充分多以后，人类是否可以获得任意点的时间信息呢？这涉及沉积成因模型的建立，需要许多观察和数据作为支撑。

从统计观点看，对数据信息进行挖掘，有经典统计和贝叶斯-拉普拉斯两类不同的思路(周永章等，2017a)。

经典统计着重频率统计，它强调：只要反复观察一个可重复的现象，直到积累了足够多的数据，就能从中推断出有意义的规律，揭示一切现象产生的原因。从理论上讲，它既不需要构建模型，也不需要默认条件，只要进行足够多次的测量，隐藏在数据背后的原因就会自动揭开面纱。如果数据量足够大，人们完全可以通过直接研究这些样本来推断总体的规律。

但当存在着大量数据，但数据又可能有各种各样的错误和遗漏的时候，如何才能从中找到真实的规律。这是贝叶斯-拉普拉斯方法关注的问题(Lake et al., 2015; 周永章等，2017a)。

贝叶斯-拉普拉斯方法则认为，在人类认识事物不全面的情况下，可以利用已有经验帮助做出大致合理的判断、决策，以后如有客观的新信息、新数据更新最初关于某个事物的信念后，就会得到一个新的、改进了的信念。这就是说，当一个人不能准确知悉一个事物的本质时，他可以依靠与事物特定本质相关的事件出现的多少去判断其本质属性的概率。用数学语言表达就是：支持某项属性的事件发生得愈多，则该属性成立的可能性就愈大。与经典统计统计学方法不同，贝叶斯-拉普拉斯方法建立在主观判断的基础上，先估计一个值，然后根据客观事实不断修正。

本文认为，可以通过现代沉积相分析获得一个定量化的分形模型(针对实际年龄)。在实际测井剖面中，可以获得层序地层和系列年龄测定(离散数据)。人们可以利用离散的年龄数据和层序地层，依据定量化的分形模型预测剖面上任意一个点的年龄值。开始时，这种预测结果是比较粗糙的，但随着实测年龄数据的增加，新信息、新数据不断更新条件下，贝叶斯-拉普拉斯方法不断迭代给出更新结果。这样，可以实现相对海平面在时间域的高精度刻画，为小尺度古环境研究和精细油藏描述提供支撑。

贝叶斯-拉普拉斯方法的数学形式如下(Lake et al., 2015；周永章等, 2017a, b)：

(6)

式中:P(A)是先验概率，P(A|B)是后验概率，表示在以后B事件发生的条件下A事件发生的条件概率。

图 8展示了贝叶斯-拉普拉斯方法分析流程图。在本文中，贝叶斯-拉普拉斯综合分析的具体过程包括(图 8)：

(1) 对古代碎屑岩沉积地层，充分利用获得的已知数据，通过岩心分析明确地层的沉积相类型，同时运用层序地层学分析的手段厘清不同井孔间地层的相对年代关系；

(2) 对上述沉积相的现代沉积对应物，以“将今论古”为原则，通过理论分析和数值验证提出适用于古代沉积相类型的分形插值函数，构建相应重建模型；

(3) 运用Fishcer图解法和INPEFA曲线等方法，通过定量研究和分析比较，得到了研究层位所对应的深度域海平面变化曲线；

(4) 将(3)中曲线离散化至可分辨的最小成因单元并编号，运用(2)建立的模型恢复各编号所对应地层间的缺失细节，依据顶底年代关系建立初步的时间域海平面曲线；

(5) 根据岩心/岩屑所处井孔的位置、地层单元的位置和岩心本身的性质(是否有沉积同生矿物、是否有古生物化石、是否有灰岩或火山碎屑岩夹层以及是否是定向岩心等)挑选样本，测定各样本的海平面信息和年龄信息；

(6) 在(1)至(5)的基础上，通过理论分析的方法明确多源间的拓扑结构，厘清各数据间的条件概率关系，建立数据同化的贝叶斯网络；

(7) 运行贝叶斯网络，调整更新海平面变化曲线；

(8) 丰富新的样本，重复(5)到(7)的过程；

(9) 完善海平面曲线，从理论分析的角度揭示其主控因素和地层记录的过程。

目前该设想还面临着两个重要难点：(1)提出适用于古代碎屑岩地层的分形插值模型；(2)明确多源数据的拓扑结构和各个变量间的条件概率关系。

这两条难点还有待进一步攻克，其解决将有助于海平面在时间域的高精度刻画，不但可以使得沉积过程“源到汇”的完整再现成为可能，还可以为地下沉积体的精确推断与表征提供重要依据。

6 结论

通过对比分析，INPEFA方法得出的曲线更能代表相对海平面变化。研究区K系列沉积时期，海平面总体上是上升的，但上升的速率在不同的位置不尽相同。总体而言，四级层序海侵体系域(TST)的上升速率要稍高于与其相邻的同级别高位体系域(HST)的速率。海平面主要以交替的快速上升(海侵)和缓慢上升(海退)为特点。Fischer方法得出的曲线在本地区不直接反映相对海平面(可容纳空间)变化，但可以指示富砂层段的位置。

实现从深度域的刻画到时间域的刻画是海平面变化研究的一项挑战，特别是在年代资料比较离散和稀疏的地区。贝叶斯-拉普拉斯方法为这项挑战的解决提供了基于大数据的新思路。开始时，人们可以利用地层的层序结构、分形特性以及离散的地层年龄等数据建立初始的较粗糙的年代关系，随着实测年龄数据的增加，新信息、新数据不断更新，这些年代信息不断细化，进而可以实现相对海平面在时间域的高精度刻画。