任意偶数阶精度有限差分法数值模拟

1、|维普贽讯 nHpy/www.cqvip.coTfi1993 年 2 月鬼方刀"卜沁兀/任意偶数阶精度有限差分法数值耨拟刘洋"李承楚电永/卜"$*>r刘洋李承蜒牟永光、任度有覘差分法数值模拟.石油地球输理勘探"卿8,33（1）：110本文航Taylor尿歎晨开戎出发.推导出了任总阶导魏的任意偶嫩阶精度差鲂格戎”文中绐出 了计算栢应差兮系菠的卧丈并计隼出了一阶、二阶、三阶和田阶导孰的若干阶精度圣令系轨。最 后对论了声波方程的任當偶數阶差令數世孵的穩定性条怜.并进行了数值棋攪试膛:黠果表期差 令精度議奇、械粧越小、戳值模拟的效果越奸.导敎 Tayhr正

2、潼模挞主题词 堆曹数蚩生里 數值棋拟 电週毛注_波动方程 構度偶數阶"ABSTRACTLiu Vang U Chengchu and Mou Yongguang, Finite-difference numerical modeling of any even-order accuracy. OGP. IPPS, 33(1)匸 110Even-order accurate difference schemes of any-order derivatives can be derived from the expansion of TayJor seThe formulae for

3、calculating relevant difference coefficients are offered hre Liu Ysr. Geophysics Research LnstLtur<!1)电砂班吟”1 of Gsosrwncei Petroleum University Reservoir Streex rChang- ping CounlVi Beijing City *102ZOO 卒文于1曲7年3 H 21 H收理Land they are used to get some-ordcr accurate difference coefficients of one(

4、two three and four-order derivatives The thesis also analyses the stable conditions of numerical solutions of any even-order differences to acoustic wave equation 自nd describes trial numerical modeling* It is shown that tht higher accuracy of difference brings less frequency difipersion and better n

5、umerical mod亡l- ing effectsSubject heading； seismic data processings numerical modeling, finite difference, wave equation * derivative* Taylor series« forward modeling, accuracy > even order地理波数值模拟是研究复杂地层资料解释的有效辅助手段。地廉波数值摸拟的方法主要 包括两类，即波动方程法和几何射线法。波动方程散值模拟方法实质上是求解地靈波波动方| 维普资讯 h ttp:/第裁雑 第1期刘洋

6、等：任意假數阶精度有限差林袪数值横拟p （工0 +人）一讥竝如2A=卫工。十 2户（工爲 + 屮（工炉一h、十os、(1-b)程'因此模根的地義波波场包含了地K波传播的所有信息，但英计算速度相对于几何肘钱法要 慢"几何射统祛也就是射线追踪苦，属于几何地廉学方法口由于它将地健浪波动理论简化为射 线理论主要考慮的是地屣浪传播的运动学特征，缺少地恳波的动力学倍息因此该方法计算 速度快。由于菠型方程模拟包含了丰富的波动信息，为研究地健波的传播机理和复杂地层的解 释提供了更參的佐证，所以諛动方程歎值模拟方法一直在地廉模拟中占有星要地位.地廉波潑动方程数疽模拟方送主要包括有限差分迭.有限

7、元法和博里叶变换法。博里叶变 换吿由于利用空闾的全部信息对波场函数进行三角画載插值*所以能更加精确地模拟地恋浚 的传播规律"同时，利用快邃傳里叶变换（FFT）进行计算*还町以提高运算效率“该方法的主要 优点是精度高，占用内存小，缺点是计算速度较慢.波动方程有限元法的主要优点是适宜于模 拟任意地质体形态可凶任意三角形逼近地层界®S保证复杂地层形态模拟的逼真性+其缺点 是占用内存和运算Jt峋较尢“有限羞分法的主要优点是计算速度快，占用内存小：缺点是精度 低，仅适合于相对较简单的地质模型.本文的目的在于建立一种能兼顾计算精度和计算速度的正演方摆便根据不同的卄算 精度和计算速度的耍

8、求，灵活运用“2K-1阶导数的2N阶精度展开式苜先分析索规有限差分方法的实现过程#有醒差分法是采用筆分的方医求取波动方程中 对空间的僞导数项"一般情况一阶和二阶空间导数采用如下公式进疔计算(1-a)上述二養分裁式均具有二阶精度f其中0（表示嘉的同阶无穷小绘。由于进行差分时、 仅利用相邻两点来计算，所以精度低。由此可见，当舉与差分的点数增加时*应当可以获痔更 高的精度"在此启发下，我们从丁矽1"级載展开式出发*推导出了任意阶导数的偶数阶精度 差分椿式，并用于声波波动方程的不同阶差分精度正演模拟，从而验证了差分格式的正确性和 有效性。设少（£雀续且具有2N

9、+ 2K-1阶导敦，则日:乂&= V S聞严 + c严"5）椚 + Oih-1）上式证明见附录 X 其中MN+K 1.6 = 0沁裤=一“刑=且=(-)；七K_mMM+l胳心J其中衷示的行列式，卫m为 1°5-ir<m + l)c(m + 2)D_ AT _F2Z 1(mu2d(m + l)35 十 2)2d* Af3V|2JC-4i丹T.CmI-l)ga烁+ 1严7SVE + 2 严T2iK (m-iyK(構+ 1严(m + 2 严gjtJC+S严+2*(mi严处Cm + 1严+ w沪7fljZiW 2d严?(m*-l)2M-=i(m + 1J3M -2f1

10、(m + 2产，a>»心巴当K=1时.即对于一阶导数的展开式系数，式C3； J4)可简化为<il)N" 1 时* Ci l/2t *) l/6i(2) N=2 时< 门=啓3, c：=l/12t = 1/30：(3) N>2 时，有=一 1尸|ffl- JNn<mi 鬥 n s - 1I.% = _ <2 + l)!舟讦eJ(4Nf g时有(-ir41ffl = ST<7>| 维普资讯 h ttp:/第裁雑 第1期刘洋等：任意假數阶精度有限差林袪数值横拟可以证明*上述一阶导数的2N阶箱度展幵系数喪达式，与文献1中的系数表达式是

11、等价的。2K阶导数的2N阶精度展开式设pGr)连续，且具有2N + 2K阶导数，则響0= 2才讹十畑户” 5少+ 0(丹+】)上式证明见附录色其中MhN + K-1* =j(一尸十气盛釘-h1K(10>m= 1且(11)| 维普资讯 h ttp:/第裁雑 第1期刘洋等：任意假數阶精度有限差林袪数值横拟| 维普资讯 h ttp:/第裁雑 第1期刘洋等：任意假數阶精度有限差林袪数值横拟(即一 (25/+ 2) E §mE，W+ 匕当K = 1时，即对于二阶导数的展开式系数，式<ll)s(12)可简化为| 维普资讯 h ttp:/4石油地球物理勘探< 1 )N = 1 时

12、* '严2 - G = 1. - /12j(2>N«2 时，応口一5/2,门=4/3 c2=-l/12t ： = l/OOi<3)JV>2 时有C 1.),+1n 产l-.r#n,o_1T<<13)5 -Ymu |rt /?! t JV )口£嗨？一产)tr - w2)p ILfm hi -|-2匸咖诃+2"(14)N (2N + 2)1J（4）押f 00时、有?S = « - I)】沙5 = 1.2宀)（% 0(16)偶数阶精度差分格式的构造单变导數的鬻樹阶糖度差分格式令H为的近似值户尸为P（R在工=朋处门阶导数的

13、近似值匸等于2K或2K-l,; = 0，±l.±2t-. h为网格大小、则根峯式嘗人（1（0可得匸阶导数的2AT阶精度有限差分格氏为理_赛一1x£, = S 学匕）E= 代一 H用差分茶数.的表达式已曲式虫、（】1）给出。前面已经证明，2K 1阶导数的2N阶精度 差分系数与持求导点呈奇对称設K阶导数的阶帝度養分系数与待求导点呈偶对称"表U 表氛表3和表屮分别给出了一阶、二阶、三阶利四阶导数的2N阶精度差分索数g的值同时 还给出了栩应阶段的谋差系数斑的值口可以看出随着精度的增加.参与计算的点也瞎加截 断课菱系数逐斷减小.壹1 一曲导數的偶養fifrfil度

14、肓限盖分格式系數惡阶导#fc d/djT的偶数阶構匿塞分累數截断逞婪抿数慚数）5J口GGg臨盼C200.500K66fi7 vlo-'40钿裁一仇1叫特3- 3333 > 10"-60 0. 7500 4O.Ulti7-7,1429 X 10_,£n r) 30(jn -V. 3000 D 03tf QF 003G1. FB73X10-J1&U D. fSS13 0. 23S1 0. 0595 -Q, 0Q99 0. 00Q3-3. 6075 X 10 _11200. 8571 D, 2579 0.0794 -0.0179 fl,. 0026 -0.0

15、002a. a25oio-sU0 U. 8750 -0, 2917 D.0972 0. 065 0. 0053 0. 007Q 0. U000-L 945a lt>E1500. S889 -0, 3111 C. J131 -0,0354 0.00F7 -0.0016 0.0002 -<i+'X<JOK*-*4. 5706 a 10f«giiift-IFVPF90 LOOOO -0, 5000 0, S333 -0. 2S00 0. 200C -0.16fi7 0. 1429 -0+1250 口第勢卷 第1期训洋等：枉盍隅数阶搐度有限绘井法觀值撲拟5*2二阶显

16、魏的幺散峭箱匱有思差分轄式累搭農度二阶显敘护的隅数阶糟度澄分專数蝦断误差廉数阶SHGctCiCtcscic；th*”2-2. 0000 I. WOD-8. 33a3X10"f1-2. 5000 1. 3333 -0. 08331, 1111X10-= 6-2. 7222 1. 5000 -0,1500 a 0111-1, 7857 X10 J8-2, U12 1, 6000 - 0. 2000 0. 024 - 0. 00183.1746 X 10102- 9272 1. 6657 0. 2381 0. 0397 0, 0050 血 0003-6. Q125XIO-* 12-2. 9

17、82S1.7K3 -0. 2679 0, 0529 -0.0089 P, 0010 0, 00011. 1893XIO-514X 023 1. 7500 0. 291?厲 064B -0.0133 0, 0021 一 0002 0, OOOQ一鱼 4381X1O_,lti一玉 0543 1. 7778 一 乱 31 H 也 0754 一 9177 0. 0035 一 0005 6 00010. OOOD*i*B4*5,07B4XL0-Tq<x>iLI¥-3. 2&9& 2. 0000 一氏 5000 此 2222 - 6 1250 0. 0300 一山毎稠

18、 0. 0408 広 031J gV*3三曲导数的数阶精度有隔淮分格式黑裁議靖度三阶导散卍ALX的fll数阶情度蹇井慕散一 -" "掘断溟差果魏（阶数）51C“以内心C门20 -1.0000 0.5000-2. 00 X 10"'40- 1. 250 1.0000 -0.12505. 8333 X IO-J60-2, 0333 L4G3S - 0. 3000 0. 0292一1 355SX 10*Ja03181 1, 7337 -D.&S30 0. <34 -6 D城3. 16&0XW3100-2. 5271 1. 9951 - 0,

19、 657J 0.1531 -0. 022S 0, D016-1. 4495 X1Q-4IE0-z. 6S69 2. 2DS1 -0.171 0. 2306 -0.0468 0. 0061 -0. 0MM1. 7620X10-'140 -2. &129 2. 3S45 -a 914 0. 3107 -&. 0777 0.0140 -0. 0016 0, 0(»1-4. 1837X 10乍弘中»Ebj.vEh»iJ.Hi-*!*!* *¥I4-V%卡«%4电10003.6809 3；81SJ 臨 54卸 L 717S 1,

20、1700 0. 7929 -0* 5Z98 0. 3467 “2. 56O9X 10w«4四脑导數的赠阶精度荷IB建分格式«第度（阶数）四阶导的度君金奉数fjrBCf26. 0000-4.00001+ CODO1.6667X 10_l449. J33J6. 50QQ2. 0000-M6&7Z.9167X 10_<&1K 3750B. 19332. B167D. 40000.0292-5.4233X IO-3R吃 7417-轧 27223. 4675-0” 64400. 084-0. 0054l O5BCX10-*1013.717410.10863. 9

21、902-0.87 &40L1531-0.01B10.0011-2.1284 XW121.14- 4479-10. 747?4.4163-1,0894L0. 2306-0. 0375!(L 004G*-D.0DO2Jk 405OX1O-51M187196>74. 723797.63521-3. 39061.1,7179n-a 936010. 5286Q. 3027 3. 0917 X 10-n混合傭导数的偶数阶精度差分格式对丁混合偏导数护】+ "/站厶皤】的高阶差分计算可腺先沿工方向求关于工的厶阶导 数，在求出这些导数的基础上然后沿工方向求关于忙的厶阶导数。根据上述思路即

22、可构造 出相应的任意偶数阶精度的差分格式。| 维普资讯 h ttp:/1998 年令打丄为的近徵值心'为心沽在卫=丙“=血St美于工的L.阶山的 厶阶馆导效的近拟值厶等于2£ 1或等于2KltLt于2K1或士1* 士魯r f = S±】，±2rr h.为沿工方向的网格大小，也为皤黑方向的网格大小，则可得混合偲导魏 的偶数阶精度有限差分搐式曲08)| 维普资讯 h ttp:/其載断溟差为严+島严其中M】 = M+KlX 购=皿+应Tm晋为厶阶导散的 2M阶稱度差分为厶阶导数的阶稱度差分系数“肾和蹲两个系数序列可分别 由式（3） All）求得口为了检测所构造的

23、望分格式的精度，我们进行了声波漩功方程的正演摸拟试验臼二维弾性 介质中的声波方程为(.19)式中关于空间的导妣，采用2N阶辎度差分*关于时同的导数，采用二阶中心差分在进行正演计算前，首托将变工和事离散化"设空阖用格离敵长度分别为 般取 2=2=&、时间步长为心“ &和 M应满足如下稔定条件（详细推导见附录C）(20)其中“为不超过“的最大奇数.表5为針对不同的糟度所廿算的穩定因子，其表达式为«5声讓連动方程礪撤阶持度基分格式的定因子阶ft24占101214161820co程定因子J 70710.6124D. $7520»55450, 541 30.

24、 E31B6 524£0. $1B90. 51436 51050,4502理论模型为一均匀介威，速度v-2 500为点健源，莅于模型中央"网格大小为山=加="=20皿时间母規为mt图1氣图丄乩图lc和图Id分别为二阶、四阶、b 阶和二十阶差分格式所模拟的i=200 ms的瞬时波场图口可以看出，当特度较低时，頻前现象 严重*随着阶数的增加，频散现象逐渐降低*模拟肢场的精度逐渐提高-显然，随着精度的增加，参与计算差分的点数也發增加，计算盘随之增大。所以，要获得较 高的稱度，必殖花费校参的丹算时间V耍想快速得列正演模拟结果丫就要以降低精度为代价°总 之，根据所

25、要隶的计算速度和精度我们总能栈劃一个阶数議佳的差分格式来进行计算.第蓟卷 第1期刘洋綁任意偶数阶轄度有限盖分医数值模曲7(d)图】旅同酚稱度有限差分医声哉数値棋披誠场点厘源倉览】对比图3】二阶八仍四阶，FBfriub二十徹車文从壮W5缓数展开式出发.推导出了任靈阶导数的偶数阶養分格式.通过声波渡动 方程的数值模姒试验，表明了塞井格式的正确性和有效性。本文推导出的任裁偶数阶精度差分格式，能兼顾计冀精度和计算速度两方面的不同要求. 当更側重于计算速度时，选用低阶精度差分格式；当追求高精度时选用高阶精度差分格式。如 果同时要考虑计算理度与精度、可选用适当精度的差分格式。1 Fornberf B” T

26、he pseuduapctra methcdfcmpsrLNoriH with finite tJifferences far the cystic wave equation, Gto physic st L&B7*52(43 ：4B35012 St金甫，关治阳橄分方程数值瞩楼，清华天学出版社,19893笔禎祥郑抽种、地倉披殖值棋拟与俯移成橡*河南科学技术岀腿社/94| 维普资讯 h ttp:/石油地球物理勘探附录A2K-1阶导数的2N阶精度展开式设护小有7N十氏一1阶导数，则戸（工）在工=斗+同处的厶V十2K- 1阶Taylor级数 展开式为CA-I)pg += pa：j + 工护

27、'g、十 o护 J* - L英中诃=M一M+lf1,M；M=N + K 表示扒Q在工=吐处的?阶导数。为茯得适=町处的2K-1阶导数的2N阶精度展开式，将上述2M个方程分别乘以“対 t- IM-1WQ -J *C| h -Cm -i J 然后相加并化简，有£煜庄_ 17?*严亠_乞 h丽* - -jMJPj如 + 乞 Jpg 郴h )軒壬M扁¥11S 2 门yy.乙 m + 1屮”.W+hCtq) +5 肝川) M2”| 维普资讯 h ttp:/<A-2J将上式写成矩昨形式，有可以看出一口严Tr 亠I）严八（.M 1）yK *-(M - i)E«2K

28、- L(-1)< -K旷 TM- 13V-J(2K严一 pJT-ip.W(iW _ Dskj<M - 1 严M1KyAFK(A-3)l(A-4)刘洋等：任意偶数阶精度肓限差分法数值模拟J 二一 c_e(rn l»2w *A/)则式可化简为-(I2)0(2V聲C胚屮n«0rc3严旷产7(AF)A-Zr r2Cn0(丈严-i(2S)K_1=(2K - l)|/2©严(2广(M2严V09i4i=rtf 1鼻8严-：I血严一;_ 0 一其系数矩阵为范德蒙矩阵。解此方稈即可得出相应的5(期=1 本附录B2K阶导数的2N阶精度展开式设仇刃具有2N + 2K阶导数侧p(Q在工=竝十泌处的2N戚K阶Taylor级数展开式 为汰+汕/, sr户5十心 =pg + 丫 竺