王之洋， 李幼銘 ， 白文磊

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

0 引言

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

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

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

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

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

1 偶應(yīng)力理論框架下的彈性波動(dòng)方程

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

1.1 平衡方程

(1)

(2)

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

(3)

(4)

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

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

(5)

(6)

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

(7)

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

(8)

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

(9)

(10)

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

eistσst=0.

(11)

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

而在偶應(yīng)力理論中，當(dāng)考慮偶應(yīng)力或體力偶時(shí)，切應(yīng)力互等定理不成立了，這導(dǎo)致應(yīng)力張量的不對(duì)稱性，這時(shí)候，應(yīng)力張量σij可以分解為對(duì)稱應(yīng)力張量τij與反對(duì)稱應(yīng)力張量rij之和：

(12)

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

(13)

μji,j+eistrst+Mi=0，

(14)

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

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

(15)

(16)

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

(17)

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

(18)

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

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

(19)

對(duì)公式(19)應(yīng)用散度定理，可得：

(20)

因?yàn)棣蘪i,j+eistσst+Mi=0，可得eistμst=0，即偶應(yīng)力張量μij是對(duì)稱的.

1.2 邊界條件

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

(21)

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

因?yàn)橛校?/p>

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

(22)

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

對(duì)公式(22)等號(hào)右邊第一項(xiàng)應(yīng)用Hamiltonian算子，可得：

(23)

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

(24)

則公式(21)可寫為：

(25)

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

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

(26)

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

(27)

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

1.3 幾何方程

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

dui=ui,jdxj，

(28)

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

(29)

(30)

(31)

(32)

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

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

dωi=ωi,jdxj，

(33)

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

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

(34)

1.4 本構(gòu)方程

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

(35)

(36)

對(duì)公式(35)和(36)應(yīng)用鏈?zhǔn)角髮?dǎo)法則以及散度定理，并相加，可得：

(37)

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

(38)

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

(39)

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

(40)

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

1.5 偶應(yīng)力理論框架下的彈性波動(dòng)方程

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

(41)

旋轉(zhuǎn)自由項(xiàng)中的η是介質(zhì)內(nèi)部微結(jié)構(gòu)相互作用導(dǎo)致的旋轉(zhuǎn)效應(yīng)的表征，如果η=0，則旋轉(zhuǎn)效應(yīng)消失，偶應(yīng)力理論框架下的彈性波動(dòng)方程就和傳統(tǒng)彈性波動(dòng)方程有一致的數(shù)學(xué)表達(dá).Teisseyre和Boratńy(tǒng)ski(2003)通過理論推導(dǎo)證明，在顆粒介質(zhì)或者有裂隙的連續(xù)體內(nèi)，由于微缺陷/微結(jié)構(gòu)的存在，應(yīng)力或者應(yīng)力場會(huì)出現(xiàn)不對(duì)稱性力學(xué)特征，由此產(chǎn)生一種新的旋轉(zhuǎn)運(yùn)動(dòng).Ba?ant(2002)認(rèn)為金屬材料的微缺陷/微結(jié)構(gòu)一般由晶體的位錯(cuò)引起，晶體位錯(cuò)在納米量級(jí)；而巖土介質(zhì)的微缺陷/微結(jié)構(gòu)一般是由微夾雜/微裂縫/微孔隙/孔洞等引起，尺度可以在微毫米或更高量級(jí).黃文雄和徐可(2014)通過實(shí)驗(yàn)證明，介質(zhì)特征尺度l與組成介質(zhì)的材料顆粒平均直徑以及幾何形狀有關(guān)，同時(shí)，也與變形局部化現(xiàn)象有關(guān).另外，與應(yīng)力集中效應(yīng)也有關(guān)(鮑亦興和毛昭宙，1993).如果忽略一切附加因素，介質(zhì)特征尺度參數(shù)可直接視為介質(zhì)微孔縫隙特征尺度參數(shù).相反，如果考慮組成介質(zhì)的材料顆粒的幾何特征、變形局部化現(xiàn)象以及應(yīng)力集中效應(yīng)，介質(zhì)特征尺度參數(shù)可視為介質(zhì)微孔縫隙特征尺度參數(shù)與系數(shù)ζ相乘的結(jié)果，系數(shù)ζ可通過實(shí)驗(yàn)獲得，王之洋等(2020)通過對(duì)高鐵實(shí)際數(shù)據(jù)和理論合成數(shù)據(jù)的對(duì)比分析，也提取了該系數(shù).當(dāng)然，為了進(jìn)一步研究介質(zhì)特征尺度參數(shù)和旋轉(zhuǎn)自由項(xiàng)的物理意義，必須結(jié)合巖石物理測量，開展對(duì)比驗(yàn)證研究和技術(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ā)，時(shí)間步長Δt=0.001 s，介質(zhì)微孔縫隙特征尺度參數(shù)為200 μm.為了排除數(shù)值頻散和邊界反射對(duì)數(shù)值模擬結(jié)果的影響，我們采用高階優(yōu)化有限差分算子以及優(yōu)化的PML(Perfectly Matched Layer，完全匹配層)邊界條件進(jìn)行數(shù)值模擬，以更好地對(duì)傳統(tǒng)彈性波動(dòng)方程以及偶應(yīng)力理論框架下的彈性波動(dòng)方程的合成記錄進(jìn)行對(duì)比分析.

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

圖1 應(yīng)用不同彈性波動(dòng)方程得到的合成波場快照(x分量)(a) 傳統(tǒng)彈性波動(dòng)方程； (b) 偶應(yīng)力理論框架下的彈性波動(dòng)方程； (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 應(yīng)用不同彈性波動(dòng)方程得到的合成波場快照(z分量)(a) 傳統(tǒng)彈性波動(dòng)方程； (b) 偶應(yīng)力理論框架下的彈性波動(dòng)方程； (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).

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

2.2 Marmousi模型

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

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

圖3 在Marmousi模型上應(yīng)用不同彈性波動(dòng)方程得到的合成地震記錄(x分量)(a) 傳統(tǒng)彈性波動(dòng)方程； (b) 偶應(yīng)力理論框架下的彈性波動(dòng)方程； (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模型上應(yīng)用不同彈性波動(dòng)方程得到的合成地震記錄(z分量)(a) 傳統(tǒng)彈性波動(dòng)方程； (b) 偶應(yīng)力理論框架下的彈性波動(dòng)方程； (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 高鐵橋墩模型

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

我們構(gòu)建了一個(gè)基于高鐵列車震源的簡化橋墩模型.速度模型大小為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個(gè)，橋墩間距為28 m，每個(gè)橋墩插入地下50 m深度直達(dá)圍巖，每個(gè)橋墩的地下部分被看作是間隔為10 m的5個(gè)點(diǎn)源，地震波在橋墩中的傳播速度為vp=4000 m·s-1.高鐵列車車廂為16節(jié)，每節(jié)車廂的長度為28 m，列車速度為300 km·h-1，視為多個(gè)移動(dòng)Ricker震源的疊加，Ricker子波主頻設(shè)置為20 Hz，時(shí)間步長Δt=0.0005 s，記錄時(shí)長為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是不同彈性波動(dòng)方程得到的合成地震記錄(z分量).從圖6可以觀察到，應(yīng)用傳統(tǒng)彈性波動(dòng)方程和偶應(yīng)力理論框架下的彈性波動(dòng)方程所得到的合成地震記錄在時(shí)域上的波形相似，我們的合成地震記錄可視為不同空間位置的以一定時(shí)間間隔激發(fā)的固定震源記錄的線性疊加，基于線性疊加的特性，可以通過設(shè)置模型參數(shù)，使模型的橋墩和橋梁長度和實(shí)際情況更加接近，從而使時(shí)域合成地震記錄與實(shí)際記錄具有更好的一致性.通過對(duì)頻譜的分析，我們發(fā)現(xiàn)，當(dāng)設(shè)置高速列車車速為300 km·h-1時(shí)，合成地震記錄的能量集中在大約10 Hz以及25 Hz左右，其中，頻率10 Hz附近的能量最大，這與實(shí)際高鐵地震記錄的頻譜能量分布基本相符.

圖6 不同彈性波動(dòng)方程得到的合成地震記錄(z分量)(a) 傳統(tǒng)彈性波動(dòng)方程； (b) 偶應(yīng)力理論框架下的彈性波動(dòng)方程； 圖(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.

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

圖7 合成地震記錄對(duì)比(z分量)(a) 時(shí)域波形對(duì)比； (b) 幅度頻譜對(duì)比.藍(lán)線為使用傳統(tǒng)彈性波動(dòng)方程得到的合成地震記錄，紅線為使用偶應(yīng)力理論框架下的彈性波動(dò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ì)力學(xué)理論，推導(dǎo)偶應(yīng)力理論框架下的彈性波動(dòng)方程；將其數(shù)學(xué)表達(dá)形式以及數(shù)值模擬結(jié)果與傳統(tǒng)彈性波動(dòng)方程進(jìn)行對(duì)比.研究探索推導(dǎo)的具非對(duì)稱性的波動(dòng)方程所具有的理論意義和應(yīng)用價(jià)值.

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

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

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