D應(yīng)變速率反演的MATLAB程序

1、1D應(yīng)變速率反演的A MATLAB程序Hai-Bin Song，Lin Chen, Jiong Zhang, Chang Yu Zhao , Chong-Zhi Dong1.導(dǎo)言 過去30多年來，裂谷和被動(dòng)大陸邊緣的演化一直受到普遍關(guān)注。構(gòu)造上活動(dòng)的裂谷、古裂谷和被動(dòng)大陸邊緣在不同類型的沉積盆地中形成了一類與伸展有關(guān)的重要盆地組合(Bally and Snelson, 1980)。伸展盆地覆蓋了全球相當(dāng)大的面積，并且蘊(yùn)藏了重要的沉積礦產(chǎn)和能源資源。大量的重要含油氣省與裂谷和被動(dòng)大陸邊緣有關(guān)(White et al., 2003;Ziegle and Cloetingh,2004)。為了進(jìn)一步認(rèn)

2、識(shí)裂谷盆地的演化，一些學(xué)者已提出了幾種模式(McKenzie, 1978; Royden and Keen,1980; Hellinger and Sclater, 1983; Rowley and Sahagian, 1986)。在模擬沉降作用時(shí)，伸展作用一般被假定為瞬時(shí)的(McKenzie, 1978)，或者是以一定的速率在有限的時(shí)期內(nèi)發(fā)生（Jarvis and McKenzie, 1980)。White(1993, 1994)提出了一種利用沉降資料反演的方法，該方法現(xiàn)已應(yīng)用于許多不同的盆地(Newman and White, 1999; Xie et al., 2006)。在本文中，我們

3、介紹的研究成果，主要是針對(duì)已有模型的重要進(jìn)展方面，為研究者研究裂谷盆地的演化提供有用的工具。該項(xiàng)工作的目的是提供一個(gè)開放的MATLAB程序，以用于根據(jù)沉降資料反演應(yīng)變速率和估算伸展系數(shù)。圖1 巖石圈隨時(shí)間伸展的物理模型示意圖。a=巖石圈厚度；T1=軟流圈溫度；u(x)=水平速度； a(z)=垂直速度。2.物理模型 我們所考慮的2D物理模型如圖1所示（其細(xì)節(jié)參見Jarvis and McKenzie (1980)資料）。地殼和巖石圈地幔以水平速度u(x)進(jìn)行伸展，與此同時(shí)，軟流圈地幔向上流動(dòng)穿過平面z=0，進(jìn)而取代向外流出的巖石圈物質(zhì)。上部表面即z=a,這里的a是巖石圈厚度，上部表面溫度為T=0

4、；平面z=0處的溫度保持為T=T1,也就是軟流圈的溫度。這種邊界條件顯然是大致的，因?yàn)閹r石圈的頂面伴隨著伸展作用的進(jìn)行一定會(huì)發(fā)生沉降，并且也會(huì)發(fā)生沉積。變形作用限定為純剪應(yīng)變方式，并且垂直速度在z=0處為G(t)a,在這里G(t)為垂直應(yīng)變速率，它是時(shí)間t的函數(shù)。垂直速度要求在z=a處消失，因此可表達(dá)為：V(z, t) （az)G(t) （1） 巖石圈隨時(shí)間變化的溫度結(jié)構(gòu)可以通過求解帶有附加對(duì)流項(xiàng)的1D熱流方程來獲得： 式中T 是溫度, t是時(shí)間, 是熱擴(kuò)散系數(shù)。物質(zhì)的平衡要求：u(x,t)=xG(t) (3) 在均勻伸展模型中并沒有包含應(yīng)變速率，因?yàn)樵谠撃Ｐ椭屑俣ㄉ煺棺饔檬撬矔r(shí)發(fā)生的。在有限

5、伸展模型中（Jarvis and McKenzie，1980），假定在巖石圈伸展作用期間應(yīng)變速率保持為一個(gè)常數(shù)(Wooler et al., 1992)。實(shí)際上，應(yīng)變速率是隨時(shí)間變化的。 理論上的沉降曲線一般是伸展系數(shù)（即總應(yīng)變）的函數(shù)。這些理論上的沉降曲線通過與實(shí)際觀測(cè)的沉降擬合可用來確定伸展系數(shù)，它與垂直應(yīng)變速率有關(guān)，即： 這里的t是伸展作用的持續(xù)時(shí)間。 如果G(t)為一常數(shù)，那么，估算的累計(jì)值為： 則理論上的水載盆地的沉降量的一般表達(dá)式為(White, 1994)： 其中： T(z,t) 是巖石圈的溫度，它是深度和時(shí)間的函數(shù)；T(z,)是巖石圈的平衡溫度結(jié)構(gòu)。其它參數(shù)的符號(hào)和數(shù)值見表1

6、(White,1994)。Q(t)是熱擾動(dòng)后的溫度結(jié)構(gòu)與平衡溫度結(jié)構(gòu)之間的差值，它是應(yīng)變速率G(t)的函數(shù)。A和B分別是地殼和巖石圈地幔的減薄因子。 根據(jù)已知的應(yīng)變速率變化來確定隨時(shí)間變化的沉降量被定義為正演問題。 反演問題可直接通過整理等式（6）獲得解決(White, 1994)，即： 給定S(t)，可以用迭代法解等式（10）得出G(t)，在迭代過程中，開始假設(shè)BQ(t)=0,即隨溫度而變化的密度變化最初被忽略，從而獲得G(t)的初始分布。然后用有限差分法，在隨后的多次迭代過程中，反復(fù)計(jì)算G(t)。表1 所有計(jì)算中共用參數(shù)的定義與取值 直接的反演需要具有差異的一系列數(shù)據(jù)，因而易受到誤差的影響

7、。對(duì)于分散和具有誤差的數(shù)據(jù)，反演問題的最好求解方法是，通過計(jì)算大量的正演模型，每次僅改變G(t)，直到計(jì)算的沉降值與觀測(cè)的沉降值之間的誤差值最小化為止。我們使用鮑威爾算法（Powell's algorithm）(Press et al., 1992)，來系統(tǒng)性的改變G(t)，以便使檢驗(yàn)函數(shù)H最小化(White, 1994)： 式中，Sio和Sic分別是觀測(cè)的沉降和計(jì)算的沉降；N是觀測(cè)沉降值的個(gè)數(shù)；M是G(t)的離散值的個(gè)數(shù)；W1、W2和W3是權(quán)重系數(shù)；G"k是估算的G(t)的二階導(dǎo)數(shù)，由三次樣條插值得出。當(dāng)計(jì)算的和觀測(cè)的所有沉降值均一致時(shí)，等式（11）的右側(cè)第一項(xiàng)為0。i為觀

8、測(cè)沉降值之間的深度間距(m)，用i除以它們之間的差值，使得求和中的每項(xiàng)成為單位方差。等式（11）中的最后3項(xiàng)是反演結(jié)果的偏差修正項(xiàng)。其中第2項(xiàng)和第3項(xiàng)可使G(t)曲線變得光滑；第4項(xiàng)隨著Gk接近于0而傾向于無窮。我們發(fā)現(xiàn)，當(dāng)權(quán)重系數(shù)取W1=0.5、W2=0.5和W3=0.05時(shí)，可以獲得比較滿意的結(jié)果。 3數(shù)值模型 Jarvis and Mckenzie (1980)用于解等式（2）的方法是一種偽解析法（pseudo -analytic method），它包含變量分離、本征函數(shù)展開式和數(shù)值積分方法，就如同四階Runge-Ikutta方法。盡管Jarvis et al. (1980)想用解析法解

9、等式（2），但他們的解是建立在數(shù)值方法之上的。由于很難找到用解析法解等式（2）的方法，我們采用有限差分法來解等式（2）。使用Crank-Nicolson的隱式有限差分格式(Lu and Guan, 2004)，等式（2）可以寫成： 式中，和h分別是時(shí)間和空間的步長(zhǎng)。設(shè)： 那么等式（12）可以寫成如下：給定初始的溫度結(jié)構(gòu)： 以及邊界條件： 等式（13）可用Thomas的代數(shù)概型法求解(Press et al., 1992)。 對(duì)于正演模擬，G(t)已被給定。因此，獲得了T(t,z)，我們就可以用等式（6）來計(jì)算沉降量，用等式（4）來計(jì)算伸展系數(shù)。 給定觀測(cè)的沉降量，并假設(shè)初始應(yīng)變速率g0(t)（

10、即：設(shè)Q(t)等于0，用等式（10）確定g0(t)），反演過程包括大量的正演模型計(jì)算，每次僅改變G(t)，直到計(jì)算的沉降值與觀測(cè)的沉降值之間的差最小化。反演的流程如圖2所示。 圖2 據(jù)沉降資料反演應(yīng)變速率的流程圖。4.假想數(shù)據(jù)將如上所述的算法應(yīng)用于一套由正演模擬所產(chǎn)生的假想數(shù)據(jù)。如下等式被用于產(chǎn)生應(yīng)變速率的分布型式，即應(yīng)變隨時(shí)間呈指數(shù)減小(White, 1994)： 式中，表示G(t)隨時(shí)間的增加而減小的快慢程度系數(shù)；G0是常數(shù)，由下式確定： 在這里，t是伸展作用的延續(xù)時(shí)間；是伸展系數(shù)。與此相似，如果G(t)隨時(shí)間呈現(xiàn)指數(shù)增加，那么： 正演模擬的應(yīng)變速率分布結(jié)果如圖3、圖4和圖5所示，它說明了

11、應(yīng)變速率G(t)的變化制約著沉降量和熱流值的變化，這里給定的總應(yīng)變即伸展系數(shù)是一定的。我們可以看到，如果在裂谷作用期間G(t)是一個(gè)常數(shù)，那么，沉降量和熱流值在裂谷作用期間對(duì)時(shí)間的導(dǎo)數(shù)大致也一個(gè)常數(shù)。當(dāng)裂谷作用停止時(shí)（即G=0），沉降量和熱流值隨時(shí)間以指數(shù)形式衰減。在第2種情況下，G(t)隨時(shí)間呈現(xiàn)指數(shù)減小，同裂谷期沉降也隨時(shí)間逐漸減小。但是，在裂谷期結(jié)束時(shí)，沉降量一定是比G(t)為常數(shù)時(shí)的沉降量要更大一些，而它們的總應(yīng)變是一定的。最后的第3種情況，如果G(t)隨時(shí)間呈現(xiàn)指數(shù)增加，同裂谷期沉降也隨時(shí)間逐漸增加，但總的同裂谷期沉降要比前兩種情況要更小一些。圖5所示的熱流變化具有與沉降量變化相似的

12、結(jié)果。沉降曲線和熱流曲線中的差異是裂谷作用期間應(yīng)變速率的不同分布所造成的結(jié)果。一旦受擾動(dòng)的溫度回歸到平衡狀態(tài)，那么，所有3種情況下的總的沉降速率都是相同的。圖3 三種應(yīng)變速率剖面，都從OMa開始，在40Ma結(jié)束，因此定義了相同的裂谷作用持續(xù)時(shí)間。在所有三種情況中，伸展作用的持續(xù)時(shí)間和總的應(yīng)變是相同的。 (a)應(yīng)變速率的分布以算術(shù)值顯示；(b)應(yīng)變速率剖面以對(duì)數(shù)值顯示。圖4 對(duì)應(yīng)于圖3中三種應(yīng)變速率分布的沉降曲線圖5 對(duì)應(yīng)于圖3中三種應(yīng)變速率分布的地表熱流變化 如上所述的反演算法對(duì)一套假想的數(shù)據(jù)進(jìn)行了應(yīng)用，這套數(shù)據(jù)是以裂谷作用期間G(t)隨時(shí)間呈指數(shù)減小的形式進(jìn)行正演模擬所產(chǎn)生的。如圖6所示，反

13、演的應(yīng)變速率在真實(shí)的值附近擺動(dòng)。計(jì)算的沉降值與理論上假設(shè)的沉降值十分吻合，除了邊界和中斷點(diǎn)以外，如圖7所示，在這里的標(biāo)準(zhǔn)偏差對(duì)所有觀測(cè)數(shù)據(jù)點(diǎn)均設(shè)定為25 m。盡管在真實(shí)值附近有較小的擺動(dòng)，但預(yù)測(cè)的沉降值與觀測(cè)的沉降值相當(dāng)吻合。如裂谷作用的延續(xù)時(shí)間已知，那么，反演的應(yīng)變速率就可以用等式（16）來估算伸展系數(shù)。根據(jù)假想數(shù)據(jù)計(jì)算的伸展系數(shù)和理論上的伸展系數(shù)變化如圖8所示。其結(jié)果表明，裂谷作用期間伸展系數(shù)的變化取決于時(shí)間。圖6 呈指數(shù)減小型的假想應(yīng)變數(shù)據(jù)所反演得到的應(yīng)變速率圖7 觀測(cè)的沉降與計(jì)算的沉降之間的關(guān)系圖8 反演計(jì)算得到的和理論上給出的真實(shí)應(yīng)變速率的累積伸展系數(shù)變化對(duì)比5.對(duì)南中國(guó)海的應(yīng)用 南

14、中國(guó)海(SCS)是西太平洋地區(qū)最大的邊緣海之一，它位于歐亞板塊、太平洋板塊和印度-澳大利亞板塊的接合部位。它是由新生代大陸裂解、海底擴(kuò)張作用形成的(Taylor and Hayes,1980; Briais et a1.,1993)。南中國(guó)海(SCS)的北部大陸邊緣經(jīng)歷了多期裂谷作用事件(Ru and Pigott,1986; Clift and Lin, 2001;He et al., 2002)。我們應(yīng)用自己研制的Matlab程序，利用鉆井獲得構(gòu)造沉降資料，反演了南中國(guó)海北部大陸邊緣的應(yīng)變速率，并估算了伸展系數(shù)。實(shí)際數(shù)據(jù)資料取自于Xieet al. (2006)。 圖9顯示了WC1411井

15、的反演應(yīng)變速率，該井位于南中國(guó)海北部大陸邊緣的珠江口盆地(PRMB)。沉降史資料來自Xie et al. (2006)。對(duì)于所有的觀測(cè)數(shù)據(jù)的標(biāo)準(zhǔn)偏差設(shè)定為25 m。圖10顯示了WC1411井的累積伸展系數(shù)，它是應(yīng)變速率的積分的指數(shù)函數(shù)。應(yīng)變速率的分布見圖9(a)，它揭示了該地區(qū)的多期裂谷作用事件。應(yīng)變速率可幫助我們了解裂谷作用過程。在WC1411井的應(yīng)變速率反演過程中，根據(jù)Yaoet al. (1994)的資料，將地殼厚度和巖石圈厚度分別設(shè)定為30km和120km，而其它所使用的參數(shù)與表1中所示相同。圖9 WC1411井的反演應(yīng)變速率。(a) WC1411井的反演應(yīng)變速率分布；(b)觀測(cè)的沉降（黑實(shí)心點(diǎn)）和預(yù)測(cè)的沉降（灰色實(shí)線）。圖10 據(jù)WC1411井的反演應(yīng)變速率估算的伸展系數(shù)。(a) WC1411井的反演應(yīng)變速率分布；(b) 由應(yīng)變速率累積推導(dǎo)出的WC1411井的累積伸展系數(shù)。6.結(jié)論 MATLAB程序的設(shè)計(jì)，主要是通過擬合伸展沉積盆地的觀測(cè)沉降史曲線，來反演隨時(shí)間變化的應(yīng)變速率的變化。正如White (1994)所指出的，該方法不需要事先提供有關(guān)裂谷作用的持續(xù)時(shí)間或者是總的應(yīng)變（即伸展系數(shù)）。用一套正演形成的假想數(shù)據(jù)對(duì)該方法進(jìn)行了測(cè)試。計(jì)算得到的沉降量與觀測(cè)的（或假想的）沉降量相當(dāng)吻合。盡管反演的應(yīng)變速率在真實(shí)值附近有所擺動(dòng)，但反映的結(jié)果反映了應(yīng)變速率變化的基本趨勢(shì)