沉积盆地按其成因可以分为伸展盆地、挠曲盆地和走滑拉分盆地^[1]。其中发育在被动陆缘之上的伸展盆地是主要的含油气盆地类型之一, 对其进行研究有助于油气资源的勘探与开发^[2]。目前针对伸展盆地开展盆地分析的方法较多, 大体可分为正演和反演两类, 反演可以恢复地质历史时期的地质状态包括沉降史、沉积史及断层发育期次等, 而正演方法则可以在反演的约束下探究盆地的演化过程及动力机制等问题。其中, 常用的反演方法包括平衡剖面方法^[3]、回剥方法^[4-5]等, 正演方法包括挠曲悬臂梁模型^[6-7]与数值模拟方法^[8-9]等。基于以上方法的盆地分析相关的商业软件也较多, 有Flex- Decomp、Stretch、Move以及Thermodel^[10]等。目前存在的问题是, 各种方法都有自身的优劣之处, 难以通过单一方法获得接近真实地质情况的所有结果。针对于此, 本研究论述了目前几个常用盆地分析软件的原理与方法, 对其各自的优势与不足进行分析, 并探讨了充分利用各自优势的软件组合使用方法, 对盆地的综合分析和油气的勘探开发具有较高的实用价值。

1 盆地分析中的反演方法

反演指的是以现今的地质状态为基础, 以古水深或断层发育期次等数据为约束恢复出初始地质状况或某个阶段地质状态的过程。通过反演可以反推盆地的构造演化过程, 对盆地的沉积、沉降速率及断层发育有更加清晰的认识。

1.1 构造沉降恢复

沉降是地壳的一种垂直运动, 是盆地发育的主导因素之一。盆地的沉降按成因可以分为构造沉降(由地壳或岩石圈动力学产生的盆地沉降)和非构造沉降(荷载导致的沉降)^{[1, 11]}。构造沉降曲线在完整的裂陷旋回上表现出两段式:裂陷阶段的快速沉降及裂后的缓慢沉降^{[1, 11]}。盆地的沉降与岩石圈的伸展减薄有关, 通过构造沉降的分析, 我们可以了解岩石圈或地壳的伸展减薄程度^[12-13]。因此, 构造沉降量和构造沉降速率反映了盆地构造演化的特征, 对盆地构造沉降史的恢复有助于了解盆地的形成机制。

总沉降量去除沉积物与水的重力作用造成的沉降量后即为构造沉降量, 只考虑沉积物的荷载而未去除水的荷载作用得到的沉降量为水载沉降, 水载沉降应用较为广泛^[14-17]。水载沉降的计算过程包括逐层去除沉积地层并进行去压实校正、沉积物重力荷载校正、水深变化校正和相对海平面变化校正^[18]。

1) 去压实校正

去压实校正是指在逐层剥去沉积层后, 其下地层在上覆荷载减小的情况下厚度恢复到该地质历史时期的真实厚度^[15]。校正过程假设沉积物在压实过程中其骨架体积和横向宽度守恒, 地层体积变化是由孔隙流体流动与孔隙体积变化引起的^[19]。地层的压实程度仅与埋深有关, 孔隙度与深度的关系通常采用指数形式^[15]:

$ \mathit{\Phi } = {\mathit{\Phi }_0}{{\rm{e}}^{ - cz}}, $

(1)

其中, Φ为孔隙度, Φ₀为初始孔隙度, c为压实系数(km^–1), z为自海底算起的深度(km)。

2) 沉积物重力载荷校正

盆地中沉积物的重力作用会促使沉降量进一步变大而产生负载沉降。在Airy均衡条件下, 沉积物造成的负载沉降效果可按如下公式进行校正:

$ {Y_{\rm{s}}} = S\left( {\frac{{{\rho _{\rm{m}}} - {\rho _{\rm{s}}}}}{{{\rho _{\rm{m}}} - {\rho _{\rm{w}}}}}} \right), $

(2)

Y_s为沉积物重力载荷校正后的基底深度, ρ_m为地幔密度, ρ_s为沉积物柱的平均密度, ρ_w为水的密度, S为去压实校正后的沉积物厚度。

3) 相对海平面变化校正和水深变化校正

海平面在地质演化历史中发生了变化, 而现今的海平面是计算构造沉降基准面, 所以水深变化数据还要考虑相对海平面变化(相对于现今海平面位置的变化)。其次, 过去地质历史时期的海平面位置与现今海平面位置之间的差值(水柱高度)会产生重力荷载作用而产生荷载沉降(S₁), 应予以剔除^[18-19]:

$ {S_{\rm{l}}} = \Delta {S_{\rm{l}}}\left( {\frac{{{\rho _{\rm{m}}}}}{{{\rho _{\rm{m}}} - {\rho _{\rm{w}}}}}} \right), $

(3)

ΔS₁为海平面相对于现今海平面的高度。

目前, 古生物资料和地球化学资料是获取古水深(W_d)的主要手段除此之外还可以依据沉积相资料来进行估计。古水深数据难以准确测定, 需要综合多种资料加以判断, 所以此项是构造沉降计算中最大的误差来源^[17]。

当岩石圈的强度非常弱, 可以认为已达到完全的重力均衡, 水载沉降量可以表示为:

$ {Y_{\rm{w}}} = {Y_{\rm{s}}} - {S_{\rm{l}}} + {W_{\rm{d}}}, $

(4)

$ {Y_{\rm{w}}} = S\left( {\frac{{{\rho _{\rm{m}}} - {\rho _{\rm{s}}}}}{{{\rho _{\rm{m}}} - {\rho _{\rm{w}}}}}} \right) - \Delta {S_{\rm{l}}}\left( {\frac{{{\rho _{\rm{m}}}}}{{{\rho _{\rm{m}}} - {\rho _{\rm{w}}}}}} \right) + {W_{\rm{d}}}, $

(5)

水载沉降量中的水柱载荷仍对地壳具有重力作用, 当水深大于1 000 m时水的载荷压力超过10 MPa, 故水的荷载作用影响较大不可以忽略, 可以进一步转化为空盆沉降(公式(6), 公式(7))并计算空盆沉降速率^[20]。

$ {Y_{\rm{E}}} = {Y_{\rm{w}}}\left( {\frac{{{\rho _{\rm{m}}} - {\rho _{\rm{w}}}}}{{{\rho _{\rm{m}}} - {\rho _{{\rm{air}}}}}}} \right), $

(6)

$ \begin{array}{c} {Y_{\rm{E}}} = {Y_{\rm{w}}}\left( {\frac{{{\rho _{\rm{m}}} - {\rho _{\rm{w}}}}}{{{\rho _{\rm{m}}}}}} \right)\\ = S\left( {\frac{{{\rho _{\rm{m}}} - {\rho _{\rm{s}}}}}{{{\rho _{\rm{m}}}}}} \right) - \Delta {S_{\rm{l}}} + {W_{\rm{d}}}\left( {\frac{{{\rho _{\rm{m}}} - {\rho _{\rm{w}}}}}{{{\rho _{\rm{m}}}}}} \right) \end{array}。$

(7)

空盆沉降排除了水的荷载作用对沉降的影响, 其沉降量完全由构造作用产生, 通过分析构造沉降量及构造沉降速率可以对演化过程中的构造作用有更加清晰的认识。

1.2 挠曲回剥

回剥是盆地分析的重要研究方法, 常用的回剥方法包括一维Airy均衡回剥与二维挠曲均衡回剥。一维回剥方法利用Airy均衡模式^[21], 认为沉积荷载造成的影响仅作用于该点的垂直方向上, 对侧向的区域没有影响。二维回剥方法^[22]认为回剥过程中应该考虑岩石圈的弹性强度, 即某一点的载荷变化不仅会影响该点处的状态还会造成岩石圈的挠曲反应, 从而对侧向的区域产生影响。

二维挠曲均衡回剥比一维Airy均衡回剥更加符合真实的地质情况。该模型由于不能恢复断层作用导致的水平方向上的变化, 所以难以恢复伸展盆地裂谷期的演化过程, 仅可以较好地恢复裂后期的演化过程。二维挠曲均衡回剥在一些盆地中得到了很好的应用^[23-25]。

Flex-Decomp软件先假定均一的拉张因子对剖面进行逐层回剥恢复, 直到恢复到裂谷作用结束时的地质状态, 后利用Stretch软件(下文将论述)正演模拟该剖面状态以获得剖面拉张因子, 后将Stretch得到的拉张因子再次带入Flex-Decomp中再次回剥直至达到最佳拟合^{[4, 26]}。Flex-Decomp的计算步骤如下^[4-5]:

1) 移除顶部的水层并进行挠曲均衡反应计算。

2) 移除剖面中的最顶部地层, 下伏地层计算去压实厚度变化。

3) 计算移除顶部地层荷载之后的挠曲均衡反应, 并考虑剩余沉积层孔隙度的增加。

4) 利用Mckenzie ^[12]的裂后热沉降模型依据估计的拉张因子与裂谷作用的年代反算热隆升。

5) 如果海平面变化过大则进行进一步的修正, 并增加海水层的荷载影响, 计算挠曲均衡反应。

6) 以上步骤建立了移除顶部地层后的挠曲均衡剖面。

7) 重复步骤1—步骤6回剥剩余层位直到裂后期的基底。

如果剖面有裂谷期的水深标志(比如暴露剥蚀等), 可以通过人为调整剖面的拉张因子以使剖面与古水深有一个很好的对应, 而此时得到的拉张因子剖面是反演得到的岩石圈的拉张因子^[27]。本软件模型认为岩石圈各层的伸展减薄程度是相同的, 在其他参数确定的情况下, 若岩石圈在伸展过程中未发生解耦, 则以此方法可以较好的反演岩石圈拉张因子。

该软件地层厚度恢复采用的方法^[15], 可以回剥得到演化过程中每一层沉积层序初始沉积时的厚度, 进而可计算每个沉积期的沉积速率。

1.3 剖面恢复

平衡剖面技术是剖面恢复的有效方法。其在发展之初是为了确定同心褶皱之下拆离断层的深度^[3], 后来这一方法被应用于逆冲推覆、伸展以及盐底辟、反转构造的研究中^{[18, 28]}。平衡剖面的恢复基于以下假设:地层的初始状态是水平的且在变形过程中遵守物质守恒的原则(厚度、面积、体积守恒), 剥蚀与压实是体积变化的唯一方式, 脆性断层是变形的主要方式且褶皱与断层有关, 压溶与构造压实对体积的影响较小^[18-19]。因此, 平衡剖面技术是以变形符合平面应变条件为前提按照面积/长度守恒的原则进行构造恢复^[29]。

2D Move是应用较广的二维平衡剖面恢复软件, 其使用的地层恢复算法有很多种, 应用较多的有单剪去褶皱法和斜剪切法。

单剪去褶皱算法以垂直或斜线方式恢复到基准面(水平面或给定基准面)以达到去褶皱的目的, 该方法认为沿剪切线方向的剪切矢量相同, 其特点是计算前后的剖面上的层面积守恒。实际的构造变形过程中地层的长度会发生变化, 所以恢复到基准面位置后地层的长度会变小^[18-19]。

斜剪切算法的基本原则是体积守恒和剪切矢量长度(沿剪切矢量方向断层面和上盘的标志层之间的距离)守恒。在恢复断层时假设断层的变形仅发生在上盘且沿系列剪切线平行发生, 上盘的变形穿透断层上盘且与断层几何状态有关。影响该算法的主要参数包括移动方向、剪切矢量和水平断距。其中剪切矢量的大小由水平断距参数定义, 斜剪切算法可以较为有效恢复犁式断层, 其倾角随深度增加而减小, 这种倾角的变化与上盘的变形是对应的^[18-19]。

平衡剖面技术在恢复构造状态的同时还要进行地层发育恢复, 地层发育史恢复需要的参数和需要进行的校正与构造沉降恢复中相同。

平衡剖面技术的优势在于可以逐步恢复各个地质时期地层发育、断层和褶皱状态, 并根据古水深数据得到一个相对可靠的古地质剖面。

2 盆地正演方法

正演即假设初始条件, 进而在反演结果及地质观测的约束下, 重现盆地演化过程的方法。正演方法可以检验反演结果并探究盆地演化中的控制因素。

2.1 挠曲悬臂梁模型

挠曲悬臂梁模型将岩石圈分为3层, 分别为上地壳、下地壳和岩石圈地幔, 其中上地壳呈脆性以简单剪切的形式伸展, 而下地壳和岩石圈地幔呈现塑性, 在张力的作用下以纯剪切形式伸展^[6]。上地壳产生的断层上下两盘由断层相连接像两个相互作用的挠曲梁, 使断层的上盘产生坳陷下盘产生隆升。上地壳的脆性伸展产生的空间由塑性的下地壳物质补充, 同样下地壳的物质上涌产生的空间由岩石圈地幔补偿, 进而会造成莫霍面的扰动^{[4, 6]}。该模型认为上地壳中的大型基底控盆断层都是平直断层, 在这里为了更好的模拟陆缘上的盆地形态采用铲式的断层来模拟上地壳的盆地形态。挠曲悬臂梁模型通过设置断层参数与反演得到的模板进行对比, 认为上地壳的脆性应变由下地壳的塑性应变补偿, 利用McKenzie的纯剪模式^[12]对下地壳的拉张因子进行计算, 并认为上下地壳与地幔的拉张因子相同, 来表征上地壳产生的伸展^[6], 当同一个点有多个断层的作用时其拉张因子等于此处各个断层产生的拉张因子的乘积^[7]。

利用Flex-Decomp得到的地质模型通过设置断层参数进行拟合, 得到岩石圈拉张因子(模型认为岩石圈各层的伸展减薄程度及拉张因子相同)。可进一步计算裂后热沉降、沉积物充填, 并可计算重力值以及热流和沉降史。软件可实现多期的裂谷作用。

软件计算的拉张因子在其他基础参数确定的情况下仅与输入断层参数有关, 而断层仅代表了脆性上地壳的伸展减薄作用。该模型考虑了上地壳的挠曲响应, 在已知断层参数的情况下可以将断层参数输入本软件中计算上地壳的伸展减薄。此方法不要求与地质模型的拟合效果, 仅通过该模型计算不同断层造成的伸展减薄(拉张因子)。

通过人为调整有效弹性厚度对比地震剖面的凹陷形态与模拟结果, 可以推测岩石圈的有效弹性厚度^[30], 有效弹性厚度越大断层影响的范围就越大, 造成的凹陷范围更大, 坡度相对平缓, 相反若有效弹性厚度较小则形成的凹陷范围较小凹陷边界坡度较大。在岩石圈解耦时代表了上地壳的有效弹性厚度而非整个岩石圈的有效弹性厚度。

2.2 二维粘弹塑性热-力学模型

有关岩石圈拉张的热-力学数值模拟程序有很多^[31-32], 但主要模拟方法差别不大。本文以Pérez- Gussinyé等开发的二维粘弹塑性热-力学有限元模拟程序^[33]为例进行阐述。该程序是在热-力学有限元计算程序MILAMIN^[34]的基础上发展而来的。在模型的左右边界给定恒定的拉张速度, 上部为自由边界, 内部单元采用拉格朗日三角网格划分, 即单元随着物质的运动而运动。三角网格分辨率可根据应用需求进行不同的划分。网格模型分为4层:上地壳、下地壳、岩石圈地幔以及软流圈^[35]。程序利用物质守恒、能量守恒和动量守恒定律相结合得到斯托克斯方程和热扩散方程。分别用来求解重力场作用下粘性流体的流动速度与壳幔变形过程中的温度场。在MILAMIN求解器的基础上增加多种物理化学过程, 可以在不同程度上模拟复杂的地质过程, 如陆缘张裂过程中的岩浆作用及蛇纹石化过程等。

采用的流变学模型为

$ \dot \varepsilon = A\sigma _{\rm{d}}^n\exp \left( { - \frac{Q}{{RT}}} \right), $

(8)

其中, $\dot \varepsilon $为应变速率, A和n为实验室测得的与物质相关的常数, σ_d为偏应力, Q为蠕变活化能, T为温度, R为气体常数。

其他热力学参数包括各层的密度、热膨胀系数、体积模量、剪切模量、热容、热导率、放射性生热率、初始摩擦系数和内聚力, 其中放射性生热主要发生在上地壳, 因此上地壳厚度的加厚会影响地壳乃至整个岩石圈的温度场分布。

为更加真实地模拟地质过程, 模型考虑粘弹塑性变形, 设定自由表面以获得稳定的动力学地貌, 利用应变软化以模拟断裂和剪切带, 并可选择是否增加地表沉积过程, 即侵蚀和沉积。岩石受力一旦达到屈服强度, 很容易沿屈服面发生变形, 因此断裂面或剪切带会产生集中形变^[36], 这个过程可用应变软化来实现。在脆性变形域, 通过降低摩擦角来实现应变软化以模拟脆性断层的发生; 在粘性变形域, 变形方式由位错蠕变转为扩散蠕变时降低了有效粘滞系数^[37], 也可以发生粘性应变软化, 岩石矿物发生结晶优选定向, 从而使得剪切带在该方向更加容易发生^[38-39]。

本数值模拟程序可以在地质观测数据的约束下对伸展盆地或者被动陆缘的形成过程进行模拟。通过对比模拟结果与地质观测的结果, 选定合适的模拟模型参数, 这些参数在一定程度上代表了真实的地质参数, 比如伸展速度, 下地壳强度等。

3 讨论

不同的盆地模拟方法或软件其研究对象及侧重点不同, 像平衡剖面恢复及回剥(Airy均衡回剥及挠曲回剥), 其研究对象为浅部的沉积盆地, 而挠曲悬臂梁则以挠曲回剥的结果作为约束研究岩石圈的伸展减薄, 但由于在真实的地质过程中这是一个综合的过程, 即当板块受到拉张应力时, 在岩石圈浅部会发育断层形成沉积盆地, 岩石圈整体会发生伸展减薄, 下部软流圈会发生上涌。所以针对不同研究对象的正反演方法得到的结果可以相互约束验证从而得到一个相对可靠的地质演化过程。综合利用每种方法或软件的优势是进行盆地模拟研究的最优选择。

在岩石圈强度较弱的情况下, 一维Airy均衡的回剥计算构造沉降的方法相对可靠^[40], 整个构造演化过程中裂后期的构造沉降是热沉降主导的, 断层作用减少, 剖面没有水平方向的变化, 仅存在垂向上的变化, 所以在其他参数可靠的情况下, 裂后期构造沉降史的恢复具有较好的可信度。但是由于裂谷期存在断层作用, 地壳上表面的变化除了垂向变化还有水平方向的变化, 是两者共同作用的结果, 而一维Airy均衡回剥构造沉降恢复方法并不能有效恢复断层作用导致的水平方向上的变化。平衡剖面技术在恢复构造演化过程中的褶皱与断层作用造成的水平方向上的变化上具有优势, 可以逐步恢复每一期的地质状态^[18], 这为裂谷期伸展盆地或被动陆缘构造沉降恢复中水平方向的变化范围提供了约束。

二维挠曲回剥模型可以对伸展盆地的裂后期地层进行恢复, 该模型假设岩石圈各层在整个伸展过程中是完全耦合的, 各层拉张减薄程度均相同。如果研究区存在与深度相关的伸展^[32], 则模型不适用, 而被动陆缘的岩石圈在伸展过程中大多数是部分解耦的^[41], 所以在研究伸展量较大的被动陆缘时该模型不适用, 而仅适用于岩石圈较冷或未发生解耦区域的伸展盆地。但同时该软件具有层拉平功能可在一定程度弥补这一劣势, 层拉平指的是可将剖面最顶部地层界面拉到统一水平面上, 这在一定意义上意味着该沉积期的沉积物供应充足, 将可容空间完全充满。在陆上伸展盆地或一直处于浅水区的沉积盆地研究中, 这一操作可以较好地恢复盆地状态, 但不适用深水盆地。此外, 层拉平技术的可行性与不整合的类型也有关。若不整合形成时是水平的则可以很好地恢复^[42], 而其他类型不整合具有大规模的侵蚀, 导致不整合界面形成时非水平状态^[42], 而强行拉平则会改变真实的断层与下伏地层的产状。该软件中的地层去压实以恢复地层厚度的算法适用于所有的地质情况, 并可以进一步计算沉积速率, 而回剥拉平恢复地质状态则仅适用于陆上裂谷盆地或水深较小的沉积盆地, 要求未发生岩石圈各层之间的解耦并且仅适用于裂后期的盆地状态恢复。平衡剖面技术同样可以进行挠曲均衡的恢复并且在有效的古水深数据的约束下可以有效的恢复盆地的演化史。

挠曲悬臂梁模型同样假定岩石圈各层伸展减薄程度相同, 并且认为断层是平直断层, 但是地震剖面显示大型的断层多呈铲式形状, 最新的挠曲悬臂梁软件添加了铲式断层的模型。在其他初始模型参数确定的情况下, 该模型计算拉张因子的过程仅与输入的断层参数有关, 而脆性断层在大多数的伸展盆地或被动陆缘中仅能代表上地壳的伸展减薄程度。进一步将平衡剖面恢复得到的每一期盆地剖面上的断层产状输入Stretch软件中利用挠曲悬臂梁模型, 可以得到每一期的上地壳拉张因子。将模拟结果的断层影响范围与真实的地质观测进行对比, 可以得到上地壳的近似有效弹性厚度, 而未发生解耦的区域可以代表整个岩石圈的有效弹性厚度^[43]。

综上, 在存在岩石圈各层解耦的情况下, 挠曲回剥与挠曲悬臂梁模型均难以恢复剖面状态, 而平衡剖面在古水深的约束下可以对盆地的剖面进行较好的恢复, 并且恢复了断层作用造成的变形。本研究探讨的优势互补的软件组合使用方法, 以平衡剖面技术恢复的剖面形态为基础, 利用挠曲悬臂梁模型, 可以正演得到脆性上地壳的拉张因子, 利用构造沉降恢复可以得到各个地质时期构造沉降曲线长度的变化与差异, 而挠曲回剥软件中的地层厚度恢复可以有效得到每个沉积期的原始地层厚度, 从而得到其沉积速率信息。典型剖面利用体积守恒原则可以得到剖面处的平均伸展速度, 进一步在热流数据的约束下对裂谷盆地和被动陆缘进行数值模拟的研究, 可以开展岩石圈伸展的动力学研究。结合沉积速率变化、运动学特征与动力学特征可以有效地确定“生储盖”组合以指导油气的勘探与开发(图 1)。

4 结语及展望

反演可恢复盆地演化的过程, 得到盆地的初始地质状态, 而正演则可以在地质观测和反演结果的约束下, 探讨成盆机制与动力学问题。伸展盆地或被动陆缘的研究方法与商业软件众多, 但其各有侧重, 本次研究论述了目前应用较多的盆地模拟软件的原理及其应用, 其中盆地模拟的反演方法包括一维Airy均衡的构造沉降计算方法、二维挠曲回剥方法、平衡剖面技术。正演方法包括挠曲悬臂梁模型和二维粘弹塑性热-力学数值模拟模型。本研究分析了不同方法及软件的优势与不足, 提出优势互补的软件组合使用方法。在岩石圈较弱情况下一维Airy均衡可以较好的进行裂后的回剥, 而岩石圈较强时则应考虑侧向的影响即挠曲回剥方法较为合适, 二者均不具备恢复裂谷期地质状态的功能。挠曲悬臂梁模型在岩石圈不发生解耦的情况下可以很好的计算岩石圈的拉张因子, 若发生解耦, 则仅能表征上地壳的拉张减薄情况。平衡剖面技术在古水深的约束下可以较好的恢复裂谷期断层状态, 其可以为构造沉降计算提供水平方向的范围约束, 同时结合挠曲悬臂梁模型可以很好地计算出每一个地质历史时期时的伸展减薄情况。二维挠曲回剥中的地层厚度恢复可以用于计算沉积速率。在以上计算结果的约束下结合热流数据通过二维粘弹塑性热-力学模型模拟盆地或被动陆缘的形成过程。综合得到的沉积速率、运动学特征、动力学特征可以为进一步为油气资源的勘探与开发提供指导。

盆地形成演化是一个复杂的地质过程, 在不同尺度具有不同的表现, 利用在不同尺度下观测得到的结果相互验证可以更好的认识这一过程并了解其机制。现今的研究大多停留在盆地的浅层构造而非整体的地质构造过程, 充分利用各个尺度方法的优势并尽量避免其存在的问题才可以更加有效地约束演化模型进而得到更加可靠的盆地演化史, 也是未来的发展方向。在获得更加准确的沉积、沉降史的基础上可以对油气的勘探具有一定的指导意义。