分數階時間導數計算方法在含黏滯流體黏彈雙相VTI介質波場模擬中的應用

胡 寧，劉 財

吉林大學地球探測科學與技術學院，長春 130026

0 引言

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

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

分數階導數實際上是整數階導數的推廣，關于分數階導數的定義，許多數學家從不同角度分別給出了不同定義，其合理性和科學性在實踐中已經得到了驗證，這個數學分支已在實際問題中得到廣泛應用。比較常見的3種分數階導數定義分別為：Grünwald-Letnikov分數階導數定義[9]、Riemann-Liouville分數階導數定義[9]和Caputo分數階導數定義[9]。Riemann-Liouville分數階導數定義是Grünwald-Letnikov分數階導數定義的擴充，Caputo分數階導數也是對Grünwald-Letnikov分數階導數定義的另一種改進[10]，在階數為負實數和正整數時三者是等價的。在實際應用中，用哪種形式的定義需依照具體情況而定。

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

1 原理

由于Grünwald-Letnikov分數階導數定義是基于Grünwald-Letnikov有限差分近似式導出的，而本文波傳播方程中整數階導數數值模擬采用高階交錯網格有限差分算法進行，故相應地本文中的分數階時間導數計算采用Grünwald-Letnikov分數階導數定義。

1.1 Grünwald-Letnikov分數階導數定義

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

f(x-mh)

(1)

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

1.2 分數階波傳播方程

本文采用基于Kjartansson常Q理論的分數階時間導數黏彈本構關系來描述雙相介質固體骨架的黏彈性效應，并引入流變學本構關系表征孔隙流體的黏滯力學行為，根據Yang等[12]提出的包含BISQ機制和固流耦合各向異性效應的孔隙彈性波動理論，給出基于分數階時間導數常Q黏彈本構關系的含黏滯流體雙相VTI介質二維波傳播方程。

固體骨架黏彈本構方程為

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

(2)

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

(3)

(4)

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

幾何方程為

(5)

式中，eij為固體骨架應變張量分量。

運動方程為

(6)

(7)

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

綜合式(3)—(7)即可導出以固相和流相位移usi和ufi為基本未知量的、基于分數階時間導數常Q黏彈本構關系的含黏滯流體雙相VTI介質波傳播方程。在二維x-z平面內，這組方程的速度-應力形式是由以下8個方程構成的含卷積運算的偏微分方程組：

(8)

其中，

2 穩定性條件和吸收邊界條件

在進行顯式有限差分數值模擬時，需要考慮計算過程的穩定性。分數階計算的穩定性較為復雜，至今缺乏系統的分析研究。本文采用與二維一階速度-應力方程時間二階、空間2N階精度交錯網格有限差分一致的穩定性條件，其近似表達式為[15]

(9)

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

為消除人工邊界反射，采用簡便實用的Cerjan衰減邊界[16]，即沿人工邊界向外擴充H個網格點，進入擴充區域內的波場通過與因子G相乘逐漸衰減為零。

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

(10)

式中，a為衰減系數，可通過多次試驗確定其最佳值。

3 三種算法及其差分格式

對于分數階波傳播方程中的整數階時間導數和空間導數，采用高階交錯網格有限差分算法進行數值計算。由于交錯網格有限差分方法將空間變量的導數在響應網格點之間的半程上進行計算，可以在不增加計算量和內存要求的前提下提高模擬精度，有效減少傳統有限差分的網格頻散。

對于分數階波傳播方程中的分數階時間導數，使用Grünwald-Letnikov分數階導數定義，并分別使用全局記憶法、短時記憶法、自適應記憶法進行數值計算。

采用Grünwald-Letnikov分數階導數定義，分數階導數的離散式為

(11)

式中，定義函數

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

特殊地，當m=0時，τ(2γ,0)=1。

當Δt→0、t=mΔt時，離散化式(11)有

(12)

式中：Δt為時間步長；(j,l)為空間離散點坐標；k為當前時間迭代次數。

3.1 全局記憶法

由于分數階時間導數具有歷史依賴性和空間全域相關特性，所以計算當前時刻的分數階導數需要用到之前所有時刻的函數值。在式(12)中不對m做任何截斷處理，即為全局記憶法，其分數階導數差分格式即為式(12)。

3.2 短時記憶法

短時記憶法通過截斷算子長度，僅考慮當前時刻t以前且對當前時刻影響較大的有限時間區間[t-L,t]。對于本文分數階波傳播方程中的分數階導數，其短時記憶法差分近似式為

(13)

3.3 自適應記憶法

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

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

(14)

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

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

ζ中每個區間內的采樣間隔d的集合為

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

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

根據上述定義，本文分數階波傳播方程中分數階導數應用自適應記憶法的差分近似式為

(15)

式中：g∈N+；mi為時間點集

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

ξ∈N1，mi≤mmax}

中的元素，mmax為當i=imax時的mi。

4 數值模擬

為了對3種方法進行對比分析，選取一個模型進行數值模擬。模型大小為3 000 m×3 000 m。介質物性參數：干燥情況下各向同性背景孔隙介質中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，ρa1=420 kg/m3，ρa3=450 kg/m3，固相體積模量Ks=80 GPa，x方向和z方向滲透率分別為k11=6×10-10m2/s和k33=1×10-9m2/s，η=0.001 Pa·s，P波和S波的品質因子分別為QP=30和QS=15，參考頻率ωr=80 Hz。其他計算參數：Δx=Δz=10 m，Δt=5×10-4s，t=0.2 s。震源子波采用Ricker子波，震源位置為(1 500 m, 1 500 m)，主頻為40 Hz，在模擬中震源分別為x、z方向集中力源。在網格點坐標為(150, 200)處設置檢波點記錄固相質點垂直速度分量的時間記錄。

4.1 全局記憶法

由于全局記憶法未對分數階微分算子進行截斷處理，在數值模擬上認為是精確的。利用全局記憶法對上述模型進行數值模擬的運行時間和所占用物理內存分別為1 395.965 962 s和5 875 MB。圖1給出了利用全局記憶法對波傳播方程進行數值模擬，固相x分量速度和流相z分量速度在t=0.15 s時刻的波場快照。

4.2 短時記憶法不同記憶長度比較

為了比較分析不同短時記憶長度對模擬結果的影響，針對上述模型，分別采用L=0.15、0.10、0.05 s三種不同短時記憶長度進行數值模擬，表1給出了所用時間和內存。從表1可以看出，隨著短時記憶長度的減小，計算所需時間和物理內存都在減??；這是因為采用的短時記憶長度越小，所計算的當前時刻之前的時間節點越少。

表1采用3種不同短時記憶長度模擬的運行時間與所占物理內存

Table1Runningtimeandoccupiedphysicalmemoryofthreedifferentshort-termmemorylengthsimulations

L/s運行時間/s物理內存/MB0.151250.06432450110.101032.77132049960.05693.6316934878

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

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

準SV波的精度越來越差。

4.3 自適應記憶法不同初始間隔比較

為了比較分析不同初始間隔對模擬結果的影響，針對上述模型，分別采用q=15、10、5三種不同初始間隔進行數值模擬，表2給出了所用時間和內存。從表2可以看出，隨著初始間隔減小，計算所需時間和物理內存都在增大；這是因為初始間隔越小，之前的時間就會被分成越多的區間，就會有更多的時間節點參與計算。

表2采用3種不同初始間隔運行時間與所占物理內存結果

Table2Runningtimeandoccupiedphysicalmemoryofthreedifferentinitialintervalsimulations

q運行時間/s物理內存/MB151243.4999644583101357.584104508351461.3141835289

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

4.4 三種不同計算方法對比分析

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

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

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

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

5 結論

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

2)全局記憶法需要較大的計算內存和計算時間，計算效率較低；短時記憶法所需計算內存和計算時間與記憶長度的設置有關，記憶長度越大，所需計算內存越大，計算時間越長，計算精度較高；自適應記憶法雖然模擬所需時間與短時記憶法相比較長，所需存儲量較大，但是精度更高。

3)本文在模擬相同模型情況下，對全局記憶法、短時記憶法、自適應記憶法所需的計算內存、模擬精度和所需時間進行了歸納總結，可為后續分數階時間導數數值模擬及新算法開發和研究提供理論參考基礎；對于分數階時間導數正演問題，需要從模擬精度、模擬所需時間和存儲量3個方面進行利弊權衡，從中選出精度較高、存儲量較小和模擬所需時間較短的數值計算方法。

(

)：

[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] 林孔容. 關于分數階導數的幾種不同定義的分析與比較[J]. 閩江學院學報, 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] 盧明輝，巴晶，楊慧珠. 含粘滯流體孔隙介質中的彈性波[J].工程力學，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] 劉財，蘭慧田，郭智奇，等. 基于改進BISQ機制的雙相HTI介質波傳播偽譜法模擬與特征分析[J].地球物理學報, 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] 董良國，馬在田，曹景忠. 一階彈性波動方程交錯網格高階差分解法穩定性研究[J]. 地球物理學報，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.