第24卷第2期
2009年4月(页码:391~397)
地球物理
IN
学进展
v01.24,No.2
Apr.2009
PROGRESSGEOPHYSICS
皮红梅,蒋先艺,刘财,等.波动方程数值模拟的三种方法及对比.地球物理学进展,2009,24(2):39l~397,DoI:10.3969/j.
issn.1004—2903.2009.02.004.PiH(in
M,JiangXY,LiuC,Pf口Z.Threemethodsofwaveequationforwardmodelingandcompari50ns.PrDg化5s锄GPo户^∥.
Chinese),2009,24(2):391~397,DoI:lO.3969/j.issn.1004—2903。2009.02.004.
波动方程数值模拟的三种方法及对比
皮红梅1’2,
蒋先艺2,
刘
财¨,
王成祥2,
姜绍辉2
(1-吉林大学地球探测科学与技术学院。长春130026}
2.中国石油东方地球物理公司物探技术研究中心,涿州072750)
摘要波动方程数值模拟方法是研究地震渡场传播的一种重要手段,本文采用交错网格高阶有限差分方法分别对
双程声波方程和双程弹性波方程进行了波场数值模拟,并且根据定位原理采用傅立叶有限差分算子进行了单程波方程数值模拟,在分析定位原理的基础上,对其计算过程稍作修改,将延拓到地面的波场直接由每个检波点接收,无需横向叠加过程,得到了单程声波方程共炮记录.基于不同波动方程的数值模拟结果表明,双程波方程结果包含直达波、多次波等干扰渡,信噪比低;单程波数值模拟结果只包含了介质分界面的一次反射波,信噪比高,但对于大角度入射渡误差较大,并且对于同一个地质模型而言,双程弹性波方程计算速度最慢,双程声波方程次之,单程声波方程计算速度最快.因此对于复杂地质模型,三种模拟方法可以取长补短,综合应用.关键词波动方程数值模拟,交错网格,高阶有限差分法,定位原理,单程波
DOI:lO.3969/j.issn.1004—2903.2009.02.004
中图分类号
P315
文献标识码
A
ThreemethodsofwaVeequationforwardmodelingandcomparisons
PIHong—meil”,
(1.CDzZFgP
JIANGXian—yi2,
o厂GP08zpZor乜£io靠Sc{P,lce
LIUCai¨,
n行d
WANG
Cheng—xian92,
JIANGShao_hui2
T■f_Il雅ozogy,.,订i以UhitJP,百f£,,(冼口”gc|Il配种130026,c^f,l口;
2.R#s∞”九&DP御ZD户m£斗£CP咒fPro,BGP,CNPC,Z^“oz^o“072750,Chi靠口)
Abstr扯tWave
equationforwardmodelingis
on
averyimportantway
to
studytheseismicpropagationunderground.In
thispaper,seismicforwardmodelingbaseddone
using
the
staggered—grid
on
two-wayacousticwaveequationandtw0-wayelasticwaVeequationis
method.made
Besides,
according
t0
high_order
differencethe
positioningprinciple,
Inthe
simulationbased
progress
one-wayacousticwaveequationis
bytheFourierfinitedifference
propagator.
to
ofcomputing,eachgeophonedirectlyreceivesthewavefieldthathasbeenextrapolated
Thethreesimulation
resultsshow
thesurfaceofthe
generate
a
modelwithoutlateralstacking.waveand
thattwo.waywa、reequations
directive
muItiplewavewhichreducetheSNRseverely,whiletheone—waywaveequationmethodgives
createscosts
higherSNR
withonlyrenectedwaye,butitelasticwaveequation
larger
errors
withwide-an一eincidentwa’re.Forthesamemodel,two—way
costs
algorithm
the
mostcomputationaltime,one_wayacousticwaveequationalgorithm
algorithminbetween.Sofor
a
theleastandtw0-wayacousticwaveequationmodelingKeywords
complex
structuremodel,
thethree
methods
can
beusedsimultaneouslywithdifferentadvantagesanddisadvantages.
high.orderdifferencemethod,positioningprinciple,
waveequationforwardmodeling,staggeredgrid,
one.waywaveequation
收稿日期作者简介
2008_06一12;
修回日期
2008一09-2z.
皮红梅,1981年生,博士研究生,2003年毕业于吉林大学地球物理专业,获理学学士学位,现在该校地球探测科学与技术学院攻读固体地球物理学专业博士学位,研究领域为复杂地震波场数值模拟.(E—mail.brook_3993@email.jlu.edu.cn)
通讯作者刘财,吉林大学教授,博士生导师,主要从事固体地球物理研究.(E-mamliucaUlu@163.com)
万方数据
392
地球物理学进展
24卷
O
引言
算效率低,正演记录的信噪比低.贺振华、熊高君(1998)等提出检波点下延原理[8 ̄儿],首次采用单程声波方程md33将炮点和检波点的地震波场同时从地面向地下进行延拓得到共炮正演记录,但该方法是从偏移的思路来进行正演,计算效率较低,于是又提出一种能快速合成共炮记录的固定检波点定位原理n引,此后,贺振华、胡光岷(2004)等人提出的数学检波器的概念[15]和贺振华、熊晓军(2005)提出的等时叠加原理[161,从理论上对单程波动方程叠前正演方法进行了归纳和总结.本文在分析固定检波点定位原理的基础上,对其计算过程稍作修改,将延拓到地面的波场直接由每个检波点接收,无需横向叠加过程,得到了单程声波方程共炮记录,对三种波动方程得到的共炮记录进行分析对比.最后进一步分析了双程声波方程和单程声波方程的自激自收记录.
1
地震数值模拟方法是地震勘探和地震学的重要基础,是根据已知的地下介质结构和相应物理参数模拟研究地震波在地下各种介质中的传播规律,并计算在地面或地下各观测点的地震记录的一种地震模拟方法.数值模拟方法主要有三大类:波动方程法,积分方程法和射线追踪法.波动方程数值解法是建立在以弹性或粘弹性理论和牛顿力学为基础的双曲型偏微分方程一地震波传播方程的理论基础上的.由于地下介质性质不同,其相应的地震波传播方程也不同.如声学介质中的声波波动方程;弹性介质中的弹性波波动方程;粘弹性介质中的粘弹性波波动方程;孔隙弹性介质(双相或多相介质)中的双相或多相介质弹性波方程;各向异性介质中的各向异性弹性波波动方程等.积分方程法是建立在以惠更斯原理为基础的波叠加原理基础上的,其数学表达形式为波动方程的格林函数域积分方程式和边界积分方程式.射线追踪法是建立在以射线理论为基础的波动方程高频近似理论基础上的,其数学表达形式为程函方程和传输方程.
波动方程法又分为双程波方程法和单程波方程法.本文采用交错网格高阶有限差分法对双程声波方程和双程弹性波方程进行数值模拟.交错网格高阶有限差分法[1 ̄3]是董良国(2000)首先提出的,他将Mad耐aga(1976)提出的交错网格法[43与Dablain(1986)提出的高阶差分方法[5]有机结合,用于求解
双程声波方程数值模拟
直角坐标系下,各向同性介质中的二维声波方
程的一阶压力一速度方程m3形式为:
警一一K(篆+碧),
垫一一上丝.
a£
p
az’
婺一一土挈.
at
p
<1)
、1
az’
式中,户一户(z,z,£)为声波波场,%=砚(z,2)为质点振动速度的水平分量,让一皱(z,名)为垂直分量,K一即2为弹性模量,口为介质的纵波速度,lD为介质密度.
采用交错网格(如图1所示)高阶有限差分方法对方程(1)进行数值离散,其2阶时间2N阶空间的离散格式为:
横向各向同性介质一阶速度一应力弹性波方程,其
后他又对该方法的稳定性问题[6](2000)和频散问题"](2004)进行了研究.双程波方程法包含了地震波场的全部信息,可以处理复杂的地质构造,但是计
1.....................................................一——..——————一..——.——..
矿%∽=一∽沪警{塾Ⅲ+竽,J)一谚(z一竿,歹)]
十塾[破(奶+学)一破(叫一竽)]},
谚
l
+
1—2
一记
.Z
+
1—2
一
妻“[户蚪}(i+咒,歹)一户外专(i+1一行,歹)],(2)
矿
l
}
1—2一
破
.£
}
1—2一
其中:∥+专(i,_『)表示户(z,z,£+譬),
万方数据
塑一
Ⅳ
惫箭一∞一∞
∑c。[户蚪}(i,歹+咒)一夕蚪}(z,歹+1一咒)].
^=1
谚(f+丢,_f)表示V。(z+等,2,t),其它与此类
2期皮红梅,等:波动方程数值模拟的三种方法及对比
393
同.厶为交错网丰子差分系识.
2双程弹性波方程数值模拟
均匀各向同性介质中的二维弹性波方程m3为:
.o等=鲁+等,P警=鲁+等,
鲁=(A+2p)卺+|;I卺,
鲁__(A+刎凳+A警,鲁=p(等+卺).
其中,让和钞。分别为质点振动速度的水平分量和垂直分量;%,%,屹为应力分量;lD为介质密度;A和口为拉梅常数.
■px
1,
●vx’,p
∥
o屹,形
图1声波方程交错网格示意图
Fig.1
Thestaggeredgriddacousticwaveequatio矗
采用交错网格(如图2所示,据Levander[193)高阶有限差分法来求解方程(3).
图2
弹性波方程交错网格示意图(据Levander1988)
F螗.Z
Thestaggeredgridofelastic
waveequation(Levander1988)
万
方数据呓{=吒吾+考惫*[蚤“《字。,一P:2孚≯
+∑c。(《.升学一≤,卜学)],
’(4)
咙U2,尸:。/2.』,s:一。,。分别为质点速度分量饥和应
力分量盯。,%的离散值,其它几个分量的离散方法与此类似.
3单程声波方程数值模拟
二维各向同性均匀介质中的单程波方程田3为:
筹一土i志庐=士i√等一磋夕,,
(5)
歹=声(走。,名,甜)为频率一波数域声波波场,“+”、“一”号分别代表上行波方程和下行波方程.则根据D.Ristow(1994)提出的傅立叶有限差分‘213(FFD)法可得用于正演的波场延拓算子为:
f声(忌。,锰l,c£,)一声(惫:,≈,(c,)e_‘√孟巧△f
{石(z,铣l,叫)一万(z,巧,(c,)fb血b∞缸
.
(6)
【户(z,z灶1,(c’)=夕(z,≈,cE,)r巩缸
其中:石=石(z,2,∞)为频率一空间域声波波场,c为
每一个延拓深度层内的背景速度,取为:
c=min(口(z,2。)),
纵…,一志一÷,
A_钟吲【蘸J’
f生竺]
6=丢[(詈)2+詈+-].
根据贺振华,熊高君提出的固定检波点定位原理[H3合成共炮记录,对其计算过程稍作修改,将延拓到地面的波场直接由每个检波点接收,无需横向叠加过程.原理表述如下:固定炮点和检波点,用下行波方程将震源波场从炮点位置延拓到反射界面,与所有反射点的位置做相关运算,相当于模拟入射波沿反射界面激发新震源或反射波,把这个新震源暂时置于所产生的位置,同时把原来的入射波继续向下延拓直到最深处,形成所有深度处的新震源,再从最深处用上行波方程将新震源向上延拓,每延拓一步,将延拓波场与该处的新震源波场叠加,将叠加后的波场继续向上延拓,直到检波器所在位置,将该延拓波场与检波器谐振因子进行相关,该相关结果就是所固定检波点的正演结果.
记:地面记录59(z,£),震源波场s(z,2,f),检
其2阶时间、2N阶空间的离散格式为:
地球物理学进展
24卷
波谐振波场g(z。,zo,t),反射系数r(z,z).其中:
s(z,州)一F五嘞一几z一五
(7)
lO,z≠z,
(z.,勐)、(z。,‰)分别为震源位置和检波点位置,本论文中震源和检波因子都用ricker子波.将震源和检波因子都变换到频率域i(z,z,cc,)、虿(以,孙,甜),并且反射系数拓展为频率不变函数尸(z,名,∞),即
尹(z,z,∞)=r(z,2).
固定检波点原理P-P波合成记录的步骤如下:(1)根据(6)式将震源波场向下延拓,直到地下最深处,并且与延拓深度上的反射系数进行相关,形成新震源:
j(z,z件1,∞)=j(z,z,,(£’)e一让。血,
(8)西(z,≈+l,叫)一j(z,乃+1,(£,)尹(z,影件1,叫).
(9)
(2)根据(6)式从地下最深处将反射界面上的所有新震源向上延拓,并与该延拓层新震源叠加,直到将波场延拓到检波器所在的深处(假设地面z=0处).
歹1(z,争1,c£,)=歹l(z,乃,倒)e-让r缸+而(z,夸l,(cJ)
(10)
其中歹l(z,名。,叫)一鄢(z,名。,c£,)
(3)将延拓到地面的波场歹,(z,z=O,∞)与检波器对应位置的值乘以频率域检波因子,即为该检波器接收到的波场值:
元(zg,∞)=歹1(zF,z—O,cc’)虿(zg,z,∞),
(11)
(4)将合成记录变换到时间域,得到以(z,£).固定新的检波点位置,重复步骤(3)和(4),可得到完整的共炮记录户(z,f).
对于转换波P—SV波的合成,只需将步骤(2)中向上延拓的波速用横波速度代替即可.而对于零偏移距正演记录,根据爆炸反射界面原理[2引,只需将步骤(1)中的新震源用地下反射系数代替,再进行步骤(2)和(3)的运算,即可获得.
4数值计算结果
图3为逆断层网格模型,介质泊松比为O.25,模型大小为5000
m×3000
m,离散网格间距为如
一如=10m,震源子波为主频25Hz的雷克子波,位于地表2500m处,采用的交错网格高阶有限差分法精度为时间2阶空间8阶,分别进行双程声波波动方程模拟和弹性波方程模拟,合成共炮记录如图4和图5所示.单程波模拟的时间采样率为lms,记录长度为2.5s,P—P波合成记录和P—SV波合成记录如图6所示.对比三种情况下的合成记录
万
方数据可以看出,弹性波波场成份丰富,较正确地反映了波在界面上的反射、透射等物理现象,符合地震波的传
播规律,图5中S波双曲线顶部出现较明显的频散现象,这主要是由于弹性波方程采用了与声波方程相同的离散参数所致,所以波场模拟需要考虑S波时,应当以S波频散最小为原则.图6单程波模拟结果不含直达波、多次波等干扰波,只含有与地质界面相对应的一次反射波,对比图5、6双曲线两翼可以发现,单程波波场能量主要集中在小炮检距范围内,这是因为单程波方程是双程波方程的近似,单程波算子存在一个最大倾角,对该倾角范围之外的波模拟误差较大.
图4双程声波共炮记录
Fig.4
Shotrecordsoftwo—wayacousticwaVeequation
图7为地层尖灭网格模型,模型大小为5000
m
×3000
m,离散网格间距为dz=如=5m,时间采样
率为1ms,记录长度为2.2s,震源子波为主频25Hz的雷克子波,双程声波方程和单程波方程自激自
收记录如图8所示.对比图8两种记录可见,单程波
2期皮红梅,等:波动方程数值模拟的j种方法及对比
395
方程记录各反射界面的反射波很明显,而双程声波受起伏界面影响,双程声波模拟产生的多次波将倾方程记录倾斜地层反射波模糊不清,这主要是由于
斜地层反射波掩盖所致.
图7地层尖灭网格模型
Fig.7
Stratapinching
model
万
方数据
396
地球物理学进展
24卷
图8自激自收记录
左:双程声波方程,右:单程声波方程
Fig.8
left:two.way
wave
Zeroo“setrecords
waVeequation
equation,right:one-way
Wang
X
M,ZhangH
in
L,WangD.Modeling
of
s幽mic
uging
a
wa、re
5
结论
propagation
order
hetemgeneous
poroelasticmedia
high-J.
基于不同波动方程数值模拟结果分析表明,双程波方程数值模拟的波场成份丰富,比较实际地反映了地震波在介质中传播的基本情况,但是当介质比较复杂时,多次波严重干扰了界面的一次反射波信息,波型识别困难.单程声波数值模拟结果可以清晰地反映地质分界面,但是受倾角限制的影响,波场能量主要集中在小炮检距范围内.因此对于复杂地质模型,三种模拟方法可以取长补短,综合应用.致谢感谢东方地球物理公司物探技术研究中心方法部田振平,卢秀丽等同志在论文完成过程中给予的热情帮助.参
考
文
[4]
sta昭eredfinitPdifferencemethod[J].
Cllinese
Geophys.(inChinese),2003,46(6),84Z~849.MadariagaR.
DynarIlicsof
an
expandingcircuIarfault[J].
BSSA,1976,66(3):639~666.
[5]DablainM
scalarwave
A.The
applicationof
high—order
dif{erencing
to
the
equation[J].Geophysics,1986,51(1)154~66.
[6]董良国,马在田,曹景忠.一阶弹性波方程交错网格高阶差分解
法稳定性研究口]..地球物理学报,2000,43(6):856~864.
Dong
L
G。MaZT,CaoJZ.A
hig卜order
study
on
stabilityoftheof
first-order
staggered—grideIasticwave
differencemethod
equation[J].chineseJ.Geophy3.(incldnese),
ZOOO。43(6):856~864.
[7]董良国,李培明.地震波传播数值模拟中的频散问题[J].天然
气工业,2004,24(6):53~56.
D0ng
L
G,LiPM.
Dispersive
problemin
sei8mic
WaVe
献(References):
[8]
propagationnumerical2004,24(6):53~56.
rIlode“ng[J].NaturalGas
Industry,
[1]董良国,马在田,曹景忠,王华忠,耿建华,雷兵,许世勇.一阶弹
性波方程交错网格高阶差分解法[J].地球物理学报,2000,43
(3):411~419.Dong
L
He
ZH,XiongG
Modeling
J,zhngYX.Nollzeroo“setSeismicBy∞e.wayAcou“cwave-equation[J].
Fo蝴8rd
SEGExpandedAbstracts,1998.
G,MaZT,CaoJZ,以nZ.Astaggered-gridhigh-
wave
[9]熊高君,贺振华,张琳,等.共炮记录正演模拟检波点下延记录
原理[J].石油物探,1999,38(2):43~49.
XiongG
orderdifferencemethodofone-orderelastic
equation[刀.
ChineseJ.Geophys.(inChinese),2000,43(3):411~419.
J,HeZH,ZhangL。甜口Z.Geophonerecordfrom
common-8hot
”cord
fo聊ard
[2]汪功礼,张庚骥,崔锋修,等.三维感应测井响应计算的交错罔
格有限差分法[J].地球物理学报,z003,46(4),561~567.Wang
G
modeling
[J].GeophysicaI
ProspectingforPetroleum,1999,38(2)143~49.
L,Zhang
g—d
GJ,CuiF
X,ct“.Applicationof
to
[10]熊高君,贺振华,黄德济.等.转换横波共炮记录单程声波方程
正演模拟[j].物探化探计算技术,1998,20(3):213~217.
Xiong
staggered
3一D
fillitedi“erencemethod
thecomputationof
induction1099ingresponse[J].ChineseJ.Geophys.(inGJ,HeZ
H,Huang
of
D
J,“以.(:ommon—shotrecords
one_way、张Ve
Chinese),2003,46(4):561~567.
forwardmodeling
convertedshea卜waVewith
techniques
[3]
王秀明。张海澜,王东.利用高阶交错网格有限差分法模拟地震波在非均匀孔隙介质中的传播[J].地球物理学报,2003,46
(6),842~849.
equation[J].(bmputiflg
for
geophysicaland
geochemicalexploration,1998,20(3):213~217.
[11]熊高君,贺振华,黄德济t等.多次反射可控制的波动方程共
万方数据
2期皮红梅,等:波动方程数值模拟的三种方法及对比
397
炮记录正演模拟口].石油地球物理勘探,1998,33(5):583~
590.
Xiong
G
J,HeZH,Huang
D
J,以以.WaVe-equation
common-shotseismogramforwardsimulation
with
mult卜
reflection
wave
constrained[J].0ilGeophysical
Pr08pecting,1998。33(5):583~590.
[12]刘札农,崔风林,张剑锋.三维复杂构造中地震波模拟的单程
波方法[J].地球物理学报,2004,47(3);514~520.
Liu
L
N.CuiF
L。ZhangJF.Seismic
nlodeling
with
one-
way、张ve
equationin
3D
complex
structures.ChineseJ.Geophys.(inChin舶e)t2004,47(3)i514~520.
[13]谢桂生,刘洪。李幼铭,等.界面起伏条件下反射/透射算子+
单程波方程的地震波模拟方法[J].地球物理学报,2005.48
(5):1172~1178.Xie
G
S.“uH,Li
Y
M,以以.Seismicmodeling
byrenection/transmissionoperator
and
one-way
wave
equation
underthecondition
of
fluctuating
reflectors[J].ChineseJ.
Geophys.(inChinese),2005.48(5):1172~1178.
[14]熊高君,贺振华,黄德济,等.改进的正演模拟定位原理[J].
石油地球物理勘探,1998,33(6):742~748.
XiongGJ,HeZH,HuangD
J,以“.1mproved
positioning
principle
of
fo唧ard
modeIing[J].oil
GeophysicalProspecting,1998,33(6):742~748.[15]贺振华。胡光岷,黄德济.数学检波器与波动方程地震叠前正演[J].成都理工大学学报(自然科学版),z004,3l(6):675~678.
HeZH,HuGM,Huang
D
J.Mathematicalgeophoneswithapplicationsto
seisITlic
prestackingforwardmodeling[J].
JouMnlof
ChengduUniversity
ofTechnology(Science&
万
方数据TechnologyEdition),2004,31(6):675~678.
[16]贺振华,熊晓军.等时叠加波动方程叠前正演[J].物探化探
计算技术,2005,27(3):194~200.
HeZH,XiongXJ.Seisnlicprestackingfo邢ard
modeling
based
on
equal-time
stack
principle
and
one-way
wave
equation[J].computingTechniquesforGeophysicaland
GeochemicalExploration,2005,27(3):194~200.
[17]牟永光,裴正林著.三维复杂介质地震数值模拟[M].北京:
石油工业出版社,2004.
Mu
Y
G,PeiZL.Seismicnumericalmodeling
for3-d
complex
media[M].BeijingfPetroleumIndustryPr船s,2004.
[18]V试euxJ.P—sVwave
propagation
in
heterogeneousmedia:
Velocity'stressfinite-differencemethod[J].
Geophysics。
1986,51(4);889~901.
[19]L_evander
A
R.
Fourtbrorderfinite_differenceILsV
seismograms[J].Geophysics,1988,53(11):1425~1436.
[20]马在田著.地震成像技术一有限差分法偏移[M].石油工业出
版社.1989.
MaZ
T.
Seismicimagingtechniques-Migration
withfinite
difference
method[M].Beijing:PetroleumIndustry
Press,
1989.
[21]
RistowD,RuhlT.Fourier“nite
difference
migration[J].
Geophysics,1994,59(12):1882~1893.
[z2]贺振华,王才经,李建朝,等.反射地震资料偏移处理与反演
方法[M].重庆大学出版社.1989.
He
Z
H,Wang
CJ,“J
C,甜“.Migrationand
conVersion
methodsfor
reflectionseismic
data[M].chongqing
University
Press,1989.
第24卷第2期
2009年4月(页码:391~397)
地球物理
IN
学进展
v01.24,No.2
Apr.2009
PROGRESSGEOPHYSICS
皮红梅,蒋先艺,刘财,等.波动方程数值模拟的三种方法及对比.地球物理学进展,2009,24(2):39l~397,DoI:10.3969/j.
issn.1004—2903.2009.02.004.PiH(in
M,JiangXY,LiuC,Pf口Z.Threemethodsofwaveequationforwardmodelingandcompari50ns.PrDg化5s锄GPo户^∥.
Chinese),2009,24(2):391~397,DoI:lO.3969/j.issn.1004—2903。2009.02.004.
波动方程数值模拟的三种方法及对比
皮红梅1’2,
蒋先艺2,
刘
财¨,
王成祥2,
姜绍辉2
(1-吉林大学地球探测科学与技术学院。长春130026}
2.中国石油东方地球物理公司物探技术研究中心,涿州072750)
摘要波动方程数值模拟方法是研究地震渡场传播的一种重要手段,本文采用交错网格高阶有限差分方法分别对
双程声波方程和双程弹性波方程进行了波场数值模拟,并且根据定位原理采用傅立叶有限差分算子进行了单程波方程数值模拟,在分析定位原理的基础上,对其计算过程稍作修改,将延拓到地面的波场直接由每个检波点接收,无需横向叠加过程,得到了单程声波方程共炮记录.基于不同波动方程的数值模拟结果表明,双程波方程结果包含直达波、多次波等干扰渡,信噪比低;单程波数值模拟结果只包含了介质分界面的一次反射波,信噪比高,但对于大角度入射渡误差较大,并且对于同一个地质模型而言,双程弹性波方程计算速度最慢,双程声波方程次之,单程声波方程计算速度最快.因此对于复杂地质模型,三种模拟方法可以取长补短,综合应用.关键词波动方程数值模拟,交错网格,高阶有限差分法,定位原理,单程波
DOI:lO.3969/j.issn.1004—2903.2009.02.004
中图分类号
P315
文献标识码
A
ThreemethodsofwaVeequationforwardmodelingandcomparisons
PIHong—meil”,
(1.CDzZFgP
JIANGXian—yi2,
o厂GP08zpZor乜£io靠Sc{P,lce
LIUCai¨,
n行d
WANG
Cheng—xian92,
JIANGShao_hui2
T■f_Il雅ozogy,.,订i以UhitJP,百f£,,(冼口”gc|Il配种130026,c^f,l口;
2.R#s∞”九&DP御ZD户m£斗£CP咒fPro,BGP,CNPC,Z^“oz^o“072750,Chi靠口)
Abstr扯tWave
equationforwardmodelingis
on
averyimportantway
to
studytheseismicpropagationunderground.In
thispaper,seismicforwardmodelingbaseddone
using
the
staggered—grid
on
two-wayacousticwaveequationandtw0-wayelasticwaVeequationis
method.made
Besides,
according
t0
high_order
differencethe
positioningprinciple,
Inthe
simulationbased
progress
one-wayacousticwaveequationis
bytheFourierfinitedifference
propagator.
to
ofcomputing,eachgeophonedirectlyreceivesthewavefieldthathasbeenextrapolated
Thethreesimulation
resultsshow
thesurfaceofthe
generate
a
modelwithoutlateralstacking.waveand
thattwo.waywa、reequations
directive
muItiplewavewhichreducetheSNRseverely,whiletheone—waywaveequationmethodgives
createscosts
higherSNR
withonlyrenectedwaye,butitelasticwaveequation
larger
errors
withwide-an一eincidentwa’re.Forthesamemodel,two—way
costs
algorithm
the
mostcomputationaltime,one_wayacousticwaveequationalgorithm
algorithminbetween.Sofor
a
theleastandtw0-wayacousticwaveequationmodelingKeywords
complex
structuremodel,
thethree
methods
can
beusedsimultaneouslywithdifferentadvantagesanddisadvantages.
high.orderdifferencemethod,positioningprinciple,
waveequationforwardmodeling,staggeredgrid,
one.waywaveequation
收稿日期作者简介
2008_06一12;
修回日期
2008一09-2z.
皮红梅,1981年生,博士研究生,2003年毕业于吉林大学地球物理专业,获理学学士学位,现在该校地球探测科学与技术学院攻读固体地球物理学专业博士学位,研究领域为复杂地震波场数值模拟.(E—mail.brook_3993@email.jlu.edu.cn)
通讯作者刘财,吉林大学教授,博士生导师,主要从事固体地球物理研究.(E-mamliucaUlu@163.com)
万方数据
392
地球物理学进展
24卷
O
引言
算效率低,正演记录的信噪比低.贺振华、熊高君(1998)等提出检波点下延原理[8 ̄儿],首次采用单程声波方程md33将炮点和检波点的地震波场同时从地面向地下进行延拓得到共炮正演记录,但该方法是从偏移的思路来进行正演,计算效率较低,于是又提出一种能快速合成共炮记录的固定检波点定位原理n引,此后,贺振华、胡光岷(2004)等人提出的数学检波器的概念[15]和贺振华、熊晓军(2005)提出的等时叠加原理[161,从理论上对单程波动方程叠前正演方法进行了归纳和总结.本文在分析固定检波点定位原理的基础上,对其计算过程稍作修改,将延拓到地面的波场直接由每个检波点接收,无需横向叠加过程,得到了单程声波方程共炮记录,对三种波动方程得到的共炮记录进行分析对比.最后进一步分析了双程声波方程和单程声波方程的自激自收记录.
1
地震数值模拟方法是地震勘探和地震学的重要基础,是根据已知的地下介质结构和相应物理参数模拟研究地震波在地下各种介质中的传播规律,并计算在地面或地下各观测点的地震记录的一种地震模拟方法.数值模拟方法主要有三大类:波动方程法,积分方程法和射线追踪法.波动方程数值解法是建立在以弹性或粘弹性理论和牛顿力学为基础的双曲型偏微分方程一地震波传播方程的理论基础上的.由于地下介质性质不同,其相应的地震波传播方程也不同.如声学介质中的声波波动方程;弹性介质中的弹性波波动方程;粘弹性介质中的粘弹性波波动方程;孔隙弹性介质(双相或多相介质)中的双相或多相介质弹性波方程;各向异性介质中的各向异性弹性波波动方程等.积分方程法是建立在以惠更斯原理为基础的波叠加原理基础上的,其数学表达形式为波动方程的格林函数域积分方程式和边界积分方程式.射线追踪法是建立在以射线理论为基础的波动方程高频近似理论基础上的,其数学表达形式为程函方程和传输方程.
波动方程法又分为双程波方程法和单程波方程法.本文采用交错网格高阶有限差分法对双程声波方程和双程弹性波方程进行数值模拟.交错网格高阶有限差分法[1 ̄3]是董良国(2000)首先提出的,他将Mad耐aga(1976)提出的交错网格法[43与Dablain(1986)提出的高阶差分方法[5]有机结合,用于求解
双程声波方程数值模拟
直角坐标系下,各向同性介质中的二维声波方
程的一阶压力一速度方程m3形式为:
警一一K(篆+碧),
垫一一上丝.
a£
p
az’
婺一一土挈.
at
p
<1)
、1
az’
式中,户一户(z,z,£)为声波波场,%=砚(z,2)为质点振动速度的水平分量,让一皱(z,名)为垂直分量,K一即2为弹性模量,口为介质的纵波速度,lD为介质密度.
采用交错网格(如图1所示)高阶有限差分方法对方程(1)进行数值离散,其2阶时间2N阶空间的离散格式为:
横向各向同性介质一阶速度一应力弹性波方程,其
后他又对该方法的稳定性问题[6](2000)和频散问题"](2004)进行了研究.双程波方程法包含了地震波场的全部信息,可以处理复杂的地质构造,但是计
1.....................................................一——..——————一..——.——..
矿%∽=一∽沪警{塾Ⅲ+竽,J)一谚(z一竿,歹)]
十塾[破(奶+学)一破(叫一竽)]},
谚
l
+
1—2
一记
.Z
+
1—2
一
妻“[户蚪}(i+咒,歹)一户外专(i+1一行,歹)],(2)
矿
l
}
1—2一
破
.£
}
1—2一
其中:∥+专(i,_『)表示户(z,z,£+譬),
万方数据
塑一
Ⅳ
惫箭一∞一∞
∑c。[户蚪}(i,歹+咒)一夕蚪}(z,歹+1一咒)].
^=1
谚(f+丢,_f)表示V。(z+等,2,t),其它与此类
2期皮红梅,等:波动方程数值模拟的三种方法及对比
393
同.厶为交错网丰子差分系识.
2双程弹性波方程数值模拟
均匀各向同性介质中的二维弹性波方程m3为:
.o等=鲁+等,P警=鲁+等,
鲁=(A+2p)卺+|;I卺,
鲁__(A+刎凳+A警,鲁=p(等+卺).
其中,让和钞。分别为质点振动速度的水平分量和垂直分量;%,%,屹为应力分量;lD为介质密度;A和口为拉梅常数.
■px
1,
●vx’,p
∥
o屹,形
图1声波方程交错网格示意图
Fig.1
Thestaggeredgriddacousticwaveequatio矗
采用交错网格(如图2所示,据Levander[193)高阶有限差分法来求解方程(3).
图2
弹性波方程交错网格示意图(据Levander1988)
F螗.Z
Thestaggeredgridofelastic
waveequation(Levander1988)
万
方数据呓{=吒吾+考惫*[蚤“《字。,一P:2孚≯
+∑c。(《.升学一≤,卜学)],
’(4)
咙U2,尸:。/2.』,s:一。,。分别为质点速度分量饥和应
力分量盯。,%的离散值,其它几个分量的离散方法与此类似.
3单程声波方程数值模拟
二维各向同性均匀介质中的单程波方程田3为:
筹一土i志庐=士i√等一磋夕,,
(5)
歹=声(走。,名,甜)为频率一波数域声波波场,“+”、“一”号分别代表上行波方程和下行波方程.则根据D.Ristow(1994)提出的傅立叶有限差分‘213(FFD)法可得用于正演的波场延拓算子为:
f声(忌。,锰l,c£,)一声(惫:,≈,(c,)e_‘√孟巧△f
{石(z,铣l,叫)一万(z,巧,(c,)fb血b∞缸
.
(6)
【户(z,z灶1,(c’)=夕(z,≈,cE,)r巩缸
其中:石=石(z,2,∞)为频率一空间域声波波场,c为
每一个延拓深度层内的背景速度,取为:
c=min(口(z,2。)),
纵…,一志一÷,
A_钟吲【蘸J’
f生竺]
6=丢[(詈)2+詈+-].
根据贺振华,熊高君提出的固定检波点定位原理[H3合成共炮记录,对其计算过程稍作修改,将延拓到地面的波场直接由每个检波点接收,无需横向叠加过程.原理表述如下:固定炮点和检波点,用下行波方程将震源波场从炮点位置延拓到反射界面,与所有反射点的位置做相关运算,相当于模拟入射波沿反射界面激发新震源或反射波,把这个新震源暂时置于所产生的位置,同时把原来的入射波继续向下延拓直到最深处,形成所有深度处的新震源,再从最深处用上行波方程将新震源向上延拓,每延拓一步,将延拓波场与该处的新震源波场叠加,将叠加后的波场继续向上延拓,直到检波器所在位置,将该延拓波场与检波器谐振因子进行相关,该相关结果就是所固定检波点的正演结果.
记:地面记录59(z,£),震源波场s(z,2,f),检
其2阶时间、2N阶空间的离散格式为:
地球物理学进展
24卷
波谐振波场g(z。,zo,t),反射系数r(z,z).其中:
s(z,州)一F五嘞一几z一五
(7)
lO,z≠z,
(z.,勐)、(z。,‰)分别为震源位置和检波点位置,本论文中震源和检波因子都用ricker子波.将震源和检波因子都变换到频率域i(z,z,cc,)、虿(以,孙,甜),并且反射系数拓展为频率不变函数尸(z,名,∞),即
尹(z,z,∞)=r(z,2).
固定检波点原理P-P波合成记录的步骤如下:(1)根据(6)式将震源波场向下延拓,直到地下最深处,并且与延拓深度上的反射系数进行相关,形成新震源:
j(z,z件1,∞)=j(z,z,,(£’)e一让。血,
(8)西(z,≈+l,叫)一j(z,乃+1,(£,)尹(z,影件1,叫).
(9)
(2)根据(6)式从地下最深处将反射界面上的所有新震源向上延拓,并与该延拓层新震源叠加,直到将波场延拓到检波器所在的深处(假设地面z=0处).
歹1(z,争1,c£,)=歹l(z,乃,倒)e-让r缸+而(z,夸l,(cJ)
(10)
其中歹l(z,名。,叫)一鄢(z,名。,c£,)
(3)将延拓到地面的波场歹,(z,z=O,∞)与检波器对应位置的值乘以频率域检波因子,即为该检波器接收到的波场值:
元(zg,∞)=歹1(zF,z—O,cc’)虿(zg,z,∞),
(11)
(4)将合成记录变换到时间域,得到以(z,£).固定新的检波点位置,重复步骤(3)和(4),可得到完整的共炮记录户(z,f).
对于转换波P—SV波的合成,只需将步骤(2)中向上延拓的波速用横波速度代替即可.而对于零偏移距正演记录,根据爆炸反射界面原理[2引,只需将步骤(1)中的新震源用地下反射系数代替,再进行步骤(2)和(3)的运算,即可获得.
4数值计算结果
图3为逆断层网格模型,介质泊松比为O.25,模型大小为5000
m×3000
m,离散网格间距为如
一如=10m,震源子波为主频25Hz的雷克子波,位于地表2500m处,采用的交错网格高阶有限差分法精度为时间2阶空间8阶,分别进行双程声波波动方程模拟和弹性波方程模拟,合成共炮记录如图4和图5所示.单程波模拟的时间采样率为lms,记录长度为2.5s,P—P波合成记录和P—SV波合成记录如图6所示.对比三种情况下的合成记录
万
方数据可以看出,弹性波波场成份丰富,较正确地反映了波在界面上的反射、透射等物理现象,符合地震波的传
播规律,图5中S波双曲线顶部出现较明显的频散现象,这主要是由于弹性波方程采用了与声波方程相同的离散参数所致,所以波场模拟需要考虑S波时,应当以S波频散最小为原则.图6单程波模拟结果不含直达波、多次波等干扰波,只含有与地质界面相对应的一次反射波,对比图5、6双曲线两翼可以发现,单程波波场能量主要集中在小炮检距范围内,这是因为单程波方程是双程波方程的近似,单程波算子存在一个最大倾角,对该倾角范围之外的波模拟误差较大.
图4双程声波共炮记录
Fig.4
Shotrecordsoftwo—wayacousticwaVeequation
图7为地层尖灭网格模型,模型大小为5000
m
×3000
m,离散网格间距为dz=如=5m,时间采样
率为1ms,记录长度为2.2s,震源子波为主频25Hz的雷克子波,双程声波方程和单程波方程自激自
收记录如图8所示.对比图8两种记录可见,单程波
2期皮红梅,等:波动方程数值模拟的j种方法及对比
395
方程记录各反射界面的反射波很明显,而双程声波受起伏界面影响,双程声波模拟产生的多次波将倾方程记录倾斜地层反射波模糊不清,这主要是由于
斜地层反射波掩盖所致.
图7地层尖灭网格模型
Fig.7
Stratapinching
model
万
方数据
396
地球物理学进展
24卷
图8自激自收记录
左:双程声波方程,右:单程声波方程
Fig.8
left:two.way
wave
Zeroo“setrecords
waVeequation
equation,right:one-way
Wang
X
M,ZhangH
in
L,WangD.Modeling
of
s幽mic
uging
a
wa、re
5
结论
propagation
order
hetemgeneous
poroelasticmedia
high-J.
基于不同波动方程数值模拟结果分析表明,双程波方程数值模拟的波场成份丰富,比较实际地反映了地震波在介质中传播的基本情况,但是当介质比较复杂时,多次波严重干扰了界面的一次反射波信息,波型识别困难.单程声波数值模拟结果可以清晰地反映地质分界面,但是受倾角限制的影响,波场能量主要集中在小炮检距范围内.因此对于复杂地质模型,三种模拟方法可以取长补短,综合应用.致谢感谢东方地球物理公司物探技术研究中心方法部田振平,卢秀丽等同志在论文完成过程中给予的热情帮助.参
考
文
[4]
sta昭eredfinitPdifferencemethod[J].
Cllinese
Geophys.(inChinese),2003,46(6),84Z~849.MadariagaR.
DynarIlicsof
an
expandingcircuIarfault[J].
BSSA,1976,66(3):639~666.
[5]DablainM
scalarwave
A.The
applicationof
high—order
dif{erencing
to
the
equation[J].Geophysics,1986,51(1)154~66.
[6]董良国,马在田,曹景忠.一阶弹性波方程交错网格高阶差分解
法稳定性研究口]..地球物理学报,2000,43(6):856~864.
Dong
L
G。MaZT,CaoJZ.A
hig卜order
study
on
stabilityoftheof
first-order
staggered—grideIasticwave
differencemethod
equation[J].chineseJ.Geophy3.(incldnese),
ZOOO。43(6):856~864.
[7]董良国,李培明.地震波传播数值模拟中的频散问题[J].天然
气工业,2004,24(6):53~56.
D0ng
L
G,LiPM.
Dispersive
problemin
sei8mic
WaVe
献(References):
[8]
propagationnumerical2004,24(6):53~56.
rIlode“ng[J].NaturalGas
Industry,
[1]董良国,马在田,曹景忠,王华忠,耿建华,雷兵,许世勇.一阶弹
性波方程交错网格高阶差分解法[J].地球物理学报,2000,43
(3):411~419.Dong
L
He
ZH,XiongG
Modeling
J,zhngYX.Nollzeroo“setSeismicBy∞e.wayAcou“cwave-equation[J].
Fo蝴8rd
SEGExpandedAbstracts,1998.
G,MaZT,CaoJZ,以nZ.Astaggered-gridhigh-
wave
[9]熊高君,贺振华,张琳,等.共炮记录正演模拟检波点下延记录
原理[J].石油物探,1999,38(2):43~49.
XiongG
orderdifferencemethodofone-orderelastic
equation[刀.
ChineseJ.Geophys.(inChinese),2000,43(3):411~419.
J,HeZH,ZhangL。甜口Z.Geophonerecordfrom
common-8hot
”cord
fo聊ard
[2]汪功礼,张庚骥,崔锋修,等.三维感应测井响应计算的交错罔
格有限差分法[J].地球物理学报,z003,46(4),561~567.Wang
G
modeling
[J].GeophysicaI
ProspectingforPetroleum,1999,38(2)143~49.
L,Zhang
g—d
GJ,CuiF
X,ct“.Applicationof
to
[10]熊高君,贺振华,黄德济.等.转换横波共炮记录单程声波方程
正演模拟[j].物探化探计算技术,1998,20(3):213~217.
Xiong
staggered
3一D
fillitedi“erencemethod
thecomputationof
induction1099ingresponse[J].ChineseJ.Geophys.(inGJ,HeZ
H,Huang
of
D
J,“以.(:ommon—shotrecords
one_way、张Ve
Chinese),2003,46(4):561~567.
forwardmodeling
convertedshea卜waVewith
techniques
[3]
王秀明。张海澜,王东.利用高阶交错网格有限差分法模拟地震波在非均匀孔隙介质中的传播[J].地球物理学报,2003,46
(6),842~849.
equation[J].(bmputiflg
for
geophysicaland
geochemicalexploration,1998,20(3):213~217.
[11]熊高君,贺振华,黄德济t等.多次反射可控制的波动方程共
万方数据
2期皮红梅,等:波动方程数值模拟的三种方法及对比
397
炮记录正演模拟口].石油地球物理勘探,1998,33(5):583~
590.
Xiong
G
J,HeZH,Huang
D
J,以以.WaVe-equation
common-shotseismogramforwardsimulation
with
mult卜
reflection
wave
constrained[J].0ilGeophysical
Pr08pecting,1998。33(5):583~590.
[12]刘札农,崔风林,张剑锋.三维复杂构造中地震波模拟的单程
波方法[J].地球物理学报,2004,47(3);514~520.
Liu
L
N.CuiF
L。ZhangJF.Seismic
nlodeling
with
one-
way、张ve
equationin
3D
complex
structures.ChineseJ.Geophys.(inChin舶e)t2004,47(3)i514~520.
[13]谢桂生,刘洪。李幼铭,等.界面起伏条件下反射/透射算子+
单程波方程的地震波模拟方法[J].地球物理学报,2005.48
(5):1172~1178.Xie
G
S.“uH,Li
Y
M,以以.Seismicmodeling
byrenection/transmissionoperator
and
one-way
wave
equation
underthecondition
of
fluctuating
reflectors[J].ChineseJ.
Geophys.(inChinese),2005.48(5):1172~1178.
[14]熊高君,贺振华,黄德济,等.改进的正演模拟定位原理[J].
石油地球物理勘探,1998,33(6):742~748.
XiongGJ,HeZH,HuangD
J,以“.1mproved
positioning
principle
of
fo唧ard
modeIing[J].oil
GeophysicalProspecting,1998,33(6):742~748.[15]贺振华。胡光岷,黄德济.数学检波器与波动方程地震叠前正演[J].成都理工大学学报(自然科学版),z004,3l(6):675~678.
HeZH,HuGM,Huang
D
J.Mathematicalgeophoneswithapplicationsto
seisITlic
prestackingforwardmodeling[J].
JouMnlof
ChengduUniversity
ofTechnology(Science&
万
方数据TechnologyEdition),2004,31(6):675~678.
[16]贺振华,熊晓军.等时叠加波动方程叠前正演[J].物探化探
计算技术,2005,27(3):194~200.
HeZH,XiongXJ.Seisnlicprestackingfo邢ard
modeling
based
on
equal-time
stack
principle
and
one-way
wave
equation[J].computingTechniquesforGeophysicaland
GeochemicalExploration,2005,27(3):194~200.
[17]牟永光,裴正林著.三维复杂介质地震数值模拟[M].北京:
石油工业出版社,2004.
Mu
Y
G,PeiZL.Seismicnumericalmodeling
for3-d
complex
media[M].BeijingfPetroleumIndustryPr船s,2004.
[18]V试euxJ.P—sVwave
propagation
in
heterogeneousmedia:
Velocity'stressfinite-differencemethod[J].
Geophysics。
1986,51(4);889~901.
[19]L_evander
A
R.
Fourtbrorderfinite_differenceILsV
seismograms[J].Geophysics,1988,53(11):1425~1436.
[20]马在田著.地震成像技术一有限差分法偏移[M].石油工业出
版社.1989.
MaZ
T.
Seismicimagingtechniques-Migration
withfinite
difference
method[M].Beijing:PetroleumIndustry
Press,
1989.
[21]
RistowD,RuhlT.Fourier“nite
difference
migration[J].
Geophysics,1994,59(12):1882~1893.
[z2]贺振华,王才经,李建朝,等.反射地震资料偏移处理与反演
方法[M].重庆大学出版社.1989.
He
Z
H,Wang
CJ,“J
C,甜“.Migrationand
conVersion
methodsfor
reflectionseismic
data[M].chongqing
University
Press,1989.