張 炯,宋海斌

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

2.中國(guó)科學(xué)院地質(zhì)與地球物理研究所 油氣資源研究院重點(diǎn)實(shí)驗(yàn)室,北京 100029)

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

張 炯1,宋海斌2

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

2.中國(guó)科學(xué)院地質(zhì)與地球物理研究所 油氣資源研究院重點(diǎn)實(shí)驗(yàn)室,北京 100029)

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

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

0 前言

隨著大量資料積累及電子技術(shù)的進(jìn)步和發(fā)展,通過(guò)地質(zhì)過(guò)程定量分析或模擬研究,加深了對(duì)地質(zhì)過(guò)程的認(rèn)識(shí),因此盆地定量模擬研究對(duì)解決地質(zhì)問(wèn)題的定量化提供了依據(jù)。其中地史和熱史的數(shù)值模擬最終目的,是為油氣生、排、運(yùn)聚、成藏模擬創(chuàng)造條件。盆地構(gòu)造～地史和熱演化史,是油氣資源評(píng)價(jià)和油氣成藏動(dòng)力學(xué)研究的重要內(nèi)容之一。

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

1 純剪切模型

M cKenzie[4]提出了巖石圈均勻拉張模型,又稱(chēng)純剪切模型。該模型描述了巖石圈對(duì)拉伸作用的基本響應(yīng),建立了拉伸系數(shù)(β)和盆地?zé)崃髦g的定量關(guān)系,奠定了拉伸盆地定量模擬的理論基礎(chǔ)。在盆地模擬和沉積盆地?zé)崾坊謴?fù)領(lǐng)域中,具有劃時(shí)代意義。構(gòu)造～熱演化模擬方法的發(fā)展,實(shí)際上即為拉伸盆地定量模擬的地質(zhì)地球物理模型的發(fā)展。M cKenzie提出的定量模型屬于瞬時(shí)均勻純剪切模型,如下頁(yè)圖1所示,即:拉伸時(shí)間短(<20M a),拉伸期間的熱擴(kuò)散忽略不計(jì),地溫梯度瞬時(shí)增高,巖石圈變形機(jī)制為純剪切變形,拉伸過(guò)程中地殼和巖石圈地幔均勻變薄。

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

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

(1)初始沉降公式:

(2)熱沉降公式:

(3)溫度異常采用傅里葉級(jí)數(shù)表達(dá):

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

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

研究分析M cKenzie均勻伸展模型的理論結(jié)果,認(rèn)為:

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

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

(3)M cKenzie模型模擬的是地表熱流,巖石圈拉張瞬間完成時(shí)熱流達(dá)到最大值,隨即呈指數(shù)衰減。伸展因子越大,初始熱流越高,在經(jīng)過(guò)200M a冷卻后,熱流趨于穩(wěn)態(tài)。

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

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

2 伸展因子估算方法

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

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

(1)已知m個(gè)時(shí)期的基底的埋藏史數(shù)據(jù),并對(duì)數(shù)據(jù)進(jìn)行轉(zhuǎn)換,得到構(gòu)造沉降史數(shù)據(jù),其中包括裂谷期沉降(Ys)和裂后期沉降(Yt(t))。

(2)已知裂后期開(kāi)始的時(shí)間,確定裂后期熱沉降史數(shù)據(jù)有n個(gè),從構(gòu)造沉降史數(shù)據(jù)中減去裂谷期沉降(Ys)的貢獻(xiàn),即得到裂后期熱沉降量(Yt(t)-Ys)。

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

(4)當(dāng)β>6.1時(shí),伸展因子取6.1。

3 多點(diǎn)伸展因子擬合界面研究

采用m atlab語(yǔ)言進(jìn)行編程,編制反演多點(diǎn)伸展因子擬合的程序,并利用m atlab軟件的GU I相關(guān)功能,實(shí)現(xiàn)PC機(jī)可視化編程,編輯獨(dú)立的計(jì)算機(jī)界面,用于模型的數(shù)值模擬工作。在PC機(jī)中有m atlab動(dòng)態(tài)鏈接庫(kù)的前提下,可以實(shí)現(xiàn)程序脫離m atlab軟件環(huán)境運(yùn)行[6～9]。

圖2為主界面,“顯示路徑”鍵:可以顯示參數(shù)輸入的路徑;“運(yùn)行”鍵:即輸出運(yùn)算結(jié)果。

圖3為時(shí)間參數(shù)文件界面,在“tim e.par”文件中,第一行中“8”代表本次參與運(yùn)算的地層時(shí)間有八個(gè);第二行分別顯示各個(gè)時(shí)間點(diǎn);第三行中“23.3”代表裂后熱沉降開(kāi)始的時(shí)間。

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

圖2 設(shè)置參數(shù)路徑Fig.2 Set the param eterspath

圖3 時(shí)間參數(shù)文件Fig.3 The file of tim e param eters

圖4 計(jì)算基本參數(shù)文件Fig.4 The file of calculating the elem entary param eters

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

圖6(見(jiàn)下頁(yè))為密度參數(shù)文件界面,在“rho_Sed.par”中,第一行數(shù)據(jù)“8”對(duì)應(yīng)八個(gè)沉積物平均密度,單位“kg/m3”,“120”代表有一百二十個(gè)點(diǎn)。

圖7和圖8為輸出文件界面,其中“heat_flux.dat”文件和“beta_inv.dat”文件是輸出結(jié)果文件?！癶eat_flux.dat”文件中,第一行數(shù)據(jù)“120”是代表本次參與擬合伸展因子的點(diǎn)數(shù)為一百二十個(gè),“6”代表是輸出六個(gè)時(shí)間點(diǎn)的熱流;“beta_inv.dat”文件中,第一行數(shù)據(jù)“120”是代表本次參與擬合伸展因子的結(jié)果有一百二十個(gè);從第二行數(shù)據(jù)開(kāi)始,第一列數(shù)據(jù)顯示為伸展因子擬合結(jié)果,第二列數(shù)據(jù)顯示為計(jì)算擬合誤差。

圖5 沉降數(shù)據(jù)參數(shù)文件Fig.5 The file of subsidence data

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

4 SO 49-25測(cè)線實(shí)例研究

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

圖7 輸出的熱流結(jié)果文件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測(cè)線北部陸架部份的伸展因子擬合的結(jié)果顯示(見(jiàn)下頁(yè)圖9),在測(cè)線0 km～20 km處伸展量相對(duì)較小,伸展因子在1.7～2.5之間,而且計(jì)算誤差多在20～90之間,分析認(rèn)為是因該區(qū)域位于珠二坳陷,裂后沉降較大,這說(shuō)明伸展因子在擬合過(guò)程中有一定的合理性。在20 km～80 km處伸展量相對(duì)較大,伸展因子變化較小導(dǎo)致,在1.3～2.0之間,而且計(jì)算誤差在20～110之間,分析認(rèn)為是因該區(qū)域位于一統(tǒng)暗沙隆起區(qū),相對(duì)裂后沉降較小導(dǎo)致。在計(jì)算西沙海槽伸展因子結(jié)果時(shí),因發(fā)生裂后沉降結(jié)果較大,導(dǎo)致伸展因子結(jié)果在2.0～6.1之間變化較大,且計(jì)算誤差在50～250之間變化,導(dǎo)致結(jié)果可信度不高。

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

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

對(duì)各個(gè)時(shí)期的熱流值分布的認(rèn)識(shí):

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

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

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

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

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

自23.8M a以來(lái),測(cè)線所在研究區(qū)處于熱沉降階段,熱流從裂陷谷期末期逐漸下降。

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

5 結(jié)論

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

致謝 作者在本文所用SO航次的地震數(shù)據(jù)是由德國(guó)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界面設(shè)計(jì)與編譯技巧[M].北京:電子工業(yè)出版社,2006.

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

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

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

TU 432

A

1001—1749(2011)01—0101—06

國(guó)家基礎(chǔ)研究發(fā)展規(guī)劃項(xiàng)目資助(2007CB411704)

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

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