反演純剪切模型中β的算法研究

張 炯,宋海斌

(1.遼寧省冶金地質勘查局 地質勘查研究院,鞍山 114038;

2.中國科學院地質與地球物理研究所 油氣資源研究院重點實驗室,北京 100029)

反演純剪切模型中β的算法研究

張 炯1,宋海斌2

(1.遼寧省冶金地質勘查局 地質勘查研究院,鞍山 114038;

2.中國科學院地質與地球物理研究所 油氣資源研究院重點實驗室,北京 100029)

以往反演伸展因子的算法多是基于單井資料進行擬合計算得到,很少把構造與熱演化史結合研究。這里利用M cKenzie在1978年提出的純剪切模型中計算沉降和熱流的公式,基于一維反演伸展因子的思想,采用M atlab軟件GU I功能,編制多點反演伸展因子界面,大大減少在對偽井數據反演伸展因子時的工作量。以南海中北部SO49-25測線為例,計算了該剖面陸架和陸坡區的伸展因子和古熱流值,并驗證了算法的可行性。

純剪切模型;伸展因子;反演

0 前言

隨著大量資料積累及電子技術的進步和發展,通過地質過程定量分析或模擬研究,加深了對地質過程的認識,因此盆地定量模擬研究對解決地質問題的定量化提供了依據。其中地史和熱史的數值模擬最終目的,是為油氣生、排、運聚、成藏模擬創造條件。盆地構造～地史和熱演化史,是油氣資源評價和油氣成藏動力學研究的重要內容之一。

二十世紀七十年代以來,許多學者對巖石圈拉張、盆地和大陸邊緣的形成演化,進行了廣泛深入的研究,提出了多種地球動力學定量模型。其中以M cKenzie[4]提出的巖石圈純剪切模型最為經典,為以后定量研究盆地與張裂大陸邊緣形成演化機制奠定了基礎。數值模擬方法必須與觀測資料相結合,只有這樣模擬結果才具有地質意義,才能為解釋地質現象提供科學依據,換言之,數值模型需要以觀測數據作為約束條件[1～4]。

1 純剪切模型

M cKenzie[4]提出了巖石圈均勻拉張模型,又稱純剪切模型。該模型描述了巖石圈對拉伸作用的基本響應,建立了拉伸系數(β)和盆地熱流之間的定量關系,奠定了拉伸盆地定量模擬的理論基礎。在盆地模擬和沉積盆地熱史恢復領域中,具有劃時代意義。構造～熱演化模擬方法的發展,實際上即為拉伸盆地定量模擬的地質地球物理模型的發展。M cKenzie提出的定量模型屬于瞬時均勻純剪切模型,如下頁圖1所示,即:拉伸時間短(<20M a),拉伸期間的熱擴散忽略不計,地溫梯度瞬時增高,巖石圈變形機制為純剪切變形,拉伸過程中地殼和巖石圈地幔均勻變薄。

M cKenzie認為,裂陷盆地總體沉降由初始沉降和熱沉降二部份組成:①伸展期間由斷層控制的沉降,取決于地殼的初始厚度和伸展系數β值的大小;②伸展之后的熱沉降,由巖石圈冷卻所引起,取決于伸展量。并認為巖石圈伸展是均勻的,斷層控制的沉降是瞬間完成的,而熱沉降速率是隨時間呈指數衰減的。M cKenzie提出的純剪切模型的基本公式以及各個參數的意義為:

圖1 純剪切模型圖(M ckenzie,1978)Fig.1 Pure shearmodel(M cKenzie,1978)

(1)初始沉降公式:

(2)熱沉降公式:

(3)溫度異常采用傅里葉級數表達:

其中 an為實常數;z為深度。(4)熱流值公式:

以上公式(1)～式(5)中的參數意義:β是伸展因子;yL是初始的巖石圈厚度;yc是初始的地殼厚度;ρm＊是0℃時地幔的密度,通常取值3 330 kg/m3;ρ＊c是0℃時地殼的密度,通常取值2 800 kg/m3;ρs是沉積物或水的密度,通常取值1 030 kg/m3;αv是熱擴散系數,通常取值3.28×10-5℃;Tm是軟流圈的溫度,通常取值1 333℃;κ是熱膨脹系數,通常取值0.007 5 cm2/s。

研究分析M cKenzie均勻伸展模型的理論結果,認為:

(1)M cKenzie均勻伸展模型是一維瞬時拉張模型,其總沉積量由二部份組成:斷層控制的初始沉降和沉降速率隨指數衰減的熱沉降。熱沉降隨伸展因子β的增大而增大,β為無窮大,代表洋中脊的情況。

(2)伸展因子β為巖石圈伸展前厚度與伸展后厚度之比,因此β越大,巖石圈伸展減薄量越大,初始沉降也就越大,因此初始熱流異常也就越大。

(3)M cKenzie模型模擬的是地表熱流,巖石圈拉張瞬間完成時熱流達到最大值,隨即呈指數衰減。伸展因子越大,初始熱流越高,在經過200M a冷卻后,熱流趨于穩態。

通過了解M cKenzie均勻伸展模型的機制,可以得出在利用M cKenzie均勻伸展模型反演計算伸展因子的計算過程中,初始沉降量取決于地殼與巖石圈厚度之比(yc/yL)和伸展因子(β)。如果對于某個盆地yc/yL已知,則理論上由斷層控制的初始沉降可用來估算伸展量。由于其它因素造成了裂谷期厚度的變化,使得這種方法可信度降低。因此在擬合伸展因子的計算過程中,很少利用考慮與初始沉降的關系。因為M cKenzie模型預測的熱沉降只取決于伸展因子,影響裂后熱沉降的因素相對較少,所以通過M cKenzie模型計算熱沉降的公式,可以估算伸展因子。擬合方法是把裂后沉積厚度經過沉積物壓實校正,在沉積物負載校正,水深校正,以及重力均衡校正之后,得到構造沉降史。使得裂后沉降部份與理論模型預測的熱沉降曲線進行對比,得到最佳擬合的伸展因子。

因此作者利用M cKenzie理論模型一維擬合伸展因子的思想,編制了多點擬合伸展因子的帶界面的模擬程序。

2 伸展因子估算方法

A llen等[1]總結了多種估算拉伸因子的方法,作者在本文中討論的是利用沉降分析估算β。首先利用地震剖面和鉆井資料,通過回剝法(Steck ler和W atts[5])可以得到水載構造沉降,即觀測沉降。從純剪切模型的初始沉降公式中可以看到,初始沉降取決于初始地殼厚度與巖石圈厚度之比和拉伸因子。因此,在已知初始地殼與巖石圈厚度的情況下,理論上可以通過斷層控制的初始沉降估算拉伸因子。但由于裂谷期厚度的劇烈變化,使得估算結果變得不可靠。常用的方法是給定不同的拉伸因子β,利用純剪切模型的熱沉降公式,計算理論的水載熱沉降曲線,通過比較模型預測的沉降與觀測沉降,以確定最佳擬合的β。

基于均勻伸展模型正演模擬研究中的計算公式,編輯擬合伸展因子的反演算法,步驟如下:

(1)已知m個時期的基底的埋藏史數據,并對數據進行轉換,得到構造沉降史數據,其中包括裂谷期沉降(Ys)和裂后期沉降(Yt(t))。

(2)已知裂后期開始的時間,確定裂后期熱沉降史數據有n個,從構造沉降史數據中減去裂谷期沉降(Ys)的貢獻,即得到裂后期熱沉降量(Yt(t)-Ys)。

(3)基于M cKenzie均勻伸展模型中的計算熱沉降量的公式：),以β=1.1為初始值,計算yt(t)。然后通過迭代運算,使得為最小,即得到擬合因子β,其中err是單個構造沉降量的擬合計算的計算誤差,單位為“m”。

(4)當β>6.1時,伸展因子取6.1。

3 多點伸展因子擬合界面研究

采用m atlab語言進行編程,編制反演多點伸展因子擬合的程序,并利用m atlab軟件的GU I相關功能,實現PC機可視化編程,編輯獨立的計算機界面,用于模型的數值模擬工作。在PC機中有m atlab動態鏈接庫的前提下,可以實現程序脫離m atlab軟件環境運行[6～9]。

圖2為主界面,“顯示路徑”鍵:可以顯示參數輸入的路徑;“運行”鍵:即輸出運算結果。

圖3為時間參數文件界面,在“tim e.par”文件中,第一行中“8”代表本次參與運算的地層時間有八個;第二行分別顯示各個時間點;第三行中“23.3”代表裂后熱沉降開始的時間。

圖4為基本參數文件界面,在“param eters.par”文件中,“125000”代表原始巖石圈厚度,單位為“m”;“31200”代表原始地殼厚度,單位為“m”;“3330”代表原始地幔密度,單位為“kg/m3”;“2800”代表原始地殼密度,單位為“kg/m3”;“1030”代表海水密度,單位為“kg/m3”;“3”代表熱導率值,單位“J/C/m/s”;“3.28E-5”代表熱擴散系數,單位“℃”;“1330”代表軟流圈頂界溫度,單位為“℃”;“1.0E-6”代表熱膨脹系數,單位“cm2/s”。

圖2 設置參數路徑Fig.2 Set the param eterspath

圖3 時間參數文件Fig.3 The file of tim e param eters

圖4 計算基本參數文件Fig.4 The file of calculating the elem entary param eters

圖5(見下頁)為沉降數據參數文件界面,在“subsidence.par”文件中,第一行數據中“8”代表有八個層位數據,單位“m”,“120”代表有一百二十個點數據。

圖6(見下頁)為密度參數文件界面,在“rho_Sed.par”中,第一行數據“8”對應八個沉積物平均密度,單位“kg/m3”,“120”代表有一百二十個點。

圖7和圖8為輸出文件界面,其中“heat_flux.dat”文件和“beta_inv.dat”文件是輸出結果文件。“heat_flux.dat”文件中,第一行數據“120”是代表本次參與擬合伸展因子的點數為一百二十個,“6”代表是輸出六個時間點的熱流;“beta_inv.dat”文件中,第一行數據“120”是代表本次參與擬合伸展因子的結果有一百二十個;從第二行數據開始,第一列數據顯示為伸展因子擬合結果,第二列數據顯示為計算擬合誤差。

圖5 沉降數據參數文件Fig.5 The file of subsidence data

圖6 沉積物平均密度參數文件Fig.6 The file of deposition average density

4 SO 49-25測線實例研究

作者利用在2D-M ove軟件中分別輸出各個時期的基底沉降,已經完成地震剖面的埋藏史恢復工作,只需提取垂直方向各個時期的沉降數據,對沉降數據做A iry均衡校正,轉化為基底的構造沉降。根據有效的多點伸展因子擬合的算法,計算不同區域的伸展因子的變化,并計算當中存在的誤差分析。圖9(見下頁)中顯示的計算誤差,是擬合伸展因子過程中的最小err[5]。

圖7 輸出的熱流結果文件Fig.7 The file of exporting the heat flow resu lt

圖8 輸出的伸展因子及擬合誤差文件Fig.8 The file of exporting the stretching factor and fitting error

25測線北部陸架部份的伸展因子擬合的結果顯示(見下頁圖9),在測線0 km～20 km處伸展量相對較小,伸展因子在1.7～2.5之間,而且計算誤差多在20～90之間,分析認為是因該區域位于珠二坳陷,裂后沉降較大,這說明伸展因子在擬合過程中有一定的合理性。在20 km～80 km處伸展量相對較大,伸展因子變化較小導致,在1.3～2.0之間,而且計算誤差在20～110之間,分析認為是因該區域位于一統暗沙隆起區,相對裂后沉降較小導致。在計算西沙海槽伸展因子結果時,因發生裂后沉降結果較大,導致伸展因子結果在2.0～6.1之間變化較大,且計算誤差在50～250之間變化,導致結果可信度不高。

SO49-25測線北部陸坡區熱流分析(見后面圖10):測線穿過北部陸架～西北海槽～北部陸坡區,西沙海槽和陸坡區熱流值較大,北部陸架區0 km～30 km熱流值總體呈增大趨勢;在30 km～40 km處,熱流有減小趨勢,分析認為該區域位于一統暗沙隆起區,使得熱流值較小;在40 km～80 km處,總體為增大趨勢,符合靠近洋盆熱流增大的趨勢。

圖9 SO49～25北部陸架區和西沙海槽～陸坡區伸展因子分布圖Fig.9 The stretching factor distribution in the northern continental shelf and X isha trough area of the SO49-25 line

對各個時期的熱流值分布的認識:

(1)23.8M a屬于裂后期開始的時間,熱流值最大,可達到80mW/m2,均在60mW/m2以上。

(2)16.5 M a熱流逐漸減小,最大值為70mW/m2,均值在58mW/m2以上。

(3)10.5 M a熱流繼續減小,最大值為68mW/m2,均在55mW/m2以上。

(4)現今熱流仍然在減小,最大值為60mW/m2,均在50mW/m2以上。

在以上熱流計算中,暫時沒有考慮放射性元素的影響。熱流Q=qh+q+Δq,q代表地殼淺層、深地殼和地幔熱源的熱流;Δq代表異常熱流;qh是來自放射性衰減的熱流。

自23.8M a以來,測線所在研究區處于熱沉降階段,熱流從裂陷谷期末期逐漸下降。

根據模擬計算,17測線陸坡區現今熱流在50mW/m2～62 mW/m2之間,加上放射性生熱貢獻(簡單取28 mW/m2),在78 mW/m2～90mW/m2之間;18測線陸架陸坡部份現今熱流主要在40mW/m2～60mW/m2之間,加上放射性生熱貢獻,在68 mW/m2～88 mW/m2之間;25測線陸坡部份在40mW/m2～58 mW/m2之間,加上放射性生熱貢獻,在68mW/m2～86mW/m2之間。自東向西,拉伸因子總體減小,熱流相應降低。

5 結論

作者在本文重點研究了M atlab軟件的GU I功能,結合M cKenzie提出的純剪切模型,主要完成以下幾點工作:①編制多點反演伸展因子的帶界面的程序,減少了在對偽井數據反演伸展因子時的工作量;②以南海中北部SO49-25測線為例,計算了該剖面陸架和陸坡區的伸展因子和古熱流值,把構造史和熱演化史有效地結合起來,為下一步討論油氣資源的分布提供參考。

致謝 作者在本文所用SO航次的地震數據是由德國BGR的Franke博士提供,在此向Franke博士表示感謝。

[1]ALLEN P A,ALLEN J R.Basin analysis,p rincip les and app lications[M].Oxford:B lackwell Scientific Publications,1990.

[2]JARV ISG T,MCKENZIED P.Sedim entary basin form ationw ith finite extension rates[J].Earth and Planetary Science Letters,1980,48(1):42.

[3]KEEN C E,DEHLER SA.Stretching and subsidence:R ifting of con jugatem argins in theNorth A tlantic region[J].Tectonics,1993,12:1219.

[4]MCKENZIE D.Som e rem arks on the developm ent of sedim entary basins[J].Earth and Planetary Science Letters,1978,40(1):25.

[5]STECKLERM S,WATTSA B.Subsidence of the atlan tic-type continentalm argin off New Yo rk[J].Earth and p lanetary sciences letters,1978,41:1.

[6]李顯宏.M atlab 7.x界面設計與編譯技巧[M].北京:電子工業出版社,2006.

[7]王正林,劉明.精通MATLAB7[M].北京:電子工業出版社,2006.

[8]楊長青,胥澤銀.基于MATLAB的面積計算方法[J].物探化探計算技術,2004,26(2):177.

[9]李曉昌.在MATLAB平臺上實現可控源音頻大地電磁反演數據三維可視化顯示[J].物探化探計算技術,2007,(增刊):69.

TU 432

A

1001—1749(2011)01—0101—06

國家基礎研究發展規劃項目資助(2007CB411704)

2010-08-09 改回日期:2010-11-02

張炯(1984-),男,碩士,從事應用地球物理研究工作。