胡 寧，劉 財(cái)

吉林大學(xué)地球探測科學(xué)與技術(shù)學(xué)院，長春 130026

0 引言

近年來，分?jǐn)?shù)階微積分被廣泛應(yīng)用于黏彈阻尼器、地震分析、反常擴(kuò)散、信號處理與控制、流體力學(xué)、圖像處理、軟物質(zhì)研究等多個(gè)領(lǐng)域。相對于整數(shù)階導(dǎo)數(shù)，分?jǐn)?shù)階微分算子可以更簡潔地描述具有歷史依賴性和空間全域相關(guān)性的復(fù)雜力學(xué)和物理過程。隨著油氣勘測對象的日趨深入與巖石模型的日趨復(fù)雜，更多學(xué)者考慮到，實(shí)際的地下介質(zhì)并不是完全彈性的，地震波在地下傳播過程中會(huì)發(fā)生能量衰減、相位畸變和頻散。引起這3種波場變化的機(jī)制為：孔隙及流體、基質(zhì)黏彈性和各向異性。為了更準(zhǔn)確地描述波在實(shí)際介質(zhì)中的傳播規(guī)律，需要考慮以上3種機(jī)制。對于雙相孔隙介質(zhì)：與Biot模型[1]相比，BISQ(Biot-squirt)模型[2]在考慮了Biot流動(dòng)的同時(shí)還考慮了噴射流動(dòng)，但其特征噴射流動(dòng)長度為虛數(shù)，在時(shí)間域處理不太合理，并且物理意義并不十分明確，本文采用不依賴于特征噴射長度的改進(jìn)的BISQ模型[3]來描述孔隙介質(zhì)；對于孔隙流體黏滯性：孔隙介質(zhì)中固、流兩相除相對平動(dòng)外，還存在相對轉(zhuǎn)動(dòng)，流體的振動(dòng)速度梯度會(huì)產(chǎn)生剪切應(yīng)變，導(dǎo)致流相中存在剪切應(yīng)力和附加法相應(yīng)力[4]，應(yīng)該考慮黏性摩擦力產(chǎn)生的應(yīng)力偏量；對于固體骨架黏彈性：含分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)的常Q(品質(zhì)因子)模型恰好可以準(zhǔn)確地描述固體骨架的黏彈性。

Kjartansson[5]提出的常Q模型認(rèn)為：衰減因子與頻率成線性反比關(guān)系，在常規(guī)地震頻帶內(nèi)，Q值與頻率基本無關(guān)?；贙jartansson常Q模型的波傳播方程是含有分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)的，Carcione[6]采用分?jǐn)?shù)階導(dǎo)數(shù)的Grünwald-Letnikov定義，用全局記憶法進(jìn)行數(shù)值模擬。因?yàn)榉謹(jǐn)?shù)階導(dǎo)數(shù)具有全局記憶特性，所以全局記憶法求解波傳播方程需要保存每一個(gè)時(shí)刻的應(yīng)力-應(yīng)變值，并且計(jì)算某個(gè)量某個(gè)時(shí)刻的時(shí)間導(dǎo)數(shù)需要之前所有時(shí)刻的所有值。即，當(dāng)前狀態(tài)與過去所有狀態(tài)有關(guān)，因而所需存儲(chǔ)量和計(jì)算量均很大，計(jì)算效率較低。Podlubny[7]提出“短時(shí)記憶法”對分?jǐn)?shù)階微分算子進(jìn)行截?cái)?，僅考慮當(dāng)前時(shí)刻t以前的有限時(shí)間區(qū)間[t-L,t]，隨著短時(shí)記憶長度L的選取不同，誤差也不同。過去不同時(shí)間點(diǎn)對當(dāng)前時(shí)間點(diǎn)的影響是不同的，接近當(dāng)前時(shí)刻時(shí)間點(diǎn)對當(dāng)前時(shí)刻的影響較大，而距離當(dāng)前時(shí)刻時(shí)間間隔較長的時(shí)間點(diǎn)對當(dāng)前時(shí)刻的影響較小。利用這一性質(zhì)，Brian等[8]于2010年提出“自適應(yīng)記憶法”。

分?jǐn)?shù)階導(dǎo)數(shù)實(shí)際上是整數(shù)階導(dǎo)數(shù)的推廣，關(guān)于分?jǐn)?shù)階導(dǎo)數(shù)的定義，許多數(shù)學(xué)家從不同角度分別給出了不同定義，其合理性和科學(xué)性在實(shí)踐中已經(jīng)得到了驗(yàn)證，這個(gè)數(shù)學(xué)分支已在實(shí)際問題中得到廣泛應(yīng)用。比較常見的3種分?jǐn)?shù)階導(dǎo)數(shù)定義分別為：Grünwald-Letnikov分?jǐn)?shù)階導(dǎo)數(shù)定義[9]、Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)定義[9]和Caputo分?jǐn)?shù)階導(dǎo)數(shù)定義[9]。Riemann-Liouville分?jǐn)?shù)階導(dǎo)數(shù)定義是Grünwald-Letnikov分?jǐn)?shù)階導(dǎo)數(shù)定義的擴(kuò)充，Caputo分?jǐn)?shù)階導(dǎo)數(shù)也是對Grünwald-Letnikov分?jǐn)?shù)階導(dǎo)數(shù)定義的另一種改進(jìn)[10]，在階數(shù)為負(fù)實(shí)數(shù)和正整數(shù)時(shí)三者是等價(jià)的。在實(shí)際應(yīng)用中，用哪種形式的定義需依照具體情況而定。

本文基于Kjartansson常Q應(yīng)力-應(yīng)變關(guān)系，將含有分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)的黏彈固體骨架各向異性本構(gòu)關(guān)系與雙相介質(zhì)理論改進(jìn)BISQ模型有機(jī)地結(jié)合起來，引入流變學(xué)本構(gòu)關(guān)系描述孔隙流體的黏滯性力學(xué)行為[11]，建立含黏滯流體黏彈雙相VTI(橫向各向同性)介質(zhì)中的分?jǐn)?shù)階波傳播方程；并在時(shí)間域進(jìn)行數(shù)值模擬，采用全局記憶法、短時(shí)記憶法、自適應(yīng)記憶法對分?jǐn)?shù)階微分算子進(jìn)行數(shù)值計(jì)算，整數(shù)階時(shí)間、空間微分算子均采用高階交錯(cuò)網(wǎng)格有限差分算法；最后對3種方法各自的特點(diǎn)進(jìn)行歸納總結(jié)，并對3種方法進(jìn)行對比分析。

1 原理

由于Grünwald-Letnikov分?jǐn)?shù)階導(dǎo)數(shù)定義是基于Grünwald-Letnikov有限差分近似式導(dǎo)出的，而本文波傳播方程中整數(shù)階導(dǎo)數(shù)數(shù)值模擬采用高階交錯(cuò)網(wǎng)格有限差分算法進(jìn)行，故相應(yīng)地本文中的分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)計(jì)算采用Grünwald-Letnikov分?jǐn)?shù)階導(dǎo)數(shù)定義。

1.1 Grünwald-Letnikov分?jǐn)?shù)階導(dǎo)數(shù)定義

基于Grünwald-Letnikov有限差分近似式，將

f(x-mh)

(1)

稱為f(x)的γ階Grünwald-Letnikov分?jǐn)?shù)階導(dǎo)數(shù)。式中：h為步長；m為節(jié)點(diǎn)數(shù)；Γ為歐拉伽馬函數(shù)；α∈R，α≠-N1，N1為除1之外的自然數(shù)集。從式(1)可以看出，求得當(dāng)前時(shí)刻分?jǐn)?shù)階導(dǎo)數(shù)需要知道從m=0到m=x/h的所有函數(shù)值，這正是分?jǐn)?shù)階導(dǎo)數(shù)的歷史依賴性和空間全域相關(guān)特性。

1.2 分?jǐn)?shù)階波傳播方程

本文采用基于Kjartansson常Q理論的分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)黏彈本構(gòu)關(guān)系來描述雙相介質(zhì)固體骨架的黏彈性效應(yīng)，并引入流變學(xué)本構(gòu)關(guān)系表征孔隙流體的黏滯力學(xué)行為，根據(jù)Yang等[12]提出的包含BISQ機(jī)制和固流耦合各向異性效應(yīng)的孔隙彈性波動(dòng)理論，給出基于分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)常Q黏彈本構(gòu)關(guān)系的含黏滯流體雙相VTI介質(zhì)二維波傳播方程。

固體骨架黏彈本構(gòu)方程為

σ=c0(t)*?te(t)-α(t)p。

(2)

式中：σ為總應(yīng)力張量；c0(t)為干燥固體骨架弛豫函數(shù)剛度矩陣；*為卷積運(yùn)算符；e(t)為固體骨架應(yīng)變張量；α(t)3×3為有效應(yīng)力之孔隙彈性系數(shù)張量；p為孔隙黏滯流體壓力[13]，可表示為

(3)

(4)

式中：η為流體黏滯系數(shù)；φ為孔隙度；vfi(vfj) (i，j=1,3)為流相質(zhì)點(diǎn)在正交方向x,z上的速度分量；xi(xj)(i，j=1,3)表示x,z軸；us=(us1，us3)T與uf=(uf1，uf3)T分別為固相和流相位移矢量；β為雙相介質(zhì)壓縮率系數(shù)；α11(t)、α33(t)為α(t)3×3的分量；t為時(shí)間；下標(biāo)i,j=(1,3)分別表示二維平面x,z方向分量，后續(xù)公式中i,j意義同此。

幾何方程為

(5)

式中，eij為固體骨架應(yīng)變張量分量。

運(yùn)動(dòng)方程為

(6)

(7)

式中：σij為總應(yīng)力張量分量；ρ1=(1-φ)ρs；ρ2=φρf；ρs為固體密度；ρf為流體密度；ρ22i=ρ2-ρ12i，ρ12i=-ρa(bǔ)i，ρa(bǔ)i為xi方向的固-流耦合附加質(zhì)量密度；rij為流體阻抗系數(shù)，rij=(kij)-1，kij為滲透率張量k的元素，對于VTI介質(zhì)，k=diag(k11,k11,k33)。

綜合式(3)—(7)即可導(dǎo)出以固相和流相位移usi和ufi為基本未知量的、基于分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)常Q黏彈本構(gòu)關(guān)系的含黏滯流體雙相VTI介質(zhì)波傳播方程。在二維x-z平面內(nèi)，這組方程的速度-應(yīng)力形式是由以下8個(gè)方程構(gòu)成的含卷積運(yùn)算的偏微分方程組：

(8)

其中，

2 穩(wěn)定性條件和吸收邊界條件

在進(jìn)行顯式有限差分?jǐn)?shù)值模擬時(shí)，需要考慮計(jì)算過程的穩(wěn)定性。分?jǐn)?shù)階計(jì)算的穩(wěn)定性較為復(fù)雜，至今缺乏系統(tǒng)的分析研究。本文采用與二維一階速度-應(yīng)力方程時(shí)間二階、空間2N階精度交錯(cuò)網(wǎng)格有限差分一致的穩(wěn)定性條件，其近似表達(dá)式為[15]

(9)

式中：Δx和Δz分別為橫向和縱向空間網(wǎng)格步長；Δt為時(shí)間間隔；v=max{vP//,vP⊥,vS}，vP//、vP⊥、vS分別為平行于地層的P波速度分量和垂直于地層的P波速度分量、S波速度分量；Cn為交錯(cuò)網(wǎng)格差分算子系數(shù)。

為消除人工邊界反射，采用簡便實(shí)用的Cerjan衰減邊界[16]，即沿人工邊界向外擴(kuò)充H個(gè)網(wǎng)格點(diǎn)，進(jìn)入擴(kuò)充區(qū)域內(nèi)的波場通過與因子G相乘逐漸衰減為零。

G=exp[-a(H-i)2]，1≤i≤H。

(10)

式中，a為衰減系數(shù)，可通過多次試驗(yàn)確定其最佳值。

3 三種算法及其差分格式

對于分?jǐn)?shù)階波傳播方程中的整數(shù)階時(shí)間導(dǎo)數(shù)和空間導(dǎo)數(shù)，采用高階交錯(cuò)網(wǎng)格有限差分算法進(jìn)行數(shù)值計(jì)算。由于交錯(cuò)網(wǎng)格有限差分方法將空間變量的導(dǎo)數(shù)在響應(yīng)網(wǎng)格點(diǎn)之間的半程上進(jìn)行計(jì)算，可以在不增加計(jì)算量和內(nèi)存要求的前提下提高模擬精度，有效減少傳統(tǒng)有限差分的網(wǎng)格頻散。

對于分?jǐn)?shù)階波傳播方程中的分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)，使用Grünwald-Letnikov分?jǐn)?shù)階導(dǎo)數(shù)定義，并分別使用全局記憶法、短時(shí)記憶法、自適應(yīng)記憶法進(jìn)行數(shù)值計(jì)算。

采用Grünwald-Letnikov分?jǐn)?shù)階導(dǎo)數(shù)定義，分?jǐn)?shù)階導(dǎo)數(shù)的離散式為

(11)

式中，定義函數(shù)

這樣，函數(shù)τ(2γ,m)有如下遞推關(guān)系：

特殊地，當(dāng)m=0時(shí)，τ(2γ,0)=1。

當(dāng)Δt→0、t=mΔt時(shí)，離散化式(11)有

(12)

式中：Δt為時(shí)間步長；(j,l)為空間離散點(diǎn)坐標(biāo)；k為當(dāng)前時(shí)間迭代次數(shù)。

3.1 全局記憶法

由于分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)具有歷史依賴性和空間全域相關(guān)特性，所以計(jì)算當(dāng)前時(shí)刻的分?jǐn)?shù)階導(dǎo)數(shù)需要用到之前所有時(shí)刻的函數(shù)值。在式(12)中不對m做任何截?cái)嗵幚?，即為全局記憶法，其分?jǐn)?shù)階導(dǎo)數(shù)差分格式即為式(12)。

3.2 短時(shí)記憶法

短時(shí)記憶法通過截?cái)嗨阕娱L度，僅考慮當(dāng)前時(shí)刻t以前且對當(dāng)前時(shí)刻影響較大的有限時(shí)間區(qū)間[t-L,t]。對于本文分?jǐn)?shù)階波傳播方程中的分?jǐn)?shù)階導(dǎo)數(shù)，其短時(shí)記憶法差分近似式為

(13)

3.3 自適應(yīng)記憶法

對于任意正實(shí)數(shù)q，設(shè)初始間隔[0,q]，這一區(qū)間內(nèi)的所有點(diǎn)都參與求和計(jì)算，過去時(shí)間歷程被分為不同長度的區(qū)間。第i區(qū)間定義為

I=[qi-1+i,qi]，i∈N+，i≠1。

(14)

對于式(14)，定義區(qū)間集合為

ζ={I=[qi-1+i,qi]:i∈N+，i≠1，i≤imax}。

ζ中每個(gè)區(qū)間內(nèi)的采樣間隔d的集合為

D={d=2i-1:i∈N+,i≠1}。

式中，imax為使得k∈I=[qi-1+i,qi]的最大i值。

根據(jù)上述定義，本文分?jǐn)?shù)階波傳播方程中分?jǐn)?shù)階導(dǎo)數(shù)應(yīng)用自適應(yīng)記憶法的差分近似式為

(15)

式中：g∈N+；mi為時(shí)間點(diǎn)集

M={mi=qi-1+(2i-1)ξ-i+1:

ξ∈N1，mi≤mmax}

中的元素，mmax為當(dāng)i=imax時(shí)的mi。

4 數(shù)值模擬

為了對3種方法進(jìn)行對比分析，選取一個(gè)模型進(jìn)行數(shù)值模擬。模型大小為3 000 m×3 000 m。介質(zhì)物性參數(shù)：干燥情況下各向同性背景孔隙介質(zhì)中vP//=5 600 m/s，vP⊥=5 167 m/s，vS=3 126 m/s，φ=0.2，ρs=2 500 kg/m3，ρf=1 040 kg/m3，ρa(bǔ)1=420 kg/m3，ρa(bǔ)3=450 kg/m3，固相體積模量Ks=80 GPa，x方向和z方向滲透率分別為k11=6×10-10m2/s和k33=1×10-9m2/s，η=0.001 Pa·s，P波和S波的品質(zhì)因子分別為QP=30和QS=15，參考頻率ωr=80 Hz。其他計(jì)算參數(shù)：Δx=Δz=10 m，Δt=5×10-4s，t=0.2 s。震源子波采用Ricker子波，震源位置為(1 500 m, 1 500 m)，主頻為40 Hz，在模擬中震源分別為x、z方向集中力源。在網(wǎng)格點(diǎn)坐標(biāo)為(150, 200)處設(shè)置檢波點(diǎn)記錄固相質(zhì)點(diǎn)垂直速度分量的時(shí)間記錄。

4.1 全局記憶法

由于全局記憶法未對分?jǐn)?shù)階微分算子進(jìn)行截?cái)嗵幚恚跀?shù)值模擬上認(rèn)為是精確的。利用全局記憶法對上述模型進(jìn)行數(shù)值模擬的運(yùn)行時(shí)間和所占用物理內(nèi)存分別為1 395.965 962 s和5 875 MB。圖1給出了利用全局記憶法對波傳播方程進(jìn)行數(shù)值模擬，固相x分量速度和流相z分量速度在t=0.15 s時(shí)刻的波場快照。

4.2 短時(shí)記憶法不同記憶長度比較

為了比較分析不同短時(shí)記憶長度對模擬結(jié)果的影響，針對上述模型，分別采用L=0.15、0.10、0.05 s三種不同短時(shí)記憶長度進(jìn)行數(shù)值模擬，表1給出了所用時(shí)間和內(nèi)存。從表1可以看出，隨著短時(shí)記憶長度的減小，計(jì)算所需時(shí)間和物理內(nèi)存都在減小；這是因?yàn)椴捎玫亩虝r(shí)記憶長度越小，所計(jì)算的當(dāng)前時(shí)刻之前的時(shí)間節(jié)點(diǎn)越少。

表1采用3種不同短時(shí)記憶長度模擬的運(yùn)行時(shí)間與所占物理內(nèi)存

Table1Runningtimeandoccupiedphysicalmemoryofthreedifferentshort-termmemorylengthsimulations

L/s運(yùn)行時(shí)間/s物理內(nèi)存/MB0.151250.06432450110.101032.77132049960.05693.6316934878

圖2給出了利用3種不同短時(shí)記憶長度對波傳播方程進(jìn)行數(shù)值模擬，固相x分量速度在t=0.15 s時(shí)刻的波場快照，可以看出，隨著短時(shí)記憶長度的減小，在慢縱波附近(如箭頭所指)波場存在明顯差異。為了更明顯地顯示差異，圖3給出了網(wǎng)格點(diǎn)坐標(biāo)為(150, 200)處不同計(jì)算方法固相質(zhì)點(diǎn)垂直速度分量的時(shí)間記錄。從圖3可以看出：隨著短時(shí)記憶長度的改變，快準(zhǔn)P波幾乎不受影響，而準(zhǔn)SV波所受影響較大(圖3b)，即：隨著采用的短時(shí)記憶長度減小，

圖1 全局記憶法黏滯相界固相速度x分量(a)和流相速度z分量(b)t=0.15 s時(shí)波場快照Fig.1 Vicious phase boundary wavefield snapshots of x-component velocity in solid phase (a) andz-component velocity fluid phase (b) at t=0.15 s by using global memory method

a.L=0.15 s；b.L=0.10 s；c.L=0.05 s。圖2 三種不同短時(shí)記憶長度黏滯相界固相速度x分量t=0.15 s時(shí)波場快照Fig.2 Vicious phase boundary wavefield snapshots of x-component velocity in solid phase at t=0.15 s by using three different short-term memory lengths

圖3 三種不同短時(shí)記憶長度固相質(zhì)點(diǎn)垂直速度分量時(shí)間記錄對比圖Fig.3 Comparison of the vertical velocity components of particles with three different short-term memory lengths in solid phase

準(zhǔn)SV波的精度越來越差。

4.3 自適應(yīng)記憶法不同初始間隔比較

為了比較分析不同初始間隔對模擬結(jié)果的影響，針對上述模型，分別采用q=15、10、5三種不同初始間隔進(jìn)行數(shù)值模擬，表2給出了所用時(shí)間和內(nèi)存。從表2可以看出，隨著初始間隔減小，計(jì)算所需時(shí)間和物理內(nèi)存都在增大；這是因?yàn)槌跏奸g隔越小，之前的時(shí)間就會(huì)被分成越多的區(qū)間，就會(huì)有更多的時(shí)間節(jié)點(diǎn)參與計(jì)算。

表2采用3種不同初始間隔運(yùn)行時(shí)間與所占物理內(nèi)存結(jié)果

Table2Runningtimeandoccupiedphysicalmemoryofthreedifferentinitialintervalsimulations

q運(yùn)行時(shí)間/s物理內(nèi)存/MB151243.4999644583101357.584104508351461.3141835289

圖4給出了利用3種初始間隔對波傳播方程進(jìn)行數(shù)值模擬，固相x分量速度在t=0.15 s時(shí)刻的波場快照，可以看出，基于自適應(yīng)記憶法、采用3種不同初始間隔得到的波場快照無明顯差異。為了更明顯地顯示差異，圖5給出了網(wǎng)格點(diǎn)坐標(biāo)為(150, 200)處不同計(jì)算方法固相質(zhì)點(diǎn)垂直速度分量的時(shí)間記錄。從圖4及圖5均可以看出：隨著初始間隔的改變，快準(zhǔn)P波與準(zhǔn)SV波所受影響均較小。因此，當(dāng)選取初始間隔在穩(wěn)定性范圍內(nèi)差異相差不懸殊時(shí)，自適應(yīng)記憶法在計(jì)算精度上是比較穩(wěn)定的。

4.4 三種不同計(jì)算方法對比分析

為了比較分析不同分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)計(jì)算方法對模擬結(jié)果的影響，針對上述模型，圖6給出了利用全局記憶法、L=0.1 s短時(shí)記憶法和q=10自適應(yīng)記憶法對波傳播方程進(jìn)行數(shù)值模擬，網(wǎng)格點(diǎn)坐標(biāo)為(150, 200)處固相質(zhì)點(diǎn)垂直速度分量的時(shí)間記錄。從圖6可以看出，當(dāng)L=0.1 s、q=10時(shí)，從計(jì)算精度上來看，自適應(yīng)記憶法精度較高，而短時(shí)記憶法精度較差。

綜上：全局記憶法計(jì)算當(dāng)前時(shí)刻的分?jǐn)?shù)階導(dǎo)數(shù)需要用到之前所有時(shí)刻的值，因此需要較大的計(jì)算內(nèi)存和計(jì)算時(shí)間，計(jì)算效率較低；短時(shí)記憶法計(jì)算當(dāng)前時(shí)刻的分?jǐn)?shù)階導(dǎo)數(shù)需要用到之前記憶長度L內(nèi)的所有函數(shù)值，所需計(jì)算內(nèi)存和計(jì)算時(shí)間與記憶長度L的設(shè)置有關(guān)，短時(shí)記憶長度L設(shè)置得越短所需計(jì)算內(nèi)存越小，計(jì)算時(shí)間越短，計(jì)算精度越差；自適應(yīng)記憶法根據(jù)對當(dāng)前時(shí)刻分?jǐn)?shù)階導(dǎo)數(shù)值的不同影響程度，在之前的不同時(shí)間段內(nèi)按照加權(quán)方式進(jìn)行抽取，來計(jì)算當(dāng)前時(shí)刻的分?jǐn)?shù)階導(dǎo)數(shù)值，隨著初始區(qū)間減小，計(jì)算所需時(shí)間和物理內(nèi)存都增大，但是，初始間隔在數(shù)值穩(wěn)定性范圍內(nèi)，自適應(yīng)記憶法計(jì)算精度較高。因此在3種計(jì)算方法中，只有用短時(shí)記憶法得出的波場快照在慢縱波附近與其他兩種方法得出的波場快照有差別。

a.q=15；b.q=10；c.q=5。圖4 三種不同初始間隔黏滯相界固相速度x分量t=0.15 s時(shí)波場快照Fig.4 Vicious phase boundary wavefield snapshots of x-component velocity in solid phase at t=0.15 s by using three different initial intervals

圖5 三種不同初始間隔固相質(zhì)點(diǎn)垂直速度分量時(shí)間記錄對比圖Fig.5 Comparison of the vertical velocity components of particles with three different initial intervals in solid phase

圖6 不同計(jì)算方法固相質(zhì)點(diǎn)垂直速度分量時(shí)間記錄對比圖Fig.6 Comparison of the vertical velocity components of particles with three different numerical methods in solid phase

5 結(jié)論

1)3種計(jì)算方法計(jì)算的波場中快準(zhǔn)P波及準(zhǔn)SV波幾乎相同，只有慢縱波附近有差別。

2)全局記憶法需要較大的計(jì)算內(nèi)存和計(jì)算時(shí)間，計(jì)算效率較低；短時(shí)記憶法所需計(jì)算內(nèi)存和計(jì)算時(shí)間與記憶長度的設(shè)置有關(guān)，記憶長度越大，所需計(jì)算內(nèi)存越大，計(jì)算時(shí)間越長，計(jì)算精度較高；自適應(yīng)記憶法雖然模擬所需時(shí)間與短時(shí)記憶法相比較長，所需存儲(chǔ)量較大，但是精度更高。

3)本文在模擬相同模型情況下，對全局記憶法、短時(shí)記憶法、自適應(yīng)記憶法所需的計(jì)算內(nèi)存、模擬精度和所需時(shí)間進(jìn)行了歸納總結(jié)，可為后續(xù)分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)數(shù)值模擬及新算法開發(fā)和研究提供理論參考基礎(chǔ)；對于分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)正演問題，需要從模擬精度、模擬所需時(shí)間和存儲(chǔ)量3個(gè)方面進(jìn)行利弊權(quán)衡，從中選出精度較高、存儲(chǔ)量較小和模擬所需時(shí)間較短的數(shù)值計(jì)算方法。

(

)：

[1] Biot M A. General Theory of Three-Dimensional Con-solidation[J]. J Appl Phys, 1941, 12: 155-164.

[2] Dvorkin J, Nur A. Dynamic Poroelasticity: A Unified Model with the Squirt and the Biot Mechanisms[J]. Geophysics, 1993，58: 524-533.

[3] Diallo M S, Appel E. Acoustic Wave Propagation in Saturated Porous Media: Reformulation of the Biot/Squirt Flow Theory[J]. J Appl Geophys, 2000, 44: 313-325.

[4]Liu Q R, Katsube N J. The Discovery of a Second Kind of Rotational Wave in Fluid Filled Porous Material[J]. The Journal of the Acoustical Society of America, 1990, 88(2): 1045-1053.

[5] Kjartansson E. Constant-Q Wave Propagation and At-tenuation[J]. Journal of Geophysical Research, 1979, 84: 4737-4748.

[6] Carcione J M. Time-Domain Modeling of Constant-Q Seismic Waves Using Fractional Derivatives[J]. Pure and Applied Geophys, 2002, 159: 1719-1736.

[7]Podlubny. Fractional Differential Equations[Z]. San Diego: Academic Press, 1999.

[8]Sprouse B P. Computational Efficiency of Fractional Diffusion Using Adaptive Time Step Memory[Z]. Dissertations & Theses-Gradworks, 2010.

[9] Caputo M, Mainardi F. Linear Models of Dissipation in Anelastic Solids[J]. La Rivista del Nuovo Cimento, 1971, 1(2): 161-198.

[10] 林孔容. 關(guān)于分?jǐn)?shù)階導(dǎo)數(shù)的幾種不同定義的分析與比較[J]. 閩江學(xué)院學(xué)報(bào), 2003, 24(5)：3-6.

Lin Kongrong, Analysis and Comparision of Different Definition About Fractional Integrals and Derivatives[J]. Journal of Minjiang University, 2003,24(5): 3-6.

[11] 盧明輝，巴晶，楊慧珠. 含粘滯流體孔隙介質(zhì)中的彈性波[J].工程力學(xué)，2009,26(5): 36-40.

Lu Minghui，Ba Jing，Yang Huizhu. Propagation of Elastic Waves in a Viscous Fluid-Saturated Prorous Solid[J]. Engineering Mechanics, 2009, 26(5): 36-40.

[12]Yang D H, Zhang Z J. Poroelastic Wave Equation Including the Biot/Squirt Mechanism and the Solid/Fluid Coupling Anisotropy[J]. Wave Motion, 2002, 35(3): 223-245.

[13] 劉財(cái)，蘭慧田，郭智奇，等. 基于改進(jìn)BISQ機(jī)制的雙相HTI介質(zhì)波傳播偽譜法模擬與特征分析[J].地球物理學(xué)報(bào), 2013, 56(10):3461-3473.

Liu Cai, Lan Huitian, Guo Zhiqi, et al. Pseudo-Spectral Modeling and Feature Analysis of Wave Propagation in Two-Phase HTI Medium Based on Reformulated BISQ Mechanism[J]. Chinese Journal of Geophysics, 2013, 56(10): 3461-3473.

[14] Qiao Z H. Theory and Modelling of Viscoelastic Ani-sotropic Media Using Fractional Time Derivative[C]// 78th EAGE Conference & Exhibition 2016. [S. l.]: EAGE, 2016.

[15] 董良國，馬在田，曹景忠. 一階彈性波動(dòng)方程交錯(cuò)網(wǎng)格高階差分解法穩(wěn)定性研究[J]. 地球物理學(xué)報(bào)，2000，43(6)：856-864.

Dong Liangguo, Ma Zaitian, Cao Jingzhong. A Study on Stability of the Staggered-grid High-Order Difference Method of First-Order Elastic Wave Equation[J]. Chinese Journal of Geophysics，2000, 43(6)：856-864.

[16] Cerjan C, Kosloff D, Kosloff R, et al. A Nonref-lecting Boundary Condition for Discrete Acoustic and Elastic Wave Equation[J]. Geophysics, 1985, 50(4): 705-708.