地震波正演模拟的任意差分精细积分方法_王润秋
2003年6月
石油地球物理勘探
第38卷 第3期
·处理方法·
地震波正演模拟的任意差分精细积分方法
王润秋
耿伟峰 王尚旭
摘 要
(石油大学资源与信息学院CNPC物探重点实验室)
王润秋,耿伟峰,王尚旭.地震波正演模拟的任意差分精细积分方法.石油地球物理勘探,2003,38(3):252～257
本文对钟万勰提出的抛物线方程单点子域积分法和强士中提出的温度场问题的任意差分法进行了改进,提出了波动方程任意差分精细积分解法,从理论上推导了任意差分精细积分公式并给出了其实现途径,对任意差分精细积分方法的空间与时间精度和差分格式的稳定性做了理论探讨,并进行了验证。理论分析和实际例证都表明本文方法精度较高。地震波正演模拟算例则表明本文方法适用于复杂地表和复杂构造地质体。关键词 正演模拟 波动方程 有限差分 任意差分 精细积分 复杂地表 地震模型
引 言
地震波数值模拟是波动方程的正演问题,是在地下介质结构和参数已知的情况下,利用数值计算的方法来研究地震波在地下介质中的传播规律,从而获得理论上地震记录的一种技术。目前,国内外在地震勘探生产实践中所采用的地震模型数值模拟方法是以射线理论(射线追踪法)和波动方程数值解法(克希霍夫积分法、有限差分法、有限单元法、边界单元法和虚谱法)为基础的。虽然当前的地震波数值模拟有许多方法,但在振幅保真、数值频散、模拟精度、边界吸收、计算速度方面还有待改进。本文所提出的波动方程任意差分精细积分解法,是根据钟万勰[1]提出的抛物线方程单点子域积分法、强士中提出的温度场问题的任意差分法进行改进而得到的。理论推导和实际算例均表明,波动方程的任意差分精细积分法(theArbi-traryDifferencePreciseIntegrationMethod,简称为ADPI法)是一种高精度地震正演方法。
[2]
=++W(x-x0,z-z0)f(t)v t x zu(x,z,0)=0
(1)
=0
t
式中:u(x,z,t)为波场函数,0
首先,对式(1)进行二维空间坐标离散(按均匀网格划分),即在任一离散点j(x,z)(点序号)处,有
222
j j j
=++v t x z(x,z)=(xj,zj) +W(xj-x0,zj-z0)f(t)
222
(2)
其次,将u在j点的x方向邻域做Taylor展开
m
ui=uj+
∑
k=1
k
Δxiuk! x
m+1
x=x
+
j
+O[(Δxi)]
上式中Δxi=xi-xj,i是j点的x邻域内的任意点,对上式移项并加权得
n
n
∑T(u
i
i=1
i
-uj)=
m
∑T×
i
i=1
x=x
方法原理
对于二维波动方程
本文于2002年9月4日收到,修改稿于2003年2月18日收到。
(02300)及 ×
∑
k=1
k
Δxiuk! x
+O[(Δxi)m+
j
1
(3)
其中:Ti是加权系数;m原则上可任意选取,在此一
第38卷 第3期 王润秋等:地震波正演模拟的任意差分精细积分方法
253
般取m=n(以m=n为例来展开讨论),则
m
c1,c2,进而可以求出tn+1时刻的值,u0的值可由初值确定。故得计算公式如下uj,n+
1
∑
m
i=1
TiΔxi=0
i
i
=2uj,ncos(aΔt)-uj,n-1+∑T(Δx
i=1m
)=2v)=0)=0
43
22
nn
×2-aa4
(7a)
∑T(Δx
i
i=1m
i
∑T(Δx
i
i=1
(4)
i
m
n2
2(Δt)a
b02d0
uj,1=uj,0cos(aΔt)+-×
aa
×[1-cos(aΔt)]+(Δt)2
a
uj,0=0
×[1-cos(aΔt)]+为初始条件。
(7b)(7c)
式(7a)为波场函数递推求解公式,式(7b)、式(7c)均
i
∑T(Δx
i
i=1m
)m-1=0)m=0
∑T(Δx
i
i=1n
i
式(4)是一个关于Ti的m维线性方程组,解之可求
2
j
出加权系数Ti,Ti(ui-uj)=。∑ xi=1
同理,可将u在j点的z方向邻域做Taylor展开。
任意差分精细积分方法的误差与稳定性分析
为了精确对比ADPI方法的误差和稳定性分析,本文用已知精确解的强迫振动问题来说明AD-PI方法的误差分析结果,我们用数值例证来分析二维稳定性。误差分析
对于强迫振动问题
=D+Ex(L-x)+Fu t x2
2
最后,式(2)右端的二阶导数可用相邻点的任意
差分来表示(暂不考虑震源),那么式(2)可表示为
2 j
= t
2m
∑
Ti(ui-uj)
2m
i=1
上式可写为
2 j
+uj∑Ti= ti=1
2m
∑Tu
i
i=1ii,n
i
(5)
u(x,0)=0=0 t
u(0,t)=0u(L,t)=0
式中:D=89.66;E=5000;F=3000;L=0.4572。该问题的精确解为
u(x,t)=v(x,t)+U(x)
式中
Z
f(a)ddZ-0a0x
Z
-00f(a)dda∞
v(x,t)=∑an+
Ln=1
+bnLL
其中
U(x)=
L
(8)
将式(5)的右端项展开到二阶,即
∑Tu
i=1
2m
ii
=
∑Tu
i
i=1
2m
i,n
+
∑Tu
i=1
2m
(t-tn)+
+
2
iui,n(t-tn)T∑2i=1
2m
这时将i点看成参数,对式(5)求解析解(即精细积分)得
uj=c1exp[ia(t-tn)]+ +c2exp[-ia(t-tn)]+
n
+a2
nnnnn
++-422
aaa
其中
2m
2m
i
n
(6)
L
a= en=
2
b∑T
i=1
2m
=
∑Tu
i=1
2m1
ii,n
∑Tu
i=1
ii,n
dn=
iui,nT2∑i=1
],n-1n a=sqrt(D)
f)=(-xFu
n-,n+
254
石油地球物理勘探2003年
L an=H(x)-U(xdx
L0LL bn=J(x)dx
nπa0L
由式(8)1;常规的有限差
PI方法准确得多,且稳定性更好;当四点二阶ADPI方法时间步长加大10倍以后,所得的结果仍比四点零阶ADPI方法准确(见图2)。
分解见表2;四点零阶ADPI方法解见表3;四点二阶ADPI方法解见表4。
表1 由式(8)解得的精确解
t0.10.30.50.70.91.0
x=110.1718330.1284330.064353
x=310.4017810.3002050.150264
x=510.3658570.2733510.1368190.0278290.0034620.34959
x=710.0878430.0656720.0329210.0067640.0009190.084028
图1 误差对比图
虚线为有限差分法与精确解的误差;
实线为四点零阶ADPI方法与精确解的误差
0.01318170.0305330.00175740.0037190.164346
0.383817
表2 常规的有限差分解(Δt=0.000005s)t0.10.30.50.70.91.0
x=110.1719480.1287080.0647730.0135010.0016630.164879
x=310.4020480.3008450.1512160.0312050.0034010.384961
x=510.3661010.2739400.1376910.0284460.0031810.350652
x=710.0879000.0658140.0331360.0069310.0008790.084309
图2 时间步长加大10倍后的误差对比图
虚线为四点零阶ADPI方法与精确解的误差;
实线为四点二阶ADPI方法与精确解的误差
ADPI稳定性分析
为了与有限差分法对比,假设Δx1=Δx2=Δx,则可得出
1=T2= T=T(Δx)
bn=T(uj-1,n+uj+1,n) a2=2T
对于有限差分格式,由文献[3]可知其稳定性条件为。Δx≤1
本文采用VonNeumann方法来分析四点零阶ADPI方法的稳定性,我们将误差的傅氏级数的一个分量表示为
Wuj,n=Vnexp(ikx)
Wuj,n-1=Vn-1exp(ikx)Wuj,n+1=Vn+1exp(ikx)Wuj-1,n=Vnexp[ik(x-Δx)]Wuj+
1,n
2
表3 四点零阶ADPI方法解(Δt=0.000005s)t0.10.30.50.70.91.0
x=110.1718340.1284380.0643620.0131890.0017540.164357
x=310.4017820.3002170.1502860.0305490.0037110.383841
x=51
x=71
0.3658590.08784330.2733630.06567460.1368390.03292540.0278450.00676760.0034550.00091780.3496120.0840339
表4 四点二阶ADPI方法解(Δt=0.00005s)t0.10.30.5
0.70.91.0
x=110.1718340.1284330.0643540.0131860.0017570.164348
x=310.4017810.300207
0.1502670.0305350.0037180.383821
x=510.365857
x=710.087843
0.2733530.06567230.1368210.03292120.0278320.00676440.0034610.00091910.3495930.0840291
(9)
=Vnexp[ik(x+Δx)]
式(7)的误差传播方程为
Wuj,n+1=2Wuj,ncos(aΔt)-Wuj,n-1+ +2
n
[1-cos(aΔt)](10)
通过表1～表4解的对比和对图1所示的误差曲线的分析可知,四点零阶ADPI方法比常规有限,
第38卷 第3期 王润秋等:地震波正演模拟的任意差分精细积分方法
255
将式(9)代入式(10)得到 Vn+
1
=2Vncos(aΔt)-Vn-1+
(11)
+2Vncos(kΔx)[1-cos(aΔt)]
假设误差传播因子为G,则Vn+1=GVn。将其代入式(11),则得到关于G的方程
G--2
+cos(kΔx)- Δx
cos(kΔx)G+1=0(12)Δx
图5 dt=0.00075s,
=1.165时,二阶ADPI方法的结果Δx
式(12)中G的模小于等于1的充要条件是
即
2nπ-π/2≤≤2nπ+
Δx2
n=0,1,2,…
因此,式(7)的稳定性条件为≤
稳定。
但是,对ADPI的二阶方法,
采用VonNeu-mann方法来分析太复杂,本文仅用数值例子说明其稳定性比其他两种差分格式更好。图3展示了强≈2
2
1.11089,与常规有限差分法相比,ADPI方法更为
≥0Δx
迫振动问题的精确解结果。当dt=0.00075s,而ADPI的二阶方法(图5)仍能保持稳定。
=Δx
1.165时,有限差分法(图4)的计算结果极不稳定,
地震模型正演实例
用ADPI方法的二阶展开差分格式进行地震资料正演,首先,给出正演所使用的震源子波与边界条件。
正演所使用的震源子波与边界条件为
f(t)=(t-t0)exp[-a(t-t0)]
通过试验选取合适的参数,具体使用的子波及
2
图3 强迫振动问题的精确解
图6 子波波形
图4 dt=0.00075s,
vΔt
=1.165时,有限差分方法的结果Δx
图7 子波的振幅谱
256
石油地球物理勘探2003年
振幅谱为图6与图7所示。
在正演过程中,使用的边界条件为单程波动方程,对二维正演,具体格式为
2 x t2 x t
x=x
22
-bv2=02
v t z
第二个模型见图9,模型为断层与背斜构造。
其观测系统为:共放36炮,双边接收,道间距为10m,炮间距为100m,正演所划分的网格为:10m×
(13)(14)(15)
=0
min
+
max
10m,第一炮在2250m处,地表放炮,地表接收,满
覆盖次数为10次。
单炮记录剖面见图10。
在Omega处理系统上对其做水平叠加(见图11)。如图11所示,叠加剖面很好地反映了地下构造情况。观察绕射波与断面波的情况,可以看出绕射波清晰,断面波发育。
x=x
min
22
-=0+bvv t z
z t
2
z t
2
z=z
+
max
-bv=0v t x
z=z
22
z=z
min
22-2+bv2v t z
(16)
式中:a=10/(3π);b=3/(3π)。
上式是在单程波动方程的基础上发展而来。亦即对变换域的平方根的多项式展开时不采用Taylor方
法(式(17)),而采用优化方法展开(式(18))。
2
1-x≈1-2
(17)
2
1-x≈-(18)
3π3c
为了验证正演方法的正确性,本文列举了三个二维模型
。
图9 断层与背斜构造的地质模型
图8 中心点震源的波场传播的时间切片
(a)中心点震源均匀模型;(b)100ms的时间切片;
(c)200ms的时间切片;(d)300ms的时间切片;
图10 地质模型的正演炮集
(a)2250m的正演炮集;(b)2750m的正演炮集;(c)3250m的正演炮集;(d)4250m的正演炮集
对匀速模型中心点震源的波场传播进行动态扫描分析(图8),可以看出,ADPI方法是正确的,子波
第三个模型见图12,模型为复杂地表与断层构造。在我国西部进行的地震勘探中,经常遇到复杂地
第38卷 第3期 王润秋等:地震波正演模拟的任意差分精细积分方法
257
的勘探是很有意义的。针对西部情况,设计一个复杂地表模型,即为高速岩层出露地表的情况
。
图11 地质模型的叠
加剖面
图13 地质模型的正演炮集
(a)2050m的正演炮集;(b)2450m的正演炮集;
(c)2850m的正演炮集;(d)3250m的正演炮集
结 束 语
通过对ADPI方法的研究可知,本文所使用的任意差分的精细积分方法比传统使用的有限差分法精度高。任意差分的精细积分方法可用于实际二维地质模型地震正演,是适合模拟复杂地质体的地震
图12 复杂地表与断层构造的地质模型
正演方法。ADPI方法可自然拓展到三维地质模型进行正演,且已经进行了初步的尝试。
参考文献
[1] 钟万勰.单点子域积分与差分.力学学报,1996,28(2):
159～163
[2] 强士中等.任意差分精细积分法及数值稳定性分析.应
用数学与力学,1999,20(3):256～262
[3] 马在田.地震成像技术—有限差分法偏移.北京:石油
工业出版社,1989
[4] 克莱鲍特著,许云译.地震成像理论及方法.北京:石油
工业出版社,1991
[5] LoewenthalDetal.Thewaveequationappliedtomi-gration.GeophysicalProspecting,1976,24:380～390
图13显示了正演所得到的单炮记录,从记录上可以看出由于地表复杂,即使地下构造简单,所得到的同相轴也大大偏离了正常形态(这种畸变既不是一般静校正问题,也不是一般动校正问题),而这种剖面又不是传统静校正方法所能对付的。这种同相轴偏离正常形态的现象称为静校不静,通过复杂地表模型正演对该问题得到验证。
通过以上二维正演的算例,证明ADPI方法对于解波动方程是可行的,并且精度较好
。
(本文编辑:张亚中)
向支持本刊的协办单位致谢!
2
OilGeophysicalProspecting
2003
LiuJun,DepartmentofEnergyResource,ChinaU-niversityofGeology,BeijingCity,100083,ChinaSelf-adaptivecoherentnoiseattenuationtechnique.
GaoShaowu,ZhaoBoandZhouXingyuan.OGP,2003,38(3):242～246
Noiseeliminationisanecessary,thereareallkindsofcoherentnoiseonseismicrecordsbecauseofmanyreasonscausedbyshooting,propagationandreceivingofseismicwave.Aimedatthecoher-entnoiseinseismicdataprocessing,thepaperpre-sentedself-adaptivecoherentnoiseattenuationtechniqueintime-spacedomain,separatingsignal/noise.Themethodlookstotheapparentvelocity,thatnoiseisonwholeprocessingplane,approxi-matelyasconstant(positiveornegativevalue),whichcaneffectivelyeliminatethelinearnoise,surfacewaveandmultipleetc.,andhasHI-FIsig-nal,thatisidealnoiseeliminationmethodforseis-micdataprocessingincomplexareaandareahav-inglowsignal-to-noiseratio,havingtheadvantagesofhighoperationspeedandgoodnoise-suppressedeffect.
Keywords:coherentnoise,seismicrecords,lineartransformation,separationofsignalandnoise,lin-earinversion,frequency-divisionprocessing
GaoShaowu,GRI,BGP,ZhuozhouCity,HebeiProvince,072750,China
Minimumdispersionalgorithmbyfinite-differencenumericmodelinganditsapplications.CaiQixin,HePeijun,QinGuangsheng,QinYaling,JiPing,MengFanbingandWangXiwen.OGP,2003,38(3):247～251,262
Thefinite-differenceapproachiscommon-usedalgorithmfornumericalsolutionofwaveequation,butordinaryfinite-differencealgorithmisdifficultinovercomingtheinterferenceofnumericdisper-sion.Onthebasisofsummarizingandanalyzingtheformermethods,thepaperintegratedandex-pandedtherelativetechniquesandgotoptimizedfi-nite-differencealgorithm,themainconceptsin-cludeshigh-orderfinite-difference,optimizeddif-ferenceparametersandFlax-CorrectedTransport(FCT)technique.Basedonthealgorithm,wewrotecorrespondingprogram,improvingprecisionofwaveequationforwardmodelingandreducingcomputationaleffortincomparisonwithordinarymethods,whichcanbewidelyusedforwavefieldmodelingandanalysisinrelieftopographyandcomplexsubsurfacegeologicstructuresandinandforwardmodelinginBMworkarea.
Keywords:waveequation,seismicmodeling,finite-differencealgorithm,dispersion,minimumdisper-sionalgorithm
CaiQixin,InstituteofGeologyandGeophysics,A-cademiaSinica,BeijingCity,100029,China
Seismicforwardmodelingbyarbitrarydifferencepreciseintegration.WangRunqiu,GengWeifeng
～257andWangShangxu.OGP,2003,38(3):252
Thepaperimprovedthefollowingmethods:single-pointsub-domainintegrationforparabolicequationpresentedbyZhongSixieandarbitrarydifferencefortemperaturefieldproblempresentedbyQiangShizhong,presentedaarbitrarydifferencepreciseintegration(ADPI)algorithmofwavee-quation.TheformulaofADPIistheoreticallyde-ducedanditsrealizationwayisgiveninthepaper;thespaceandtimeresolutionandthestabilityofdifferenceformatarediscussedtheoreticallyandprovedinrealcase.Thetheoreticalanalysisandre-alcasesshowedthehighprecisionofthepresentedmethodinthepaper.Thepresentedseismicfor-wardmodelingcasesshowedthatthemethodissuitableforcomplexnear-surfaceandstructuralgeologicbody.
Keywords:forwardmodeling,waveequation,finitedifference,arbitrarydifference,preciseintegra-tion,complexnear-surface,seismicmodel
WangRunqiu,SchoolofResourceandinformationtechnology,UniversityofPetroleum,BeijingCity,100083,China
Usingmultivariateseismicattributionstopredictreservoirinformation.HeBizhu,ZhouJieand
～262WangGonghuai.OGP,2003,38(3):258
Seismicattributionanalysistechniqueplaysamoreandmoreimportantroleinlateralreservoirprediction.Usingmultivariateseismicattributionanalystechniqueforpredictionofreservoirinfor-mationismainlybasedonspatialvariationofrockcharactersandfluidpropertiesinreservoir,result-edinvariationofseismicattributionssuchasseis-micreflectedwaveform,amplitudeandfrequencyetc..Thepapermainlydescribedthatusingclassi-fication(cluster)analysisandmultivariateseismicattributionregressiontechniquetoanalyzethever-ticalandlateralreservoirvariationcharactersinspecialareas,andputsstressonthepredictedreservoirdistributionscopeandphysicalpropertiesofreservoirinamodelofregionalgeologiccon-
2003年6月
石油地球物理勘探
第38卷 第3期
作者简介
杨举勇 高级工程师,1963年生;1984年毕业于江汉石油学
院地震解释专业。现在塔里木油田分公司从事地震勘探工程管理和技术方法研究工作。
韩立强 工程师,1972年生;1992年毕业于西北大学地质专
业。曾长期从事地震采集方法研究,现主要从事海底电缆地震采集方法研究。现任东方地球物理公司海上事业部海龙公司副总工程师。
王俊茹 副教授,1953年生;一直从事物探教学、科研及开
发工作。现在石家庄经济学院资源环境与工程系从事教学与科研工作。
刘 俊 工程师,1970年生;1992年毕业于吉林工业大学计
算机科学与工程系,2000年毕业于中国地质大学(北京)地矿专业,获硕士学位。曾参加并主持多项部、局级科研课题,获得过局级科技进步一等奖、部级科技进步一等奖。现在中国地质大学(北京)能源系攻读博士学位,研究方向为地震资料处理、综合解释。
高少武 高级工程师,1965年生;1986年毕业于西安地质学
院物探专业,1995年获成都理工学院硕士学位。曾从事工程物探、解释、方法研究、储层物性参数反演和储层预测、高分辨率地震数据处理以及软件设计与研究工作。在国内外有关刊物上发表论文十余篇。现在东方地球物理公司从事地震数据处理、方法研究以及计算机软件设计与研发等工作。
蔡其新 高级工程师,1965年生;1988年毕业于江汉石油学
院,现为中国科学院地质和地球物理研究所博士研究生,主要从事石油地球物理勘探研究工作。已发表论文近20篇,现任中原油田物探研究院副院长。
王润秋 副教授,1957年生;1982年毕业于南京大学数学系
计算数学专业,现为中国科学院地质与地球物理研究所在职博士研究生。现在石油大学资源与信息学院从事地球物理勘探基础理论、地震资料处理方法与软件研制方面的研究与教学工作。发表论文10余篇,获省部级科技进步奖4项。
何碧竹 高级工程师,1965年生;现为中国地质大学能源系
博士研究生。现在河南濮阳中原油田勘探开发科学研究院工作。
王西文 高级工程师,1956年生;1982年毕业于西安地质学
院,1987年在该校硕士研究生毕业后留校任教,1989年调入西北石油地质研究所,1997年进入中国科学院攻读博士学位;现在中国石油勘探开发研究院西北分院地球物理研究所工作。
宋延杰 副教授,1963年生;1988年毕业于大庆石油学院,
获硕士学位,现在大庆石油学院攻读博士学位。主要从事测井解释模型、剩余油评价技术、低孔渗储层解释方法、气层解释方法、声波测井理论及应用等研究工作。曾获中国石油天然气集团公司科技进步二等奖一项,大庆石油管理局科技进步二等奖、三等奖各一项,黑龙江省教育委员会科技进步二等奖、三等奖各一项;获中国石油天然气集团公司优秀教学成果一等奖、三等奖各一项。发表研究论文10余篇。
张军华 副教授,1965年生;1987年毕业于华东石油学院物
探专业,获学士学位,1995年毕业于石油大学应用地球物理专业,获硕士学位,2002年毕业于石油大学地球探测与信息技术专业,获博士学位。致力于地震资料处理和解释方法研究,曾在国内核心刊物上发表论文10余篇。
迟新刚 1976年生;现为成都理工大学信息工程学院应用
地球物理专业硕士研究生,从事地震属性参数提取和资料处理研究。
高楚桥 副教授,1966年生;1987年毕业于江汉石油学院,
1998年在石油勘探开发科学研究院获博士学位。现在江汉石油学院物探系从事测井教学与科研工作。
任广慧 高级工程师,1966年生;1988年毕业于石油大学应
用物理专业。现在中原石油勘探局测井公司资料解释中心从事测井资料评价等工作。
卢刚臣 工程师,1972年生;1995年毕业于石油大学综合勘
探专业,获学士学位。现在大港油田物探公司研究所从事地震资料综合解释工作,目前在石油大学(北京)攻读综合勘探专业硕士学位。
孙志华 高级地质师,1966年生;1988年毕业于江汉石油学
院石油地质专业,1991年获中国地质大学(北京)石油地质专业硕士学位,2001年获石油大学(北京)地球科学系矿产资源与勘探专业博士学位。现在中国地质大学从事科研工作。
王增明 高级工程师,1963年生;1983年毕业于华东石油学
院物探专业。一直从事地震勘探方法研究工作,发表过多篇论文。现任胜利油田物探公司西部公司经理。张丽莉 硕士研究生,1975年生;2000年毕业于江汉石油学
院。主要研究方向为成像测井资料的解释与图像处理以及综合物探研究工作。
陈清礼 副教授,1965年生;分别于1987年、1990年和2001
年在中国地质大学获学士、硕士和博士学位。多年来一直从事电法勘探的方法研究、试验和生产工作,曾主持和参与30多项科研课题,发表论文20余篇。现为电子科技大学博士后。
刘堂宴 副教授,1962年生;主要研究方向是岩石物理学,
测井储层描述与评价及核磁共振机理研究。
熊 翥 教授级高工;1965年毕业于北京石油学院地球物
理勘探专业。长期从事地球物理勘探方法研究与应用,出版过四本专著,并在不同期刊、学报、文集上发表论文40余篇,十分注重理论与实际的结合以及应用技术的研究。现任本刊主编。