【doc】 VTI介质准P波频率空间域组合边界条件研究

【doc】 VTI介质准P波频率空间域组合边界条件研究 VTI介质准P波频率空间域组合边界条件研 究 , { Z-. 』?{ ?? 谣睁糖; 第44卷第4期 2005年7月 石油物探 GE0PHYSICAIPR0SPECTINGFORPETR0LEUM VoI.44,No.4 Ju1.,2005 文章编号:1000—1441(2005)04—0301—07 VTI介质准P波频率空间域组合边界条件研究 吴国忱,梁锴 (1.同济大学教育部海洋地质重点实验室,上海200092;2.石油大学(华 东)地球资源与信息学 院,山东东营257061) 摘要:讨论了vTI介质中准P波频率空间域的边界条件.首先给出了vTI介质准P波波动方程,阐述了特征分 析方法的基本原理,讨论了边界上的反射系数与入射角和度量准纵,横波各向异性强度因子的关系;然后利用特 征分析法结合Kelvin-Christoffel方程,构造了准P波波动方程在不同边界和角点的频率域吸收边界条件,利用 最佳匹配层法构造了衰减边界条件;最后利用数值模拟对构造的边界条件进行了验证.为了得到好的吸收效 果,将吸收边界条件和衰减边界条件有机地结合起来,即先用最佳匹配层法衰减传播到边界的入射波能量,然后 再用吸收边界条件吸收边界反射,最终使边界反射降低到可以忽略不计.数值模拟的共炮点记录说明了组合边 界条件的良好效果. 关键词:边界条件;频率空间域;vTI介质;最佳匹配层;波动方程;组合 中图分类号:P631.4文献标识码:A 利用波动方程进行正演模拟,边界条件是关键 问题之一.典型的边界条件有2种,即单程波方程 构成的吸收边界和沿传播方向逐渐衰减的衰减边 界[1].构造边界条件的方法很多.Reynolds利用 波动方程分解法得到了透明边界条件,其特点是 O.角入射时反射系数为0[2].Clayton和Engquist 通过旁轴近似提出了吸收边界[3];利用特征方程分 析法[4],组合方向法[5]和优化系数法等也可以构造 吸收边界条件[6q].Cerjan提出了直接衰减法构 造衰减边界条件[9;Marfurt提出用海绵边界法构 造衰减边界[1..,Shin将其运用到频率域中…];Be— renger提出用最佳匹配层(PML)法来构建衰减边 界,PML法已在有限差分和有限元法正演模拟中 得到了广泛应用[1.Chew证明最佳匹配层可以 用于二维P波和S波耦合时的波场模拟[1.;Rap— papport将PML法引入到各向异性介质的边界条 件中[. 我们采用特征分析方法[2分离单程波,然后根 据不同边界区域有选择性地压制边界反射波,从而 得到VTI介质中准P波波动方程的时间域和频率 域吸收边界条件;采用最佳匹配层法[15,16作为衰 减边界条件,将这2种边界条件有机地组合在一 起,很好地消除了边界反射. 1VTI介质准P波波动方程 方程完全描述.由本构方程,运动微分方程和几何 方程可以建立一般各向异性介质中弹性波波动方 程.结合VTI介质的弹性矩阵,可得到VTI介质 中的弹性波波动方程.将平面波方程代入VTI介 质的弹性波波动方程中并忽略体力项,可推导出 Kelvin-Christoffel方程.该方程有非零解,其系数 矩阵的行列式为零,通过求解可获得VTI介质纵, 横波耦合的频散关系方程. 弹性矩阵中弹性参数C(,J—l,6)的物理意 义很不直观,Thomsen提出了一套表征TI介质弹 性性质的参数[1,即 一 ?警 一 ?警 ,:=鱼(1)百u y一 „ 一 (c13+c44)一(c33一c44) o一——=厂 式中,和VS0分别为准P波和准S波垂直TI介 质各向同性面的相速度;e和),是度量准P波和准 S波各向异性强度的参数;是影响垂直TI介质 收稿日期:2004—09—27;改回日期:2004—11—14. 假设地球介质为各向异性弹性介质,地震波在曩乏工19作65--.男剧 教授?博士在读,现从事地震波传 各向异性弹性介质中传播的波动特征可以由波动基金项目:国家 863”k[“~|(2002AA614010—3)资助. 羹?.{ ?302?石油物探第44卷 对称轴方向附近准P波速度大小的参数. 根据Alkhalifah声波假设思想q..,令VTI 介质纵,横波耦合的频散关系方程中横波速度为0 (一0),即可得到准P波的频散关系方程.将其 两边进行反傅里叶变换,并乘以准P波波场 p(x,Y,,),则得到VTI介质准P波波动方程,即 a?paz L[(1+2,)谝(耄+)+墙卜 2(~--编导(象+)(2) (2)式忽略了横波,比弱各向异性近似或小角 度近似更精确地描述了准P波在VTI介质中传播 的运动学与动力学特征,能够用来模拟各向异性介 质的零偏移距准P波剖面.(2)式是三维VTI介 质准P波波动方程,二维VTI介质准P波波动方 程为 -(1+2,)谝+墙一 2(~--)谝(3) 对(3)式做傅里叶变换,将p(z,X,,)变换到频 率域F(z,,co),得 (1+2,,?2_23ZF墙等+ 2(?一.~-,m4O?/”+F一0(4) 2特征分析方法的基本原理 二维VTI介质弹性波波动方程为 f:+!!?丛,生l OtPOx.Pa.P3x3z 1:+棼+詈势 式中,和分别为方向和方向的位移;C C13,c33,C4t为弹性参数;p为介质密度.将弹性参 数用Thomsen参数表征,同时令一0,再引人 l一J—l+—2e 一 干(6) 则弹性参数可表示为 fl1一lD(1+2e)墙一ID祈 c13一确厕一棚 c33一一ID旗 f44—0 将(6)式和(7)式代人(5)式,得 f=貉+1一 遁+诟 令和73分别表示沿X方向和2方向的速 度,和分别表示方向和方向的正应力,用 速度和应力作为中间变量将(8)式的方程降阶,得 到谏度一廊力场方程 式为 式中 一一 1 OtPOx 一 OtPa(9) 鲁一+棚3zVz 鲁一棚差+卺 令U一(,Vz,a.z:z,),则(9)式的矩阵形 A== B== oU —A. OU十,_oU(1O) 00 00 0 0 00 00 0 0 通过特征分析可以将A分解成左行波和右 行波.矩阵A的特征值从大到小分别为.=,, 2—0,3—0,4一一1.1一1,4一一1分别表 示沿z轴负方向和正方向传播的准P波的特征速 度.设所对应的特征向量为Z(—Z),A 是矩阵A的相似对角矩阵(其对角线上的元素为 矩阵A的特征值),矩阵P是由特征向量z组成 的矩阵,则 A—JAAPA 同理有 B一ABP口 式中 AA—A三_A(11)A口一A+A OOOOO1一POO 1一POOOOOOO 第4期吴国忱等.VTI介质准P波频率空间域组合边界条件研究?303? . 0堕 V3 0lD0V3.J 0000 0--V30..—1— lD o.一 00--V3 对于角边界,A和B均会产生边界反射. 在左上角边界,边界反射波相当于沿z轴正 方向传播的右行准P波和沿轴正方向传播的下 行准P波,因此 AAP U一 一0 AiPB3U 一0 可以吸收左上角的边界反射波,左上角吸收边界条 件为 a d U , A+a d U+B+au a(16) 同理可得,右上角,左下角,右下角的吸收边界 条件分别为 a d U , A—a d U z +B+3U a d U , A+a d U z +B一瓦3U a d U , a A— d U z +B一3U 3反射系数分析 (右上角) (左下角)(17) (右下角) 图1是入射波和反射波示意图.设入射角为 口,传播角为传播方向与轴正方向的夹角0,则入 射波和反射波的传播角分别为0.一一兀/2+口和 0z:~/2--a.传播角为0的平面简谐波相速度为 砌)一{(n升COS2叶 ? 304?石油物探第44卷 (zsin20+cosz0)z_8r/sin20cosz01 左边界 入射波 . x 反射波 图1边界反射葸图 已知 (一号+a)(号一a) 引入变量 “oh(号一a) 一 {1[;(120S20~+sin2a)+ J(Zcos口+sin口)--8r/sinZacos口 设左边界处波场为 l2 lRl= (18) F=exp()+ Rhexp(一j)(2o) 将波场表达式代人左边界吸收边界条件,整理可得 I1一?s口胡sin口I2胡一sinz口c.s口1l雾藩l (21) 右边界反射系数的模与左边界相同.同理可 求得顶,底边界的反射系数 (22) 式中,一(口) 一 {(nZa+cos2o~)+ (zsinZa+cosZa)--8rlsinzacos口 可见反射系数与入射角a,Thomsen参数,和 有关,而与准P波速度无关. 设计2组均匀介质,A组,不变,B组不变 (表1).反射系数曲线如图2所示,虚线为左,右 边界反射系数,实线为顶,底边界反射系数.对于 各向同性介质,左,右,顶,底边界反射系数相同(图 2a中的?曲线).对于VT1介质,如不变,则随 着,的增大反射系数整体趋势增大,且极小值对应 的入射角增加(图2a中的?和?曲线);如,不变, 则随着的增大反射系数整体趋势减小,且极小值 对应的入射角减小(图2b中的?和?曲线).对于 椭圆各向异性介质,其反射系数的趋势与各向同性 的相近,对于相同的入射角,左,右边界的反射系数 比各向同性的稍小,而顶,底边界的反射系数比各 向同性的稍大,且随着,和的增大差异增大(图 2b中的?曲线).在0%60.时,各介质反射系数均 较小,吸收效果明显. 表1模型参数 图2反射系数分析曲线 a,=0,0.2,—0;be=0.2.—0,0.2 第4期吴国忱等.VTI介质准P波频率空间域组合边界条件研究?305? 4吸收边界条件 f+1— 02Uz 3t一普3zOto33z3l.2 IloJ墨kx--:.u-~走k—x—k.}[.A.x:l==.24I耘礤z 0) 4~7) 10)3+?z走一走(25) +媚一 2v3xOzz3t一0(26). uu/ 一口++ 2妻v—3xOzz3t一0(右边界).,u,I./ 一 +一 !二丛曼:空一n Z3x3zOt 一 a4p一 + 呈 2v 二 3 一 3x23zO-t一0u +一一十蒜一 — 4~— --v] 一n 41V33xOzOt 一++ 二丛一n (顶边界) (底边界) o4p一 +a,4一口一十 二丛:一n 4vl口33xOzOt 04p 一一口. 0 z 4 a p一04 a p 3t t 4733一一一口?za一a一 二丛:一n 4v1hazazaf2 (左下角) (右下角) (27) 糈【27)瓦快到频率圾,就得到频率至1日J域?饭 动方程的边界条件 ?Fq-j.+警+ j?塑一0(左边界)3J?——xO—z2一u开 ?FJ?.+32F— j?霎一0(右边界)J?——u自夼 ?F+祈?aF十J. ? .O2F+ 一0(顶边界) ?F一筹a2F— j?塑一0(底边界) 3J?—-—x2—c3z—uJ氐开 Fq-j口l?.aF十J..OF一 2 壑笔晏一0(左上角)_一u上用 F—jl~FITJ.?33 d F+ ? 2— 4~--v432F一 0(右上角)?—一u自上用, F-k-j.aF—. 蓑+ 2— 4~--v]— 32 — F 一0(左下角)3_—x—3z—u『?用 ?F—j.aF一.篓+ ? 2— 4~--v]— 32 — F一 0(右下角)(28)3?——x—Oz—u自『,用Lz (左上角)5衰减边界条件 (右上角) 我们采用最佳匹配层法来构建衰减边界条件, 即在一定宽度的匹配层内给波动方程的导数项前 乘以一个衰减项(频率的函数),使传播到边界的入 ? . (j6?汕物{{l{弟巷 射波能龋很低.边界反射波能量弱要得到妤的 衰减敫粜.对匹配层的宽l望有一定的要求. 频率域的十算是按频率片埘空问网格进行整 体求解方程绀.1墨i此数据量很大.计于20021c} 的网格.采用压缩仔储技术方程【系数矩阵的数 据量约为?,J87MB..了得到较好的故果.最佳匹 配层宽嘘应为20,{『I1,网格点如采用1个网 格点宽瘦来计算200_l的有效网格.则需计算 包含匹配层的2t2?的网格.此时仍采压辅 存储技木.方程纽系数矩阡的数据醚约为844MB. 是原来的2.【9倍.. 我仃J采Hf个网格点宽度的匹配层.使匹配 层内波的能量有一定的莨减.然J再采用?盟收边界 条忭对边界反射进行吸收这样采用组合边界计 算上述横掣.{需汁算220,:2(1的网格.数据量约 为j11MI{.是来的J.3倍. 6数值算俐 6.1均匀vTI介质 模型参数为.一270,1)m/s.,一0.189.一 0.20I{:模型网格为】?】x【41.阐格间距为10m.采 样间隔为2his;震源为主频Hz的雷克子波.位 ]二模型中心录用频率域有限差分法进行数值模 拟.得到寸,葶剪率成分的渡场.图3为40Hz频 率成分的波场 b 3均匀v11介质I『I?z颁率战丹{噍场照 未界条件:}吸收边界杀怍;【边晔条件 rhl嗣:可见术加边条件叫?.边界反射很强. 入射波边界反射渡1fH:干涉;JJH上吸收边界条件 以后.边界压射波能量大大触弱,人射波的波动特 征清晰.但有抖动现象.说明边抖反射没有完全衰 减;采川组合边界条件[】11个网格点宽度的最佳fJ[ 层法与公式(18))后,人射陂的抖动现象基本消 失+边界反射雉弱. 6.2水平vri介质 我”J设计r2层水-vr】介质蟆型.:第【层 的参数为.一2768ms,E—LJ.【8,3-nL】{;第2 层的参数为f一32】m.e一一1.1.8一(】.2O4.= 模型有效网格为…I1(l】.网格和时间步同 上;震源为主频Hz的雷克子波,4和图5为 如Hz频率成分的披场快照手u共地点录 由4可见.波动特能清晰体现.只采用哦乏 收边界条f?l:叫还仔在少最边界眨目十;采用组合边界 条件后.边反射几被啵收lf净.人射波的扰动 是实际界面反射波与入射液f涉的结果 由圈5可弛.只采川吸收边条件时还存在少 量边界反射波.大饱榆距直达波与顶边界反射波锕 nb 网1水屠状Vr]介质如Hz频率成抒?照 ?峨收边条什;h毒【舟兼什 圈j水平层状vTI介质共炮点记录 a疆牝边扞条件Ih纽爵边条件 第4期吴国忱等.VTI介质准P波频率空间域组合边界条件研究?307? 互干涉,产生多个同相轴,直达波和实际界面反射 波在左,右边界产生了边界反射;采用组合边界条 件后,几乎没有边界反射,直达波同相轴很干净, 左,右边界的边界反射很小. 7结束语 本文所讨论的边界条件可以很好地消除边界 反射,使数值模拟的精度得到了保证. 1)利用特征分析法结合Kelvin—Christoffel方 程构造的时间域和频率域吸收边界条件,能吸收大 量的边界反射,且计算量增加不大; 2)利用最佳匹配层法构造的衰减边界,可以大 大衰减边界反射波的能量,但其要求有较宽的匹配 层宽度,这会使计算量大大增加; 3)对各个边界反射系数的分析表明,本文所讨 论的吸收边界条件对各向同性介质和VTI介质均 有较好的效果,但在高角度时反射系数较大,仍有 少量边界反射波没有被吸收; 4)将吸收边界条件和衰减边界条件组合在 一 起,先采用宽度较小的最佳匹配层使传播到边 界的入射波能量得到衰减,然后再采用吸收边界 条件对边界反射波进行吸收,这样不仅消除了边 界反射,而且由于匹配层宽度较窄,计算量增加 也不大. 参考文献 1张宝金.地震波参数反演及其可信度分析[D]:[学位论 文].上海:同济大学,2003 2ReynoldsACBoundaryconditionsforthenumerical solutionofwavepropagationproblems[J].Geophysics, 1978,43(6):1o99,1110 3ClaytonR,EngquistBAbsorbingboundaryconditions foracousticandelasticwaveequations[J].BullSeisSoc Am,1977,67:1529~1540 4董良国.地震波数值模拟与反演中的几个关键问题研 究[D]:[学位论文].上海:同济大学,2003 5KeysRGAbsorbingboundaryconditionsforacoustic media[J].Geophysics,1985,50(4):705~708 6崔兴福,张关泉.地震波方程人工边界的两种处理方法 EJ].石油物探,2003,42(4):452~455 7罗大清,宋炜,吴律,等.一种有效的处理模型角点反射 的方法[J].石油物探,2000,39(4):26~31 8余德平.数值模拟技术的发展现状[J].勘探地球物理进 展,2003,26(5—6):407~412 9CerjanC,KosloffD,KosloffR,eta1.Anonreflecting boundaryconditionfordiscreteacousticandelasticwave equations[J].Geophysics,1985,50(4):705~708 10MarfurtKJ.Accuracyoffinite-differenceandfinite- elementmodelingofthescalarandelasticwaveequa— tionsEJ].Geophysics,1984,49(5):533~549 11ShinCSpongeboundaryconditionforfrequency-domain modeling[J].Geophysics,1995,60(6):1870~1874 12BerengerJP.Aperfectlymatchedlayerforabsorption ofelectromagneticwaves[J].JComputPhys,1994, 114:185,200 13ChewWC,LiuQH.Perfectlymatchedlayersfor elastodynamics:Anewabsorbingboundarycondition [J].JCompAcoust,1996,4(4):341~359 14RappaportC~LPerfectlymatchedabsorbingboundary conditionsbasedonanisotropiclossymappingofspace [J].IEEEMicrowaveGuidedWaveLett,1995,5(3): 90,92 15ZengYQ,HeJQ,IiuQH.Theapplicationofthe perfectlymatchedlayerinnumericalmodelingofwave propagationinporoelasticmedia[J].Geophysics,2001, 66(4):1258,1266 16WangTI.TangX~LFinite-differencemodelingofe— lasticwavepropagation:Anonsplittingperfectly matchedlayerapproach[J].Geophysics,2003,68(5): 1749--1755 17ThomsenL.Weakelasticanisotropy[J].Geophysics. 1986,51(10):1954,1966 18AlkhalifahT.Acousticapproximationforprocessingin transverselyisotropicmedia[J].Geophysics,1998,63 (2):623,631 19AlkhalifahT.AnacousticwaveequationforanisotrO picmedia[J].Geophysics,2000,65(4):1239~1250 20AlkhalifahT.Anacousticwaveequationforortho— rhombicmedia[J].Geophysics,2003,68(4):1169, 1172 (本文编辑:戴春秋)