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

时间：2024-05-22

张炯,宋海斌

(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

1001—1749(2011)01—0101—06