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

1D应变速率反演的A MATLAB程序Hai-Bin Song，Lin Chen, Jiong Zhang, Chang Yu Zhao , Chong-Zhi Dong1.导言 过去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.物理模型 我们所考虑的2D物理模型如图1所示（其细节参见Jarvis and McKenzie (1980)资料）。地壳和岩石圈地幔以水平速度u(x)进行伸展，与此同时，软流圈地幔向上流动穿过平面z=0，进而取代向外流出的岩石圈物质。上部表面即z=a,这里的a是岩石圈厚度，上部表面温度为T=0℃；平面z=0处的温度保持为T=T1,也就是软流圈的温度。这种边界条件显然是大致的，因为岩石圈的顶面伴随着伸展作用的进行一定会发生沉降，并且也会发生沉积。变形作用限定为纯剪应变方式，并且垂直速度在z=0处为G(t)a,在这里G(t)为垂直应变速率，它是时间t的函数。垂直速度要求在z=a处消失，因此可表达为：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)。表1 所有计算中共用参数的定义与取值 直接的反演需要具有差异的一系列数据，因而易受到误差的影响。对于分散和具有误差的数据，反演问题的最好求解方法是，通过计算大量的正演模型，每次仅改变G(t)，直到计算的沉降值与观测的沉降值之间的误差值最小化为止。我们使用鲍威尔算法（Powells algorithm）(Press et al., 1992)，来系统性的改变G(t)，以便使检验函数