1D应变速率反演的MATLAB程序

1D1D 应变速率反演的应变速率反演的 A MATLABA MATLAB程序程序 Hai-Bin Song，Lin Chen, Jiong Zhang, Chang Yu Zhao , Chong-Zhi Dong 1. 1.导言导言 过去 30 多年来，裂谷和被动大陆边缘的演化一直受到普遍关注。构造上活动的裂谷、 古裂谷和被动大陆边缘在不同类型的沉积盆地中形成了一类与伸展有关的重要盆地组合 (Bally and Snelson, 1980)。伸展盆地覆盖了全球相当大的面积，并且蕴藏了重要的沉积矿产 和能源资源。大量的重要含油气省与裂谷和被动大陆边缘有关(White et al., 2003;Ziegle and Cloetingh,2004)。为了进一步认识裂谷盆地的演化，一些学者已提出了几种模式 (McKenzie, 1978; Royden and Keen,1980; Hellinger and Sclater, 1983; Rowley and Sahagian, 1986) 。 在模拟沉 降作用时， 伸展作用一般被假定为瞬时的(McKenzie, 1978)， 或者是以一定的速率在有限的时 期内发生（Jarvis and McKenzie, 1980)。White(1993, 1994)提出了一种利用沉降资料反演的方 法，该方法现已应用于许多不同的盆地(Newman and White, 1999; Xie et al., 2006)。 在本文中，我们介绍的研究成果， 主要是针对已有模型的重要进展方面， 为研究者研究 裂谷盆地的演化提供有用的工具。 该项工作的目的是提供一个开放的MATLAB程序，以用于 根据沉降资料反演应变速率和估算伸展系数。 图 1岩石圈随时间伸展的物理模型示意图。 a=岩石圈厚度；T1=软流圈温度；u(x)=水平速度； a(z)=垂直速度。 2. 2.物理模型物理模型 我们所考虑的2D 物理模型如图 1 所示 （其细节参见Jarvis andMcKenzie (1980)资料） 。 地壳和岩石圈地幔以水平速度u(x)进行伸展， 与此同时， 软流圈地幔向上流动穿过平面z=0， 进而取代向外流出的岩石圈物质。 上部表面即 z=a,这里的 a 是岩石圈厚度，上部表面温度为 T=0℃；平面 z=0 处的温度保持为 T=T1,也就是软流圈的温度。这种边界条件显然是大致的， 因为岩石圈的顶面伴随着伸展作用的进行一定会发生沉降， 并且也会发生沉积。 变形作用限 定为纯剪应变方式，并且垂直速度在 z=0 处为 G(t)a,在这里 G(t)为垂直应变速率，它是时间 t 的函数。垂直速度要求在z=a 处消失，因此可表达为： 1 V(z, t) ＝（a－z)G(t)（1） 岩石圈随时间变化的温度结构可以通过求解带有附加对流项的1D 热流方程来获得： 式中 T 是温度, t 是时间, κ 是热扩散系数。物质的平衡要求： u(x,t)=xG(t)(3) 在均匀伸展模型中并没有包含应变速率，因为在该模型中假定伸展作用是瞬时发生的。 在有限伸展模型中（Jarvis and McKenzie，1980） ，假定在岩石圈伸展作用期间应变速率保持 为一个常数(Wooler et al., 1992)。实际上，应变速率是随时间变化的。 理论上的沉降曲线一般是伸展系数β （即总应变）的函数。这些理论上的沉降曲线通过 与实际观测的沉降拟合可用来确定伸展系数β ，它与垂直应变速率有关，即： 这里的Δ t 是伸展作用的持续时间。 如果 G(t)为一常数，那么，估算的累计值为： 则理论上的水载盆地的沉降量的一般表达式为(White, 1994)： 其中： T(z,t) 是岩石圈的温度，它是深度和时间的函数；T(z,∞)是岩石圈的平衡温度结构。其 它参数的符号和数值见表1 (White,1994)。 Q(t)是热扰动后的温度结构与平衡温度结构之间的 差值，它是应变速率 G(t)的函数。A 和 B 分别是地壳和岩石圈地幔的减薄因子。 根据已知的应变速率变化来确定随时间变化的沉降量被定义为正演问题。 反演问题可直接通过整理等式（6）获得解决(White, 1994)，即： 给定 S(t)，可以用迭代法解等式（10）得出 G(t)，在迭代过程中，开始假设 BQ(t)=0,即 随温度而变化的密度变化最初被忽略，从而获得G(t)的初始分布。然后用有限差分法，在随 后的多次迭代过程中，反复计算G(t)。 2 表 1所有计算中共用参数的定义与取值 直接的反演需要具有差异的一系列数据， 因而易受到误差的影响。 对于分散和具有误差 的数据，反演问题的最好求解方法是，通过计算大量的正演模型，每次仅改变G(t)，直到计 算的沉降值与观测的沉降值之间的误差值最小化为止。我们使用鲍威尔算法（Powell s algorithm） (Press et al., 1992)， 来系统性的改变 G(t)， 以便使检验函数 H 最小化(White, 1994)： 式中，Sio和 Sic分别是观测的沉降和计算的沉降；N 是观测沉降值的个数；M 是 G(t)的 离散值的个数；W1、W2 和 W3 是权重系数；G“k 是估算的 G(t)的二阶导数，由三次样条插 值得出。当计算的和观测的所有沉降值均一致时，等式（11）的右侧第一项为 0。σ i 为观 测沉降值之间的深度间距(m)， 用σ i 除以它们之间的差值， 使得求和中的每项成为单位方差。 等式（11）中的最后3 项是反演结果的偏差修正项。其中第2 项和第 3 项可使 G(t)曲线变得 光滑；第 4 项随着 Gk 接近于 0 而倾向于无穷。我们发现，当权重系数取 W1=0.5、W2=0.5 和 W3=0.05 时，可以获得比较满意的结果。 3 3．数值模型．数值模型 Jarvis and Mckenzie (1980)用于解等式（2）的方法是一种伪解析法（pseudo -analytic ） ，它包含变量分离、本征函数展开式和数值积分方法，就如同四阶Runge-Ikutta 方 法。尽管 Jarvis et al. (1980)想用解析法解等式（2） ，但他们的解是建立在数值方法之上的。 由于很难找到用解析法解等式（ 2）的方法，我们采用有限差分法来解等式（2） 。使用 Crank-Nicolson的隐式有限差分格式(Lu and Guan, 2004)，等式（2）可以写成： 3 式中，τ和 h 分别是时间和空间的步长。设： 那么等式（12）可以写成如下： 给定初始的温度结构： 以及边界条件： 等式（13）可用 Thomas 的代数概型法求解(Press et al., 1992)。 对于正演模拟，G(t)已被给定。因此，获得了 T(t,z)，我们就可以用等式（6）来计算沉 降量，用等式（4）来计算伸展系数。 给定观测的沉降量，并假设初始应变速率g0(t)（即：设 Q(t)等于 0，用等式（10）确定 g0(t)） ，反演过程包括大量的正演模型计算，每次仅改变 G(t)，直到计算的沉降值与观测的 沉降值之间的差最小化。反演的流程如图2 所示。 图 2据沉降资料反演应变速率的流程