南海北部似海底反射层速度结构全波形反演

霍元媛,张 明

(1.中石化华北分公司勘探开发研究院,郑州 450006;2.广州海洋地质调查局,广州 510760)

0 前言

天然气水合物是一种新型能源,并作为一种潜在的石油天然气替代品得到了广泛的研究[1～3]。其中,似海底反射(BSR)是识别天然气水合物最直观的证据,它是一个从较高速度降至较低速度的界面反射标志,一般与现代海底近于平行,并且多与海底沉积层反射斜交。相对于海底反射,具有较强的反射振幅和极性反转特征[4～8]。

为了正确识别BSR,利用BSR上方、下方的弹性参数结构,来识别并分析天然气水合物及其下方游离气的分布特征和形成机制,许多学者对BSR层上方、下方的地震速度结构进行了研究[9～13]。通常层速度是基于常规叠加速度分析,由Dix公式计算均方根速度得到的,但是利用地震叠加速度进行水合物沉积层的层速度预测,存在较大的误差,其主要原因是在现有的地震勘探频率下,地震叠加速度的变化不能由地震叠加振幅敏感地反映出来,使得地震叠加速度的选取精度难以保证,造成由此求取的层速度分辨率低,误差较大[14]。Singh[15]在1993年首次提出了全波形反演的方法,并将这一方法应用于温哥华岸外的天然气水合物反射层速度结构的研究中,此后,全波形反演的方法得到了广泛的应用和发展[16～20]。在我国天然气水合物的勘探中,全波形反演方法还未发挥过作用。因此,作者在前人研究的基础上,提出了一种新的全波形反演方法,并将其应用到对我国南海北部海域天然气水合物BSR速度结构的研究中,为我国天然气水合物首钻目标优选提供了技术支持。

1 全波形反演

波形反演是一个强非线性问题,目标函数存在众多极值,传统的迭代方法使用的局部最优化法,易于收敛于局部极值,往往得不到全局最优解。全局优化算法对整个模型空间进行搜索,通常能够寻找到全局最优值。而采用全局搜索的方法寻找最优解时,当模型参数很多时,运算量太大。因此作者在本文采用了一种折中的办法,即先用全局搜索方法来解决非线性程度较高的旅行时反演,再用局部搜索方法来解决非线性程度较低的波形反演问题。

目前,比较理想的用于全局优化问题求解的算法主要是:Monte Carlo算法、模拟退火算法以及遗传算法。Sambridge[21]通过一个多元二次函数的最优化问题,阐述了当未知数个数增加时,遗传算法的效率明显高于Monte Carlo算法。因此作者利用遗传算法来进行旅行时反演,获得背景速度模型,以进一步提高传统波形反演方法的计算效率。作者在本文建立的全波形反演方法主要分二个步骤:①用全局搜索的遗传算法构建基于旅行时反演的低频背景速度模型;②用步骤①得到的模型作为初始模型,利用波形资料应用共轭梯度算法,来求取高频速度结构。这样最后得到的模型中就既有低频(长波长)信息,又有高频(短波长)信息。

波形反演的目标函数,是使实际的观测记录波形与计算合成记录的波形差异最小。反演所使用的地震数据,经过严格的拉冬变换,转化为截距时间～慢度域(τ-P域)数据。在τ-P域中,反射以椭圆形轨迹为特征。全波形反演所使用的二种算法如下所述。

1.1 遗传算法

遗传算法是一种应用广泛的随机搜索方法,它是模仿生物进化的过程,按照适者生存的法则建立的一种最优化方法。遗传算法与传统反演方法相比,具有以下一些特点[22]:

(1)遗传算法对模型参数的编码进行操作,而不对参数进行直接操作。

(2)遗传算法同时对多个模型进行操作,而不是对单独的一个模型进行操作。

(3)遗传算法使用随机的优化策略,而不是确定的优化规则。

作者在本文主要通过以下几个步骤来实现旅行时反演:

(1)确定旅行时反演参数,并指定模型参数的搜索范围。在旅行时,反演的参数为层速度,搜索范围为1 400 m/s到2 600 m/s,允许变化范围为±10 m/s。

(2)随机产生初始模型种群。通过随机选择产生的初始模型种群。群体的维数越大,其代表性越广泛,最终进化到最优解的可能性越大。但维数大的群体,势必会造成计算时间的增加。通过对不同模型规模的试验分析,当模型个数为30时,基本能得到满意的反演结果,因此可将模型规模设置为30。

(3)确定目标函数。将目标函数定义为实测的旅行时数据与理论计算旅行时数据的误差平方和:

其中 tmi和tni分别为第i个接收点上实测的和理论计算的旅行时;n为数据个数。

(4)模型参数编码。编码是应用遗传算法时要解决的首要问题,采用格雷码对参数进行编码,这样有利于提高遗传算法的局部搜索能力,并使得交叉、变异等遗传操作便于实现。

(5)目标函数值转换为适应度值其中 σ是目标函数的标准偏差。

(6)算子描述。首先计算个体的相对适应度pi=fi/∑fi,然后根据概率分布,从初始种群中随机选一些染色体,构成一个新的种群。为了避免少数个体霸占整个群体,这里采用了从子代和父代中挑选出的最好的n个可以直接进入子代的策略,这种选择方式收敛速度也比较快。

交叉运算决定了遗传算法的全局搜索能力。为了加速收敛,作者在本文中选择了多点交叉,随机选择交叉位置,交叉算子设为0.9。

变异运算是指将个体染色体编码串中的某些基因,用其它等位基因来替换,从而形成一个新的个体。经过试验,可将变异算子设定为0.05。

(7)终止条件。设置最大遗传代数K,若已完成K次计算,则停止。在本文中K=50。

1.2 共轭梯度算法

采用共轭梯度法[23]得到短波长分量,目标函数为

其中 dcal为由模型M计算出的波场;dobs为观测到的波场;m0为初始模型;模为用数据协方差矩阵cD加权的L2模;为用模型协方差矩阵cM加权的L2模。

作者在本文测试了几种不同的cD和cM结合方式,最后设定为1.0 I和2.0 I,该模型参数为纵波速度。此外在制作合成地震记录时,反射体处出现的不连续,会引入强振幅和高频能量,使得反演不稳定。因此由遗传算法得到的背景速度模型必须进行平滑,在本文中采用三角平均滤波器对模型进行平滑。

1.3 正演算法

进行波形反演,首先要选择一种合适的正演模拟算法。正演模拟有多种方法,其中射线追踪法计算成本比较低,但精确度低并且难以处理复杂介质中不同类型波之间的耦合[24];差分法应用范围比较广泛,能够较精确地模拟任意非均匀介质中的地震波场,但计算成本较高,且需要采用吸收边界条件。作者在本文中采用的是慢度法[25]进行正演,这种方法能够提供全波场解,包括地震波在地层内传播时形成的一次波、多次波和转换波。它的实现步骤为:

(1)通过傅立叶变换和汉克尔变换,将本构方程和弹性动力学方程变换到频率～慢度域。

(2)在频率～慢度域得到反射,透射系数矩阵递推公式。

(3)分别对频率和慢度进行积分变换,得到时～空域的地震道集。

图1 某线偏移剖面(BSR位于1.65 s处,具有强振幅、反极性并与周围地层斜交的特征)Fig.1 The migrated section(BSR is shown clearly on the migrated section at 1.65s,which has characteristics of strong amplitude,reverse polarity and crossing obliquely with surrounding sediments)

2 南海北部海域天然气水合物BSR速度结构反演

作者在本文采用了以上方法,对南海海域的地震数据进行了反演。预处理包括观测系统定义,速度分析,动校正与叠前时间偏移,得到了偏移剖面(见图1)。从偏移剖面上可以看到,箭头所指处的BSR具有与海底近似平行,并与海底沉积层反射斜交,以及强振幅和负极性的特征。作者抽取四个CMP道集形成一个超道集,以增强信噪比。并且此处的BSR与海底几乎都是水平的,满足了在反演方法中,假设速度仅随深度变化的条件。

在进行正演计算时,震源子波由距震源一定距离点上直接观测得到,再通过慢度法得到合成记录。考虑到可能有游离气的存在,使得纵波速度显著下降,而横波和密度受到的影响较小,为了使反演结果更加稳定,可以将横波速度和密度参数固定不变,而只将纵波速度作为自由变量。在反演中,根据Hamilton[26]的经验公式,计算出横波速度和密度。纵波速度搜索范围为1 400 m/s到2 600 m/s。

经过大量试验结果对比,作者认为将遗传算法的参数中,交叉概率设为0.9;变异概率设为0.05比较合适。在实际应用时,可在此基础上,根据具体情况进行调整。模型群体的大小,可视问题的复杂程度,计算机容量,以及计算速度来进行综合考虑。太大占用内存多且计算速度慢,太小则搜索不够彻底,应根据具体情况权衡选择,在本文中将初始种群个数设为30。

图2是遗传算法迭代过程中的平均目标函数的统计图。由图2中可以看到,在迭代的初始阶段,群体具有较快的收敛速度,在迭代十八次后,收敛速度趋于平缓。可以认为,此时种群中的个体已经达到真实数值附近。

图2 目标函数平均值统计图Fig.2 Average value of objective function changes with iteration

在反演时,非常容易受到噪声的影响。当数据含噪声水平过高时,遗传算法会出现种群中差别很大的模型的目标函数值都相等,由此造成自动停止迭代;或者收敛速度非常缓慢,造成已达到最大迭代次数却未得到满意解的情况。因此,通过对不同信噪比的数据进行试验发现,信噪比大于10的数据,其反演得到的结果较为可靠。

图3是波形反演的最终结果。由图3可知,一个速度明显高于周围速度的高速带(箭头1所指处)和一个低速带(箭头2所指处)。在300 m处高速向低速转换,由2 000 m/s骤降至1 450 m/s,这个速度转换边界对应着地震剖面上的BSR,其上方大约为15 m厚的高速带,这表明此处存在天然气水合物,并暗示水合物层的厚度可能在15 m左右。在高速带的下方,可以观测到一个1 450 m/s的速度极小值,远低于沉积背景速度,这有力地证明了游离气的存在。

图3 反演结果图(箭头1指示可能的水合物层,箭头2指示可能的游离气层)Fig.3 Inversion result(Arrow 1 implies gas hydrate sediments,while Arrow 2 implies free gas sediments)

由于纵波速度与水饱和度之间的高度非线性关系,直接根据纵波速度预测游离气的含量是不可靠的。而常用于估算沉积层中游离气含量的Domenico[27]的方法,对于我国南海海域的深部沉积层也不一定适用。因此,我们将纵波速度作为孔隙空间中水饱和度的函数,采用Biot-Geerts ma[28,29]和Gass mann[30]方程计算其值。在计算时,我们采用的模型参数,包括了全部含水饱和时的纵波速度(2 000 m/s)、全部含水饱和时的横波速度(600 m/s)、孔隙度(0.55)、水的体积模量(2.1×109Pa)、气体的体积模量(0.096×109Pa)、固体物质的密度(2.6 g/cm3)、水的密度(1.035 g/cm3)以及气体的密度(0.3 g/cm3)。沉积物种类包括方解石33.06%、石英和斜长石31.75%、云母类34.53%,沉积物骨架的有效体积模量为53×109Pa。其中,气体的属性是根据BSR处的静岩压力计算得到,并考虑了与理想气体的偏差。模型数值参考了前人的研究成果[27、31～34]。我们给定一个纵波速度(1 450 m/s),计算出孔隙空间中的水饱和度大约低于99%,这说明气体含量至少为1%孔隙体积。

至于甲烷的来源,主要有三种可能的来源:①地球深部非生物成因(无机成因)的甲烷;②海底沉积岩中有机质生成的甲烷;③综合以上二种成因的混合来源。我们的研究区位于南海北部陆缘,珠江口盆地珠二坳陷白云凹陷南侧,处在陆架到深海的陆坡位置,其南与南部神狐一统暗沙隆起区相接,海底地形总体呈东北高、西南低的斜坡形态。根据古气候条件,古地理环境及地震相特征综合分析认为:白云凹陷在始新世～早渐新世(断陷期),具有形成大型中深湖相沉积的条件——潮湿的气候环境,全封闭的深洼陷,高的沉积速率[36]。在该时期形成了巨厚的文昌～恩平组烃源岩,目前已处于产生裂解气阶段。近年来,在白云凹陷北坡～番禺低隆起的一系列天然气发现及钻探结果表明,白云北坡气藏区气藏充满度高,气源来自深部的恩平～文昌组烃源岩。通过对区内地质调查站位资料的分析,在区内浅表层沉积物中,普遍存在游离气。甲烷碳同位素δ13C测试结果显示,研究区海域δ13C(‰)值在-46.2‰～-74.3‰(PDB)之间,平均为-60.9‰(PDB),除二个样品的δ13C值为-46.2‰和-51‰(PDB)外,大多数样品的δ13C值小于-57‰(PDB),证实了这一区域浅表层沉积物天然气为生物气。许多站位顶空气甲烷含量在垂向上保持相对较高的丰度,特别是白云凹陷内的甲烷,含量接近120μl/kg和200μl/kg,暗示其深部可能有持续稳定的游离甲烷供应,来源于深部的热解气。由此推测,研究区浅部地层中的天然气可能兼有生物气和热解气。

3 结论

作者在本文中将基于遗传算法的全波形反演方法,应用于我国南海北部海域的地震资料,反演得到了高分辨率的BSR速度结构,清晰的识别出了BSR的速度反转特征。波形反演对于天然气水合物的识别具有非常重要的作用,如果可以在钻井目标区进行有针对性的波形反演,从而得到该区域精细的速度结构,就可以有效地指导井位的确定。在利用遗传算法进行反演时,正演过程非常耗时,而种群优化的过程又需要反复进行正演计算。遗传算法非常适合于并行计算,为了加快计算速度,我们正在进行对遗传并行算法的研究。

通过对反演结果的分析发现,研究区处的BSR,是由其上的高速层和下方的低速层共同作用所造成的。这个高速层可以解释为含天然气水合物的沉积层,而低速层可以解释为含至少1%游离气的薄层。并且,通过与此处AVO属性剖面,瞬时属性剖面,以及波阻抗反演剖面的对比,可以得出相同的结论。

另外,通过对研究区域附近钻探结果,及区内地质调查站位资料的分析,我们认为研究区浅部地层中的天然气,可能兼有生物气和热解气。致 谢:感谢雷新华教授对我的指导和帮助;

感谢广州海洋地质调查局各位老师的大力支持。

[1] LEE J H,BAEK Y S,RYU B J,et al.A seismic survey to detect natural gas hydrate in the East Sea of Korea[J].Mar Geophys Res,2005,3,26(1):51.

[2] TREHU A M,TORRESM E,LONG P E,et al.Three-dimensional distribution of gas hydrate beneath southern Hydrate Ridge:constraints from ODP Leg 204[J].Earth Planet SciLet,2004,222:845.

[3] TATSUO SAEKI,MASAO HAYASH I,TAKAO INAMOR I,et al.Velocity structure of the Kumano basin in the Nankai trough[M].Proceedings of the Fifth International Conference on Gas Hydrates,2005.

[4] SHIPLEY T H,HOUSTON M H,BULLER R T,et al.Seismic evidence for wide-spread possible gas hydrate horizons on continental slopes and margins[J].Am Assoc Pet GeolBull,1979,63:2204.

[5] MACKAYM E,JARRARD R D,WESTBROOK G K,et al.Origin of bottom-simulating reflectors;geophysical evidence from the Cascadian accretinary prism[J].Geology,1994,22:459.

[6] HOLBROOKW S,HOSK INS H,WOOD W T,et al.Methane hydrate and free gas on the Blake Ridge from vertical seismic profiling[J].Science,1996,273:1840.

[7] RE MPEL A W,BUFFETT B A.Formation and accumulation of gas hydrate in porous media[J].Geophys Res,1997,102:10151.

[8] HYNDMAN R,SPENCE G A.Seismic study of methane hydrates marine bottom simulating reflectors[J].Geophys Res,1992,97:6638.

[9] CHENGW B,LEE C S,L IU C S.Velocity structure in marine sediments with gas hydrate reflectors in offshore S W Taiwan,from OBS data tomography[J].Ter Atm Oce,2006,17(4):739.

[10]ZILLMER M,FLUEH E R,PETERSEN J.seismic investigation of a bottom simulating reflector and quantification of gas hydrate in the Black Sea[J].Geophys J Int,2005,161(3):662.

[11]JAISWAL P,ZELT C A,PECHER IA.seismic characterization of a gas hydrate system in the Gulf ofMexico using wide-aperture data[J].Geophy J Int,2006,165(1):108.

[12]NETZEBAND GL,HUBSCHER C P,GAJEWSKID J,et al.Seismic velocities from the Yaquina forearc basinoff Peru:evidence for free gas within the gas hydrate stability zone[J].Int J Earth Sci,2005,94(3):420.

[13]ECKER C,DVORK IN J,NUR A M.Est imating the amount of hydrate and free gas from marine seismic data[J].Geophysics,2000,65(2):565.

[14]陈建文,吴志强,龚建明.天然气水合物的地球物理识别技术[J].海洋地质动态,2004,20(6):1.

[15]S INGH S C,M INSHULL TA,SPENCE G D.Velocity structure of gas hydrate reflector[J].Science,1993,260:204.

[16]M INSHULL T A,SINGH S C,WESTBROOK G K.seismic velocity structure at a gas hydrate reflector,offshore western Columbia,from full waveform inversion[J].Geophys Res,1994,99:4715.

[17]PECHER I A,MINSHULL T A,SINGH S C,et al.Velocity structure of a bottom s imulating reflector offshore Peru:Results from full waveform and travel time inversion of wide-angle seismic data[J].Geophys,1997,102:15345.

[18]PECHER I A,RANERO C R,HUENE R,et al.The nature and distribution of bottom simulating reflectors at the Costa Rican convergent margin[J].Geophys J Int,1998,133:219.

[19]YUAN T,SPENCE GD,HYNDMAN R D,et al.Seismic velocity studies of a gas hydrate simulating reflector on the northern Cascadia continent margin:Amplitude modeling and full wave for minversion[J].Geophys Res.,1999,104:1179.

[20]宋海斌,OSAMU M,KURAMOTO S.天然气水合物似海底反射层的全波形反演[J].地球物理学报,2003,46(1):42.

[21]SAMBRIDGE M,DRIJKONINGEN G.Genetic algorithms in seismic waveform inversion[J].Geophys,JInt.,1992,109:323.

[22]石琳珂,孙铭心,王广国,等.地球物理遗传反演方法[M].北京:地震出版社,2000.

[23]KORMEND I F,D IETR ICH M.Non-linear waveform inversion of plane wave seismograms in stratified elastic media[J].Geophysics,1991,56:664.

[24]KENNETT B L N,KERRY N J.seismic waves in a stratified half-space[J].Geophys J R Astron Soc,1979,57:557.

[25]张繁昌,印兴耀,赵剑.慢度法全波场模拟及其水平慢度积分的实现[J].石油地球物理勘探,2003,38(6):597.

[26]HAM ILTON E L.Vp/Vs and poisson’s ratios in marine sediments and rocks[J].Acoust Soc Am,1979,66(4):1093.

[27]DOMEN ICO S N.Elastic properties of unconsolidated porous sand reservoirs[J].Geophysics,1977,42:1339.

[28]B IOTM A.Theory of propagation of elastic waves in a fluid-saturated porous solid[J].Acoustic Soc Am,1956,28:168.

[29]GEERTS MA J.Velocity-log interpretation:The effect of rock bulk compressibility[J].Soc Pet Eng J,1961,1:235.

[30]GASS MANN F.Uber die Elastizitat poroser Medien[J].Veirteljahrsschr Natur forsche Gesellschaft Ges Zurich,1951,96:1.

[31]ANGENHEISTER G,ED.,LANNDOLT-BOERNSTE IN.Numerical data and functional relationships in science and techno ledge[M].Physical properties of rocks.Springer,Berlin,1982.

[32]ATK INS P W.Physical chemistry[M].Oxford Univ Press,Oxford,1978.

[33]WY MAN R E.Petrogeo physics,the interrelationships of petrophysics,geology and geophysics,AAPG continuing education-Stratigraphic interpretation of seismic data[M].Am Assoc Pet Geol,1982.

[34]房殿勇,王汝建,邵磊,等.南海ODP184站深海相渐新统硅质成岩作用[J].海洋地质与第四纪地质,2002,22(2):75.

[35]CASTAGNA J P,BATZLE M L,EAST WOOD R L.Relationships be tween compress ional-wave and shear-wave velocities in elastic silicate rocks[J].Geophysics,1985,50:571.

[36]张光学,黄永样,陈邦彦.海域天然气水合物地震学[M].北京:海洋出版社,2003.