水庫(kù)蓄放水對(duì)庫(kù)底巖石介質(zhì)彈性波速影響的數(shù)值模擬

周 斌 孫 峰 薛世峰 蔣海昆 史水平

1)廣西壯族自治區(qū)地震局，南寧 530022

2)中國(guó)石油大學(xué)(華東)，青島 266580

3)中國(guó)地震臺(tái)網(wǎng)中心，北京 100045

0 引言

彈性波速是組成地下介質(zhì)的巖石性質(zhì)、巖體結(jié)構(gòu)及其賦存熱應(yīng)力環(huán)境的綜合反映。盡管影響巖石介質(zhì)彈性波速的因素較多，如巖性、孔隙幾何形態(tài)、孔隙度、膠結(jié)和壓實(shí)情況、原地應(yīng)力、溫度、孔隙流體壓力和飽和度等，但有效應(yīng)力是最重要的影響因素之一(葛洪魁等，2001)，在一定地質(zhì)條件下，各種尺度介質(zhì)的彈性波速都與有效應(yīng)力存在著依存關(guān)系(尹尚先等，2003)。在水庫(kù)蓄、放水過(guò)程中，由于地表水體的加、卸載作用及孔隙壓力的擴(kuò)散作用，在一定程度上改變了水體下方巖石介質(zhì)賦存的應(yīng)力環(huán)境，可能會(huì)導(dǎo)致巖石介質(zhì)彈性波速發(fā)生響應(yīng)性的變化。大量觀測(cè)事實(shí)(馮德益等，1993;常寶琦，1997;陳俊華等，2007;盧顯等，2013)已證實(shí)了這種波速異常變化的存在，陳蒙等(2013)的研究還發(fā)現(xiàn)，云南省賓川縣大銀甸水庫(kù)周圍地區(qū)地下波速變化與水庫(kù)水位變化之間存在著較好的一致性。那么，在水庫(kù)地表水體加、卸載及滲透過(guò)程中，庫(kù)底巖石介質(zhì)有效應(yīng)力的改變到底對(duì)彈性波速的影響程度有多大?影響范圍有多廣?彈性波速變化與孔隙壓力、有效應(yīng)力變化之間到底存在何種定量關(guān)系?對(duì)于這些問(wèn)題，目前國(guó)內(nèi)外尚少見(jiàn)相關(guān)的報(bào)道。本文以四川省紫坪鋪水庫(kù)為例，在前期建立該水庫(kù)及鄰近地區(qū)地質(zhì)構(gòu)造與水文地質(zhì)結(jié)構(gòu)模型，并推導(dǎo)出巖石介質(zhì)變形與流體滲流耦合數(shù)學(xué)模型的基礎(chǔ)上(周斌，2010;周斌等，2010)，進(jìn)一步建立了有效應(yīng)力對(duì)巖石彈性波速影響的數(shù)學(xué)模型，利用有限元方法定量計(jì)算了庫(kù)底巖石介質(zhì)的孔隙壓力、附加有效應(yīng)力和彈性波速的變化，討論了彈性波速隨孔隙壓力、附加有效應(yīng)力變化的規(guī)律。此研究是一項(xiàng)嘗試性的工作，對(duì)進(jìn)一步探討水庫(kù)誘發(fā)地震成因機(jī)制及水庫(kù)地震可能的前兆具有一定的借鑒意義。

1 有效應(yīng)力對(duì)彈性波速影響的數(shù)學(xué)模型

1.1 巖石彈性波速隨有效應(yīng)力的變化規(guī)律

在均勻、各向同性的理想彈性體中，往往存在著2種彈性波，其波速分別為

式(1)中:VP，VS分別代表巖石彈性縱、橫波速度，K，G分別為巖石體積模量和剪切模量，ρ為巖石密度。

巖石受到載荷作用時(shí)圍壓上升，巖石內(nèi)部孔隙、裂隙的變形使巖石體積模量和剪切模量均增大，波速隨之上升;孔隙流體壓力的作用與圍壓相反，抵抗巖石孔隙、裂隙的變形。根據(jù)巖石孔彈性理論(Biot，1956)，在外部荷載和孔隙流體壓力共同作用下巖石介質(zhì)的宏觀力學(xué)響應(yīng)由有效應(yīng)力來(lái)控制(壓為正)，即

式(2)中:σij和σ'ij分別為全應(yīng)力張量和有效應(yīng)力張量，P為孔隙流體壓力，δij為Kroneker符號(hào)，α為有效應(yīng)力系數(shù)，或稱為Biot系數(shù)，其大小取決于巖石的孔隙、裂隙發(fā)育程度。

在巖石介質(zhì)彈性波速度影響因素的研究中，其所處的應(yīng)力狀態(tài)一直作為主要的影響因素而受到關(guān)注，當(dāng)前室內(nèi)試驗(yàn)研究的主要認(rèn)識(shí)包括:隨著有效應(yīng)力的增大，VP、VS均增大，低應(yīng)力下彈性波速隨有效圍壓增大迅速，增速隨有效圍壓的增大而減小，在高應(yīng)力下趨于一常數(shù)值;孔隙幾何形態(tài)對(duì)巖石波速的應(yīng)力敏感性有重要影響，裂隙類巖石的彈性波速比孔隙類巖石隨應(yīng)力增大而增加迅速;巖石彈性波速敏感性與孔隙流體類型、孔隙幾何形態(tài)和巖石的構(gòu)造位置有關(guān);水飽和可大大增加巖石的縱波波速，但對(duì)橫波波速影響較小。

1.2 有效應(yīng)力對(duì)巖石彈性波速影響的數(shù)學(xué)模型

為了方便應(yīng)用，需要建立有效應(yīng)力對(duì)巖石彈性波速影響的數(shù)學(xué)模型，波速模型應(yīng)能很好地描述波速隨應(yīng)力的變化規(guī)律，同時(shí)有明確的物理含義，并不失一般性。表1列出有關(guān)文獻(xiàn)提出的有效應(yīng)力對(duì)波速影響的模型。

由前面分析可知，在有效應(yīng)力作用下，巖石波速隨有效應(yīng)力的增大而增大，最后趨于定值。表1中的模型①限定采用對(duì)數(shù)函數(shù)進(jìn)行擬合;依據(jù)模型②—④，高圍壓下巖石波速隨應(yīng)力的增大呈近似線性增大，而非趨于常數(shù);模型⑤比較復(fù)雜，各參數(shù)對(duì)波速的影響不獨(dú)立;模型⑥中，常數(shù)a代表零圍壓下的波速，但是，即使同一巖性的巖石，零圍壓下的波速也非常離散，遠(yuǎn)非常數(shù)?；谝陨戏治?，選擇⑦、⑧式作為彈性波速－有效應(yīng)力數(shù)學(xué)模型的基本形式，兩式可統(tǒng)一表達(dá)為

表1 現(xiàn)有的波速－應(yīng)力模型Table1 Existing velocity-stress relationships

式(3)中:A，B，C為回歸系數(shù)，σveff為有效圍壓。當(dāng)σveff→∞ 時(shí)，V=A，表明A代表高圍壓下的巖石波速;當(dāng)σveff→0時(shí)，V=A－B，表明A－B為有效圍壓為零時(shí)的波速，B值代表有效圍壓引起的波速變化量，C值代表波速隨有效圍壓增加的梯度。

將式(3)在σv0處泰勒展開(kāi)并略去高階項(xiàng)可得:

式(4)中:σv0為初始狀態(tài)的有效應(yīng)力，單位MPa，B，C為回歸系數(shù)。

至此就建立了有效應(yīng)力變化量與彈性波速度變化量直接聯(lián)系的數(shù)學(xué)模型。

2 數(shù)值模擬及結(jié)果分析

2.1 有限元模型、參數(shù)及邊界條件

紫坪鋪水庫(kù)是岷江上游的一座大型水利樞紐，位于四川省都江堰市麻溪鄉(xiāng)。在地質(zhì)構(gòu)造上，該水庫(kù)位于龍門山造山帶中段的前緣拆離帶內(nèi)。龍門山地區(qū)深部結(jié)構(gòu)具層圈性，最重要的滑脫界面有15～20km深處的低速層、寒武系泥質(zhì)粉砂巖層、志留系千枚巖和泥頁(yè)巖層以及中下三疊統(tǒng)富膏鹽巖層，這一構(gòu)造特征決定了龍門山造山帶深部的多級(jí)滑脫和淺層的多層次推覆(湯良杰等，2008)。庫(kù)區(qū)及鄰近地區(qū)主要發(fā)育了茂縣－汶川斷裂(WMF)、北川－映秀斷裂(BYF)、通濟(jì)場(chǎng)斷裂(TJCF)、安縣－灌縣斷裂(AGF)和廣元－大邑?cái)嗔?GDF)5條主干斷裂。北川－映秀斷裂以東的沉積蓋層構(gòu)造在垂向上被三疊系雷口坡組和嘉陵江組主滑脫面隔開(kāi)，其下屬于印支期前沉積的地層，主要發(fā)育早期的正斷層;后者為印支期間及以后沉積的地層，構(gòu)造變形樣式更為復(fù)雜，除發(fā)育一系列傾向NW的疊瓦狀逆沖推覆構(gòu)造外，還發(fā)育了雙重構(gòu)造、反沖斷層、三角帶、滑脫構(gòu)造等多種構(gòu)造樣式(見(jiàn)周斌等，2010中圖1，2)。地層可劃分為基底和蓋層兩大部分，基底地層主要為一套海底噴發(fā)火山巖建造和淺海泥頁(yè)巖、碎屑巖建造，蓋層主要發(fā)育震旦系—中三疊統(tǒng)海相碎屑巖及碳酸鹽巖序列和上三疊統(tǒng)—第四系陸相碎屑巖序列(金文正等，2008)。

考察水庫(kù)地表水體加、卸載及滲透過(guò)程中庫(kù)底巖石介質(zhì)有效應(yīng)力與彈性波速的關(guān)系，本應(yīng)考慮實(shí)際的地形及地質(zhì)構(gòu)造情況建立三維模型。但由于受研究區(qū)當(dāng)前地質(zhì)構(gòu)造與水文地質(zhì)結(jié)構(gòu)的研究程度及我們所掌握資料的限制，本文將其簡(jiǎn)化成平面應(yīng)變問(wèn)題。利用前期工作建立的紫坪鋪水庫(kù)庫(kù)區(qū)及鄰近地區(qū)地質(zhì)構(gòu)造與水文地質(zhì)結(jié)構(gòu)模型(周斌，2010;周斌等，2010)，進(jìn)行了有限單元?jiǎng)澐?，?2513個(gè)三角形單元，6434個(gè)節(jié)點(diǎn)，模型的底邊界設(shè)置為水平滑動(dòng)，為了避免剛體位移，兩端點(diǎn)固定。為考察水庫(kù)蓄、放水過(guò)程中，地表水體下方巖石介質(zhì)的孔隙壓力、附加有效應(yīng)力(附加平均有效應(yīng)力)變化及其與彈性波速變化的關(guān)系，分別在不同深度、不同構(gòu)造部位設(shè)置了6個(gè)觀測(cè)點(diǎn)(圖1)。

圖1 有限單元?jiǎng)澐旨坝^察點(diǎn)位置示意圖Fig.1 The sketch map of finite element division and location of observation points.

各地質(zhì)體力學(xué)參數(shù)與物性參數(shù)采用了周斌等(2010)文獻(xiàn)中的取值。為獲得彈性波速與有效應(yīng)力關(guān)系的回歸系數(shù)，利用前人文獻(xiàn)(Nur et al.，1969;Han，et al.，1986;劉永祥等，1995;Prasad et al.，1997;葛洪魁等，2001;馬中高等，2006)中的實(shí)驗(yàn)數(shù)據(jù)，并結(jié)合中國(guó)石化股份有限公司南方勘探開(kāi)發(fā)分公司在川東北地區(qū)的大量測(cè)試數(shù)據(jù)，分砂巖、泥(頁(yè))巖、石灰?guī)r、火成巖4類進(jìn)行了統(tǒng)計(jì)回歸，并按照巖性相近、賦值相同的原則，大致確定了不同地質(zhì)體的回歸系數(shù) B，C(表2)。

依據(jù)2005年3月31日至2008年4月14日期間紫坪鋪水庫(kù)壩前水位變化資料，將相對(duì)于2005年3月31日基準(zhǔn)水位(海拔高度約754.5m)的水位變化用分段直線加載函數(shù)在庫(kù)體處施加動(dòng)水頭邊界和動(dòng)荷載邊界(圖2)，計(jì)算時(shí)共設(shè)置了38個(gè)計(jì)算步，每步歷時(shí)30天。

2.2 數(shù)值模擬結(jié)果分析

2.2.1 附加有效應(yīng)力變化

圖3為紫坪鋪水庫(kù)地表水體加、卸載及滲透過(guò)程中，不同計(jì)算時(shí)步庫(kù)底附加有效應(yīng)力變化的云圖。由該圖可以看出，庫(kù)底正下方黑色虛線所圍的區(qū)域主要受庫(kù)體荷載作用的影響，附加有效應(yīng)力呈壓應(yīng)力狀態(tài)，在2005年3月31日水位開(kāi)始抬升后的第570天(2006年10月22日)和990天(2007年12月16日)，水庫(kù)水位分別相對(duì)于基準(zhǔn)水位抬升了約120m和118.5m，附加有效應(yīng)力也分別達(dá)到最大變幅0.59MPa和0.51MPa。圖中紅色虛線所圍的區(qū)域主要發(fā)育泥盆系—下三疊統(tǒng)飛仙關(guān)組海相碳酸鹽巖地層，巖性堅(jiān)硬，裂隙相對(duì)發(fā)育，滲透條件好，且可以通過(guò)通濟(jì)場(chǎng)斷裂與庫(kù)底直接連通。雖然該區(qū)域距離庫(kù)底較遠(yuǎn)，但隨著蓄水時(shí)間的延長(zhǎng)和附加水頭壓力的擴(kuò)散，附加有效應(yīng)力逐漸轉(zhuǎn)變?yōu)閺垜?yīng)力狀態(tài)，最大變幅達(dá)到－0.5MPa左右。

表2 不同地質(zhì)體彈性波速－有效應(yīng)力關(guān)系的回歸系數(shù)Table2 The velocity-stress regression coefficient of geologic bodies

圖2 動(dòng)水頭邊界和動(dòng)荷載邊界加載方式Fig.2 Dynamic water level and loading boundary.

2.2.2 彈性波速變化

圖3 附加有效應(yīng)力變化云圖Fig.3 The cloud map of additional effective stress.

圖4，5分別為不同計(jì)算時(shí)步彈性縱波波速增量值ΔVP和彈性橫波波速增量值ΔVS變化的云圖。由圖4可以看出，在地表水體加、卸載及滲透過(guò)程中，由附加有效應(yīng)力變化造成ΔVP變化的范圍，主要分布在地下5km以上的區(qū)域(圖4中紅色虛線以上的區(qū)域)，特別是圖4中黑色虛線所圍的區(qū)域變化較大。在2005年3月31日水位開(kāi)始抬升后的第570天(2006年10月22日)，附加有效應(yīng)力最大變幅達(dá)到±0.51MPa左右，ΔVP的最大變幅達(dá)到±0.013km/s左右。隨著水位的不斷抬升及蓄水時(shí)間的延長(zhǎng)，距離庫(kù)底較遠(yuǎn)的泥盆系—下三疊統(tǒng)飛仙關(guān)組海相碳酸鹽巖地層中(圖4中紅色虛線所圍的區(qū)域)，受地表水體附加水頭壓力擴(kuò)散作用的影響，附加有效應(yīng)力逐漸轉(zhuǎn)變?yōu)閺垜?yīng)力狀態(tài)，ΔVP顯著降低。由圖5可以看出，ΔVS變化與ΔVP變化的圖像類似，其區(qū)別在于ΔVS變化的最大幅度為±0.01km/s，小于ΔVP的最大變幅。

2.2.3 彈性波速與孔隙壓力、附加有效應(yīng)力的關(guān)系

圖6為各觀測(cè)點(diǎn)孔隙壓力、附加有效應(yīng)力和彈性波速增量值變化的曲線，其中S1～S6中的綠色曲線表示庫(kù)水位變化(相對(duì)于2005年3月31日基準(zhǔn)水位)在庫(kù)體處施加的動(dòng)荷載，也大體表示了庫(kù)水位變化的過(guò)程。

觀測(cè)點(diǎn)S1(埋深約1.5km)位于北川－映秀斷裂與通濟(jì)場(chǎng)斷裂間的近地表逆沖推覆體中，巖性為三疊系須家河組陸相碎屑巖。由于該觀測(cè)點(diǎn)巖體滲透性較差，且與水庫(kù)無(wú)直接的水力聯(lián)系，故孔隙壓力變化幅度很小，在蓄水1140天后(第38個(gè)計(jì)算步，對(duì)應(yīng)2008年4月14日)僅達(dá)到－0.016MPa。附加有效應(yīng)力主要受控于庫(kù)水位動(dòng)荷載作用，并表現(xiàn)出與庫(kù)水位變化的一致性，在蓄水990天后(第33個(gè)計(jì)算步，對(duì)應(yīng)2007年12月16日)達(dá)到最大值0.407MPa。彈性波速增量值的變化與附加有效應(yīng)力變化表現(xiàn)出較好的一致性，ΔVP的變化幅度高于ΔVS，在蓄水990天后分別達(dá)到最大值0.0084km/s和0.0064km/s。

圖4 彈性縱波波速增量值變化云圖Fig.4 The cloud map of increased elastic P-wave velocity.

圖5 彈性橫波波速增量值變化云圖Fig.5 The cloud map of increased elastic S-wave velocity.

圖6 各觀測(cè)點(diǎn)孔隙壓力、附加有效應(yīng)力及彈性波速增量變化曲線Fig.6 The change of pore pressure，additional effective stress and increased elastic wave velocity at each observation point.

觀測(cè)點(diǎn)S2(埋深約1km)位于庫(kù)體正下方通濟(jì)場(chǎng)斷裂滑脫面以上的碎屑巖中，巖體破碎程度較高、次級(jí)斷裂發(fā)育且與水庫(kù)有直接的水力聯(lián)系。隨著地表水體附加水頭壓力的擴(kuò)散，該觀測(cè)點(diǎn)孔隙壓力持續(xù)升高，在蓄水1140天后達(dá)到最大變幅－0.59MPa。附加有效應(yīng)力變化與庫(kù)水位變化大體相似，但由于受孔隙壓力持續(xù)升高的影響，附加有效應(yīng)力曲線呈趨勢(shì)性的下降，在蓄水570天后(第19個(gè)計(jì)算步，對(duì)應(yīng)2006年10月22日)上升到最大值0.50MPa，蓄水1140天后降低到最小值－0.067 MPa，表明該觀測(cè)點(diǎn)附加有效應(yīng)力變化同時(shí)受控于庫(kù)水位動(dòng)荷載作用和孔隙壓力擴(kuò)散作用。彈性波速變化與附加有效應(yīng)力變化同步，ΔVP的變化幅度稍高于ΔVS，在蓄水570天后分別上升到最大值0.0120km/s和0.0089km/s，蓄水1140天后分別降低到最小值－0.0016km/s和－0.0012km/s。

觀測(cè)點(diǎn)S3(埋深約4.5km)位于通濟(jì)場(chǎng)斷裂帶上，滲透條件好，與水庫(kù)有直接的水力聯(lián)系。隨著地表水體沿?cái)嗔淹ǖ老蛏畈康臐B流，該觀測(cè)點(diǎn)的孔隙壓力響應(yīng)略滯后于地表水體附加水頭壓力變化，在蓄水1020天后(第34個(gè)計(jì)算步，對(duì)應(yīng)2008年1月15日)達(dá)到最大變幅－0.72MPa。由于該觀測(cè)點(diǎn)距離庫(kù)底較遠(yuǎn)，且不在庫(kù)底正下方，受庫(kù)水位動(dòng)荷載作用的影響較弱，因此附加有效應(yīng)力變化主要受控于孔隙壓力擴(kuò)散作用，隨孔隙壓力的持續(xù)升高而同步降低，在蓄水1020天后達(dá)到最小值－0.32MPa。ΔVP和ΔVS隨附加有效應(yīng)力的降低而同步降低，ΔVP降低的幅度稍高于 ΔVS，在蓄水 1020天后分別降低到最小值－0.00070 km/s和－0.00068km/s。

觀測(cè)點(diǎn)S4(埋深約4.5km)位于庫(kù)體正下方安縣－灌縣斷裂上盤的三疊系須家河組陸相碎屑巖地層中，巖性軟弱，滲透條件差。該觀測(cè)點(diǎn)與觀測(cè)點(diǎn)S1的孔隙壓力、附加有效應(yīng)力和彈性波速增量值曲線的變化形態(tài)相似，區(qū)別是該觀測(cè)點(diǎn)彈性波速增量值的變化幅度更小，ΔVP和ΔVS的最大值僅為0.00099km/s和0.00095km/s。

觀測(cè)點(diǎn)S5(埋深約9km)位于安縣－灌縣斷裂與通濟(jì)場(chǎng)斷裂在地腹的交會(huì)處的下方，主要發(fā)育泥盆系—下三疊統(tǒng)飛仙關(guān)組海相碳酸鹽巖，巖性堅(jiān)硬，裂隙相對(duì)發(fā)育，滲透條件好。該觀測(cè)點(diǎn)距離庫(kù)底較遠(yuǎn)，在庫(kù)體荷載作用下產(chǎn)生的彈性附加應(yīng)力場(chǎng)較弱，附加水頭壓力擴(kuò)散作用對(duì)附加有效應(yīng)力變化起主導(dǎo)性的作用，孔隙壓力、附加有效應(yīng)力和彈性波速增量值曲線的變化形態(tài)與觀測(cè)點(diǎn)S3相似，區(qū)別是該觀測(cè)點(diǎn)ΔVP和ΔVS的變化幅度比觀測(cè)點(diǎn)S3降低了1個(gè)數(shù)量級(jí)。

觀測(cè)點(diǎn)S6位于庫(kù)體正下方約10km處的泥盆系—下三疊統(tǒng)飛仙關(guān)組海相碳酸鹽巖中。該觀測(cè)點(diǎn)孔隙壓力、附加有效應(yīng)力和彈性波速變化曲線形態(tài)與觀測(cè)點(diǎn)S2相似，區(qū)別是該觀測(cè)點(diǎn)ΔVP和ΔVS的變化幅度比觀測(cè)點(diǎn)S2至少降低了2個(gè)數(shù)量級(jí)。

3 討論與結(jié)論

(1)本文以紫坪鋪水庫(kù)為例，利用有限元方法數(shù)值模擬了水庫(kù)水體加、卸載及滲透過(guò)程中，庫(kù)底巖石介質(zhì)孔隙壓力、附加有效應(yīng)力和彈性波速增量值的變化過(guò)程。模擬結(jié)果顯示:1)隨觀測(cè)點(diǎn)埋藏深度的不同及與水庫(kù)水力聯(lián)系強(qiáng)弱的差異，孔隙壓力對(duì)地表水體附加水頭壓力變化的響應(yīng)表現(xiàn)出不同的特征。在與水庫(kù)有直接水力聯(lián)系的淺部地層中(如S2，S3觀測(cè)點(diǎn))，孔隙壓力響應(yīng)迅速、變化幅度高;在與水庫(kù)無(wú)明顯水力聯(lián)系的地層或有直接水力聯(lián)系的深部地層中(如S1，S4，S5，S6觀測(cè)點(diǎn))，孔隙壓力響應(yīng)不明顯;孔隙壓力響應(yīng)的另一個(gè)特征是滯后于附加水頭壓力的變化，距離庫(kù)底越遠(yuǎn)，滯后現(xiàn)象就越為突出。2)在與水庫(kù)無(wú)直接水力聯(lián)系且滲透性較差的地層中(如S1，S4觀測(cè)點(diǎn))，附加有效應(yīng)力主要受庫(kù)體動(dòng)荷載作用的影響;在與水庫(kù)有直接水力聯(lián)系的斷裂帶上或其附近滲透性較好的地層中，附加有效應(yīng)力同時(shí)受控于庫(kù)體動(dòng)荷載作用和附加水頭壓力擴(kuò)散作用(如S2，S3，S5觀測(cè)點(diǎn))。3)彈性波速增量值變化同步于附加有效應(yīng)力變化，并表現(xiàn)出與庫(kù)水位變化形態(tài)的相似性;不同觀測(cè)點(diǎn)ΔVP和ΔVS的大小與埋藏位置、深度及附加有效應(yīng)力的變化幅度有關(guān);同一觀測(cè)點(diǎn)，ΔVP變化的幅度大于ΔVS。

(2)室內(nèi)巖石力學(xué)測(cè)試表明，隨著有效應(yīng)力的增大，VP、VS均增大，低應(yīng)力下彈性波速隨有效圍壓增大迅速，增速隨有效圍壓的增大而減小，在高應(yīng)力下趨于一常數(shù)值。由式(4)也可以看出，在有效圍壓變化不大時(shí)，彈性波速增量值隨地層初始圍壓的升高呈指數(shù)衰減。下面用一個(gè)簡(jiǎn)單的模型作進(jìn)一步的分析。若假設(shè)庫(kù)底為均質(zhì)、各向同性的巖石，密度為2.3g/cm3，彈性波速－有效應(yīng)力關(guān)系的回歸系數(shù)B為1.5，C為0.3，且將庫(kù)底附加有效應(yīng)力變化設(shè)定為2種理想情況:1)庫(kù)底以下15km范圍內(nèi)的附加有效應(yīng)力變化均為0.5MPa，2)庫(kù)底以下15km范圍內(nèi)的附加有效應(yīng)力變化均為－0.5MPa，利用式(4)計(jì)算得到庫(kù)底彈性波速增量值隨深度變化的曲線(圖7)。可以看出，無(wú)論附加有效應(yīng)力的變化量是0.5MPa，還是－0.5MPa，隨深度的增加彈性波速增量值都表現(xiàn)出快速的指數(shù)衰減，在5km以下彈性波速增量值的變化已經(jīng)很小。數(shù)值模擬得到同樣的結(jié)果，在水庫(kù)水體加、卸載及滲透過(guò)程中，地下巖石介質(zhì)彈性波速增量值的變化并不明顯，最大變幅僅為±0.013km/s，變化范圍主要集中在庫(kù)底5km以上的局部區(qū)域。

(3)利用地震臺(tái)陣(臺(tái)網(wǎng))觀測(cè)資料反演地殼速度結(jié)構(gòu)，是監(jiān)測(cè)地下介質(zhì)屬性和應(yīng)力狀態(tài)變化最直接的方法。2008年5月12日汶川8.0級(jí)地震后，一些研究者(Chen et al.，2010;劉志坤等，2010;趙盼盼等，2012)利用四川數(shù)字地震臺(tái)網(wǎng)和川西流動(dòng)地震觀測(cè)臺(tái)陣的寬頻帶連續(xù)波形資料，通過(guò)環(huán)境噪聲方法研究了汶川地震前后地殼速度的變化特征，劉志坤等(2010)對(duì)橫跨紫坪鋪水庫(kù)的MXI－YZP臺(tái)站對(duì)相對(duì)波速變化的反演結(jié)果及Chen等(2010)對(duì)鄰近庫(kù)區(qū)的龍門山斷裂帶相對(duì)波速變化的反演結(jié)果，在汶川地震之前的時(shí)段與本文數(shù)值模擬的結(jié)果在曲線變化形態(tài)上大體一致，但變化幅度略有差異。周龍泉(2009)依據(jù)紫坪鋪水庫(kù)數(shù)字地震監(jiān)測(cè)臺(tái)網(wǎng)記錄資料，利用震源位置和速度結(jié)構(gòu)聯(lián)合反演的方法獲得了水庫(kù)蓄水后庫(kù)區(qū)彈性波速變化的圖像，結(jié)果顯示波速異常主要集中在庫(kù)底5km以上的區(qū)域，與本文數(shù)值模擬的結(jié)果一致，但異常幅度比本文結(jié)果高出1個(gè)數(shù)量級(jí)。造成以上反演結(jié)果與本文數(shù)值模擬結(jié)果差異的原因可能有3個(gè):1)由于紫坪鋪水庫(kù)距2008年5月12日汶川8.0級(jí)地震震中的距離大約6km，因此庫(kù)底巖石介質(zhì)的彈性波速不但受庫(kù)水加、卸載及滲透過(guò)程中附加有效應(yīng)力變化的影響，還可能與汶川8.0級(jí)地震成核過(guò)程中應(yīng)力場(chǎng)的變化有關(guān)，而本文的數(shù)值模擬工作僅考慮了前一方面。2)在不同的地應(yīng)力環(huán)境及埋深條件下，初始狀態(tài)的有效應(yīng)力σv0存在著差異。由于紫坪鋪庫(kù)區(qū)處于以擠壓為主的構(gòu)造應(yīng)力環(huán)境中，依據(jù)斷層形成機(jī)制的Anderson理論，大地靜力場(chǎng)的垂向分量(ρgh)應(yīng)小于σv0，本文在數(shù)值模擬時(shí)采用了大地靜力場(chǎng)的垂向分量替代σv0進(jìn)行估算，這在一定程度上會(huì)減弱有效應(yīng)力對(duì)巖石彈性波速的影響。3)在水庫(kù)水體加、卸載及滲透過(guò)程中，庫(kù)底巖石的孔隙幾何形態(tài)、孔隙度及孔隙流體飽和度等可能會(huì)發(fā)生變化，從而對(duì)巖石介質(zhì)彈性波速產(chǎn)生影響，本文在數(shù)值模擬時(shí)僅考慮了有效應(yīng)力單一因素的影響，可能會(huì)造成計(jì)算結(jié)果與真實(shí)情況的誤差。

(4)在擁有更多實(shí)驗(yàn)數(shù)據(jù)的基礎(chǔ)上建立更為合理的統(tǒng)計(jì)回歸公式，是利用數(shù)值模擬方法探索庫(kù)水加、卸載及滲透過(guò)程中庫(kù)區(qū)地球物理場(chǎng)變化規(guī)律的關(guān)鍵。但由于實(shí)驗(yàn)設(shè)備、研究目的和區(qū)域的限制，不同學(xué)者在實(shí)驗(yàn)條件、研究?jī)?nèi)容及巖石類型等方面存在很太差別。因此，系統(tǒng)而廣泛地收集整理、綜合這些資料，考慮巖石與巖體的區(qū)別，是進(jìn)一步研究中應(yīng)開(kāi)展的工作。

圖7 彈性波速增量值隨深度變化的曲線Fig.7 Change of elastic wave velocity with depths.

常寶琦，沈立英，梁紀(jì)彬.1997.水庫(kù)地震主震發(fā)生時(shí)間的預(yù)測(cè)［J］.華南地震，17(2):76—82.

CHANG Bao-qi，SHEN Li-ying，LIANG Ji-bin.1997.Prediction about the main shock origin time of reservoir earthquake［J］.South China Journal of Seismology，17(2):76—82(in Chinese).

陳俊華，甘家思，李普麗，等.2007.三峽水庫(kù)蓄水后巴東地區(qū)波速比(VP/VS)研究［J］.大地測(cè)量與地球動(dòng)力學(xué)，27(?？?:12—14，42.

CHEN Jun-hua，GAN Jia-si，LI Pu-li，et al.2007.On wave velocity ratio(VP/VS)in Badong area after impoundment of Three Gorges Reservoir［J］.Journal of Geodesy and Geodynamics，27(Suppl):12—14，42(in Chinese).

陳蒙，Gregor H，王寶善，等.2013.基于地震背景噪聲的與水庫(kù)水位變化相關(guān)的地下波速變化的研究［A］.見(jiàn):中國(guó)地球物理學(xué)會(huì)編.中國(guó)地球物理2013.419—420.

CHEN Meng，Gregor H，WANG Bao-shan，et al.2013.Study on underground rock wave velocity changes associated with reservoir water level based on the seismic background noise［A］.In:China Geophysical Society(ed).Geophysics in China 2013.419—420(in Chinese).

馮德益，虞雪君，盛國(guó)英.1993.波速異常的進(jìn)一步研究與討論(三):水庫(kù)誘發(fā)地震前的VP/VS異常［J］.西北地震學(xué)報(bào)，15(3):38—43.

FENG De-yi，YU Xue-jun，SHENG Guo-ying.1993.Further study and discussion on the anomalies of seismic wave velocity(Ⅲ):VP/VSvaration before reservoir-induced seismic events［J］.Northwestern Seismological Journal，15(3):38—43(in Chinese).

葛洪魁，陳颙，韓德華.2001.有效應(yīng)力對(duì)巖石彈性波速的影響［J］.地球物理學(xué)報(bào)，44(增刊):152—160.

GE Hong-kui，CHEN Yong，HAN De-hua.2001.The effect of effective stress on rock elastic wave velocities［J］.Chinese J Geophys，44(suppl):152—160(in Chinese).

季鐘霖，李建林，熊舜華，等.1993.含油頁(yè)巖彈性特征及其與油氣的關(guān)系［J］.地球物理學(xué)報(bào)，36(2):242—255.

JI Zhong-lin，LI Jian-lin，XIONG Shun-hua，et al.1993.Elastic character of oil-rock and relation with oil-gas formation［J］.Chinese J Geophys，36(2):242—255(in Chinese).

金文正，湯良杰，楊克明，等.2008.龍門山?jīng)_斷帶構(gòu)造特征研究主要進(jìn)展及存在問(wèn)題探討［J］.地質(zhì)論評(píng)，54(1):37—46.

JIN Wen-zheng，TANG Liang-jie，YANG Ke-ming，et al.2008.Progress and problem of study on characters of the Longmen Mountain thrust belt［J］.Geological Review，54(1):37—46(in Chinese).

劉永祥，吳福元，張世紅.1995.高溫高壓下巖石彈性波速特征及其在深部地質(zhì)研究中的意義［J］.地球物理學(xué)進(jìn)展，10(3):82—94.

LIU Yong-xiang，WU Fu-yuan，ZHANG Shi-hong.1995.The characteristics of elastic wave velocity of rocks at high T and P:implications for the study of deepseated geology［J］.Progress in Geophysics，10(3):82—94(in Chinese).

劉志坤，黃金莉.2010.利用背景噪聲互相關(guān)研究汶川地震震源區(qū)地震波速度變化［J］.地球物理學(xué)報(bào)，53(4):853—863.doi:10.3969/j.issn.0001－5733.2010.04.010.

LIU Zhi-kun，HUANG Jin-li.2010.Temporal changes of seismic velocity around the Wenchuan earthquake fault zone from ambient seismic noise correlation［J］.Chinese J Geophys，53(4):853—863(in Chinese).

盧顯，張曉東，周龍泉，等.2013.紫坪鋪水庫(kù)區(qū)域地震波速比計(jì)算及研究［J］.中國(guó)地震，29(2):236—245.

LU Xian，ZHANG Xiao-dong，ZHOU Long-quan，et al.2013.Calculation and research on VP/VSratio in the Zipingpu reservoir area［J］.Earthquake Research in China，29(2):236—245(in Chinese).

馬中高，伍向陽(yáng)，王中海.2006.有效壓力對(duì)巖石縱橫波速度的影響［J］.勘探地球物理進(jìn)展，29(3):183—186.

MA Zhong-gao，WU Xiang-yang，WANG Zhong-hai.2006.The effect of effective pressure on velocity of rock longitudinal and transverse wave［J］.Progress in Exploration Geophysics，29(3):183—186(in Chinese).

湯良杰，楊克明，金文正，等.2008.龍門山?jīng)_斷帶多層次滑脫帶與滑脫構(gòu)造變形［J］.中國(guó)科學(xué)(D輯)，38(增刊Ⅰ):30—40.

TANG Liang-jie，YANG Ke-ming，JIN Wen-zheng，et al.2008.Multi-level decollement zones and detachment deformation of Longmenshan thrust belt，Sichuan Basin，southwest China［J］.Science in China(Ser D)，38(suppl):30—40(in Chinese).

尹尚先，王尚旭.2003.彈性模量、波速與應(yīng)力的關(guān)系及其應(yīng)用［J］.巖土力學(xué)，24(增刊):597—601.

YIN Shang-xian，WANG Shang-xu.2003.Relation of stresses with elastic modulus and velocities and its application［J］.Rock and Soil Mechanics，24(suppl):597—601(in Chinese).

趙盼盼，陳九輝，劉啟元，等.2012.汶川地震區(qū)地殼速度相對(duì)變化的環(huán)境噪聲自相關(guān)研究［J］.地球物理學(xué)報(bào)，55(1):137—145.doi:10.6038/j.issn.0001－5733.2012.01.013.

ZHAO Pan-pan，CHEN Jiu-hui，LIU Qi-yuan，et al.2012.Crustal velocity changes associated with the Wenchuan M8.0 earthquake by auto-correlation function analysis of seismic ambient noise［J］.Chinese J Geophys，55(1):137—145(in Chinese).

周斌.2010.水庫(kù)誘發(fā)地震時(shí)空演化特征及其動(dòng)態(tài)響應(yīng)機(jī)制研究:以紫坪鋪水庫(kù)為例［D］:［學(xué)位論文］.北京:中國(guó)地震局地質(zhì)研究所.

ZHOU Bin.2010.Evolution characteristics of reservoir-induced seismicity in space-time and its dynamic response mechanism－A case study of the Zipingpu reservoir［D］.Dissertation.Institute of Geology，China Earthquake Administration，Beijing(in Chinese).

周斌，薛世峰，鄧志輝，等.2010.水庫(kù)誘發(fā)地震時(shí)空演化與庫(kù)水加卸載及滲透過(guò)程的關(guān)系:以紫坪鋪水庫(kù)為例［J］.地球物理學(xué)報(bào)，53(11):2651—2670.doi:10.3969/j.issn.0001－5733.2010.11.013.

ZHOU Bin，XUE Shi-feng，DENG Zhi-hui，et al.2010.Relationship between the evolution of reservoir induced seismicity in spac-time and the process of reservoir water body load-unloading and water infiltration:A case study of Zipingpu reservoir［J］.Chinese J Geophys，53(11):2651—2670(in Chinese).

周龍泉.2009.紫坪鋪水庫(kù)庫(kù)區(qū)三維速度結(jié)構(gòu)［J］.國(guó)際地震動(dòng)態(tài)，4:5—6.

ZHOU Long-quan.2009.The study on 3D velocity structure of Zipingpu reservoir［J］.Recent Developments in World Seismology，4:5—6(in Chinese).

Biot M A.1956.General solution of the equation of elasticity and consolidation for porous material［J］.Jour Appl Mech，78:91—96.

Chen J H，F(xiàn)roment B，Liu Q Y，et al.2010.Distribution of seismic wave speed changes associated with the 12 May 2008 MW7.9 Wenchuan earthquake［J］.Geophys Res Lett，L18302.doi:10.1029/2010GL044582.

Desbrandes R，Bassiouni Z.1991.Understanding the effect of pore pressure on acoustic velocity measurements，F(xiàn)inal Report［R］.GRI 5087－215－1648.Baton Rouge:Louisiana State university.

Eberhart-Phillips D M.1989.Investigation of crustal structure and active tectonic processes in the Coast Ranges，central Califonria［D］.Ph D thesis.Palo Alto:Stanford University.

Greenfield R J，Graham E K.1994.A convenient functional form consistent with elasticity theory for fitting velocities of rocks containing cracks［J］.EOS，75(16):337.

Han D H，Nur A，Morgan，D.1986.Effect of porosity and clay content on wave velocities in sandstones［J］.Geophysics，51:2093—2107.

Khaksa A，Greiffiths C M.1996.Influence of effective stress on the acoustic velocity and Log derived porosity［J］.SPE 36981.

Nur A，Simmons G.1969.The effect of saturation on velocity in low porosity rocks［J］.Earth Planet Sci Letters，2:183.

Prasad M，Manghnani M H.1997.Effects of pore and differential pressure on compressional wave velocity and quality factor in Berea and Michigan sand stones［J］.Geophysics，62(4):1163—1176.

Weper W W，Christensen N J.1991.A seismic velocity-confining pressure relation with application［J］.Int J Rock Mech Min Sci＆ Geomech Abstr，28(5):451—456.