天然⽓⽔合物似海底反射层的全波形反演
天然⽓⽔合物似海底反射层的全波形反演
宋海斌1,2 Matsubayashi Osamu 2 K uram oto Shin ’ichi 2
1中国科学院地质与地球物理研究所，北京 1001012⽇本产业技术综合研究所，筑波 3058567
摘 要 建⽴了天然⽓⽔合物似海底反射层(BSR )研究的全波形反演⽅法.这是⼀种将⽔平层状弹性介质的反射
共中⼼点道集转换为截距时间-⽔平慢度域的反演⽅法.反演过程中采⽤了全局搜索⽅法与⾮线性局部搜索⽅法.分两步进⾏.第⼀步是根据⾛时数据应⽤⾮常快速模拟算法求得速度结构的长波长分量.第⼆步，利⽤波形资料⽤共轭梯度法求得速度的短波长扰动分量.这样，最后反演得到的速度结构模型包含了长波长与短波长分量.反演中利⽤了多⽹格参数化技术.⽇本东南海海槽双BSR 的速度结构的反演表明，全波形反演是天然⽓⽔合物BSR 研究的重要⼿段之⼀.
关键词 全波形反演 天然⽓⽔合物 似海底反射层 东南海海槽沈阳旅游攻略必去景点
⽂章编号 0001-5733(2003)01-42-05 中图分类号 P739 收稿⽇期 2001-12-01,2002-09-05收修定稿.
FU LL WAVEFORM INVERSION OF GAS H YDRATE 2RE LATE D
BOTTOM SIMU LATING REFLECTORS
S ONG H AI BI N 1,2 MATS UBAY ASHI O S AM U 2 K URAMOT O S HI N ’ICHI 2
1Institute o f G eology and G eophysics ,C AS ,Beijing 100101,China
2National Institute o f Advanced Industrial Science and Technology ,T sukuba 3058567,Japan
Abstract A full waveform inversion strategy is propod for gas hydrate 2related bottom simulating reflector.This method is ud for stratified elastic media and is im plemented in intercept time 2slowness domain.It consists of tw o steps with a combination of a global arch method and a nonlinear local arch method.In the first step ,long wavelength parts of velocity structure are estimated by very fast simulating alg orithm ,when travelling time in formation is ud.In the cond step ,short wavelength velocity variations are included into the m odel by conjugate gradient method ,and waveform information in intercept time 2slowness domain is ud.Multi 2grid parameter technique is ud in our full waveform inversion.The method is applied to study a double BSR in the eastern Nankai accretionary wedge ,and the results show that the full waveform inversion is an im portant tool for gas hydrate rearch.
K ey w ords Full waveform inversion ,G as hydrates ,Bottom 2simulating reflector ,Eastern Nankai trough.
1 引 ⾔
⼤陆边缘天然⽓⽔合物稳定带的底界通常由似
海底反射层BSR (Bottom Simulating Reflector )识别[1].与天然⽓⽔合物有关的BSR 在世界⼤陆边缘的地震剖⾯上⼴泛分布[1～3].BSR 与海底起伏相似，通常与地层层理斜交，并表现为⼀强振幅的负极性反射.BSR 通常被认为是⾼速含天然⽓⽔合物沉积层与下
伏低速含游离⽓沉积层的反射边界[4].
1993年Singh 等在Science 上发表了应⽤全波形反演⽅法研究温哥华岸外的BSR 上、下⽅速度结构的论⽂[4].他们的研究表明,BSR 下⽅存在着30m 厚的低速游离⽓带.⼤洋钻探146航次[5,6]基本上证实了这⼀结果.这⼀波形反演⽅法被⼴泛应⽤于世界⼤陆边缘的天然⽓⽔合物研究中，如温哥华岸外的Cascadia ⼤陆边缘[4,6～8]、哥伦⽐亚西海岸[9]、秘鲁岸外[10]，布莱克海台[11]、哥斯达黎加岸外[12]、阿
基⾦项⽬ 国家⾃然科学基⾦(49904007)、ST A Fellowship 、中国科学院全国优秀博⼠学位论⽂专项
资⾦、国家重点基础研究项⽬(G 20000467).作者简介 宋海斌，男,1968年⽣,1990年毕业于同济⼤学,1996年获博⼠学位,1998年中国科学院地球物理研究所博⼠后出站，副研究员.主
要从事海洋地球物理、反射地震研究.E2mail:*****************:///doc/d7233f5f312b3169a451a4f0.html
第46卷第1期
2003年1⽉
地球物理学报
CHI NESE
JOURNA L
OF
GE OPHY SICS
V ol.46,N o.1
Jan.,2003
拉伯海的Makran增⽣楔[13]等.这些波形反演得到的速度结构被⽤于研究⽔合物与游离⽓的分布，包括厚度、饱和度、甲烷的数量.在天然⽓⽔合物的地震研究中，全波形反演被证实是⾮常有⽤的研究⽅法.本⽂建⽴⼀种⽤于天然⽓⽔合物BSR研究的全波形反演⽅法，并应⽤于⽇本东南海海槽增⽣楔的双BSR的速度结构的研究.
2 全波形反演
本⽂建⽴的全波形反演思路与Singh等[4]相近，是⼀种⽔平层状弹性介质的截距时间(τ)2⽔平慢度(p)域反演⽅法，分两步.第⼀步是根据τ2p域⾛时数据应⽤⾮常快速模拟退⽕算法求得速度结构的长波长部分.第⼆步，基于第⼀步的速度模型作为初始模型，利⽤波形资料⽤共轭梯度法求得速度的短波长扰动部分.这样，最后反演得到的模型包含了速度结构的长波长与短波长部分.反演结合了全局搜索⽅法与⾮线性局部搜索⽅法.
这⼀反演⽅法包含了⾛时反演与波形反演，全局反演与局部优化反演.⾛时反演⽤于求取长波长速度信息，波形反演求取短波长速度信息.前者⾮线性程度较⾼，应⽤了全局搜索⽅法,⽽后者⾮线性程度较低，应⽤了局部搜索⽅法.这也是当前实现地震波形反演的现实途径.全局反演相当耗时,⼤量模型参数的全局反演不太现实，因此采取了折衷的⽅法.正演是⽤K ennett[14]的反射透射矩阵⽅法实现的.反演
的数据是共中⼼点道集经过平⾯波分解得到τ2p域地震数据.并需考虑震源⼦波估算、地震数据标定与组合校正等各种预处理.反演中应⽤的快速模拟退⽕算法与共轭梯度算法如下所述. 2.1 ⾮常快速模拟退⽕算法
会议决议模拟退⽕算法⾃1983年K irkpatrick等[15]提出以来，已成为⼀种重要的⾮线性优化问题的计算⽅法，其最⼤的优点是可以得到问题的全局最优解.常规的模拟退⽕法即K irkpatric等提出的Metropo1is 算法，包括以下步骤:(1)给定模型每⼀参数变化范围，在这个范围内随机选择⼀个初始模型m0，并计算相应的⽬标函数值E(m0).(2)对当前模型m0进⾏扰动产⽣⼀个新模型m，计算相应的⽬标函数值E(m)，得到ΔE=E(m)-E(m0).(3)若ΔE< 0，则新模型m被接收；若ΔE>0，则新模型m按概率:P=exp(-ΔE/T)进⾏接收,T为温度.当模型被接收时，置m0=m,E(m0)=E(m).(4)在温度T下，重复⼀定次数的扰动和接收过程，即重复(2)、(3)步骤.(5)缓慢地降低温度T.
(6)重复(2)、(5)步骤，直⾄收敛条件满⾜为⽌.上述常规的模拟退⽕法包括两个步骤：对模型随机扰动并计算能量变化(即⽬标函数变化)；模型按概率P进⾏接收.在温度很⾼时，使⽬标函数增⼤的差模型能被接收，因此它具有跳出局部极值的能⼒,⽽当温度较⼩时，只接收好模型.算法只要从⾜够⾼的温度开始并缓慢降温，则最终得到全局最优解.
Ingber在其提出的⾮常快速模拟算法[16]中，采⽤了依赖于温度似Cauchy分布产⽣新模型，即:
m i=m0i+y i(B i-A i),
y i=T sign(u-0.5)[(1+1/T)2u-1-1],
式中,m0i为当前模型的第i个变量,u为0～1之间均匀分布的随机数,[A i,B i]为m i的取值范围，且要求扰动后的m i∈[A i,B i],sign()为符号函数.采⽤依赖于温度似Cauchy分布产⽣新模型改进加快了模拟退⽕算法的收敛速度.在我们的反演中，模型参数为层速度,⽬标函数为τ2p域⾛时计算值与观测值的均⽅差.
2.2 共轭梯度算法
建⽴在最⼩⼆乘法基础上的最优化⽅法为⼤多数反问题的计算提出了统⼀的框架.与频率域-慢度域波形反演相似[17]，截距时间2慢度域波形反演可通过观测波形数据与计算波形数据的极⼩化,⽤局部⾮线性优化⽅法(共轭梯度法)迭代修改模型完成.截距时间2慢度域正演问题可以简化表⽰为
d=f(m),
其中,d为数据向量，即τ2p域计算波场;m为模型向量,⽤等厚度的⽔平层状介质描述地下模型，模型参数即为各层的弹性参数；如果反演中保持横波速度、密度不变，则模型参数为纵波速度;m0为先验模型，即初始的弹性参数值;d obs为τ2p域观测值.最⼩⼆乘意义下的⽬标函数为
S(m)=1
2
‖d-d obs‖2D+‖m-m0‖2M,
其中,‖d‖2D=d T C-1D d,‖m‖2M=m T C-1M m,C D为数据协⽅差矩阵,C M为模型协⽅差矩阵.
⾮线性迭代反演就是寻找模型使⽬标泛函极⼩.在梯度的导引下求得⽬标函数的极⼩值，有了梯度，可以利⽤最优化技术中的很多⽅法求解，如最速下降法、拟⽜顿法、变尺度法、共轭梯度法等.共轭梯度迭代法步骤[17]如下:
(1)由第n次迭代得到的模型m n计算合成数据
d n=f(m n),n=0,1,2,…;
34
1期宋海斌等：天然⽓⽔合物似海底反射层的全波形反演
(2)计算合成数据与观测值的差值
Δd
n
=d n-d obs
及当前模型与先验模型的差Δm n=m n-m0;
(3)计算⽬标函数值
S(m n)=Δd T n C-1DΔd n+Δm T n C-1MΔm n
并检验它是否满⾜终⽌条件;
(4)计算最速上升⽅向(梯度)
Γ
n
=C M F T n C-1DΔd n+Δm n,
其中,F n=5f/5m m=m
n关于除夕的作文
为Frechet导数，可以解析求得.为程序编制⽅便，我们采⽤了数值解.最速上升⽅向(梯度)的计算式中取掉第2项，即不考虑先验模型在梯度计算中的影响，这是可取的，本⽂计算中采⽤了这⼀⽅案.
(5)计算共轭⽅向
Φ
n
=Γn+σnΦn-1,Φ0=Г0,
σ
n =
(Γn-Γn-1)Γn
ΓT
n-1
Γ
n-1
;
(6)计算使⽬标函数极⼩的最佳步长µn;
(7)计算新的模型
m n+1=m n-µnΦn.
3 ⽇本东南海海槽双BSR数据的反演
2000年夏天SF J2000(法⽇合作地震带探测Seism ogenic zone experiment by French2Japane cooperate study)航次在东南海海槽地区进⾏了地震数据采集.⽤道间距为12.5m的360道⽔听器记录地震数据.采样间隔为1ms，记录长度为6s.近炮检距是275m.5个不同体积的⽓枪组合成2.47×10-3 m3⽓枪阵作为震源.我们对其中的HR09测线(图1)进⾏了旨在获取海底浅部结构及BSR分布的清晰图像的地震数据处理.地震处理流程与⽂献[3]相似，包括观测系统定义、共中⼼点道集选排、速度分析、动校正与偏移.图2是HR09测线⼀段偏移剖⾯.第⼀BSR与第⼆BSR均具有强振幅、负极性、切穿地层层理反射的特点.
保用
C MP2140处第⼀BSR与第⼆BSR清晰，第⼀BSR的坡度约为4°，海底与第⼆BSR基本与它平⾏，因此可以近似为⼀维结构.C MP2140附近的4个
C MP道集组成道间距为12.5m的超道集，并转换成τ2 p域数据(图3).在τ2p域拾取海底、第⼀BSR、第⼆BSR及第⼆BSR下⽅两个反射层的⾛时(图4).⽤⾮常快速模拟退⽕算法通过⾛时观测值与计算值的极⼩化求得层速度，反演利⽤了10个随机初始模型，模型的垂直⾛时可以在拾取的垂直⾛时附近32ms的范围内变动.最终反演的速度结构见图5中的曲线1,⾛时观测值与计算值见图4.对曲线1进⾏平滑得速度趋势曲线2(图5)，并作为波形反演的初始模型.
震源⼦波从海底反射估算.在反演中进⾏了组合校正.由于震源排列垂直于测线，因此只对检波器排列进⾏了组合校正.横波速度通过Hamilton经验公式[18]根据纵波速度计算，密度则通过G ardner 公式计算.横波速度、密度在反演中保持不变.这时假设BSR两侧横波速度与密度变化较⼩，这样的假设在BSR的全波形反演中被经常采⽤.
初始平滑速度模型离散成每层厚度为10m的模型，进⾏第⼀轮波形反演，结果见图5.然后把这个反演结果离散成每层厚度5m的模型，并进⾏⼜⼀轮的波形反演，最终结果见图6，因此，我们的反演中利⽤了多⽹格参数化技术.
图1 东南海海槽地区⽔深与HR09测线位置Fig.1 Bathymetry and location of ismic reflection profi
le HR09in the area of
eastern Nankai T rough.
图2 HR09偏移剖⾯段
Fig.2 Section of M igration profile of HR09
44地球物理学报(Chine J.G eophys.) 46卷
图3 C MP2140～2137转换得到的截距时间-慢度域地震数据
Fig.3 T rans formed ismic data of C MP2140～
恐怖故事大全
2137in
intercept time 2slowness domain
图4 截距时间-慢度域5个反射的
⾛时观测值与计算值
遗闻逸事Fig.4 The obrved and calculated travelling time of five
reflectors in intercept time 2slowness domain
图5 速度模型⽐较
青少年活动1 模拟退⽕算法反演的速度模型;2 基于曲线1的速度趋势
模型;3 以曲线2为初始模型，全波形反演的10m 厚度速度模型;4 以曲线1为初始模型，全波形反演的10m 厚度速度模型.
Fig.5 C om paris on of velocity m odels 图6 速度模型⽐较
1 以图5曲线3为初始模型，全波形反演的5m 厚度速度模型;2 基于图5曲线1的速度趋势模型;3 以图5曲线4为初始
模型，全波形反演的5m 厚度速度模型.
Fig.6 C om paris on of velocity m odels
为了考察反演的灵敏度，我们⽐较了不同初始模型的反演结果.以阶梯速度模型为初始模型的10m 厚度模型反演结果见图5中的曲线4.曲线4与基于平滑初始模型的曲线3相似.以图5曲线4为初始模型,5m 厚度模型的反演结果见图6曲线3，它与曲线2⾮常接近.不同初始模型的反演结果⽐较表明，反演结果是可靠的.
图6中的反演结果表明，有两个速度⾼于趋势速度的⾼速带，它们的下⽅则存在低速带.上、下两个⾼速到低速的转换边界分别对应第⼀BSR 与第⼆BSR.第⼀BSR 上⽅与第⼆BSR 上⽅的⾼速带与沉
积物中含天然⽓⽔合物有关,⽽第⼆BSR 下⽅的低速带与沉积物中含游离⽓有关.第⼀BSR 下⽅的低速带可以解释成不含天然⽓⽔合物的饱和⽔沉积物或者是含极少游离⽓的沉积物.
4 ⼩结什么是企业愿景
常规速度分析与振幅分析很难区分天然⽓⽔合物楔模型与游离⽓带模型[9]，难以揭⽰似海底反射层的成因与性质.全波形反演则正确地解释了双BSR 的反极性特征，给出了⾼分辨率的速度结构，可
5
41期宋海斌等：天然⽓⽔合物似海底反射层的全波形反演
以⽤于天然⽓⽔合物与游离⽓分布的研究.
双BSR可能揭⽰着⼀个动态演化的天然⽓⽔合物系统，因此对它的研究备受关注.对第⼆BSR通常有两种解释[19].第⼀种观点认为第⼆BSR是古BSR，对应着古温压条件下的天然⽓⽔合物稳定带的底界反射.第⼆种观点认为第⼆BSR对应不同于甲烷⽔合物,⽽由多种⽓体成分形成的天然⽓⽔合物地层的底界反射.反演给出的天然⽓⽔合物层与游离⽓层的分布与其他地质、地球化学资料的结合则可望给出双BSR⽐较合理的解释.
本⽂给出的是单分量地震数据全波形反演，它可以揭⽰单参数———纵波速度的分布特征.⽽多分量地震数据的反演能较好地分辨横波速度与密度[20～22],⼀些⼤陆边缘已⽤或正在⽤海底地震仪采集三分量⼴⾓数据，可以预见，多参数的全波形反演是新世纪天然⽓⽔合物研究的重要⼿段.本⽂对天然⽓⽔合物BSR的全波形反演的⽅法研究与应⽤研究均是初步的，进⼀步⼯作正在进⾏之中.
感谢参加SF J2000调查的所有科研⼈员与船员，感谢⽇本产业技术综合研究所Okuda博⼠、T anahashi博⼠、K ano博⼠的热情帮助.感谢袁天声博⼠的有益建议.感谢⽇本振兴集团JST与国际交流中⼼J ISTEC在第⼀作者作为ST A Fellow期间的⼤⼒⽀持.
参考⽂献
[1] Shipley T H,H ouston M H,Buller R T,et al.Seism ic evidence for
wide2spread possible gas hydrate horiz ons on continental slopes and
margins.Am.Assoc.Pet.G eol.Bull.,1979,63:2204～2213 [2] 宋海斌，耿建华,W ong H owkin等.南海北部东沙海域天然⽓⽔
合物的初步研究.地球物理学报,2001,44(5):684～691
S ONG Haibin,GE NGJianhua,W ONG H owkin,et al.A prelim inary
study of gas hydrates in D ongsha region north of S outh China Sea.
Chine J.G eophys.,2001,44(5):684～691
[3] 宋海斌，松林修，仓本真⼀.西南海海槽地震资料处理及其
BSR特征.地球物理学报,2001,44(6):785～790
S ONG Haibin,M atsubayashi O,K uram oto S.Seism ic data
processing of western Nankai T rough and the character of its BSR.
Chine J.G eophys.,2001,44(6):785～790
[4] S ingh S C,M inshull T A,S pence G D.Velocity structure of a gas
