胡 寧，劉 財(cái)

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

0 引言

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

Kjartansson[5]提出的常Q模型認(rèn)為：衰減因子與頻率成線(xiàn)性反比關(guān)系，在常規(guī)地震頻帶內(nèi)，Q值與頻率基本無(wú)關(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)與過(guò)去所有狀態(tài)有關(guān)，因而所需存儲(chǔ)量和計(jì)算量均很大，計(jì)算效率較低。Podlubny[7]提出“短時(shí)記憶法”對(duì)分?jǐn)?shù)階微分算子進(jìn)行截?cái)?，僅考慮當(dāng)前時(shí)刻t以前的有限時(shí)間區(qū)間[t-L,t]，隨著短時(shí)記憶長(zhǎng)度L的選取不同，誤差也不同。過(guò)去不同時(shí)間點(diǎn)對(duì)當(dāng)前時(shí)間點(diǎn)的影響是不同的，接近當(dāng)前時(shí)刻時(shí)間點(diǎn)對(duì)當(dāng)前時(shí)刻的影響較大，而距離當(dāng)前時(shí)刻時(shí)間間隔較長(zhǎng)的時(shí)間點(diǎn)對(duì)當(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í)際問(wèn)題中得到廣泛應(yīng)用。比較常見(jiàn)的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ù)也是對(duì)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é)合起來(lái)，引入流變學(xué)本構(gòu)關(guān)系描述孔隙流體的黏滯性力學(xué)行為[11]，建立含黏滯流體黏彈雙相VTI(橫向各向同性)介質(zhì)中的分?jǐn)?shù)階波傳播方程；并在時(shí)間域進(jìn)行數(shù)值模擬，采用全局記憶法、短時(shí)記憶法、自適應(yīng)記憶法對(duì)分?jǐn)?shù)階微分算子進(jìn)行數(shù)值計(jì)算，整數(shù)階時(shí)間、空間微分算子均采用高階交錯(cuò)網(wǎng)格有限差分算法；最后對(duì)3種方法各自的特點(diǎn)進(jìn)行歸納總結(jié)，并對(duì)3種方法進(jìn)行對(duì)比分析。

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)

稱(chēng)為f(x)的γ階Grünwald-Letnikov分?jǐn)?shù)階導(dǎo)數(shù)。式中：h為步長(zhǎng)；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ù)的歷史依賴(lài)性和空間全域相關(guān)特性。

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

本文采用基于Kjartansson常Q理論的分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)黏彈本構(gòu)關(guān)系來(lái)描述雙相介質(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的元素，對(duì)于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ì)算過(guò)程的穩(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)格步長(zhǎng)；Δt為時(shí)間間隔；v=max{vP//,vP⊥,vS}，vP//、vP⊥、vS分別為平行于地層的P波速度分量和垂直于地層的P波速度分量、S波速度分量；Cn為交錯(cuò)網(wǎng)格差分算子系數(shù)。

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

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

(10)

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

3 三種算法及其差分格式

對(duì)于分?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)格頻散。

對(duì)于分?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í)間步長(zhǎng)；(j,l)為空間離散點(diǎn)坐標(biāo)；k為當(dāng)前時(shí)間迭代次數(shù)。

3.1 全局記憶法

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

3.2 短時(shí)記憶法

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

(13)

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

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

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

(14)

對(duì)于式(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ù)值模擬

為了對(duì)3種方法進(jìn)行對(duì)比分析，選取一個(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 全局記憶法

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

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

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

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

圖1 全局記憶法黏滯相界固相速度x分量(a)和流相速度z分量(b)t=0.15 s時(shí)波場(chǎng)快照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í)記憶長(zhǎng)度黏滯相界固相速度x分量t=0.15 s時(shí)波場(chǎng)快照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ǎng)度固相質(zhì)點(diǎn)垂直速度分量時(shí)間記錄對(duì)比圖Fig.3 Comparison of the vertical velocity components of particles with three different short-term memory lengths in solid phase

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

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

為了比較分析不同初始間隔對(duì)模擬結(jié)果的影響，針對(duì)上述模型，分別采用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種初始間隔對(duì)波傳播方程進(jìn)行數(shù)值模擬，固相x分量速度在t=0.15 s時(shí)刻的波場(chǎng)快照，可以看出，基于自適應(yīng)記憶法、采用3種不同初始間隔得到的波場(chǎng)快照無(wú)明顯差異。為了更明顯地顯示差異，圖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ì)算方法對(duì)比分析

為了比較分析不同分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)計(jì)算方法對(duì)模擬結(jié)果的影響，針對(duì)上述模型，圖6給出了利用全局記憶法、L=0.1 s短時(shí)記憶法和q=10自適應(yīng)記憶法對(duì)波傳播方程進(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ì)算精度上來(lái)看，自適應(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ù)需要用到之前記憶長(zhǎng)度L內(nèi)的所有函數(shù)值，所需計(jì)算內(nèi)存和計(jì)算時(shí)間與記憶長(zhǎng)度L的設(shè)置有關(guān)，短時(shí)記憶長(zhǎng)度L設(shè)置得越短所需計(jì)算內(nèi)存越小，計(jì)算時(shí)間越短，計(jì)算精度越差；自適應(yīng)記憶法根據(jù)對(duì)當(dāng)前時(shí)刻分?jǐn)?shù)階導(dǎo)數(shù)值的不同影響程度，在之前的不同時(shí)間段內(nèi)按照加權(quán)方式進(jìn)行抽取，來(lái)計(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í)記憶法得出的波場(chǎng)快照在慢縱波附近與其他兩種方法得出的波場(chǎng)快照有差別。

a.q=15；b.q=10；c.q=5。圖4 三種不同初始間隔黏滯相界固相速度x分量t=0.15 s時(shí)波場(chǎng)快照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í)間記錄對(duì)比圖Fig.5 Comparison of the vertical velocity components of particles with three different initial intervals in solid phase

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

5 結(jié)論

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

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

3)本文在模擬相同模型情況下，對(duì)全局記憶法、短時(shí)記憶法、自適應(yīng)記憶法所需的計(jì)算內(nèi)存、模擬精度和所需時(shí)間進(jìn)行了歸納總結(jié)，可為后續(xù)分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)數(shù)值模擬及新算法開(kāi)發(fā)和研究提供理論參考基礎(chǔ)；對(duì)于分?jǐn)?shù)階時(shí)間導(dǎo)數(shù)正演問(wèn)題，需要從模擬精度、模擬所需時(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] 董良國(guó)，馬在田，曹景忠. 一階彈性波動(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.