有限差分法地震波传播数值模拟

有限差分法有限差分法地震波传播地震波传播数值模拟数值模拟主要内容主要内容一、一、有限差分法地震波传播数值模拟方法有限差分法地震波传播数值模拟方法二、网格频散问题二、网格频散问题三、算法稳定性问题三、算法稳定性问题四、边界问题（自由边界、吸收边界）四、边界问题（自由边界、吸收边界）五、震源模拟问题五、震源模拟问题一、地震波传播数值模拟技术地震波模拟就是模拟地震波在地下的传播，模拟的目的就是在给定地质模型情况下，预测可能接收的地震记录。“Seismicmodelingisoneofthecornerstonesofgeophysicaldataprocessing.”-----------J.M.Carcione正演模拟地震波传播理论岩石物理数据采集理论指导物性参数研究传播规律指导设计观测系统地震解释验证数据处理参数反演提供理论数据试验处理流程提供正演方法地震波传播数值模拟应用领域T=2000ms炮点断层下覆界面反射能量强T=2300ms炮点炮点位于11km处的单炮记录断层下部界面反射较右边弱炮点位于17km处的单炮记录断层下部界面反射与右边相当二、地质模型的描述建模方法要和模拟方法相适应。1、介质用层状或块状描述，输入各块边界信息（折线顶点）以及介质参数。三维时按片输入。------FE方法、射线方法。2、建立模型的物性参数场。------FD方法。•Ray-tracingmethods(Asymptoticmethods)------高频近似，适应简单模型和特定波型，振幅计算困难，效率高三、地震波场数值模拟方法三、地震波场数值模拟方法综述综述•Integral-equationmethods(BasedHuygen’sprinciple)•Volumeintegralequations•Boundaryintegralequations------具有解析特征，只适用特殊模型，如除局部外介质均匀声波方程：声波方程：)(),(1)1(022rrtrftPKPrr−+∂∂=Δ⋅∇δρ应力运动方程：形变位移方程：)(21,,ijjiijuue+=本构方程：klijklijec=σiijijuf&&ρσ=+,弹性波方程：弹性波方程：•Directmethods(GridmethodsorFull-waveequationmethods)------没有限制，网格小时精度很高，效率是问题11、有限元法（、有限元法（8080年代初）年代初）不受边界几何形状的限制，具有灵活的剖分方法，效率低，不受边界几何形状的限制，具有灵活的剖分方法，效率低，模拟精度有待提高。模拟精度有待提高。22、傅立叶变换法（、傅立叶变换法（8080年代中期）年代中期）利用傅立叶变换求波场的空间微分，空间精度高，边利用傅立叶变换求波场的空间微分，空间精度高，边界处理困难，效率不高。界处理困难，效率不高。33、有限体法（、有限体法（9595年）年）介于有限元和有限差分之间，具有有限元的剖分灵活介于有限元和有限差分之间，具有有限元的剖分灵活性，又具有差分的效率。但处理边界困难。性，又具有差分的效率。但处理边界困难。44、有限差分法（、有限差分法（6868年）年）--------兼顾精度和效率，使用方便灵活，兼顾精度和效率，使用方便灵活，MarmousiMarmousi和盐丘模型均使用和盐丘模型均使用FDFD方法方法Directmethods：采取的措施zz高阶差分解法高阶差分解法----提高计算精度，减小数值频散提高计算精度，减小数值频散zz采用基于特征分析方法得到的吸收边界条件采用基于特征分析方法得到的吸收边界条件zz数值频散数值频散影响模拟精度，严重时造成计算溢出影响模拟精度，严重时造成计算溢出zz人为边界反射人为边界反射处理不当影响模拟效果及模拟效率处理不当影响模拟效果及模拟效率常规差分法存在的问题四、有限差分法地震波传播数值模拟方法数值频散问题数值频散问题------------高阶差分解法高阶差分解法)()]()(2)([121)(222MMmMmxoxmxfxfxmxfCxxfΔ+Δ−+−Δ+Δ=∂∂∑=()())(122122121)(NNnNnxOnxxfnxxfxxfCΔ+⎭⎬⎫⎩⎨⎧⎥⎦⎤⎢⎣⎡−Δ−−⎥⎦⎤⎢⎣⎡−Δ+Δ=∂∂∑=声波：弹性波：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡MMMMMMMM2222666644442222321321321321LMOMMMLLL⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡)()(3)(2)(1MMMMMCCCCM⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡0001M=高阶差分系数的确定：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡−−−−−−−−0001)12(531)12(531)12(53112531)()(3)(2)(11212121255553333MMLMOMMMLLLCCCCNNNNNNNNNNNNN)(2xoΔ)(10xoΔ数值频散试验dx=dz=10m,dt=1ms高阶差分为何会消除数值频散？高阶差分为何会消除数值频散？2M阶空间差分数值频散关系：其中，xxkΔλπ=Δ=φ2)()]()(2)([121)(222MMmMmxoxmxfxfxmxfCxxfΔ+Δ−+−Δ+Δ=∂∂∑=MMmMMmMMconstmCMcc2122)(22021)!22()1(21φφ⋅+=+−+=∑=+zz当当22M=2M=2时，时，zz当当22M=4M=4时，时，zz当当22M=6M=6时，时，zz当当22M=8M=8时，时，zz……......1212202φ−≈cc9014202φ−≈cc56016202φ−≈cc315018202φ−≈ccxxkΔλπ=Δ=φ2其中，数值频散理论分析小结数值频散理论分析小结zz，即空间离散造成的数值频散出现在正常波形之，即空间离散造成的数值频散出现在正常波形之后，作为一个后，作为一个““尾巴尾巴””出现。出现。zz对于相同频率子波，空间网格间距越大，数值频散越对于相同频率子波，空间网格间距越大，数值频散越严重；对于相同的空间网格间距，子波频率越高，数严重；对于相同的空间网格间距，子波频率越高，数值频散越严重。值频散越严重。zz只要只要，即只要一个波长包含几个空间步，即只要一个波长包含几个空间步长，随着差分精度长，随着差分精度22MM的提高，上述高阶差分解法产生的提高，上述高阶差分解法产生的数值频散会逐渐减小的数值频散会逐渐减小。。0cc12≤Δ=Δ=xxkλπφ不同差分精度空间频散曲线不同差分精度时间频散曲线五、边界问题五、边界问题吸收边界条件内部边界条件自由边界条件zz方程方程++边界条件数学上是非病态的边界条件数学上是非病态的zz方程方程++边界条件可以近似描述无限介质中的物理过程边界条件可以近似描述无限介质中的物理过程zz边界条件和内部点的计算方式是相容、不冲突的边界条件和内部点的计算方式是相容、不冲突的zz边界条件的实施是稳定的边界条件的实施是稳定的zz吸收效率要高吸收效率要高zz计算效率要高，不能因边界条件而引入过量的计算，失计算效率要高，不能因边界条件而引入过量的计算，失去采用边界条件的意义去采用边界条件的意义连续问题离散问题设计吸收边界条件的目标：zz吸收程度不高吸收程度不高zz边角（边角（33DD））及四个角（及四个角（22DD））处的人为边界反射难以处理处的人为边界反射难以处理zz存在计算稳定性问题存在计算稳定性问题一般吸收边界条件存在的问题：一般吸收边界条件存在的问题：吸收边界条件类型：吸收边界条件类型：zzABCsBasedonParaxialApproximationABCsBasedonParaxialApproximationzzABCsBasedonDampingMethodABCsBasedonDampingMethodzzPMLPML（（PerfectMatchedLayerPerfectMatchedLayer））absorbingboundaryconditionsabsorbingboundaryconditionsABCsBasedonParaxialApproximation012=−+ttzzxxPvPPAcousticWave:(BottomBoundary)222xzkvk−±=ω2221ωωxzkvvk−=⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧+⎟⎟⎠⎞⎜⎜⎝⎛−−=LL222222281211ωωωxxzkvkvvkvkzω≈01=∂∂+∂∂tPvzPFirst-orderApproximation:Third-orderApproximation:(45Degree)⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧⎟⎟⎠⎞⎜⎜⎝⎛−−=222222281211ωωωxxzkvkvvk01434332323223=∂∂−∂∂∂+∂∂∂−∂∂∂tPvtxPvzxPvtzPSecond-orderApproximation:(15Degree)⎭⎬⎫⎩⎨⎧−≈222211ωωxzkvvk02122222=∂∂−∂∂+∂∂∂xPvtPvtzPElasticWave:(BottomBoundary)xzzzxxttWUUU)(2222βαβα−++=xzzzxxttU)(2222βααβ−++=ABCsBasedonDampingMethod)exp(*),(),(2djiPjiPoldnewσ−=d•ABCsBasedonCharacteristicAnalysiszUMMxUMMtULLLLLLLLBBBBBBAAAAAA∂∂∂∂∂∂)()()(1)(1)(1)(1+−−−+−−−+++=如2D弹性波方程可以分解为：xUMLLAAA∂∂)(1+−xUMLLAAA∂∂)(1+−zUMLLBBB∂∂)(1+−左边界右边界上边界下边界左上角右上角左下角右下角zUxUtUBA∂∂+∂∂=∂∂1zUxUtUBA∂∂+∂∂=∂∂2zUxUtUBA∂∂+∂∂=∂∂1zUxUtUBA∂∂+∂∂=∂∂2zUxUtUBA∂∂+∂∂=∂∂11zUxUtUBA∂∂+∂∂=∂∂12zUxUtUBA∂∂+∂∂=∂∂21zUxUtUBA∂∂+∂∂=∂∂22⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛=00000211ccccA其中：⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛−−=00000212ccccA⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛=ccccB00000211⎟⎟⎟⎠⎞⎜⎜⎜⎝⎛−−=ccccB00000212我们采用的吸收边界条件如下（2D）：边界吸收前后波前快照X分量Z分量边界吸收前边界吸收后吸收边界条件试验（各向同性介质）X分量Z分量边界吸收前边界吸收后吸收边界条件试验（TI介质）3D吸收边界条件三维声波方程转化为一阶方程组(简写为ABC）。222222222tuc1zuyuxu∂∂=∂∂+∂∂+∂∂zUCyUBxUAtU∂∂+∂∂+∂∂=∂∂左面A1BC右面A2BC前面AB1C后面AB2C底面ABC1顶面ABC2上左棱A1BC1上右棱A2BC1上前棱AB1C1上后棱AB2C1下左棱A1BC2下右左棱A2BC2下前棱AB1C2下后棱AB2C2左前棱A1B1C左后棱A1B2C右前棱A2B1C右后棱A2B2C上左前角A1B1C1上左后角A1B2C1上右前角A2B1C1上右后角A2B2C1下左前角A1B1C2下左后角A1B2C2下右前角A2B1C2下右后角A2B2C23D吸收边界条件对三维地震波数值模拟显得尤为重要！前人所用的一种弹性波吸收边界条件底边界吸收边界条件：()()⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧+=βτρ25.0oldxzoldxxvv()()⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧+=βτρ15.0oldzzoldzzvv()()⎪⎭⎪⎬⎫⎪⎩⎪⎨⎧+=vCCColdZoldzzxxβττ11333135.0()(){}voldzoldzzzzβττρ15.0+=()(){}voldxoldxzxzβττρ25.0+=其中，1986年，Baylissetal.在使用MacCormack分裂算法模拟二维各向同性介质中的