基于緊致交錯差分格式的二維聲波及黏滯聲波方程數(shù)值模擬

鄧光校, 鮑羽, 汪勇*

(1.中國地質(zhì)大學(武漢)地球物理與空間信息學院, 武漢 430074; 2.中國石化碳酸鹽巖縫洞型油藏提高采收率重點實驗室,烏魯木齊 830011; 3.油氣資源與勘探技術教育部重點實驗室(長江大學), 武漢 430100)

地震波場正演模擬是針對地下介質(zhì)運用數(shù)值計算方法去模擬地震波在地下的傳播過程,進而得到位于地面或地下觀測點的地震記錄,是地震勘探以及地震學的重要基礎。基于有限差分數(shù)值模擬的發(fā)展以及實際生產(chǎn)的需求,針對如何去提高有限差分的計算時的效率[1]、模擬精度[2]、算法穩(wěn)定性[3-5]、處理復雜介質(zhì)[6-7]、吸收邊界條件[8-9]和壓制數(shù)值頻散[10-11]等研究方向上,研究人員陸續(xù)提出了多種方法,也做了大量研究,也取得較為重大的研究成果。

差分計算時增加網(wǎng)格點數(shù)是提高數(shù)值,模擬精度最直接的方法,但網(wǎng)格點增加必然導致計算時計算量以及存儲量的增加。然而緊致差分就可以很好地解決這個問題,緊致差分是一種隱式差分格式,穩(wěn)定性也較好,也正是具備這些優(yōu)勢也使它成為目前研究較多的有限差分方法之一。Du等[12]為了降低頻散將交錯網(wǎng)格與緊致有限差分結合,并將其應用于垂直橫向各向同性(vertical transversely isotropic, VTI)介質(zhì)的一階速度應力方程的數(shù)值模擬。Liu等[13-14]證明了在不考慮訪問數(shù)組的額外成本時,高階顯式差分能夠被相同階的隱式差分所替代。并且將其運用到了聲波以及彈性波的數(shù)值模擬進行了對比。楊寬德等[15]把緊致差分運用于包含流體的多孔隙各向異性介質(zhì)中,并且發(fā)現(xiàn)通量校正傳輸(flux-corrected transport,FCT)緊致差分方法能有效壓制粗網(wǎng)格條件下模擬彈性波場所引起的數(shù)值頻散以及源噪聲。周成堯等[16]基于五點八階超緊致差分格式在黏彈性介質(zhì)中實現(xiàn)了聲波方程的正演模擬。緊致差分格式后也逐漸被研究人員依次應用于聲波、彈性波和復雜介質(zhì)等地震波場數(shù)值模擬中。交錯網(wǎng)格差分格式最早是由Madariaga[17]提出的,交錯差分提高了數(shù)模模擬時的局部精度,也使得收斂速度更快。因此將交錯網(wǎng)格技術引入緊致差分,可以進一步提高數(shù)值模擬的精度,更好地壓制數(shù)值頻散。

由于地層的黏彈性所導致的地層對地震波產(chǎn)生的吸收、衰減作用,使得地震波在地下的傳播規(guī)律變得復雜,所以研究地震波在黏彈性介質(zhì)中的傳播規(guī)律對地震勘探有著重要的意義。研究人員也基于不同方法對不同的黏彈介質(zhì)模型進行了數(shù)值模擬研究,如李曉波等[18]基于斑狀飽和介質(zhì)的黏彈特性進行了地震波模擬;羅文山等[19]針對黏滯波動方程推導了分數(shù)階拉普拉斯算子黏滯波動方程的一階速度-壓力形式,在時間域進行了正演模擬。姚振岸等[20-21]在黏彈各向異性介質(zhì)中進行了微地震波場模擬。

關于聲波方程在緊致交錯差分格式下的穩(wěn)定性條件以及該差分格式運用于黏滯聲波方程的數(shù)值模擬尚未有相關方向的文獻發(fā)表。現(xiàn)分析和對比它與傳統(tǒng)的中心差分和交錯差分格式,在模擬精度、數(shù)值頻散和穩(wěn)定性條件3個方面的區(qū)別,并將其應用于對一階速度-應力聲波方程組和黏滯聲波方程組的數(shù)值模擬。

1 二維黏滯聲波方程

在研究介質(zhì)的非完全彈性時,地震勘探中常采用開爾文-佛各特(Kelvin-Voigt)體的模型,認為地震波吸收系數(shù)與地震信號頻率的平方成正比,這也與實際的地質(zhì)情況更加貼近。該介質(zhì)中黏滯聲波方程的一階速度-應力波動方程組為

(1)

2 緊致交錯有限差分方法

Madariaga[17]最早提出波動方程交錯網(wǎng)格數(shù)值模擬方法,差分精度為O(Δt2+Δx2),進而Levander[22]將交錯網(wǎng)格方法應用于求解P-SV波方程,差分精度也被提高到O(Δt2+Δx4),董良國等[4]進一步提出了一階彈性波方程的交錯網(wǎng)格高階差分解法,使得差分精度達到O(Δt2M+Δx2N)。現(xiàn)如今,交錯有限差分方法已經(jīng)發(fā)展成一種高效且被廣泛應用的數(shù)值模擬方法。

假設函數(shù)f(x)連續(xù),則一階導數(shù)的2N階精度交錯網(wǎng)格差分可以表示為

(2)

式(2)中:fi為節(jié)點函數(shù)值;f′i為一階導數(shù)值;Δx為網(wǎng)格大小;cn為差分系數(shù),如表1所示。

表1 一階導數(shù)的差分系數(shù)表

對式(1)中的vx(t)分別在t+Δt/2時刻和t-Δt/2時刻通過泰勒公式展開相減,略去高次項,可得vx(t)的時間二階精度近似式為

(3)

代入黏滯聲波方程[式(1)],將vx對時間的導數(shù)轉換為P對空間的導數(shù),可以得到vx的二階精度兩層顯示差分格式,同理也可以得到vz和P的時間二階精度遞推式為

(4)

利用式(4)就可以進行黏滯聲波方程的地震波場時間層的推進。

由式(2)可以得到,常規(guī)交錯差分格式達到2N階空間差分精度,需要用到2N個節(jié)點的函數(shù)值。如果想要提高差分格式數(shù)值的計算精度,最直接有效的方法就是在差分計算時增加網(wǎng)格點的數(shù)目,但這也使得計算量以及所需的存儲空間增加。緊致差分格式很好地解決了這個問題,早期的緊致差分格式是建立在Hermite多項式構造的基礎之上。

將緊致差分格式與交錯網(wǎng)格技術結合,最早由Nagarajan等[23]提出,并且將其用于大渦模擬(large eddy simulation)問題,緊致交錯差分僅用4個網(wǎng)格節(jié)點得到6階空間精度的模擬效果。隨后,Boersma[24]提出了最高到12階空間精度的緊致交錯差分格式,用于可壓縮流體的Navier-Stokes方程的數(shù)值模擬。

若函數(shù)f(x)是連續(xù)的,則一階導數(shù)的2N(N=2～5)階精度緊致交錯差分格式表示為

(5)

式(5)中:a、bn為差分系數(shù),如表1所示。

從式(2)和式(5)可以看出,二者均是利用半網(wǎng)格點上的變量值計算整網(wǎng)格點上的一階導。區(qū)別在于,常規(guī)交錯差分格式(staggered finite difference,SD)在求取中心點處的導數(shù)值時,僅需要知道周圍網(wǎng)格點的函數(shù)值,而緊致交錯差分格式(staggered compact finite difference,SCD)還需使用相鄰點的導數(shù)值,所以可以認為SD格式是SCD格式的特例,即式(5)中的a=0。由于緊致差分格式是空間隱式差分,若要達到2N階空間差分精度,僅需2N-2個節(jié)點函數(shù)值,能夠節(jié)省計算量和存儲空間。

根據(jù)式(5)所表示的SCD格式可以求得式(4)中的各偏導數(shù)的值。對二維地震波場離散化后假設在x方向有n個節(jié)點,z方向有m個節(jié)點。緊致交錯差分格式在進行一階導數(shù)的求取時,是把變量的值和導數(shù)值分別置于兩套網(wǎng)格中,因此可以把差分格式定義成向前和向后差分兩種類型。以變量P為例,將導數(shù)值固定在整網(wǎng)格上,變量P的值固定在半網(wǎng)格上。然后根據(jù)式(5)求?P/?x和?P/?z的值,則差分格式可以寫成矩陣形式為

(6)

式(6)中:

向前差分時:

向后差分時:

根據(jù)式(6),由應力場P求它的空間一階導數(shù)?P/?x和?P/?z可表示為

(7)

表1給出了矩陣中的差分系數(shù)a和bn(n=1～4),可以求得4～10階的空間差分精度近似值。

緊致交錯差分格式與常規(guī)交錯差分格式定義變量類似,利用相似的方法,將應力P固定于時間整網(wǎng)格上,速度vx和vz固定于時間半網(wǎng)格,可得到時間為二階精度的一階速度-應力黏滯聲波方程的緊致交錯差分格式表示為

(8)

為了說明SD和SCD格式的區(qū)別,這里直接給出一階速度-應力黏滯聲波方程交錯網(wǎng)格時間二階精度差分格式為

(9)

式(9)中:E1和E2分別為x方向的向前和向后差分系數(shù)矩陣;F1和F2分別為z方向的向前和向后差分系數(shù)矩陣。E1和E2的表達式如下,其中的ci(i=1～5)為表1中的差分系數(shù)。

2.1 精度分析

當時間差分精度相同時進行數(shù)值模擬,不論是利用常規(guī)交錯網(wǎng)格差分(SD)格式,還是利用緊致交錯網(wǎng)格差分(SCD)格式,它們在時間層的推進方式是相同的。唯一區(qū)別在于SD和SCD空間偏導數(shù)求取時的差分格式不同,例如SCD格式表示的向前差分為?P/?x=PB1A-1/Δx,SD格式則為?P/?x=PE1/Δx。在進行二者近似精度的比較時,需要獲取它們一階導數(shù)的截斷誤差去進行比較,如表2所示。不難看出,在同樣空間近似精度的條件下,兩種格式的截斷誤差有較大的差別。SD格式計算一階偏導數(shù)4～10階差分精度的截斷誤差約是SCD格式的1.46、3.18、9.71和8.7倍,數(shù)據(jù)表明SCD與SD相比,其截斷誤差更小和差分精度也更高。

表2 一階導數(shù)SD和SCD方法的截斷誤差主項系數(shù)比較

為了比較SD、SCD格式和常規(guī)中心差分格式(finite difference,FD)的數(shù)值模擬精度,以一維平面簡諧波初值問題進行試算。一維平面諧波初值問題數(shù)學表達式為

(10)

上述偏微分方程的達朗貝爾解,即精確的解析解為

(11)

設置一維介質(zhì)模型長度200 m,v=50 m/s,時間步長為0.1 ms,空間網(wǎng)格大小設置為1 m。3種差分格式的空間差分固定為四階精度進行數(shù)值模擬。可得圖1所示的1 s時刻的?u/?x分量波場快照,圖1中黑色的實線部分是套用一維平面簡諧波的達朗貝爾解計算得到的精確解析解。不難看出,緊致交錯差分的數(shù)值模擬結果更加接近精確的解析解,說明其具有更好的模擬精度,常規(guī)交錯網(wǎng)格次之,有限差分格式的精度最差。由定量計算可以得知,該時刻FD、SD和SCD的數(shù)值解與解析解之間的相對誤差依次為24%,8.5%和5.2%,導致這一誤差的結果正是由于它們在計算空間導數(shù)時存在不同截斷誤差大小,這也與表2的理論分析結果相符。需要說明的是,這里為了比較精度差異,選取的網(wǎng)格較大,導致波形方波化,若減小網(wǎng)格大小,能夠使波形足夠光滑,模擬精度也更高。

圖1 3種差分格式的一維聲波方程數(shù)值模擬

2.2 頻散分析

如果在進行數(shù)值模擬的時候采用的空間網(wǎng)格過大,就會得到較大的求解誤差,產(chǎn)生數(shù)值頻散,因此在進行判別數(shù)值模擬方法的優(yōu)劣時,頻散分析是一種重要的分析手段,也是進行數(shù)值模擬確定網(wǎng)格大小的重要依據(jù)。利用一階導數(shù)的緊致交錯差分格式[式(5)]進行頻散分析,依據(jù)數(shù)值波數(shù)與真波數(shù)的比值,分析該方法的適用條件。令

(12)

將式(12)代入式(5)中可得

(13)

利用歐拉公式可得

(14)

將B=(Ik′)A代入式(14)可得

(15)

假設φ=kh,φ′=k′h,那么修正波數(shù)與真波數(shù)之比就可以定義為

(16)

在理想情況下,如果不存在數(shù)值頻散則波數(shù)比R恒等于1。R偏離1越大,則說明該方法的數(shù)值頻散越嚴重,反之則說明該方法能更好地壓制數(shù)值頻散。同樣的方法可以得到一階導數(shù)的常規(guī)交錯差分(SD)和中心差分(FD)格式的波數(shù)比,圖2為3種格式在不同差分精度情況下的波數(shù)比曲線。

圖2 3種差分格式的波數(shù)比曲線

將φ∈[0,π] 作為坐標系的橫軸,可由波數(shù)與空間步長相乘得到,單位波長中采樣點個數(shù)N=2π/φ,因此橫坐標亦可以看作N由∞逐漸減小至2。根據(jù)圖2中的速度比曲線不難得出,交錯差分以及緊致交錯差分的數(shù)值波數(shù)相較于常規(guī)的中心差分更接近于真波數(shù),直接表明這兩種方法在壓制數(shù)值頻散上的效果更好;緊致交錯差分格式與常規(guī)交錯差分格式相比,壓制數(shù)值頻散的特征更加明顯,例如圖2中6階SCD曲線接近10階SD,所以SCD方法能使用更少的網(wǎng)格點,可以節(jié)省計算資源;3種方法均在10階精度的前提下,為保證無數(shù)值頻散,要求波數(shù)比大于99.99%,則FD、SD和SCD方法每個波長內(nèi)需要采樣11.6、5.2和4.2個點。所以緊致交錯差分更能適用于粗網(wǎng)格時的數(shù)值模擬,能夠更好地壓制數(shù)值頻散。

前人研究表明FD格式不利于壓制數(shù)值頻散,這里只比較SD和SCD兩種差分方法的差異。模型大小設置為4 000 m×4 000 m,縱波速度給定為4 000 m/s,在模型的中心加載雷克子波震源,峰值頻率25 Hz。時間步長2 ms,空間網(wǎng)格20 m,再針對一階速度-應力聲波方程采用SD和SCD格式在不同空間差分精度進行波場模擬,進而得到同一時刻的應力P分量的波場快照(若未加說明,波場快照和地震記錄均為應力P分量),如圖3所示。

圖3 SD和SCD不同差分精度計算的400 ms時刻波場快照

從圖3可以看出,同前面的理論分析結果一致,SCD格式在壓制數(shù)值頻散方面比SD格式更具有優(yōu)勢,圖3(d)為6階SCD格式,其波場快照與圖3(c)所示的10階SD格式接近,要遠遠好于圖3(a)所示的6階SD格式。為了更清楚地觀察差異,提取模型x=2 000 m,z=2 000 m處質(zhì)點的振動曲線,見圖4。可以看出,6階SD的振動曲線存在明顯的尾波,即數(shù)值頻散,而10階SD、6階和10階SCD的振動曲線則基本重合,波形平滑,無尾波現(xiàn)象,這表明SCD格式能用較低的差分階數(shù),能夠使用粗網(wǎng)格計算,提高了運算效率。

圖4 SD和SCD不同差分精度計算的波場記錄

2.3 穩(wěn)定性分析

有限差分數(shù)值模擬中另外一個非常重要的問題就是穩(wěn)定性條件,也是影響差分方法計算效率的重要因素之一。這里參照董良國等[4]和杜啟振等[5]使用的Fourier方法對差分格式[式(8)]進行穩(wěn)定性分析,省略掉推導過程,直接利用式(17)表示的二維聲波方程(2M,2N)階差分精度的緊致交錯網(wǎng)格差分格式的穩(wěn)定性條件。

(17)

定義α=vΔt/h為差分格式的Courant數(shù),并與常規(guī)交錯網(wǎng)格差分格式的穩(wěn)定性條件進行對比,對比結果如表3所示。從表3來看,差分格式的時間差分精度越高,差分格式也越穩(wěn)定,但是更高階的時間精度在計算時需要轉變?yōu)楦唠A的空間導數(shù),會相應地增加計算量。當差分精度相同時,SCD格式的穩(wěn)定性條件比SD格式要略微嚴格,也就是說空間步長相同時,允許采納的時間步長要略小。

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

2.4 PML邊界條件

數(shù)值模擬過程模型大小的限制在計算時的設置都會導致人工邊界的存在,倘若不對其進行處理,就會產(chǎn)生邊界反射,從而對正常的地震波場造成干擾,因此處理人工邊界的好壞也直接影響到了數(shù)值模擬的精度和效率。對于邊界的處理,采用完全匹配層(perfectly matched layer,PML),按照王守東[25]和劉有山等[26]推導聲波PML控制方程的思路和方法,略去推導過程,這里直接給出一階速度-應力黏滯聲波方程的PML邊界條件的控制方程組,當η=0時,即可得到聲波方程的PML控制方程,即

(18)

式(18)中:P1和P2為引入的中間變量;d(x)和d(z)分別為x和z方向的衰減系數(shù),在這兩個方向上起到衰減的作用。高剛等[27]的研究指出,d(x)和d(z)的衰減函數(shù)類型和層數(shù)會影響邊界反射的衰減效果,認為余弦型吸收衰減函數(shù)效果較好,所以在模型試算中采用該衰減函數(shù),表示為

(19)

式(19)中:αi為邊界層的吸收衰減因子;B為衰減幅度因子,計算時取500;PML為整個吸收邊界層的厚度,在進行實際的數(shù)值模擬時,可令其為20。為了驗證所使用差分格式在PML條件下對邊界反射的吸收效果,在2.2節(jié)中的模型之上采用SCD(2,10)(時間2階,空間10階精度)格式進行數(shù)值模擬,空間網(wǎng)格10 m,時間步長1 ms,對比未采用PML邊界條件和采用PML邊界條件時在600 ms時刻的波場快照,如圖5(a)和圖5(b)所示,并將二者做差,得到圖5(c)只包含邊界反射的波場快照。不難看出,SCD格式在PML控制方程作用下,人工邊界反射得到了有效吸收,并且有效波沒有發(fā)生明顯改變,處理結果十分可靠。

圖5 600 ms時刻波場快照

3 模型試算

第2節(jié)通過理論分析和數(shù)值模擬討論了二維聲波方程的緊致交錯差分格式的計算精度、頻散關系、穩(wěn)定性和邊界條件,然后基于SCD(2,10)差分格式對式(1)表示的二維黏滯聲波方程進行試算,以驗證該方法是否適用于黏滯聲波方程的數(shù)值模擬。

3.1 均勻介質(zhì)模型

圖6 350 ms時刻波場快照

提取數(shù)值模擬時得到的波場快照中x=2 000 m的波形記錄,如圖7所示。不難看出,波形曲線平滑,并且不存在數(shù)值頻散,且隨著品質(zhì)因子的減小,能看到地震波振幅發(fā)生顯著衰減,頻率降低,波形展寬,相應的波長也增加了。

圖7 不同品質(zhì)因子時的波形曲線

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

設置4層層狀介質(zhì)模型,利用SCD(2,10)格式進行數(shù)值模擬,模型參數(shù)如表4所示。

表4 多層層狀介質(zhì)模型參數(shù)

模型深度和長度均設置成2 000 m,時間步長Δt=0.5 ms,空間網(wǎng)格大小h=10 m,在地表(x=1 000 m,z=0 m)位置加載主頻30 Hz雷克子波震源,由地表接收的單炮記錄(記錄長度1.5 s)如圖8所示,為了方便比較顯示,對地震記錄進行了增益處理,提高深層反射波振幅強度。

圖8 地面地震記錄

利用四層水平層狀介質(zhì)模型對地面接收到的地震波場進行模擬,發(fā)現(xiàn)地震記錄清晰可見,不存在數(shù)值頻散和邊界反射情況,反射波顯示清晰,充分說明本文算法可以有效地模擬黏滯聲波在多層介質(zhì)模型中的傳播。

圖9給出檢波點(x=1 500 m,z=0 m)接收到的振動曲線,可以看出,由于介質(zhì)的吸收衰減作用,黏滯聲波中的反射波振幅比聲波模擬要小,且隨著傳播時間和距離的增加,反射波信號變寬,頻率變低,這也更加符合實際地震信號特征。

圖9 檢波器(x=800 m,z=0 m)接收到的波形曲線

3.3 Marmousi模型

為了檢驗緊致交錯差分格式是否適用于復雜的介質(zhì)的數(shù)值模擬,采用經(jīng)典的二維Marmousi縱波速度模型進行模型試算。模型如圖10所示,當速度v的范圍是1 729～5 500 m/s,品質(zhì)因子Q通過經(jīng)驗公式Q=14v2.2計算得到,其中v的單位是km/s,可以算出品質(zhì)因子范圍為46.7～595.6。模型大小設置為501×501個網(wǎng)格點,空間網(wǎng)格大小固定為8 m,時間步長0.5 ms,縱波震源位于(x=2 000 m,z=0 m),加載30 Hz的Ricker子波,采樣時長3 s。

圖10 Marmousi模型

用Marmousi模型對黏滯聲波方程和聲波方程進行數(shù)值模擬得到了如圖11所示的波場快照,可以看到,波場特征清晰,無明顯數(shù)值頻散以及邊界反射,邊界的吸收效果良好。闡明了本文方法對復雜介質(zhì)的適用性。圖12分別是由聲波方程以及黏滯聲波方程得到的地震記錄,由于波前擴散導致深部的反射振幅明顯減弱,因此對地震記錄進行了自動增益控制處理,加強了深層的能量使顯示更為清晰,時窗設置為500 ms,通過對比可以發(fā)現(xiàn):黏滯聲波的數(shù)值模擬結果除了波前擴散造成的振幅衰減,地震波傳播過程中還受到了介質(zhì)的吸收衰減作用,從而反射波的振幅衰減幅度要大于聲波方程中的模擬結果,而且地震波的頻率也會隨之變化,頻率會隨著地震波的傳播距離的增加而減小。

圖11 Marmousi模型波場快照及黏滯聲波模擬

圖12 Marmousi 模型的地面地震記錄(z=0 m)

對于進行AGC處理后的地面地震記錄,提取單道的地震記錄波形圖(x=1 000 m)如圖13所示。可以看到,黏滯聲波記錄的地震波振幅明顯小于聲波,子波的延續(xù)時間大于聲波,反射波個數(shù)也更少,并且不難發(fā)現(xiàn),傳播的時間越長這種差異越明顯,進而導致了黏滯聲波模擬結果從垂向分辨率上小于聲波。

圖13 檢波器(x=1 000 m)接收的地震記錄

對模擬記錄進行廣義S變化,再進行時頻分析,得到的結果如圖14所示,可以看到兩道模擬記錄在頻譜上的差別。從結果來看,聲波的模擬記錄[圖14(a)]的地震波振幅隨著時間的增加而減小,導致這種變化的原因是波前擴散。但是振幅雖然衰減了,反射波的主頻卻未發(fā)生明顯改變,幾乎與激發(fā)主頻平直,也是30 Hz左右。而黏滯聲波記錄[圖14(b)]除了振幅衰減,反射波主頻也會降低到10 Hz左右。通過對比發(fā)現(xiàn),兩個記錄的時頻譜,聲波的數(shù)值模擬結果其能量的峰值相對集中也更加清晰,與時間域的主要反射波對應;而黏滯聲波記錄中,1.4 s和2.5 s位置的能量峰值相對模糊一些,這也從頻域說明了地震垂向分辨率的降低。

圖14 x=1 000 m處單道地震記錄時頻分析結果

4 結論與建議

利用聲波方程結合緊致交錯差分格式,在模擬精度、頻散曲線和穩(wěn)定性分析3個方面,研究認為SCD格式僅需使用2N-2個節(jié)點,就能使一階導數(shù)的差分精度達到2N階精度,少于FD格式(2N+1個)和SD格式(2N個);與SD格式相比,SCD格式具有更高的模擬精度和略微嚴格的穩(wěn)定性條件;此外,6階SCD格式的頻散關系幾乎可以等同于10階SD格式,這也說明相同精度下的SCD格式其數(shù)值頻散更小,更加適用于粗網(wǎng)格下的大尺度的地震波場數(shù)值模擬。

其次,基于黏滯聲波的PML控制方程,在均勻介質(zhì)模型、水平層狀介質(zhì)模型和Marmousi模型上利用SCD格式,進行了數(shù)值模擬。最終可見地震波場特征清晰,表明本文方法對于復雜介質(zhì)的適用性,也進一步說明該差分格式不僅可以用于黏滯聲波的波場模擬,并且還能夠推廣到二維或三維的各向異性介質(zhì)以及雙相介質(zhì)等復雜介質(zhì)的聲波或彈性波數(shù)值模擬中,同時為研究地震波傳播規(guī)律、逆時偏移以及全波形反演等工作提供了一種行之有效的方法。

為減少計算內(nèi)存和時間,所用的差分格式時間為二階,在數(shù)值模擬時可以通過提高時間差分的階數(shù),從而達到增加數(shù)值模擬計算精度的目的。此外,也可以通過優(yōu)化差分系數(shù)的方法進一步提升該方法的優(yōu)勢。