緊致交錯(cuò)網(wǎng)格優(yōu)化差分系數(shù)二維聲波方程數(shù)值模擬

汪 勇 王 鵬 蔡文杰 桂志先

(油氣資源與勘探技術(shù)教育部重點(diǎn)實(shí)驗(yàn)室(長(zhǎng)江大學(xué))，湖北武漢 430100)

0 引言

地震波場(chǎng)數(shù)值模擬是應(yīng)用數(shù)值方法模擬地震波在地下介質(zhì)中的傳播過(guò)程，得到在地面或地下各觀測(cè)點(diǎn)地震記錄的一種方法，是地球物理勘探和地震學(xué)的重要基礎(chǔ)。隨著有限差分?jǐn)?shù)值模擬技術(shù)的發(fā)展，針對(duì)實(shí)際應(yīng)用中提高有限差分計(jì)算效率[1]、模擬精度[2]、算法穩(wěn)定性[3-4]，處理復(fù)雜介質(zhì)[5-7]、吸收邊界條件[8-10]和壓制數(shù)值頻散[11-13]等多方面需求，已研發(fā)了多種實(shí)用技術(shù)，并取得相應(yīng)研究成果。

提高差分格式數(shù)值計(jì)算精度最直接的方法就是在模擬時(shí)增加(加密)網(wǎng)格節(jié)點(diǎn)數(shù)，但隨之也增加了運(yùn)算量和存儲(chǔ)空間。緊致差分方法能較好地解決這個(gè)矛盾； 同時(shí)緊致差分是一種隱式差分格式，穩(wěn)定性較強(qiáng)。這些優(yōu)勢(shì)使其成為目前頗受關(guān)注的有限差分方法之一，并被廣泛應(yīng)用于聲波、彈性波和復(fù)雜介質(zhì)等地震波場(chǎng)數(shù)值模擬中[14-22]，提高了地震波場(chǎng)數(shù)值模擬的精度和計(jì)算效率。

而Madariaga[23]提出的交錯(cuò)網(wǎng)格差分格式既可提高數(shù)值模擬的局部精度，還能加快收斂速度。將交錯(cuò)網(wǎng)格與緊致差分技術(shù)結(jié)合，可進(jìn)一步提高數(shù)值模擬的精度和壓制數(shù)值頻散的能力[24-26]。

為使差分格式在較大波數(shù)范圍內(nèi)減少頻散誤差，Kim等[27]基于頻散關(guān)系保持(Dispersion relation preserving，DRP)的思路優(yōu)化了緊致差分格式，指明采用這種優(yōu)化差分系數(shù)方法可使差分算子最大程度地逼近空間偏導(dǎo)數(shù)算子。Liu等[28]基于頻散關(guān)系提出了優(yōu)化的時(shí)空域有限差分系數(shù)，在不增加計(jì)算成本的情況下可顯著提高模擬精度。Zhang等[29-30]應(yīng)用最大范數(shù)的目標(biāo)函數(shù)及模擬退火算法求解目標(biāo)函數(shù)，優(yōu)化了一階和二階常規(guī)中心差分格式的差分系數(shù)。Liu[31-32]利用最小平方法優(yōu)化了二階導(dǎo)數(shù)中心差分和一階導(dǎo)數(shù)交錯(cuò)差分系數(shù)。Yu等[33]基于DRP思想，優(yōu)化得到5階精度的組合型緊致差分系數(shù)。Ren等[34]以優(yōu)化后時(shí)空域差分格式進(jìn)行聲波和彈性波方程的數(shù)值模擬。Yang等[35-36]利用極值逼近(Minimax approximation)算法優(yōu)化了交錯(cuò)網(wǎng)格差分系數(shù)，當(dāng)波數(shù)較大時(shí)可提高模擬精度。

本文首先討論了一階導(dǎo)數(shù)的2N階緊致交錯(cuò)差分格式的建立方法； 基于頻散關(guān)系保持的思想，通過(guò)最小平方法建立了最小數(shù)值波數(shù)誤差的目標(biāo)函數(shù)，利用拉格朗日乘數(shù)法求解目標(biāo)函數(shù)，得到優(yōu)化后的4～10階精度差分系數(shù)； 然后分析和對(duì)比優(yōu)化前、后的緊致交錯(cuò)網(wǎng)格差分格式的模擬精度、數(shù)值頻散和聲波方程的穩(wěn)定性條件，并用優(yōu)化后的緊致交錯(cuò)差分格式對(duì)一階速度—應(yīng)力聲波方程組進(jìn)行數(shù)值模擬。研究結(jié)果表明，采用優(yōu)化差分算子能有效壓制數(shù)值頻散，提高差分近似導(dǎo)數(shù)的精度。

1 二維聲波方程的高階差分近似

Madariaga[23]早年提出的波動(dòng)方程交錯(cuò)網(wǎng)格數(shù)值模擬方法的差分精度為O(Δt2+Δx2)； Levander[37]應(yīng)用交錯(cuò)網(wǎng)格方法求解P-SV波方程，將差分精度提高到O(Δt2+Δx4)； 董良國(guó)等[3]提出了一階彈性波方程交錯(cuò)網(wǎng)格高階差分解法，其差分精度達(dá)到O(Δt2M+Δx2N)。現(xiàn)今，交錯(cuò)網(wǎng)格有限差分方法已成為一種高效且應(yīng)用廣泛的數(shù)值模擬方法。

根據(jù)彈性力學(xué)分析，二維情況下一階應(yīng)力—速度聲波方程組(假設(shè)體力為零)可表示為

(1)

式中：P為應(yīng)力(標(biāo)量)；Vx和Vz分別表示x和z方向質(zhì)點(diǎn)振動(dòng)速度分量；ρ為介質(zhì)密度；v為地震波速度。

1.1 時(shí)間2M階近似

用交錯(cuò)網(wǎng)格法進(jìn)行聲波方程數(shù)值模擬時(shí)，速度和應(yīng)力分量分別是在t和t+Δt/2時(shí)刻計(jì)算的，其中Δt為時(shí)間步長(zhǎng)。用泰勒公式可得速度分量Vx的2M階精度時(shí)間遞推式

(2)

當(dāng)M=1時(shí)，即可得到時(shí)間二階精度差分近似

(3)

同理可得速度分量Vz和應(yīng)力分量P的時(shí)間二階精度差分近似

(4)

(5)

利用應(yīng)力—速度方程組，將式(3)～式(5)中的變量對(duì)時(shí)間的導(dǎo)數(shù)轉(zhuǎn)換為對(duì)空間的導(dǎo)數(shù)，則可得

(6)

1.2 空間2N階緊致交錯(cuò)差分近似

緊致交錯(cuò)差分格式最早由Nagarajan等[38]提出，用于大渦模擬(Large eddy simulation)問(wèn)題，該差分格式利用4個(gè)網(wǎng)格節(jié)點(diǎn)可得到6階空間精度。隨后，Boersma[39]提出最高可達(dá)12階空間精度的緊致交錯(cuò)差分格式，用于可壓縮流體的Navier-Stokes方程的數(shù)值模擬。類似于常規(guī)交錯(cuò)網(wǎng)格數(shù)值模擬，緊致交錯(cuò)網(wǎng)格法求取變量的導(dǎo)數(shù)時(shí)，也是在相應(yīng)的變量網(wǎng)格點(diǎn)之間的半程上進(jìn)行的。

一階導(dǎo)數(shù)的2N階精度緊致交錯(cuò)差分格式為

(7)

式中：h為空間網(wǎng)格尺寸；a和bn是差分系數(shù)。利用式(7)近似計(jì)算式(6)中變量對(duì)空間的一階導(dǎo)數(shù)，確保2N階差分精度的關(guān)鍵是確定差分系數(shù)a和bn。利用泰勒級(jí)數(shù)展開和待定系數(shù)法，即用

(8)

可求得緊致交錯(cuò)網(wǎng)格2N階空間差分精度的差分系數(shù)。如當(dāng)差分精度2N=10時(shí)，上述方程組共有5個(gè)方程，可確定a和bn(n=1～4)共5個(gè)差分系數(shù)。

此外，從上述公式可看出，常規(guī)交錯(cuò)差分格式(Staggered finite difference，SD)求取中心點(diǎn)處的導(dǎo)數(shù)值時(shí)，僅用到周圍網(wǎng)格點(diǎn)的函數(shù)值，而緊致交錯(cuò)差分格式(Staggered compact finite difference，SCD)卻還用到相鄰點(diǎn)導(dǎo)數(shù)值。緊致交錯(cuò)網(wǎng)格法只需利用2N-2個(gè)節(jié)點(diǎn)即可達(dá)到2N階空間差分精度，而常規(guī)交錯(cuò)網(wǎng)格法則須使用2N個(gè)節(jié)點(diǎn)，因此SCD方法可節(jié)省存儲(chǔ)空間。表1列出4～10階泰勒級(jí)數(shù)展開得到的差分系數(shù)。

表1 緊致交錯(cuò)差分格式一階導(dǎo)數(shù)的差分系數(shù)

1.3 差分格式的矩陣形式

(9)

其中

向前差分為

向后差分為

據(jù)式(9)，由應(yīng)力(標(biāo)量)場(chǎng)P求取其空間一階導(dǎo)數(shù)

(10)

上述各式中差分系數(shù)a和bn(n=1～4)由表1給出，可求得一階導(dǎo)數(shù)的4～10階空間差分精度近似值。

2 緊致交錯(cuò)格式差分系數(shù)優(yōu)化及分析

在地震波場(chǎng)數(shù)值模擬中，為得到清晰的波場(chǎng)記錄，需要高精度的空間離散差分算子。差分算子的精度取決于差分系數(shù)和差分階數(shù)，差分階數(shù)越高，差分算子越精確，但其運(yùn)算量也會(huì)隨之增加。優(yōu)化差分系數(shù)可使差分算子最大程度地逼近空間偏導(dǎo)數(shù)算子，本節(jié)對(duì)一階導(dǎo)數(shù)的緊致交錯(cuò)差分系數(shù)進(jìn)行優(yōu)化，進(jìn)一步提高緊致交錯(cuò)差分格式的分辨率，并分析優(yōu)化后的差分格式的精度、頻散和穩(wěn)定性。

2.1 優(yōu)化差分系數(shù)

首先對(duì)一階導(dǎo)數(shù)的緊致交錯(cuò)差分格式(式(7))進(jìn)行數(shù)值頻散分析。令

(11)

式中：h為網(wǎng)格尺寸；k為真波數(shù)。若數(shù)值模擬時(shí)不存在誤差，則B=(ik)A。但在差分計(jì)算中，可能會(huì)產(chǎn)生不同結(jié)果，即產(chǎn)生數(shù)值波數(shù)(又稱修正波數(shù))。可令B=(ik′)A,k′為修正波數(shù)。

將式(11)代入式(7)，并用歐拉公式化簡(jiǎn)可得

(12)

式中：ω=kh;ω′=k′h。

在理想情況下，若不存在數(shù)值頻散，則ω′=ω。它們的差異越大，說(shuō)明該方法的數(shù)值頻散越嚴(yán)重；反之則說(shuō)明該方法能更好地壓制數(shù)值頻散。依據(jù)修正波數(shù)盡可能在大范圍內(nèi)接近真波數(shù)，以4點(diǎn)緊致交錯(cuò)差分格式為例說(shuō)明優(yōu)化的思路和方法。

由式(7)可得

(13)

為了滿足2階和4階精度泰勒公式截?cái)嗾`差的要求，式中差分系數(shù)a、b1和b2須滿足下列方程組

(14)

即為式(8)表示的前兩個(gè)方程。

為了確定上述方程中的未知差分系數(shù)，按照Tam等[40]的思路和方法，在某個(gè)選定波數(shù)范圍內(nèi)，確定式(13)中3個(gè)未知差分系數(shù)，使得修正波數(shù)盡可能地接近真波數(shù)，故定義誤差函數(shù)為

(15)

(16)

其中

(17)

由式(14)和式(16)中的3個(gè)方程，可確定3個(gè)優(yōu)化后差分系數(shù)，即a=0.185715、b1=0.985714、b2=0.128572。將優(yōu)化后的差分系數(shù)代入式(13)，則得到優(yōu)化后的緊致交錯(cuò)差分格式(Optimized staggered compact finite difference，OSCD)。

需要說(shuō)明的是： 4階精度OSCD格式在式(14)基礎(chǔ)上，需加入一個(gè)表示DRP要求的積分方程，所以不能嚴(yán)格滿足空間6階精度泰勒公式截?cái)嗾`差的要求，只能達(dá)到空間4階精度近似； 而常規(guī)SCD格式的3個(gè)方程即能滿足空間6階精度截?cái)嗾`差的要求。因此，若要達(dá)到2N階空間差分精度，SCD方法只需利用2N-2個(gè)節(jié)點(diǎn)，而OSCD方法需使用2N個(gè)節(jié)點(diǎn)，與SD方法使用的節(jié)點(diǎn)數(shù)相同。

采用同樣方法，對(duì)式(7)表示的6～10階OSCD格式做差分系數(shù)優(yōu)化，優(yōu)化后差分系數(shù)列于表2。

表2 優(yōu)化的緊致交錯(cuò)網(wǎng)格一階導(dǎo)數(shù)的差分系數(shù)

本文方法與文獻(xiàn)[33]方法基本一樣，不同之處在于： 文獻(xiàn)中原差分格式是3點(diǎn)6階組合型緊致差分格式，為了加入代表DRP思路的積分方程，在原差分格式基礎(chǔ)上構(gòu)造了新格式，增加了1個(gè)網(wǎng)格節(jié)點(diǎn)，降低了1階精度，得到的是4點(diǎn)5階迎風(fēng)型組合緊致差分格式； 本文的優(yōu)化格式與原緊致交錯(cuò)差分格式相比，在使用同樣網(wǎng)格節(jié)點(diǎn)進(jìn)行差分計(jì)算時(shí)，降低了2階精度。

本文求取優(yōu)化差分系數(shù)的方法本質(zhì)上是根據(jù)泰勒級(jí)數(shù)展開式，而文獻(xiàn)[35]和[36]則是根據(jù)函數(shù)逼近的思路。文獻(xiàn)[35]首先根據(jù)緊致交錯(cuò)差分格式的頻散關(guān)系和極小極大原則建立了誤差函數(shù)，在所有逼近格式中尋找誤差最小且滿足誤差標(biāo)準(zhǔn)的三角函數(shù)多項(xiàng)式，求取過(guò)程中使用的是列梅茲(Remez)迭代算法。文獻(xiàn)[35]指出該方法求得的差分系數(shù)在中等和大波數(shù)條件下比泰勒級(jí)數(shù)得到的要好，但算法復(fù)雜，實(shí)現(xiàn)難度較大，這也限制了該方法的應(yīng)用。此外，該方法得到的差分格式是在給定的誤差標(biāo)準(zhǔn)條件下得到的，但不能明確差分格式的近似精度，即通常意義上的2N階。

2.2 頻散分析

對(duì)優(yōu)化前、后一階導(dǎo)數(shù)的緊致交錯(cuò)差分格式進(jìn)行數(shù)值頻散分析，通過(guò)數(shù)值波數(shù)與真波數(shù)的比值評(píng)判優(yōu)化效果并確定適用的空間網(wǎng)格尺寸。

據(jù)式(12)，修正波數(shù)與真波數(shù)之比可定義為

(18)

在理想情況下，若不存在數(shù)值頻散則波數(shù)比R恒等于1。R偏離1越大，說(shuō)明該方法的數(shù)值頻散越嚴(yán)重； 反之則說(shuō)明該方法能較好地壓制數(shù)值頻散。利用表1和表2中不同精度的差分系數(shù)可求取優(yōu)化前后的波數(shù)比曲線。為了比較，計(jì)算得到常規(guī)交錯(cuò)差分格式(SD)的波數(shù)比，計(jì)算公式見文獻(xiàn)[35]，不再贅述公式。

圖1為3種差分格式在不同差分網(wǎng)格節(jié)點(diǎn)個(gè)數(shù)(即差分算子長(zhǎng)度)情況下的波數(shù)比曲線，如圖1a為4個(gè)差分節(jié)點(diǎn)，根據(jù)前文分析，SD、SCD和OSCD格式分別可達(dá)到4階、6階和4階差分精度。

此外，從圖中可見：低階差分精度時(shí)，三者波數(shù)比曲線差別最大； 隨差分精度增加，差異變小。顯然，OSCD格式在壓制數(shù)值頻散方面的優(yōu)勢(shì)主要體現(xiàn)在低階差分精度時(shí)。因此，可用較低(4或6)階OSCD格式達(dá)到高階SCD或SD格式的效果。

圖1 不同網(wǎng)格節(jié)點(diǎn)數(shù)三種差分格式的波數(shù)比曲線

2.3 精度分析

為了比較SD、SCD和OSCD三種格式的近似精度，對(duì)比分析其截?cái)嗾`差。從表3可見： 在達(dá)到相同空間近似精度的條件下，4和6階精度的OSCD格式具有更小的截?cái)嗾`差，8和10階時(shí)誤差介于SD和SCD之間，所以低階差分時(shí)OSCD格式的精度優(yōu)勢(shì)更明顯。

表3 不同差分格式的一階導(dǎo)數(shù)截?cái)嗾`差主項(xiàng)系數(shù)比較

利用一維平面簡(jiǎn)諧波初值問(wèn)題，比較SD、SCD和OSCD格式的數(shù)值模擬精度。一維平面諧波初值問(wèn)題可表示為

(19)

該偏微分方程的達(dá)朗貝爾解，即精確的解析解為

(20)

圖2 三種差分格式的一維聲波方程數(shù)值模擬快照

通過(guò)定量計(jì)算，可得到該時(shí)刻、4000～5000m之間，OSCD、SCD和SD格式數(shù)值解與解析解之間的相對(duì)誤差依次為0.33%、0.77%和0.97%，這正是由于它們?cè)谟?jì)算空間導(dǎo)數(shù)時(shí)存在不同截?cái)嗾`差造成的，這也與表3的理論分析結(jié)果相一致。相對(duì)誤差定義為

(21)

2.4 穩(wěn)定性分析

按照董良國(guó)等[3]和杜啟振等[4]使用的Fourier方法對(duì)差分格式進(jìn)行穩(wěn)定性分析，略去推導(dǎo)過(guò)程，直接給出式(22)表示的二維聲波一階速度—應(yīng)力方程組(2M，2N)階差分精度的緊致交錯(cuò)網(wǎng)格差分格式的穩(wěn)定性條件。

(22)

其中

表4 二維聲波方程不同差分格式的Courant數(shù)

3 模型試算

采用OSCD差分格式對(duì)簡(jiǎn)單模型和復(fù)雜的Marmousi模型進(jìn)行二維聲波方程進(jìn)行數(shù)值模擬，檢驗(yàn)該方法數(shù)值模擬的適用性。

3.1 均勻介質(zhì)模型

為比較SD、SCD和OSCD三種差分方法的差異，設(shè)置模型尺寸為3000m×3000m，縱波速度為2000m/s，模型中心加載正弦衰減子波，峰值頻率為30Hz，空間網(wǎng)格為12m×12m，時(shí)間步長(zhǎng)為1ms。采用三種格式的時(shí)間2階、空間4～8階差分精度對(duì)一階速度—應(yīng)力聲波方程進(jìn)行波場(chǎng)模擬，得到同一時(shí)刻的應(yīng)力P分量的波場(chǎng)快照(圖3)。

在時(shí)間和空間網(wǎng)格尺寸相同的情況下，圖3a所示的4～8階SD格式的波場(chǎng)快照頻散非常嚴(yán)重；圖3b的4～8階SCD格式波場(chǎng)快照好于SD格式，但4階SCD格式頻散依然明顯，6階在0°和90°方向上也存在較明顯頻散，8階精度時(shí)波場(chǎng)快照較清晰，頻散得到較好壓制； 圖3c的4～8階OSCD格式波場(chǎng)快照均很清晰，幾乎看不到數(shù)值頻散。與圖1的理論分析結(jié)果一致，4階OSCD格式(圖3c左)不僅明顯好于4階SCD格式(圖3b左)，也要好于使用了相同節(jié)點(diǎn)數(shù)的6階SCD格式(圖3b中)，達(dá)到了8階SCD格式(圖3b右)的效果。若要SD和SCD格式達(dá)到與OSCD相同的數(shù)值模擬效果，在相同網(wǎng)格步長(zhǎng)條件下，必須增加差分節(jié)點(diǎn)提高差分精度，所以O(shè)SCD格式能使用較少的差分節(jié)點(diǎn)，從而減少了運(yùn)算量，提高了計(jì)算效率。

3.2 水平層狀介質(zhì)模型

設(shè)置四層水平層狀介質(zhì)模型，模型長(zhǎng)度和深度均為2000m，各層厚度均為500m，縱波速度以500m/s間隔，從3000m/s增至4500m/s，密度由Gardner經(jīng)驗(yàn)公式得到。在地表(1000m，0)處激發(fā)主頻為20Hz的雷克子波(震源)，時(shí)間步長(zhǎng)Δt=1ms，記錄時(shí)間長(zhǎng)度為1s。邊界處理采用分裂的PML邊界條件對(duì)邊界反射進(jìn)行吸收衰減，其控制方程見文獻(xiàn)[7]，此處不贅述。

理論分析和均勻介質(zhì)模擬都說(shuō)明了常規(guī)SD格式與SCD和OSCD格式相比，不利于壓制數(shù)值頻散，所以下文中不再對(duì)常規(guī)SD格式進(jìn)行分析，只比較SCD和OSCD的數(shù)值模擬結(jié)果。在都使用4個(gè)網(wǎng)格節(jié)點(diǎn)情況下，也就是利用4階OSCD格式和6階SCD格式進(jìn)行數(shù)值模擬，并采用不同的網(wǎng)格尺寸進(jìn)行計(jì)算，數(shù)值模擬地表接收的單炮記錄(長(zhǎng)度為1s)如圖4所示。所有地震記錄使用同樣的參數(shù)進(jìn)行增益處理，并提取x=1000m處質(zhì)點(diǎn)模擬得到的不含直達(dá)波單道地震記錄(圖5)。

從圖4和圖5可看出，隨著網(wǎng)格尺寸h的增加，數(shù)值頻散總體上逐漸變得嚴(yán)重。圖4顯示6階SCD格式和4階OSCD格式分別最大在15m和17m時(shí)地震記錄清晰，無(wú)明顯數(shù)值頻散和邊界反射； 圖5中單道記錄在波尾處也無(wú)抖動(dòng)，SCD和OSCD分別從16m和18m開始出現(xiàn)“拖尾”現(xiàn)象。僅以該模型的單個(gè)變量而言，15m網(wǎng)格SCD格式的網(wǎng)格數(shù)約為133×133，17m網(wǎng)格OSCD格式的網(wǎng)格數(shù)為117×117，相對(duì)減少內(nèi)存23%，也減少了運(yùn)算時(shí)間，這對(duì)于大尺度模型的數(shù)值模擬或逆時(shí)偏移還是比較有意義的。

圖3 4階(左)、6階(中)、8階(右)三種格式模擬的600ms時(shí)刻波場(chǎng)快照

3.3 Marmousi模型

為了驗(yàn)證OSCD格式對(duì)復(fù)雜介質(zhì)的適用性，用經(jīng)典二維Marmousi縱波速度模型(圖6)做數(shù)值模擬。模型的速度范圍是1729～5500m/s，根據(jù)Gardner公式ρ=0.31v0.25，得到密度范圍是2.00～2.67g/cm3。模型范圍是501×501個(gè)網(wǎng)格點(diǎn)，空間網(wǎng)格尺寸取12m，時(shí)間步長(zhǎng)取0.8ms，縱波震源位于地表中心位置，激發(fā)30Hz的Ricker子波，采樣時(shí)間3s，采用20層PML吸收邊界。

在上節(jié)主要對(duì)比了4階OSCD與6階SCD格式，本節(jié)對(duì)比6階OSCD與8階SCD格式數(shù)值模擬的結(jié)果。一方面驗(yàn)證OSCD格式對(duì)復(fù)雜介質(zhì)的適用性； 另一方面驗(yàn)證本文方法對(duì)其他精度SCD格式的優(yōu)化效果。圖7是分別使用優(yōu)化前8階SCD和優(yōu)化后6階OSCD格式得到的地震記錄，時(shí)間精度均為2階。為了顯示清晰，對(duì)地震記錄進(jìn)行了瞬時(shí)自動(dòng)增益控制(AGC)處理，時(shí)窗為500ms。

圖4 水平層狀介質(zhì)的不同網(wǎng)格尺寸時(shí)6階SCD (a)和4階OSCD (b)格式模擬地面地震記錄

圖5 不同網(wǎng)格尺寸的6階SCD(上)、4階OSCD(下)格式單道地震記錄對(duì)比

從圖7中虛線方框所圈定范圍的地震波場(chǎng)看，優(yōu)化前的8階SCD格式存在較明顯的數(shù)值頻散，不利于波場(chǎng)特征分析或偏移成像處理； 而優(yōu)化差分系數(shù)后的6階OSCD格式數(shù)值頻散減弱，直達(dá)波、反射波及繞射波等波型清晰可見，且邊界吸收效果較好，驗(yàn)證了本文算法對(duì)復(fù)雜模型的適用性。

4 結(jié)論

基于頻散關(guān)系保持的思想，本文利用最小平方法和拉格朗日乘數(shù)法對(duì)一階導(dǎo)數(shù)的緊致交錯(cuò)差分格式的差分系數(shù)進(jìn)行了優(yōu)化，對(duì)比了優(yōu)化前、后差分格式的頻散關(guān)系和模擬精度，以及一階速度—應(yīng)力聲波方程組優(yōu)化前、后緊致交錯(cuò)差分格式的穩(wěn)定性條件，獲得以下認(rèn)識(shí)。

(1)相同的差分精度時(shí)，與常規(guī)SD和SCD格式相比，OSCD格式具有更小的數(shù)值頻散、更小的截?cái)嗾`差、更高的計(jì)算精度。

(2)相同的網(wǎng)格參數(shù)時(shí)，優(yōu)化后4階OSCD格式在壓制數(shù)值頻散方面優(yōu)于6階SCD格式，能達(dá)到8階SCD格式的效果，即可使用更少的計(jì)算節(jié)點(diǎn)(算子長(zhǎng)度)，提高了計(jì)算效率。

(3)相同的計(jì)算節(jié)點(diǎn)時(shí)，優(yōu)化后的OSCD格式雖然比優(yōu)化前的SCD格式精度降低2階，但能更好地壓制數(shù)值頻散，即OSCD格式能使用更粗的空間網(wǎng)格進(jìn)行計(jì)算，因而節(jié)省了內(nèi)存容量，更適用于粗網(wǎng)格下的大尺度的地震波場(chǎng)數(shù)值模擬。

(4)優(yōu)化差分系數(shù)后，聲波方程的穩(wěn)定性條件比優(yōu)化前略微嚴(yán)格，但總體上差別不大。Courant數(shù)隨空間差分精度的增加而減小，隨時(shí)間差分精度的增加而增加。高階的時(shí)間精度會(huì)轉(zhuǎn)換為更高階的空間導(dǎo)數(shù)，勢(shì)必會(huì)增加內(nèi)存和運(yùn)算時(shí)間，所以數(shù)值模擬中，應(yīng)兼顧計(jì)算精度與效率，選取適宜的差分階數(shù)、空間和時(shí)間步長(zhǎng)。

(5)本文從積分方程表示的修正波數(shù)與真波數(shù)之間的誤差出發(fā)，利用拉格朗日乘數(shù)法求取優(yōu)化差分系數(shù)，積分上限是采用多次試驗(yàn)得到的3π/4。后續(xù)研究中，擬從方程最優(yōu)解角度選定適用于不同差分階數(shù)的積分上限，從而得到不同的優(yōu)化差分系數(shù)，并進(jìn)行對(duì)比和分析，找到最佳優(yōu)化差分系數(shù)。