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

11D应变速率反演的AMATLAB程序Hai-BinSong，LinChen,JiongZhang,ChangYuZhao,Chong-ZhiDong1.导言过去30多年来，裂谷和被动大陆边缘的演化一直受到普遍关注。构造上活动的裂谷、古裂谷和被动大陆边缘在不同类型的沉积盆地中形成了一类与伸展有关的重要盆地组合(BallyandSnelson,1980)。伸展盆地覆盖了全球相当大的面积，并且蕴藏了重要的沉积矿产和能源资源。大量的重要含油气省与裂谷和被动大陆边缘有关(Whiteetal.,2003;ZiegleandCloetingh,2004)。为了进一步认识裂谷盆地的演化，一些学者已提出了几种模式(McKenzie,1978;RoydenandKeen,1980;HellingerandSclater,1983;RowleyandSahagian,1986)。在模拟沉降作用时，伸展作用一般被假定为瞬时的(McKenzie,1978)，或者是以一定的速率在有限的时期内发生（JarvisandMcKenzie,1980)。White(1993,1994)提出了一种利用沉降资料反演的方法，该方法现已应用于许多不同的盆地(NewmanandWhite,1999;Xieetal.,2006)。在本文中，我们介绍的研究成果，主要是针对已有模型的重要进展方面，为研究者研究裂谷盆地的演化提供有用的工具。该项工作的目的是提供一个开放的MATLAB程序，以用于根据沉降资料反演应变速率和估算伸展系数。2.物理模型我们所考虑的2D物理模型如图1所示（其细节参见JarvisandMcKenzie(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岩石圈随时间伸展的物理模型示意图。a=岩石圈厚度；T1=软流圈温度；u(x)=水平速度；a(z)=垂直速度。2V(z,t)＝（a－z)G(t)（1）岩石圈随时间变化的温度结构可以通过求解带有附加对流项的1D热流方程来获得：式中T是温度,t是时间,κ是热扩散系数。物质的平衡要求：u(x,t)=xG(t)(3)在均匀伸展模型中并没有包含应变速率，因为在该模型中假定伸展作用是瞬时发生的。在有限伸展模型中（JarvisandMcKenzie，1980），假定在岩石圈伸展作用期间应变速率保持为一个常数(Wooleretal.,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)。3直接的反演需要具有差异的一系列数据，因而易受到误差的影响。对于分散和具有误差的数据，反演问题的最好求解方法是，通过计算大量的正演模型，每次仅改变G(t)，直到计算的沉降值与观测的沉降值之间的误差值最小化为止。我们使用鲍威尔算法（Powell'salgorithm）(Pressetal.,1992)，来系统性的改变G(t)，以便使检验函数H最小化(White,1994)：式中，Sio和Sic分别是观测的沉降和计算的沉降；N是观测沉降值的个数；M是G(t)的离散值的个数；W1、W2和W3是权重系数；Gk是估算的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．数值模型JarvisandMckenzie(1980)用于解等式（2）的方法是一种伪解析法（pseudo-analyticmethod），它包含变量分离、本征函数展开式和数值积分方法，就如同四阶Runge-Ikutta方法。尽管Jarvisetal.(1980)想用解析法解等式（2），但他们的解是建立在数值方法之上的。由于很难找到用解析法解等式（2）的方法，我们采用有限差分法来解等式（2）。使用Crank-Nicolson的隐式有限差分格式(LuandGuan,2004)，等式（2）可以写成：表1所有计算中共用参数的定义与取值4式中，τ和h分别是时间和空间的步长。设：那么等式（12）可以写成如下：给定初始的温度结构：以及边界条件：等式（13）可用Thomas的代数概型法求解(Pressetal.,1992)。对于正演模拟，G(t)已被给定。因此，获得了T(t,z)，我们就可以用等式（6）来计算沉降量，用等式（4）来计算伸展系数。给定观测的沉降量，并假设初始应变速率g0(t)（即：设Q(t)等于0，用等式（10）确定g0(t)），反演过程包括大量的正演模型计算，每次仅改变G(t)，直到计算的沉降值与观测的沉降值之间的差最小化。反演的流程如图2所示。图2据沉降资料反演应变速率的流程图。54.假想数据将如上所述的算法应用于一套由正演模拟所产生的假想数据。如下等式被用于产生应变速率的分布型式，即应变随时间呈指数减小(White,1994)：式中，τ表示G(t)随时间的增加而减小的快慢程度系数；G0是常数，由下式确定：在这里，Δt是伸展作用的延续时间；β是伸展系数。与此相似，如果G(t)随时间呈现指数增加，那么：正演模拟的应变速率分布结果如图3、图4和图5所示，它说明了应变速率G(t)的变化制约着沉降量和热流值的变化，这里给定的总应变即伸展系数β是一定的。我们可以看到，如果在裂谷作用期间G(t)是一个常数，那么，沉降量和热流值在裂谷作用期间对时间的导数大致也一个常数。当裂谷作用停止时（即G=0），沉降量和热流值随时间以指数形式衰减。在第2种情况下，G(t)随时间呈现指数减小，同裂谷期沉降也随时间逐渐减小。但是，在裂谷期结束时，沉降量一定是比G(t)为常数时的沉降量要更大一些，而它们的总应变β是一定的。最后的第3种情况，如果G(t)随时间呈现指数增加，同裂谷期沉降也随时间逐渐增加，但总的同裂谷期沉降要比前两种情况要更小一些。图5所示的热流变化具有与沉降量变化相似的结果。沉降曲线和热流曲线中的差异是裂谷作用期间应变速率的不同分布所造成的结果。一旦受扰动的温度回归到平衡状态，那么，所有3种情况下的总的沉降速率都是相同的。图3三种应变速率剖面，都从OMa开始，在40Ma结束，因此定义了相同的裂谷作用持续时间。在所有三种情况中，伸展作用的持续时间和总的应变是相同的。(a)应变速率的分布以算术值显示；(b)应变速率剖面以对数值显示。6如上所述的反演算法对一套假想的数据进行了应用，这套数据是以裂谷作用期间G(t)随时间呈指数减小的形式进行正演模拟所产生的。如图6所示，反演的应变速率在真实的值附近摆动。计算的沉降值与理论上假设的沉降值十分吻合，除了边界和中断点以外，如图7所示，在这里的标准偏差对所有观测数据点均设定为25m。尽管在真实值附近有较小的摆动，但预测的沉降值与观测的沉降值相当吻合。如裂谷作用的延续时间已知，那么，反演的应变速率就可以用等式（16）来估算伸展系数。根据假想数据计算的伸展系数和理论上的伸展系数变化如图8所示。其结果表明，裂谷作用期间伸展系数的变化取决于时间。图4对应于图3中三种应变速率分布的沉降曲线图5对应于图3中三种应变速率分布的地表热流变化7图6呈指数减小型的假想应变数据所反演得到的应变速率图7观测的沉降与计算的沉降之间的关系图8反演计算得到的和理论上给出的真实应变速率的累积伸展系数变化对比85.对南中国海的应用南中国海(SCS)是西太平洋地区最大的边缘海之一，它位于欧亚板块、太平洋板块和印度-澳大利亚板块的接合部位。它是由新生代大陆裂解、海底扩张作用形成的(TaylorandHayes,1980;Briaiseta1.,1993)。南中国海(SCS)的北部大陆边缘经历了多期裂谷作用事件(RuandPigott,1986;CliftandLin,2001;Heetal.,2002)。我们应用自己研制的Matlab程序，利用钻井获得构造沉降资料，反演了南中国海北部大陆边缘的应变速率，并估算了伸展系数。实际数据资料取自于Xieetal.(2006)。图9显示了WC1411井的反演应变速率，该井位于南中国海北部大陆边缘的珠江口盆地(PRMB)。沉降史资料来自Xieetal.(2006)。对于所有的观测数据的标准偏差设定为25m。图10显示了WC1411井的累积伸展系数，它是应变速率的积分的指数函数。应变速率的分布见图9(a)，它揭示了该地区的多期裂谷作用事件。应变速率可帮助我们了解裂谷作用过程。在WC1411井的应变速率反演过程中，根据Yaoetal.(1994)的资料，将地壳厚度和岩石圈厚度分别设定为30km和120km，而其它所使用的参数与表1中所示相同。图9WC1411井的反演应变速率。(a)WC1411井的反演应变速率分布；(b)观测的沉降（黑实心点）和预测的沉降（灰色实线）。96.结论MATLAB程序的设计，主要是通过拟合伸展沉积盆地的观测沉降史曲线，来反演随时间变化的应变速率的变化。正如White(1994)所指出的，该方法不需要事先提供有关裂谷作用的持续时间或者是总的应变（即伸展系数）。用一套正演形成的假想数据对该方法进行了测试。计算得到的沉降量与观测的（或假想的）沉降量相当吻合。尽管反演的应变速率在真实值附近有所摆动，但反映的结果反映了应变速率变化的基本趋势。而且，反演的应变速率可以用来估算伸展系数。这些模拟所得到的有关应变速率在时间上的变化资料，为岩石圈变形的动力学模型提供了所需要的约束。应变速率不仅在时间上是变化的，而且在空间上也是变化的。最近，WhiteandBellingham(2002)将1D应变速率反演推进到了2D的情况；Shillingtonetal.(2008)描述了一种类似的程序算法，该算法反演2D随深度变化的应变速率。为了处理真正的实际问题，将来的发展要在该程序和White最近研究的基础上，实施2D的应变速率反演程序。图10据WC1411井的反演应变速率估算的伸展系数。(a)WC1411井的反演应变速率分布；(b)由应变速率累积推导出的WC1411井的累积伸展系数。