偶應力理論框架下的彈性波數(shù)值模擬與分析

王之洋， 李幼銘 ， 白文磊

1 北京化工大學信息科學與技術(shù)學院， 北京 100029 2 中國科學院地質(zhì)與地球物理研究所， 中國科學院油氣資源研究重點實驗室， 北京 100029 3 非對稱性彈性波動方程聯(lián)合研究組， 北京 100029 4 高鐵地震學聯(lián)合研究組， 北京 100029

0 引言

地球介質(zhì)是一種非均勻、非完全彈性、各向異性、多相態(tài)的介質(zhì)(傅承義等，1985)，在經(jīng)典地震學理論框架下，先人發(fā)展了數(shù)不勝數(shù)的地震技術(shù)，為社會發(fā)展做出了巨大貢獻；可是，當前技術(shù)仍有難以逾越的障礙，亟需破局.因此，尋求更接近真實地球介質(zhì)的波動理論，構(gòu)建更加符合物理實際的地球介質(zhì)模型，是地球物理學發(fā)展的一個趨勢.

2019年8月底，中國科學院地質(zhì)與地球物理研究所、北京大學以及西安交通大學等單位組成的高鐵地震學聯(lián)合研究組，在河北省定興縣野外采集資料，檢波器分別放置于高鐵的高架橋沿線以及橋墩處.通過對中國科學院地質(zhì)與地球物理研究所“高鐵地震學研究協(xié)調(diào)辦公室”提供的數(shù)據(jù)的分析，可以觀察到大量的旋轉(zhuǎn)運動分量.由于基于經(jīng)典連續(xù)介質(zhì)力學推導的彈性波動方程，從理論出發(fā)點上就去除了獨立的旋轉(zhuǎn)自由項，且在其理論框架內(nèi)，介質(zhì)被視為一個連續(xù)的質(zhì)量體(Yang et al.，2002)，事實上，不論是人造還是天然的介質(zhì)都存在著復雜的內(nèi)部結(jié)構(gòu)(Zhu et al.，2019)，沒有介質(zhì)是理想的連續(xù)統(tǒng)，廣義連續(xù)介質(zhì)力學更適合描述具有更加復雜內(nèi)部結(jié)構(gòu)的情況(Jirásek，2004).于是，我們嘗試啟用廣義連續(xù)介質(zhì)力學理論，推導偶應力理論框架下的彈性波動方程，將其表達形式以及數(shù)值模擬結(jié)果與傳統(tǒng)彈性波動方程進行對比，研究發(fā)現(xiàn)推導的具非對稱性的波動方程所具有的理論意義和應用價值.

傳統(tǒng)地震學是基于經(jīng)典的連續(xù)介質(zhì)力學原理建立的，在該理論中，應用角動量守恒定律，我們可以得到應力張量是對稱的結(jié)論.而對于廣義連續(xù)介質(zhì)力學，組成介質(zhì)的每個點都被視為一個體積不為零的剛性微元體，微元體上可以承受均布的體力偶和面力偶以產(chǎn)生旋轉(zhuǎn)的作用，在應用角動量守恒定律時，必然導致應力張量是非對稱的.應力張量是否對稱，或者說力學上是否具有對稱性，是連續(xù)介質(zhì)力學和廣義連續(xù)介質(zhì)力學的本質(zhì)區(qū)別和重要差異.經(jīng)典連續(xù)介質(zhì)力學理論下的旋轉(zhuǎn)運動，在小變形假設(shè)下，其不造成形變(Aki and Richards，2002)，所以在小變形假設(shè)下，位移梯度全部用對稱的應變張量表示是合理的.對于各向同性介質(zhì)，旋轉(zhuǎn)造成的位移擾動包含在橫波分量里，可以通過求取位移旋度的二分之一獲得旋轉(zhuǎn)分量.而廣義連續(xù)介質(zhì)力學理論下的旋轉(zhuǎn)運動則不完全相同，其是由于介質(zhì)內(nèi)部存在的微孔縫隙結(jié)構(gòu)所導致的力學不對稱性所產(chǎn)生的.

廣義和經(jīng)典連續(xù)介質(zhì)力學理論是并行發(fā)展的兩個分支，廣義連續(xù)介質(zhì)力學(Cosserat and Cosserat，1909；Mindlin and Tiersten，1962；Koiter，1964；Eringen and Suhubi，1964；Mindlin，1964，1965；Eringen，1966，1990，1999；Aifantis，1999；Yang et al.，2002；Lam et al.，2003；Askes and Aifantis，2011；Hadjesfandiari and Dargush，2011；Chen et al.，2011；Polizzotto，2013；Li et al.，2015；Gopalakrishnan，2016；De Domenico et al.，2019；Zhu et al.，2019)通過對經(jīng)典連續(xù)介質(zhì)力學基本假設(shè)、原理或者限制條件進一步放松，引入狀態(tài)變量的高階導數(shù)項，描述由假設(shè)偶應力存在而導致的介質(zhì)內(nèi)部微結(jié)構(gòu)相互作用所產(chǎn)生的介質(zhì)非均勻性，以豐富經(jīng)典連續(xù)介質(zhì)力學理論的內(nèi)容.廣義連續(xù)介質(zhì)力學理論可追溯至19世紀末，其最基本的理論方法就是偶應力理論，其他理論方法分支就是在偶應力理論的基礎(chǔ)上被建立起來的.Voigt(1887)指出關(guān)于物體的一部分對其鄰近部分的作用可能引起體力偶和面力偶的猜想，并認為應力張量可能是非對稱的.Cosserat和Cosserat(1909)驗證了這個猜想，提出了Cosserat連續(xù)介質(zhì)理論.在該理論中，每個組成介質(zhì)的點都被看成是一個非零體積的剛性微元體，具有包括旋轉(zhuǎn)和位移等6個自由度.Toupin(1962)對該理論進行了拓展，建立了完全彈性固體有限變形的本構(gòu)關(guān)系.Mindlin和Tiersten(1962)將Toupin(1962)的理論進行了簡化，其假定微元體與周圍宏觀介質(zhì)的相對旋轉(zhuǎn)為零，曲率張量表示為位移的二階梯度，建立了約束轉(zhuǎn)動的偶應力理論，至此，偶應力理論的框架被基本建立起來.Yang等(2002)引入高階力矩平衡理論，提出僅具有一個材料尺度參數(shù)的改進Cosserat連續(xù)介質(zhì)理論，這一理論也被稱為修正的偶應力理論(Park and Gao，2006；Ma et al.，2008；Kong et al.，2008).Hadjesfandiari和Dargush(2011)提出了與Yang等(2002)完全相反的偶應力理論，他們認為介質(zhì)的應變能密度函數(shù)與經(jīng)典應變張量以及曲率張量的反對稱部分有關(guān)，而與曲率張量的對稱部分無關(guān).通過對實驗數(shù)據(jù)和實際應用的分析，Yang等(2002)的理論更加符合物理實際(Reddy and Kim，2012；Jung et al.，2014；Shaat et al.，2014；王之洋等，2020).Chen等(2011)提出了一種新的修正的偶應力理論，并將其推廣到各向異性介質(zhì).

本文研究推導偶應力理論框架下的具非對稱性的彈性波動方程，并應用該方程進行數(shù)值模擬試驗，將合成記錄與傳統(tǒng)彈性波動方程的合成結(jié)果以及實際數(shù)據(jù)進行對比分析，研究包含介質(zhì)特征尺度參數(shù)的具非對稱性的彈性波動方程對地震波傳播的影響與規(guī)律.本文首先從經(jīng)典平衡方程、幾何方程和本構(gòu)方程出發(fā)，在經(jīng)典平衡方程中加入體力偶和面力偶的成分，在幾何方程中新增加對稱曲率張量和位移高階導數(shù)的關(guān)系，在本構(gòu)方程加入對稱曲率張量和對稱偶應力張量的關(guān)系，構(gòu)建新的平衡方程、幾何方程和本構(gòu)方程，并推導偶應力理論框架下包含介質(zhì)特征尺度參數(shù)的非對稱彈性波動方程的數(shù)學表達式，同時，描述引入的旋轉(zhuǎn)自由項以及介質(zhì)特征尺度參數(shù)的物理意義.其次，在均勻介質(zhì)模型和Marmousi模型上進行數(shù)值模擬測試，并將數(shù)值模擬結(jié)果與傳統(tǒng)彈性波動方程的數(shù)值模擬結(jié)果進行對比，分析非對稱彈性波動方程所增加的旋轉(zhuǎn)自由項對地震波傳播的影響.最后，在高鐵橋墩模型上進行數(shù)值模擬試驗，并與實際數(shù)據(jù)以及傳統(tǒng)彈性波動方程的合成記錄進行對比，進一步研究分析介質(zhì)特征尺度參數(shù)和旋轉(zhuǎn)自由項對地震波傳播的影響.

1 偶應力理論框架下的彈性波動方程

相比于傳統(tǒng)彈性波動方程的推導，平衡、幾何、本構(gòu)三組方程的基本形式?jīng)]有改變，只是對其進行了拓展和補充.基于此邏輯，推導得到的偶應力理論框架下的彈性波動方程只會在傳統(tǒng)方程的基礎(chǔ)上增加了旋轉(zhuǎn)自由項，其通過位移的高階導數(shù)描述由力學不對稱特征所導致的旋轉(zhuǎn)運動.如果去掉此旋轉(zhuǎn)自由項，偶應力理論框架下的彈性波動方程就和傳統(tǒng)彈性波動方程有一致的數(shù)學表達.

1.1 平衡方程

(1)

(2)

應力張量σij和偶應力張量μij均滿足線動量平衡方程(公式(3))和角動量平衡方程(公式(4))：

(3)

(4)

其中，Vi,Ωi分別為線動量和角動量.Fi,Mi分別是體力和體力偶.eijk為置換符號.rj是位矢.ui是位移矢量，ωi是旋轉(zhuǎn)矢量.

對公式(3)和(4)分別應用散度定理，由表面積分變換為體積分，則公式(3)和(4)可改寫為：

(5)

(6)

從公式(6)中可略去高階微量，由公式(5)和(6)可以分別得到微分形式的平衡方程：

(7)

μji,j+eistσst+Mi=0.

(8)

公式(7)和(8)可分別展開為：

(9)

(10)

在經(jīng)典連續(xù)介質(zhì)力學中，沒有考慮偶應力與體力偶的作用，也即μij=0,Mi=0，因此，可由公式(8)，得到：

eistσst=0.

(11)

公式(11)表明，經(jīng)典連續(xù)介質(zhì)力學理論框架下，應力張量是對稱的.這就意味著應力張量有6個獨立量，這里有3個平衡方程，另外的3個方程為本構(gòu)方程，6個方程可以定解.

而在偶應力理論中，當考慮偶應力或體力偶時，切應力互等定理不成立了，這導致應力張量的不對稱性，這時候，應力張量σij可以分解為對稱應力張量τij與反對稱應力張量rij之和：

(12)

如果將公式(12)代入公式(7)和(8)，將應力張量σij用對稱應力張量τij與反對稱應力張量rij表示，則公式(7)和(8)可以改寫為：

(13)

μji,j+eistrst+Mi=0，

(14)

其中，s,t=1,2,3.

同時基于公式(8)，建立反對稱應力張量rij與偶應力張量μij的關(guān)系式：

(15)

(16)

將公式(15)等號兩邊求散度，并將公式(16)代入，可得：

(17)

將公式(12)和(17)一齊代入公式(13)，可消除反對稱應力張量rij，得到只包含對稱應力張量τij以及偶應力張量的偏斜部分mij的平衡方程：

(18)

公式(18)即為偶應力理論框架下的平衡方程，其構(gòu)建了偶應力張量、應力張量、體力偶和體力的關(guān)系，如果去除了偶應力張量和體力偶，該平衡方程就為經(jīng)典理論的平衡方程.

進一步探討偶應力張量μij的對稱性.Yang等(2002)基于連續(xù)體內(nèi)偶力矩為零的假設(shè)，引入了高階力矩平衡關(guān)系，對偶應力張量μij施加了約束，見公式(19)：

(19)

對公式(19)應用散度定理，可得：

(20)

因為μji,j+eistσst+Mi=0，可得eistμst=0，即偶應力張量μij是對稱的.

1.2 邊界條件

考慮不存在體力和體力偶的單位體積元能量守恒(Hadjesfandiari and Dargush，2011)：

(21)

其中，n是方向矢量，δu為虛位移矢量，δω為虛旋轉(zhuǎn)矢量，δw為能量改變量.

因為有：

n·μ·δω=n·μ·nn·δω+n·μ·(I-nn)·δω

(22)

其中，nn，(I-nn)分別代表法向和切向方向.

對公式(22)等號右邊第一項應用Hamiltonian算子，可得：

(23)

并考慮表面Sa是光滑的，即有：

(24)

則公式(21)可寫為：

(25)

公式(25)表明，因為表面Sa上偶應力張量的法向分量μn n與應力張量σij的組合作為δu的系數(shù)，因此在偶應力介質(zhì)中僅僅包含5個邊界條件，即應力矢量p(n)的3個分量以及力偶矢量m(n)的2個切向分量.

由公式(25)，可得到偶應力理論框架下的邊界條件：

(26)

m(n)=n·μ·(I-nn).

(27)

真實的邊界上一般沒有力矩作用，這意味著真實邊界上自然邊界條件處處為零，即m(n)=0.但偶應力張量μij在物體內(nèi)部卻會產(chǎn)生，在任意體積(包括了微元體積)的表面上卻存在著非零的m(n).

1.3 幾何方程

這里僅僅考慮微小變形|ui,j|?1的運動學(幾何方程).在笛卡爾坐標系中，參考構(gòu)型中的兩點A(位置為xi)和B(位置為xi+dxi)之間的相對位移為(Aki and Richards，2002)：

dui=ui,jdxj，

(28)

位移梯度張量ui,j可分解為對稱的無窮小應變張量εij和反對稱的無窮小旋轉(zhuǎn)張量ωij之和：

(29)

(30)

(31)

(32)

因此，旋轉(zhuǎn)矢量ωi可以通過求取位移場旋度的二分之一得到.

在微小變形|ui,j|?1假設(shè)下，經(jīng)典連續(xù)介質(zhì)力學理論下的旋轉(zhuǎn)運動，是剛性旋轉(zhuǎn)，其不造成形變(Aki and Richards，2002)，位移梯度可以全部用對稱的應變張量表示.偶應力理論下，引入了一種假設(shè)偶應力存在而出現(xiàn)的新的旋轉(zhuǎn)運動，使微元線段dxi產(chǎn)生了扭轉(zhuǎn)和彎曲形變，其由相鄰兩點A,B之間的相對轉(zhuǎn)動dωi描述，并引入對稱的曲率張量ij描述這種形變(Yang et al.，2002)：

dωi=ωi,jdxj，

(33)

這樣，在微小變形中，連續(xù)體A點臨近B點的形變總結(jié)起來有4種：隨著A點的平動，相對A點的伸縮，相對A點的轉(zhuǎn)動以及微元線段AB的扭轉(zhuǎn)彎曲.

通過構(gòu)建曲率張量以及應變張量與位移的關(guān)系，我們可以推導得到幾何方程，相比于經(jīng)典理論的幾何關(guān)系，偶應力理論框架下增加了曲率張量與位移之間的關(guān)系，幾何方程見公式(34)：

(34)

1.4 本構(gòu)方程

由虛功原理，在公式(7)兩端乘以虛位移δui,在公式(8)兩端乘以虛旋轉(zhuǎn)δωi，并在整個體積Va上積分，可得(Koiter，1964)：

(35)

(36)

對公式(35)和(36)應用鏈式求導法則以及散度定理，并相加，可得：

(37)

將公式(1)和(2)代入公式(37)，并應用散度定理，將公式(37)等號右邊第一項的面積分變?yōu)轶w積分，可得：

(38)

將平衡方程(公式(18))、幾何方程(公式(34))以及公式(12)，代入公式(38)并化簡，可得功共軛關(guān)系：

(39)

對于各向同性介質(zhì)，應用廣義胡克定律以及功共軛條件(公式(39))，可得偶應力理論框架下的本構(gòu)方程：

(40)

其中，λ,μ為Lamé常數(shù)，η為反映介質(zhì)微旋轉(zhuǎn)運動特性的參數(shù)，η=μl2，l為介質(zhì)特征尺度參數(shù)，其是為了平衡無量綱的應變和有量綱的曲率張量(其量綱為長度的倒數(shù))兩項間的量綱，也為了描述偶應力張量和曲率張量之間的本構(gòu)關(guān)系.

1.5 偶應力理論框架下的彈性波動方程

對于各向同性介質(zhì)，將幾何方程(公式(34))和本構(gòu)方程(公式(40))以及平衡方程(公式(18))聯(lián)立，并忽略體力以及體力偶，可推導得到偶應力理論框架下的彈性波動方程：

(41)

旋轉(zhuǎn)自由項中的η是介質(zhì)內(nèi)部微結(jié)構(gòu)相互作用導致的旋轉(zhuǎn)效應的表征，如果η=0，則旋轉(zhuǎn)效應消失，偶應力理論框架下的彈性波動方程就和傳統(tǒng)彈性波動方程有一致的數(shù)學表達.Teisseyre和Boratńy(tǒng)ski(2003)通過理論推導證明，在顆粒介質(zhì)或者有裂隙的連續(xù)體內(nèi)，由于微缺陷/微結(jié)構(gòu)的存在，應力或者應力場會出現(xiàn)不對稱性力學特征，由此產(chǎn)生一種新的旋轉(zhuǎn)運動.Ba?ant(2002)認為金屬材料的微缺陷/微結(jié)構(gòu)一般由晶體的位錯引起，晶體位錯在納米量級；而巖土介質(zhì)的微缺陷/微結(jié)構(gòu)一般是由微夾雜/微裂縫/微孔隙/孔洞等引起，尺度可以在微毫米或更高量級.黃文雄和徐可(2014)通過實驗證明，介質(zhì)特征尺度l與組成介質(zhì)的材料顆粒平均直徑以及幾何形狀有關(guān)，同時，也與變形局部化現(xiàn)象有關(guān).另外，與應力集中效應也有關(guān)(鮑亦興和毛昭宙，1993).如果忽略一切附加因素，介質(zhì)特征尺度參數(shù)可直接視為介質(zhì)微孔縫隙特征尺度參數(shù).相反，如果考慮組成介質(zhì)的材料顆粒的幾何特征、變形局部化現(xiàn)象以及應力集中效應，介質(zhì)特征尺度參數(shù)可視為介質(zhì)微孔縫隙特征尺度參數(shù)與系數(shù)ζ相乘的結(jié)果，系數(shù)ζ可通過實驗獲得，王之洋等(2020)通過對高鐵實際數(shù)據(jù)和理論合成數(shù)據(jù)的對比分析，也提取了該系數(shù).當然，為了進一步研究介質(zhì)特征尺度參數(shù)和旋轉(zhuǎn)自由項的物理意義，必須結(jié)合巖石物理測量，開展對比驗證研究和技術(shù)生成.

2 數(shù)值模擬與分析

2.1 均勻各向同性介質(zhì)模型

首先，我們構(gòu)建均勻各向同性介質(zhì)模型，模型大小為nx=nz=400，網(wǎng)格間距為dx=dz=8 m，縱波速度為2.0 km·s-1，橫波速度為1.1547 km·s-1，密度為2000 kg·m-3，采用主頻為25 Hz的Ricker子波，并在模型中間激發(fā)，時間步長Δt=0.001 s，介質(zhì)微孔縫隙特征尺度參數(shù)為200 μm.為了排除數(shù)值頻散和邊界反射對數(shù)值模擬結(jié)果的影響，我們采用高階優(yōu)化有限差分算子以及優(yōu)化的PML(Perfectly Matched Layer，完全匹配層)邊界條件進行數(shù)值模擬，以更好地對傳統(tǒng)彈性波動方程以及偶應力理論框架下的彈性波動方程的合成記錄進行對比分析.

圖1和圖2分別是應用傳統(tǒng)彈性波動方程和偶應力理論框架下的彈性波動方程進行數(shù)值模擬所得到的x分量和z分量的合成波場快照.從圖1和圖2中可以看出，不論是x分量，還是z分量，偶應力理論框架下的彈性波動方程增加的旋轉(zhuǎn)自由項所帶來的位移擾動都對地震波傳播帶來了明顯的影響，出現(xiàn)了新的波場分量.比較圖1a和圖1b的x分量波場快照，偶應力理論框架下的旋轉(zhuǎn)運動，對縱波沒有影響，但卻使得橫波在傳播過程中出現(xiàn)了物理頻散，如果將圖1a和圖1b相減，可以更加清晰地觀察到橫波波場的變化.圖2中z分量波場快照的比較，同樣可以得到此結(jié)論.

圖1 應用不同彈性波動方程得到的合成波場快照(x分量)(a) 傳統(tǒng)彈性波動方程； (b) 偶應力理論框架下的彈性波動方程； (c) 圖(a)和圖(b)的差值.Fig.1 Synthetic wavefield snapshots (x component) using different elastic wave equations(a) The conventional elastic wave equations; (b) Elastic wave equations in the frame of the couple stress theory; (c) The difference between figure (a) and figure (b).

圖2 應用不同彈性波動方程得到的合成波場快照(z分量)(a) 傳統(tǒng)彈性波動方程； (b) 偶應力理論框架下的彈性波動方程； (c) 圖(a)和圖(b)的差值.Fig.2 Synthetic wavefield snapshots (z component) using different elastic wave equations(a) The conventional elastic wave equations; (b) Elastic wave equations in the frame of the couple stress theory; (c) The difference between figure (a) and figure (b).

偶應力理論框架下的彈性波動方程包含獨立的旋轉(zhuǎn)自由項，可描述由考慮介質(zhì)內(nèi)部微孔縫隙結(jié)構(gòu)相互作用所帶來的旋轉(zhuǎn)運動以及其位移擾動的傳播，這種位移擾動可以在合成記錄中被清晰地觀察到，同時，該旋轉(zhuǎn)運動對縱波沒有影響，而使橫波的波場出現(xiàn)了新的成分，導致其以頻散的方式傳播.

2.2 Marmousi模型

為了進一步對比應用傳統(tǒng)彈性波動方程和偶應力理論框架下的彈性波動方程得到的合成地震記錄的差異，我們在Marmousi模型上進行測試.速度模型大小為737×750，橫向采樣間隔為12.5 m，縱向采樣間隔為4 m，密度ρ=1500 kg·m-3，橫波速度和縱波速度的比值為vp/vs=1.7.采用主頻為25 Hz的Ricker子波，震源位置設(shè)置在模型頂層中間處，時間步長Δt=0.0005 s，記錄時長為t=10 s，同樣，設(shè)置介質(zhì)微孔縫隙特征尺度參數(shù)也為200 μm.

圖3和圖4分別是應用傳統(tǒng)彈性波動方程和偶應力理論框架下的彈性波動方程對Marmousi模型進行數(shù)值模擬所得到的x分量和z分量的合成地震記錄.對于較復雜的Marmousi模型，介質(zhì)內(nèi)部孔縫隙所帶來的位移擾動的影響依然存在，從圖3和圖4中可以清晰地觀察到，相比于應用傳統(tǒng)彈性波動方程得到的合成地震記錄，偶應力理論框架下的彈性波動方程下合成地震記錄，在圖中紅色箭頭標示區(qū)域，存在新的波場擾動，即使我們設(shè)置的介質(zhì)微孔縫隙特征尺度參數(shù)在微米量級，相比于子波波長，量級上的差距巨大，但其所帶來的對地震波傳播的影響仍然能夠在地震記錄上被觀察到.由于采用高階優(yōu)化的差分算子進行數(shù)值模擬，數(shù)值頻散已經(jīng)得到了有效的壓制，因此，圖3c和圖4c中觀察到的地震波場差異，有理由相信，是地震波傳播中出現(xiàn)的物理頻散.

圖3 在Marmousi模型上應用不同彈性波動方程得到的合成地震記錄(x分量)(a) 傳統(tǒng)彈性波動方程； (b) 偶應力理論框架下的彈性波動方程； (c) 圖(a)和圖(b)的差值.Fig.3 Synthetic shot records (x component) for Marmousi model using different elastic wave equations(a) The conventional elastic wave equations; (b) Elastic wave equations in the frame of the couple stress theory; (c) The difference between figure (a) and figure (b).

圖4 在Marmousi模型上應用不同彈性波動方程得到的合成地震記錄(z分量)(a) 傳統(tǒng)彈性波動方程； (b) 偶應力理論框架下的彈性波動方程； (c) 圖(a)和圖(b)的差值.Fig.4 Synthetic shot records (z component) for Marmousi model using different elastic wave equations(a) The conventional elastic wave equations; (b) Elastic wave equations in the frame of the couple stress theory; (c) The difference between figure (a) and figure (b).

2.3 高鐵橋墩模型

最后，我們應用偶應力理論框架下的彈性波動方程，在高鐵橋墩模型上進行數(shù)值模擬試驗，并將合成地震記錄與應用傳統(tǒng)彈性波動方程得到的合成地震記錄進行對比.

我們構(gòu)建了一個基于高鐵列車震源的簡化橋墩模型.速度模型大小為400×100，網(wǎng)格間距為1 m×1 m，分為兩層，第一層是厚度為10 m的低速土壤層，密度為400 kg·m-3，縱波速度為0.5 km·s-1；第二層是厚度為90 m的高速圍巖層，密度為1400 kg·m-3，縱波速度為1.6 km·s-1，橫波速度和縱波速度的比值為vp/vs=1.7.橋墩數(shù)量為3個，橋墩間距為28 m，每個橋墩插入地下50 m深度直達圍巖，每個橋墩的地下部分被看作是間隔為10 m的5個點源，地震波在橋墩中的傳播速度為vp=4000 m·s-1.高鐵列車車廂為16節(jié)，每節(jié)車廂的長度為28 m，列車速度為300 km·h-1，視為多個移動Ricker震源的疊加，Ricker子波主頻設(shè)置為20 Hz，時間步長Δt=0.0005 s，記錄時長為t=10 s，設(shè)置表層土壤和深層圍巖的微孔縫隙特征尺度參數(shù)均為700 μm.基于高鐵列車震源的簡化橋墩模型如圖5所示.

圖5 基于高鐵列車震源的簡化橋墩模型Fig.5 A simplified bridge pier model of the source based on high-speed trains

圖6是不同彈性波動方程得到的合成地震記錄(z分量).從圖6可以觀察到，應用傳統(tǒng)彈性波動方程和偶應力理論框架下的彈性波動方程所得到的合成地震記錄在時域上的波形相似，我們的合成地震記錄可視為不同空間位置的以一定時間間隔激發(fā)的固定震源記錄的線性疊加，基于線性疊加的特性，可以通過設(shè)置模型參數(shù)，使模型的橋墩和橋梁長度和實際情況更加接近，從而使時域合成地震記錄與實際記錄具有更好的一致性.通過對頻譜的分析，我們發(fā)現(xiàn)，當設(shè)置高速列車車速為300 km·h-1時，合成地震記錄的能量集中在大約10 Hz以及25 Hz左右，其中，頻率10 Hz附近的能量最大，這與實際高鐵地震記錄的頻譜能量分布基本相符.

圖6 不同彈性波動方程得到的合成地震記錄(z分量)(a) 傳統(tǒng)彈性波動方程； (b) 偶應力理論框架下的彈性波動方程； 圖(c)，(d)分別是圖(a)，(b)的幅度頻譜.Fig.6 Synthetic shot records using different elastic wave equations (z component)(a) The conventional elastic wave equations; (b) Elastic wave equations in the frame of the couple stress theory; Figures (c) and (d) are the amplitude spectrum of figures (a) and (b), respectively.

進一步分析，將傳統(tǒng)彈性波動方程和偶應力理論框架下的彈性波動方程所得到的合成地震記錄在時頻域進行對比，見圖7.通過對圖7的觀察，我們發(fā)現(xiàn)，由于引入了包含介質(zhì)特征尺度參數(shù)的旋轉(zhuǎn)自由項，應用偶應力理論框架下的彈性波動方程所得到的合成地震記錄在時間域波形上出現(xiàn)了衰減，其幅度頻譜相較于傳統(tǒng)彈性波動方程的合成記錄，大于20 Hz的頻譜能量都出現(xiàn)了小幅衰減，其中，20 Hz到30 Hz頻段的能量衰減相對最大，且整體能量向低頻10 Hz處移動，更趨近于實際高鐵地震記錄的頻譜.我們認為，由于考慮了介質(zhì)內(nèi)部的微孔縫隙結(jié)構(gòu)相互作用，應用偶應力理論框架下的彈性波動方程進行數(shù)值模擬，可以進一步增強地層的濾波效應.

圖7 合成地震記錄對比(z分量)(a) 時域波形對比； (b) 幅度頻譜對比.藍線為使用傳統(tǒng)彈性波動方程得到的合成地震記錄，紅線為使用偶應力理論框架下的彈性波動方程得到的合成地震記錄.Fig.7 Comparison of synthetic shot records (z component)(a) The waveform comparison in time domain; (b) The comparison of the amplitude spectrum. The blue line is the synthetic shot record using the conventional elastic wave equations, and the red line is the synthetic shot record using elastic wave equations in the frame of the couple stress theory.

3 結(jié)論

本文嘗試啟用廣義連續(xù)介質(zhì)力學理論，推導偶應力理論框架下的彈性波動方程；將其數(shù)學表達形式以及數(shù)值模擬結(jié)果與傳統(tǒng)彈性波動方程進行對比.研究探索推導的具非對稱性的波動方程所具有的理論意義和應用價值.

本文推導的偶應力理論框架下的彈性波動方程，相比傳統(tǒng)彈性波動方程，既增加了獨立的旋轉(zhuǎn)自由項，又引入了介質(zhì)特征尺度參數(shù)，描述由介質(zhì)內(nèi)部微孔縫隙結(jié)構(gòu)相互作用產(chǎn)生的不均勻性所導致的旋轉(zhuǎn)運動以及其位移擾動的傳播.旋轉(zhuǎn)自由項中的參數(shù)η是反映介質(zhì)微旋轉(zhuǎn)運動特性的參數(shù).如果η=0，則旋轉(zhuǎn)效應消失，偶應力理論框架下的彈性波動方程就和傳統(tǒng)彈性波動方程有一致的數(shù)學表達.

通過對均勻模型和Marmousi模型上的數(shù)值模擬結(jié)果對比分析，可以發(fā)現(xiàn)，考慮介質(zhì)內(nèi)部微孔縫隙相互作用所帶來的位移擾動對地震波傳播的影響是存在的，介質(zhì)微孔縫隙結(jié)構(gòu)的特征尺度，相比于子波波長，量級上的差距是巨大的，但其所帶來的對地震波傳播的影響仍然能夠在地震記錄上被觀察到.另外，偶應力理論框架下的旋轉(zhuǎn)運動，對縱波沒有影響，但卻使得橫波在傳播過程中出現(xiàn)了物理頻散.為了進一步研究分析介質(zhì)微孔縫隙相互作用對地震波傳播的影響，我們在高鐵橋墩模型上應用不同的彈性波動方程進行數(shù)值模擬試驗，并將合成地震記錄與高鐵實際記錄進行對比，結(jié)果表明，即使介質(zhì)微孔縫隙結(jié)構(gòu)的特征尺度為微米量級，合成的地震波場仍然出現(xiàn)了衰減，且其頻譜能量往低頻端移動.這表明介質(zhì)內(nèi)部的微孔縫隙結(jié)構(gòu)相互作用可以增強地層的濾波效應，同時，也說明了地震記錄對介質(zhì)微孔縫隙結(jié)構(gòu)的敏感性，可以啟發(fā)后續(xù)的介質(zhì)內(nèi)部微孔縫隙特征尺度參數(shù)反演的研究.

致謝感謝中國科學院地質(zhì)與地球物理研究所、北京大學和西安交通大學提供的數(shù)據(jù)支持，感謝中國科學院地質(zhì)與地球物理研究所提供的計算資源.