1、第四章第四章 弹性波场数值模拟弹性波场数值模拟内容提要内容提要第一节 波动方程的交错网格有限差分模拟第二节 二阶波动方程的有限差分模拟引言引言 地震数学模型分类：1)根据地质构造形态、特征的不同假设分：一维地质模型：地层是水平的，沿横向地层参数是不变的，也即水平层状介质模型。在地面上任何一点得到的地震记录是一样的，是一道地震记录。二维地质模型：地质体是一个二度体，沿某个方向无限延伸，在垂直于延伸方向得到的地震剖面是一样的。三维地质模型：真正的地质体一般都是三维的。这时沿地面上不同方向的测线得到的地震响应一般是不同的。引言引言地震数学模型分类：2）根据对地层物性参数的不同假设，可以是只考虑各层在
2、速度上的差别，设每层速度是常数或连续变化；考虑各层速度和密度都不一样；给出各层的弹性参数，甚至吸收系数等参数,即构造模型和岩性模型。3)按照所依据的地震波理论，可以是运用运动学理论（射线理论）计算波的旅行时；考虑波的动力学特点，把界面的反射系数、透射系数、各种多次波、吸收损失、波前扩散等的影响都考虑进去，在计算出的地震响应中反映波的传播时间、振幅、相位和方向性等特点。在实际的计算方法上分别用射线法、绕射叠加（即物理地震学）法和解波动方程法。引言引言 主要优点：改变模型参数很方便，也可以灵活地按具体要求选用不同的理论和公式，特别适用于人机联作解释过程中反复修改模型和计算模型的地震响应。在理论研究
3、中也很有用处。缺 点：由于各种地质现象往往是非常复杂的，理论计算难以精确地反映真实情况，不可避免地要忽略了许多因素。特别是要研究各种构造形态复杂的地层圈闭、岩性的变化都有不少困难。在解决这些问题时，地震物理模型技术有其独特的优点，有时是很好的配合和补充。地震数学模型技术的特点三维声波方程111yxzxyzvPvvKtxyzvPtxvPtyvPtz vx,vy,vz为质点的振动速度，P为压力，K为体积模量，为密度，v为地震波速度。2Kv1.1 二阶精度交错网格法各个物理量的空间配置关系各个物理量的空间配置关系按一定的规律按一定的规律三维情况下，各物理量在节点上和网格边上的配置情况如上右图所示。三
4、维情况下，各物理量在节点上和网格边上的配置情况如上右图所示。编编号号1 1为压力和体积模量的位置为压力和体积模量的位置，编号编号2 2为密度和为密度和x方向速度分量的位置方向速度分量的位置，编编号号3 3为密度和为密度和y方向速度分量的位置方向速度分量的位置，编号编号4 4为密度和为密度和z方向速度分量的位置方向速度分量的位置。三维规则网格三维规则网格 二维交错网格二维交错网格 三维交错网格三维交错网格,p i k1,p ik,1p i k 1,1p ik12,zvi k 12,xvik12,1xvik12,zvi k 1.1 二阶精度交错网格法交错网格解法的优点：l精度高l能求解非均匀介质中
5、的波场，而不需要特别处理两种不同介质之间的边界条件(应力、位移连续)111yxzxyzvPvvKtxyzvPtxvPtyvPtz 1.1 二阶精度交错网格法交错网格差分格式：用中心差商近似声波方程中的一阶偏导数112,1122,1122,1122,2,yxzi j kxxxi j k nyyyi j k nzzzi j knnvPvvKtxyzp i j kp i j kpttvij k nvij k nvxxvvi jk nvi jk nyyvi j knvi j knnvznz 1.1 二阶精度交错网格法121212122121,112211211222,111,11,1yxzxxxj k
6、yyykikzzznni jnijvvPPvPtxtytzvj kvj knnnnvttvvikvikttvi jvi jvttjijiknkn ,1.1 二阶精度交错网格法1122112211221122,1122,1122,1,1,1,ij k ni jk ni j knp ij k np i j k npxxp i jk np i j k npyyp i j knp i j k npzz1.1 二阶精度交错网格法将差商代入声波方程1122111122221122,xxyyzzvij k nvij k nxvi jk nvi jk np i j k np i k j nKtyvi j kn
7、vi j knz 1.1 二阶精度交错网格法11112222,1,1,1xxvij k nvij k np ij k np i j k ntx 11112222,1,1,1yyvi jk nvi jk np i jk np i j k nty 11112222,1,1,1zzvi j knvi j knp i j knp i j k ntz 1.1 二阶精度交错网格法整理上述方程得各量的更新式1122112211221122,xxyyzzvij k nvij k nxvi jk nvi jk np i j k np i k j nK tyvi j knvi j knz1.1 二阶精度交错网格法
8、11112222,1,1,xxtvij k nvij k np ij k np i j k nx11112222,1,1,yytvi jk nvi jk np i jk np i j k ny11112222,1,1,zztvi j knvi j knp i j knp i j k nz1.1 二阶精度交错网格法加上源项后11112222,1,1,xxtvij k nvij k np ij k np i j k nx11112222,1,1,yytvi jk nvi jk np i jk np i j k ny11112222,1,1,zztvi j knvi j knp i j knp i
9、j k nz1122112211221122,xxyyzzvij k nvij k nxvi jk nvi jk np i j k np i k j nK ts i k j nyvi j knvi j knz1.1 二阶精度交错网格法或1122112211221122,xxyyzzvij k nvij k nxvi jk nvi jk np i j k np i k j nK tyvi j knvi j knz111111222222,1,1,xxxtv ij k nv ij k np ij k np i j k nf ij k nx111111222222,1,1,yyytv i jk nv
10、 i jk np i jk np i j k nf i jk ny111111222222,1,1,zzztv i j knv i j knp i j knp i j k nf i j knz1.1 二阶精度交错网格法速度模型速度模型1.1 二阶精度交错网格法模拟实例模拟实例一个共炮集模拟记录（已消除了直达波）32546879 00.20.40.60.81.01.22.41.61.82.0Time(s)Distance(km)1.1 二阶精度交错网格法1.2 弹性波动方程的交错网格有限差分模拟 交错网格法适用于求解一阶偏微分方程组。将弹性介质中的波动方程组变交错网格法适用于求解一阶偏微分方程组
11、。将弹性介质中的波动方程组变形为一阶偏微分方程组。形为一阶偏微分方程组。弹性介质的弹性介质的运动平衡微分方程运动平衡微分方程为为 222222xyxxxzxyxyyyzyzyzxzzzuftxyzvftxyzwftxyz 式中，式中，,u v w为质点振动的位移分量，为质点振动的位移分量，,iiijif 分别为正应力、剪应力和外力分别为正应力、剪应力和外力分量。分量。1.2 弹性波动方程的交错网格有限差分模拟 以质点振动的速度分量以质点振动的速度分量,xyzvut vvt vwt 来替换运动平来替换运动平衡微分方程中的位移分量，得衡微分方程中的位移分量，得 xyxxxxzxyyxyyyzyzy
12、zxzzzzvftxyzvftxyzvftxyz 1.2 弹性波动方程的交错网格有限差分模拟 将将应变应变位移关系位移关系(*)(*)代入应力代入应力-应变满足的虎克定律应变满足的虎克定律（*）(*)(*)(*)(*)uvwxyz 1,21,21,2xxxyyyxzzzyzuuveexyxvuweeyzxwwveezyz2,22,22,2xxxxxzxzyyyyxyxyzzzzyzyzeeeeee1.2 弹性波动方程的交错网格有限差分模拟 有有 222xxxzyyxyzzyzuvwuuwxyzxzxuvwvuvxyzyyxuvwwwvxyzzyz 1.2 弹性波动方程的交错网格有限差分模拟 求
13、其对求其对 t t 的的一次导数一次导数，并将并将,xyzvut vvt vwt 代入其中，代入其中，得得 2,2,2,yxxxxzxzxzyyyyxyyxxzyyzyxzzzzzvvvvvvtxyzxtzxvvvvvvtxyzytyxvvvvvvtxyzztyz 1.2 弹性波动方程的交错网格有限差分模拟 未知量为未知量为,xyzxxyyzzxzxyyzv v v的一阶弹性波为的一阶弹性波为 2,2,2,yxxxxzxzxzyyyyxyyxxzyyzyxzzzzzvvvvvvtxyzxtzxvvvvvvtxyzytyxvvvvvvtxyzztyz,xyxxxxzxyyxyyyzyzyzxzz
14、zzvftxyzvftxyzvftxyz ,xxyyzz,xyxzyz,xyzv v v在节点上，在节点上，在面的中心点上，在面的中心点上，在楞上。在楞上。1.2 弹性波动方程的交错网格有限差分模拟 112211112222112211122111122221122111112222,2,1,nnnnxxxxxxnnnnyyzznnxxnnxzxznzi j ki j kvij kvij ktxvi jkvi jkvi j kvi j kyzvij kvij kij kij kztv11221122112211221211111111222222221,1,nznnnnxxxxxxxnnnnx
15、zxzxzxzij kvi j kxvij kvij kij ki j kfij ktxijkijkij kij kyz 1.3 高阶交错网格法 地震波方程的离散化必将涉及到地震波场的数值逼近问题。地震波场的数值模拟精度一方面依赖于剖分网格的形状和大小，另一方面取决于离散波场的时间微分和空间微分的逼近误差。这里主要讨论规则网格和交错网格上的差分算子的高阶近似，截断误差，差分系数的收敛速度以及与虚谱差分算子精度的对比。u x推导函数 一阶导数的6阶精度差分系数。设 有7阶导数，则 在 和 处的7阶泰勒展开式为 u x u x1xix1xix1.3 高阶交错网格法高阶规则网格法2i 1i i2i
16、1i 3i 3i 35735791135722221!3!5!7!iiiiiixxxuuuuxuuOxxxxx357357922357222222221!3!5!7!iiiiiixxxuuuuxuuOxxxxx357357933357333322221!3!5!7!iiiiiixxxuuuuxuuOxxxxx同理得：1.3 高阶交错网格法高阶规则网格法由于一阶导数6阶精度中心差分近似式可表示为33333333312355555555777777232221!1!1!232223!3!3!232225!5!5!2227!7!iiiiiiiiiiiuuuxxxxxxxxxuuuxxxaaaxxxu
17、uuxxxxxuuxx577700327!iLiuxxERxux则有LERLe777iuxx为误差项，其系数为（即不含导数项和x）的系数，项的系数。例如该例中的1.3 高阶交错网格法高阶规则网格法化简可得系数方程 1112225531231 212301230aaa 777321123112322227!7!7!21!LLLmmeaaamaL求解系数方程得 1.3 高阶交错网格法高阶规则网格法 u x u xximx任意2L阶精度中心有限差分系数计算公式推导如下。设有2L+1阶导数，则在处的2L+1阶泰勒又有由于一阶导数2L阶精度中心差分近似式可表示为展开式为1.3 高阶交错网格法高阶规则网格
18、法将上述L个方程代入、化简，有式中，差分系数由以下方程确定解得1.3 高阶交错网格法高阶规则网格法L 中心差分近似的截断误差系数为中心差分近似的极限，即时，有于是有其中，一阶导数的中心差分算子长度为2L。1.3 高阶交错网格法高阶规则网格法2i u x21123,iiiiiiuuu uuu推导函数交错网格一阶导数的6阶精度差分系数。一阶导数的6阶精度差分系数计算由离散点确定。交错网格上1i 1.3 高阶交错网格法高阶交错网格法i2i 1i 3i 12i12ji 2yx 要确定网格点 上的一阶导数，取，u x u x12ji 531135,jjjjjjuuuuuu于是有在的一阶导数由离散点确定。
19、根据泰勒展开式，有1.3 高阶交错网格法高阶交错网格法2i 1i i2i 1i 3i jj+1j+2j-2j-1j+3j-3j-5j-4j+4j+5系数方程求解该线性方程可得于是若算子的对称点为i，则上式可改写为2i 1i i2i 1i 3i 1.3 高阶交错网格法高阶交错网格法 u x u x0212xxmx同理，可以推导出交错网格任意2L阶精度有限差分系数计算公式。设有2L+1阶导数，则，在处的2L+1阶泰勒展开式为由于交错网格一阶导数2L阶精度差分近似式可表示为1.3 高阶交错网格法高阶交错网格法其中，差分系数由以下方程确定其中，差分系数由以下方程确定将上述将上述L个方程代入、化简，有个
20、方程代入、化简，有1.3 高阶交错网格法高阶交错网格法当当截断误差系数为截断误差系数为1.3 高阶交错网格法高阶交错网格法当当于是于是1.3 高阶交错网格法高阶交错网格法交错网格法交错网格法规则网格法规则网格法211221!LLLmmemaL截断误差系数比较截断误差系数比较2i 1i i2i 1i 3i 2i 1i i2i 1i 3i 3i 1.3 高阶交错网格法高阶交错网格法截断误差系数比较1.3 高阶交错网格法高阶交错网格法阶数a1a2a3a4截断误差系数25.000000010-11.666666710-146.666666710-1-8.333333310-2-3.333333310-
21、267.500000010-1-1.500000010-11.666666710-27.142857110-388.000000010-1-2.000000010-13.809523810-2-3.571428610-3-1.587301510-3阶数a1a2a3a4截断误差系数21.00000004.166666710-241.1250000-4.166666710-2-4.687500010-361.1718750-6.510416710-24.687500010-36.975446410-481.1962891-7.975260410-29.570312510-3-6.975446410
22、-4-1.186794710-4规则网格一阶导数的偶数阶精度有限差分格式系数表交错网格一阶导数的偶数阶精度有限差分格式系数表00.050.10.150.20.250.30.350.40.450.500.10.20.30.40.50.60.70.80.91设设有有一一函函数数 i tu te，其其理理论论导导数数为为 i tu ti e。其其规规则则网网格格中中心心差差分分商商为为 11Lmmutau tm tu tm tt 12sinLmmiu tu tamtt，微微分分算算子子的的振振幅幅谱谱为为 111122sinsin2222sin2sinLLmmmmLLmmmmAamtamftttka
23、mk ftamttN N N 为为离离散散样样点点数数，k k=1 1，2 2，N N/2 2。以以 k k/N N 为为横横坐坐标标，以以振振幅幅谱谱为为纵纵坐坐标标画画图图(时时间间采采样样间间隔隔设设为为 1 1)，并并对对振振幅幅谱谱归归一一化化。2L=6KA1.3 高阶交错网格法高阶交错网格法设有一函数设有一函数 i tu te，其理论导数为，其理论导数为 i tu ti e。其其交错交错网格中心差分商为网格中心差分商为 11212122Lmmmmu tau ttu ttt 1221sin2Lmmimu tu tatt，微分算子的振幅谱为微分算子的振幅谱为 11111112212si
24、nsin21222sin21sin21sincos,sin,22LLmmmmLLmmmmLLLmmmmmmmAatamftttkamk ftamttNNamamamk 2L=6KA1.3 高阶交错网格法高阶交错网格法不同长度的中心差分算子（a）和交错网格差分算子（b）的振幅谱与虚谱差分算子的振幅谱的对比。图中的线条从下到上对应L=1-10。误差比较#1.3 高阶交错网格法高阶交错网格法对对 111yxzxyzvPvvKtxyzvPtxvPtyvPtz 变形变形 第一式对时间求导，第第一式对时间求导，第 2 2-4 4 式分别对式分别对 x，y，z 求导求导。2.1 三维声波方程规则网格有限差分
25、模拟 2222222222111111yxzxyzvPvvKtx ty tz tvPx tPPPPKxxvPy tyyvPz tzztxxyyzz 2.1 三维声波方程规则网格有限差分模拟 假设介质均匀，假设介质均匀，22222222PKPPPtxyz ，2Kv,v为声波速度。为声波速度。2.1 三维声波方程规则网格有限差分模拟 二阶精度的二阶中心差商代替二阶偏导数，得二阶精度的二阶中心差商代替二阶偏导数，得 2222,12,11,2,1,1,2,1,1,2,1,P i j k nP i j k nP i j k ntP ij k nP i j k nP ij k nxK i j kP i j
26、k nP i j k nP i jk ni j kyP i j knP i j k nP i j knz 2.1 三维声波方程规则网格有限差分模拟 更新格式更新格式 2222,12,11,2,1,1,2,1,1,2,1,P i j k nP i j k nP i j k nP ij k nP i j k nP ij k nxP i jk nP i j k nP i jk nK tyP i j knP i j k nP i j knz 2.1 三维声波方程规则网格有限差分模拟 更新格式更新格式的稳定性条件：的稳定性条件：对离散格式作时间和空间的四重对离散格式作时间和空间的四重 FourierFo
27、urier 变换，变换，000,1,1,xxyzitxyzikxxyzP i j k nP k k kP i j k nP k k keP ij k nP k k ke 00022000000200022222xxyxzzikxikxikxikxititikxikxPePPexPePPeK tPePPeyPePPez 2.1 三维声波方程规则网格有限差分模拟 22222222yxxxzzititikxikxikxikxikxikxeeK teeeeeexyz 22222222sinsinsin222sin2yxzkykxkytK txyz 根据根据 Von NeummanVon Neumman
28、 稳定性分析法知，上式右端的取值必须满足左稳定性分析法知，上式右端的取值必须满足左侧的取值范围，即侧的取值范围，即 2222222sinsinsin22201yxzkykxkyK txyz 2.1 三维声波方程规则网格有限差分模拟 2222222sinsinsin22201yxzkykxkyK txyz 取取,xyzkkk的最大值，即的最大值，即,xyzkkkxyz，则有，则有 222211101K txyz 或或 2221111v txyz 2.1 三维声波方程规则网格有限差分模拟 另一种求稳定性条件的方法 222222,12,1,1,2,1,xyzxP i j k nP i j k nP
29、i j k nK tL P i j k nL P i j k nL P i j k nP ij k nP i j k nP ij k nL P i j k nx 对上式作空间变量的 Fourier 变换，时间导数用二阶中心差商代替 11222222nnnnxyzPPPvkkkPt 222224sin2xxxkxikxikxxeekxx 其中 2.1 三维声波方程规则网格有限差分模拟 12222212nnnxyzPvtkkkPP 22222112110nnxyznnvtkkkPPPP 上式稳定的条件是增长矩阵的谱半径上式稳定的条件是增长矩阵的谱半径 1A。A 的特征值方程的特征值方程 22222
30、2210 xyzvtkkk 222222240 xyzvtkkk （此时（此时1）22222222xyzvtkkk 2222204xyzvtkkk 2.1 三维声波方程规则网格有限差分模拟 2222204xyzvtkkk 将将2222222222224sin4sin4sin,yxzkykxkzxyzkkkxyz代入上式代入上式，得，得 22222222222sinsinsin01yxzkykxkzvtxyz 同样，取同样，取,xyzk kk的最大值，即的最大值，即,xyzkkkxyz，则有，则有 2222211101vtxyz，2221111v txyz 2.1 三维声波方程规则网格有限差分模
31、拟 2222222222222xuuvwvwtxxyzxyzuvuwfyyxzzxvuvuvwtxyxyxyzuyx 22222222yzwwvfzzyzwuwwvtxzxyyzuvwuvfzxyzzxy 2.2 三维弹性波方程规则网格有限差分模拟 均匀介质均匀介质 222222222222222222222222222222222222222222222xyzuuuuvwftxyzx yx zvvvvuwftyxzy xy zwwwwuvftxyzz xz y 若直接离散化，求出的是均匀介质中的波场。在遇到介质的分界面时，需若直接离散化，求出的是均匀介质中的波场。在遇到介质的分界面时，需要根
32、据位移、应力连续性条件对界面处的波场作衔接处理，即根据连续性条件，要根据位移、应力连续性条件对界面处的波场作衔接处理，即根据连续性条件，给出边界处位移分量满足的更新格式。给出边界处位移分量满足的更新格式。2.2 三维弹性波方程规则网格有限差分模拟 直接离散化非均匀介质的波动方程直接离散化非均匀介质的波动方程 2222222222222xuuvwvwtxxyzxyzuvuwfyyxzzxvuvuvwtxyxyxyzuyx 22222222yzwwvfzzyzwuwwvtxzxyyzuvwuvfzxyzzxy 2.2 三维弹性波方程规则网格有限差分模拟 为为说明问题的方便说明问题的方便，将，将三维
33、方程三维方程退化为二维方程，即介质在退化为二维方程，即介质在 y 方向是均匀方向是均匀不变的，于是不变的，于是 2222222222222422xuuuwwtxxxx zxxzuwufzzzxz 2.2 三维弹性波方程规则网格有限差分模拟 以中心差商代替一阶微商，以中心差商代替一阶微商，,1,1,2i kikikxx ,1,12u i k nu i k nu i k ntt ,1,1,2u i k nu ik nu ik nxx 以二阶中心差商代替二阶微商，以二阶中心差商代替二阶微商，222,1,2,1,u i k nu ik nu i k nu ik nxx 222,12,1u i k nu
34、 i k nu i k nu i k ntt 2.2 三维弹性波方程规则网格有限差分模拟 2222222222222422xuuuwwtxxxx zxxzuwufzzzxz 112,1,1,1,1,1,1,22222,1,11,11,11,1,1,1,1,1,1222422nnni kikiki ki ki ki knnnnnikikiki kiki ki knnnnikikikiki kikiki kikikni kuuuuuuuutxxwwwwx zwwx z ,1,1,1,1,1,1,11,1,22,1,1,2222ni ki ki ki ki ki ki knnnni ki kikik
35、i knnnni ki ki kx i kuuwwzx zuuufz 2.2 三维弹性波方程规则网格有限差分模拟 2,1,1,11,1,1,222,1,1,2222,1,11,11,11,12,1,1,1,1,222422i kikiknnnnni ki ki kikiki knnniki kiki ki knnnnikikikiki kikiki kikiktuuuuuxtuuuxtwwwwx zt,1,122,1,1,1,1,1,11,1,222,1,1,2222nni ki ki ki ki ki ki ki knnnni ki kikiki knnnni ki ki kx i kwwx
36、 zttuuwwzx ztuuufz 2.2 三维弹性波方程规则网格有限差分模拟 前面讲述的有限差分方法中，都以泰勒展开的方式推导出微商的差商形式，以差商代替微商，从而实现了偏微分方程的离散化，得到计算波场的差分格式。这里讨论一种利用 Fourier 变换求空间导数的方法，通常称为伪谱法或虚谱法。声波方程 22222222PKPPPtxyz 中的空间偏导数用 Fourier 方法求取。,xik xxP ky z tP x y z t edx ,xik xxP x y z tP ky z t edx ,xxik xik xxxxP x y z tP ky z t edxik P ky z ted
37、xxx 2222,xxik xxik xxxP x y z tP ky z t edxxxikP ky z tedx ,P x y z tx ,xxik P ky z t 22,P x y z tx 2,xxikP ky z t ,P x y z t对 x 的一阶导数可以通过以下方法求取：(1)将,P x y z t对 x 作 Fourier 变换，得到,xP ky z t；(2)以xik乘,xP ky z t，得到,xxik P ky z t；(3)对,xxik P ky z t作反 Fourier 变换，得到,P x y z tx。高阶导数可通过类似的方法求导，例如，二阶偏导数22,P x
38、 y z tx可由2,xxikP ky z t的反 Fourier 变换得到。用 Fourier 变换求空间导数的优点是精度高。通常的规则网格二阶差分格式求解波动方程时，需要每个波长采8-10个点才能避免网格色散（波的传播速度与传播方向和频率有关），而伪谱法则每个波场采个点即可。以一阶声波方程为例来说明产生网格色散的原因和克服办法。111yxzxyzvPvvKtxyzvPtxvPtyvPtz 差分格式：1122112211221122,xxyyzzp i j k np i k j nvij k nvij k nxvi jk nvi jk nK tyvi j knvi j knz 差分格式：11
39、221122,1,1,xxvij k nvij k ntp ij k np i j k nx 11221122,1,1,yyvi jk nvi jk ntp i jk np i j k ny 11221122,1,1,zzvi j knvi j kntp i j knp i j k nz 设差分格式的解为设差分格式的解为 0000,exp,exp,exp,expxyzxxxyzyyxyzzzxyzp i j k nPiin tk i xk j yk k zvi j k nviin tk i xk j yk k zvi j k nviin tk i xk j yk k zvi j k nviin
40、 tk i xk j yk k z 其中其中ii为虚数单位，为虚数单位，0P，0 xv，0yv，0zv为振幅，为振幅，,xyzk kk为数值为数值波数。波数。将它们代入差分格式中，有 110022011220112201122expexpexpexpexpexpexpexpxxxyyyzzzPiitPiitvii kxii kxxvK tii kyii kyyvii kzii kzz 将将它们它们代入代入差分格式中差分格式中，有，有 11002211022expexpexpexpxxxxviitviittPii kxii kxx 11002211022expexpexpexpyyyyviitv
41、iittPii kyii kyy 11002211022expexpexpexpzzzzviitviittPii kzii kzy 进一步整理进一步整理四个等式四个等式后后得得 000111102222sinsinsinsinyxzxyzvvvPtK tkxkykzxyz 120012sinsinxxkxtPvxt 120012sinsinyykytPvxt 120012sinsinzzkztPvxt 后三式代入第一式，后三式代入第一式，得得 222112222222112222sinsinsinsinxyzK ttkxxK tK tkykzyz 以以2Kv代替代替K，则，则 222211si
42、nsin2211sinsin22xyztkxv txkykzyz 这是数值计算时的色散关系式，而理论色散关系式为 2222xyzvkkk (*)(*)当,0txyz 时，222211112222sinsinsinsinxyzkxkykztv txyz(*)(*)的极限就是（*）式。但在实际计算中,0txyz 不成立，因此，数值色散关系式（*）与理论色散关系式（*）之间存在误差，这就引起数值计算时波动的色散，即波的传播速度不为常数，而是与频率、传播方向有关的变量。为考虑问题的简便，假设介质为二维介质，剖分区域时用正方形网格，边长为，Sv t ，N为一个波长所采样点数。这时，数值频散关系式变为 2
43、2221cossinsinsinsin22SkkSN 为波的传播方向与x轴的夹角。求解非线性方程：2222cossin1sinsinsin022kkSf kSN 可以求得k与的关系。现考察0,90,180,270时的数值相速度v。以这些角度代入22221cossinsinsinsin22SkkSN，可得数值波数 121sinsinSkSN。数值相速度为 11211sinsinsinsinvvkSSNSNSN v为实际地震波传播速度。也可直接求出沿45,135,225,315方向数值波数和数值相速度 12 21sinsin2SkNS，112sinsin2vvSNNS。数值相速度的数值。设0.5S
44、，则由 11sinsinvvkSNSN（网格方向）和 112sinsin2vvSNNS（对角方向）20N时算出的数值相速度分别为0.996892v，0.998968v；10N时算出的数值相速度分别为0.987264v，0.995817v。可见，这两个数值相速度都小于真实速度，波在不同方向的传播速度不同，这就是数值计算中波的色散现象。数值速度的值与真实速度值的差异还与N有关，N越大，数值速度的值与真实速度值的差异越小。色散程度与空间采样点密度有关。用色散程度与空间采样点密度有关。用 22221cossinsinsinsin22SkkSN 计算出数值波数后再计算相速度并以理论相速度归一化后画出数值
45、计算出数值波数后再计算相速度并以理论相速度归一化后画出数值相速度与传播角度和每个波长的采样点数之间的关系曲线相速度与传播角度和每个波长的采样点数之间的关系曲线。在获得在获得该曲线时，该曲线时，0.5S。归一化数值相速度随传播角度和采样密度之间的关系曲线 网格长度固定，对不同的频率，低频采样率高，高频采样率低。在同一方向，高频的传播速度低，低频的传播速度高。100Hz5025750h=4 mv=2000m/sf=80Hz=25mttx 150Hz6030900120h=4 mv=2000 m/sf=160Hz=12.5mttx 90Hz主频的子波主频的子波h=4 mv=2000 m/stx 30Hz主频的子波主频的子波,t0=0.01 s 220000exp1 2ff ttf ttt 由图可见，每个波长采 5 个样点时，数值相速度与传播方向有很强的关系，而每个波长采 10 个样点时，数值相速度与传播方向的关系不是很强了。因此，每个波长至少采 10 个点以上，网格色散才不太明显。这是指导数值模拟中如何选取空间采样间隔的原则之一。
