利用背景噪聲研究紫坪鋪水庫水位加卸載對庫區(qū)地下介質(zhì)波速的影響

安艷茹 張曉東

1）中國地震局地震預(yù)測研究所，北京市復(fù)興路63號 100036

2）中國地震臺(tái)網(wǎng)中心，北京市西城區(qū)三里河南橫街5號 100045

0 引言

中國是世界上水庫數(shù)量最多的國家，截至2011年底，全國已建成水庫97246座，總庫容8104.1億m3，在建水庫756座，總庫容 1219.0億 m3（孫振剛等，2013）。中國也是世界上水庫地震多發(fā)地之一，尤其是近些年在地震多發(fā)的川滇地區(qū)建設(shè)了大量高壩、大庫容的梯級電站，水庫地震時(shí)有發(fā)生，大型水庫的水庫地震問題亦引起了地震學(xué)者的高度重視。汶川地震后，由于紫坪鋪水庫與其特殊的時(shí)空關(guān)系，更是引發(fā)廣泛關(guān)注（雷興林等，2008；Shemin et al，2009；馬文濤等，2011；Deng et al，2010；盧顯等，2010）。雷興林等（2008）認(rèn)為，紫坪鋪水庫在蓄水過程中對其地下的龍門山中央斷層和山前斷層有明顯的作用；Shemin等（2009）認(rèn)為蓄水后的庫侖應(yīng)力增量足以引發(fā)地震；馬文濤等（2011）認(rèn)為汶川地震與紫坪鋪水庫在時(shí)空上均有一定的聯(lián)系；Deng等（2010）認(rèn)為紫坪鋪水庫周邊的M-t圖像沒有增強(qiáng)，庫水不能滲透到汶川地震深度14km的震源處；盧顯等（2010）認(rèn)為都江堰震群與紫坪鋪水庫無關(guān)。無論如何，水庫的蓄水對地下介質(zhì)產(chǎn)生影響都是個(gè)不爭的事實(shí)，而利用數(shù)字地震資料研究庫區(qū)介質(zhì)和應(yīng)力場變化則是十分有效的方法。目前有關(guān)庫區(qū)水位變化與庫區(qū)地震的關(guān)系、庫區(qū)介質(zhì)變化與庫區(qū)地震的關(guān)系等方面的研究較多（周斌等，2010；王惠琳等，2012a、2012b；盧顯等，2013），但有關(guān)庫區(qū)水位變化與庫區(qū)介質(zhì)變化的研究還相對較少。本文基于紫坪鋪水庫數(shù)字地震臺(tái)網(wǎng)的連續(xù)波形數(shù)據(jù)，利用噪聲互相關(guān)技術(shù)研究水位的加卸載和滲透過程對庫區(qū)地下介質(zhì)波速的影響，以求深化對庫區(qū)地下結(jié)構(gòu)、狀態(tài)和動(dòng)力學(xué)過程的認(rèn)識。

1 紫坪鋪水庫及地震臺(tái)網(wǎng)概況

紫坪鋪水庫位于四川省都江堰市區(qū)西北方向9km處的岷江上游麻溪鄉(xiāng)，距成都市區(qū)約60km，距汶川8.0級大地震微觀震中最近處僅6km左右。紫坪鋪水庫于2001年3月29日正式動(dòng)工興建，2006年12月竣工，最大壩高156m，總庫容11.12×108m3。本文利用四川省地震局水庫地震研究所提供的水位數(shù)據(jù)，繪制了紫坪鋪水庫水位變化圖（圖1）。由圖1可見，2005年9月30日下閘蓄水，數(shù)天內(nèi)水位迅速升高至死水位817m以上，蓄水后最高水位為875.18m。2005年9月30日～2007年底，紫坪鋪水庫共經(jīng)歷了3次大規(guī)模蓄水、2次泄水過程。2005年12月5日水位達(dá)835.91m，2006年10月14日水位上升到最高值875.18m，2007年12月12日水位再次上升到873.39m。

圖1 紫坪鋪水庫水位變化

紫坪鋪水庫位于青藏高原東緣的龍門山造山帶的中段。庫區(qū)及鄰區(qū)被主干斷裂分隔，可劃分為茂汶韌性剪切帶、中央推覆構(gòu)造帶、龍門山前緣拆離帶、前陸擴(kuò)展變形帶和川西前陸盆地等5個(gè)地質(zhì)構(gòu)造單元，它們在物質(zhì)組成、構(gòu)造層次和變形樣式等方面具有明顯的差異。紫坪鋪水庫庫體坐落于龍門山前緣拆離帶內(nèi)，鄰近中央推覆構(gòu)造帶和前陸擴(kuò)展變形帶（周斌等，2010）。庫區(qū)地下水分為碎屑巖類裂隙孔隙水和碳酸鹽巖巖溶水等2大類（圖2），如圖2所示，庫區(qū)西南主要巖性為碳酸鹽巖，東北主要巖性為碎屑巖和部分碳酸鹽巖。

紫坪鋪水庫數(shù)字遙測地震臺(tái)網(wǎng)于2004年8月16日正式采集地震信息，2005年6月27日通過驗(yàn)收（胡先明等，2006）。臺(tái)網(wǎng)包括7個(gè)數(shù)字遙測地震臺(tái)，使用短周期地震儀，頻帶寬度1～40Hz。此外，該臺(tái)網(wǎng)共享四川區(qū)域數(shù)字臺(tái)網(wǎng)中YZP（油榨坪）臺(tái)的數(shù)據(jù)。8個(gè)地震臺(tái)大致均勻地分布在庫區(qū)周邊，平均臺(tái)間距約為14km，臺(tái)站分布如圖3所示，臺(tái)站及儀器參數(shù)見表1。

圖2 庫壩區(qū)地下水類型（據(jù)王云基（2001））

圖3 紫坪鋪水庫臺(tái)網(wǎng)臺(tái)站分布

表1 紫坪鋪水庫臺(tái)網(wǎng)臺(tái)站及儀器參數(shù)表

2 觀測資料

本研究使用紫坪鋪臺(tái)網(wǎng)及YZP臺(tái)垂直分量的連續(xù)波形數(shù)據(jù)，時(shí)間范圍2005年1月1日～2008年1月1日。圖4為各臺(tái)站數(shù)據(jù)連續(xù)記錄情況，可用數(shù)據(jù)的天數(shù)為1003天，最長連續(xù)空缺時(shí)間為27天。

圖4 各臺(tái)數(shù)據(jù)連續(xù)情況

3 理論基礎(chǔ)與數(shù)據(jù)處理

利用背景噪聲測量地下介質(zhì)波速連續(xù)變化的主要技術(shù)路徑是：比較經(jīng)驗(yàn)格林函數(shù)（CCFs）與參考格林函數(shù)（REF）。其中，CCFs表征的是一段時(shí)間內(nèi)地下介質(zhì)的狀態(tài)，而 REF代表的是地下介質(zhì)的背景狀態(tài)。計(jì)算過程分為4步：①計(jì)算臺(tái)站對間不同時(shí)間段的噪聲互相關(guān)函數(shù)，即經(jīng)驗(yàn)格林函數(shù)CCFs；②獲取參考格林函數(shù)REF，并計(jì)算每個(gè)CCFs與REF的走時(shí)延遲Δt；③對不同臺(tái)站對間不同時(shí)間錯(cuò)動(dòng)的走時(shí)延遲進(jìn)行平均，并使用一個(gè)相對速度均勻一致變化（Δv/v=常數(shù)）的簡單模型（Lecocq et al，2014）來表征庫區(qū)地下介質(zhì)波速的變化。

本文使用MSNoise軟件（Lecocq et al，2014）處理數(shù)據(jù)。該軟件已應(yīng)用于奧克蘭火山地區(qū)（Kasper，2014）及地震前后（Berkeley Seismological Laboratory，2014）的分析研究中，并得到了穩(wěn)定可靠的結(jié)果。

3.1 噪聲互相關(guān)提取經(jīng)驗(yàn)格林函數(shù)CCFs

3.1.1 單臺(tái)數(shù)據(jù)預(yù)處理

本文研究的是庫區(qū)地下介質(zhì)微小的波速變化，因此能否得到穩(wěn)定可靠的CCFs非常關(guān)鍵。首先，我們必須通過數(shù)據(jù)預(yù)處理，獲得較為純凈的背景噪聲，以保證CCFs的質(zhì)量。將每個(gè)臺(tái)站的連續(xù)波形數(shù)據(jù)合并成每天一個(gè)文件，然后進(jìn)行預(yù)處理，包括：去平均、尖滅處理，高通、低通濾波，降采樣，時(shí)間域歸一化，譜白化處理等。通常，波形數(shù)據(jù)總會(huì)存在一個(gè)非零的均值，這會(huì)影響到數(shù)據(jù)的分析；為了使數(shù)據(jù)標(biāo)準(zhǔn)化，適用于各種標(biāo)準(zhǔn)公式，必須在數(shù)據(jù)分析前去平均值（鄭治真，1979）。另外，在對數(shù)據(jù)進(jìn)行譜域操作，如FFT、濾波時(shí)，若數(shù)據(jù)的兩端不為零，則會(huì)出現(xiàn)譜域假象，因而實(shí)際數(shù)據(jù)經(jīng)常需要做尖滅處理，以使數(shù)據(jù)兩端在短時(shí)間窗內(nèi)逐漸變成零值。本文研究的是庫區(qū)淺層介質(zhì)波速的變化，關(guān)注的是相對高頻的信息，所以我們選擇了0.1～2.0Hz的濾波參數(shù)。為了減少計(jì)算量和存儲(chǔ)量，將數(shù)據(jù)降采樣到10Hz。為了減少地震信號、儀器異常信號和臺(tái)站附近非穩(wěn)定噪聲源的影響，我們對數(shù)據(jù)進(jìn)行了時(shí)間域歸一化處理（Tukey，1962）。為了使信號在不同頻率上的能量更為均衡，我們對數(shù)據(jù)進(jìn)行了譜白化處理。

3.1.2 互相關(guān)計(jì)算提取經(jīng)驗(yàn)格林函數(shù)

利用背景噪聲提取地下結(jié)構(gòu)信息設(shè)想的探索最初可追溯到20世紀(jì)50、60年代，Aki（1957）認(rèn)為可以用背景噪聲提取地下結(jié)構(gòu)中的面波頻散性質(zhì)；Claerbout（1968）提出可以利用背景噪聲恢復(fù)一維層狀介質(zhì)的反射響應(yīng)。類似想法的首次成功應(yīng)用是在太陽地震學(xué)研究中，即通過對太陽表面噪聲進(jìn)行互相關(guān)計(jì)算，成功地提取出了聲波時(shí)距曲線（Duvall et al，1993）。隨后，噪聲互相關(guān)技術(shù)在超聲學(xué)領(lǐng)域獲得重大進(jìn)展。Weaver等（2001）在聲學(xué)上給出了隨機(jī)波場的互相關(guān)特性，即通過從發(fā)散的或隨機(jī)的波場中提取格林函數(shù)來計(jì)算地球的彈性響應(yīng)。證明該特性的實(shí)例就是基于一個(gè)彈性體中有發(fā)散場的模型

式中，x為點(diǎn)的位置；t為時(shí)間；i為虛數(shù)單位；an為模態(tài)激發(fā)函數(shù)；un和ωn為地球的本征方程和本征頻率。發(fā)散場的一個(gè)重要性質(zhì)就是模態(tài)振幅是不相關(guān)的隨機(jī)變量

式中，δnm為克羅內(nèi)克函數(shù)；F為頻譜能量密度。因式（2）的交叉項(xiàng)在取平均時(shí)消失，所以場中x、y點(diǎn)的相互關(guān)系簡化為

式中，τ為時(shí)間錯(cuò)動(dòng)。比較發(fā)現(xiàn)，式（3）與x、y點(diǎn)間真實(shí)的格林函數(shù)相比，僅差一個(gè)振幅因子F。因此，根據(jù)足夠長時(shí)間內(nèi)得到的場和場的關(guān)系，可從發(fā)散場中提取2點(diǎn)的格林函數(shù)，這正是噪聲成像的理論基礎(chǔ)。之后，該方法在海洋聲學(xué)（Roux et al，2004）和地震學(xué)（Shapiro et al，2004、2005）等領(lǐng)域得到了迅速的發(fā)展。

數(shù)學(xué)上，頻率域與時(shí)間域的互相關(guān)計(jì)算是等價(jià)的。若2個(gè)序列x（t）、y（t）的傅里葉變換分別是X（f）、Y（f），則頻率域的互相關(guān)公式為

其中，X*（f）為X（f）的復(fù)共軛；而C（f）的傅里葉反變換正是時(shí)間域的互相關(guān)c（t）。

本文將1天的數(shù)據(jù)分成48段×30min（50%重疊），在頻率域?qū)γ總€(gè)臺(tái)站對的每段數(shù)據(jù)做互相關(guān)運(yùn)算，并疊加形成每天的經(jīng)驗(yàn)格林函數(shù)CCF。在測量波速時(shí)，會(huì)因?yàn)榕cREF相關(guān)度較差的CCF而產(chǎn)生強(qiáng)烈的變化；因此我們分別以10、30、50天作為滑動(dòng)窗，疊加獲得 CCFs，以提高相關(guān)度，獲得更為可靠的波速變化。圖5中，以BAJ-GHS臺(tái)站對為例，展示了互相關(guān)提取的經(jīng)驗(yàn)格林函數(shù)以50天為滑動(dòng)窗疊加后的能量干涉圖，并繪出了不同長度的滑動(dòng)窗獲得的CCFs與REF的互相關(guān)系數(shù)曲線。由圖5可見，10天滑動(dòng)窗的互相關(guān)系數(shù)曲線抖動(dòng)較明顯，30、50天滑動(dòng)窗的互相關(guān)系數(shù)曲線較為平滑，能夠突出顯著的變化。

圖5 BAJ-GHS臺(tái)站對50天疊加的CCFs的能量干涉圖（a）、紫坪鋪水庫水位變化（b）、CCFs與REF互相關(guān)系數(shù)的變化（c）

3.2 計(jì)算走時(shí)延遲Δt

為了研究庫區(qū)水位對地下不同深度介質(zhì)波速的影響，分析水位壓力及滲透的作用過程，我們分別對 1～2s（0.5～1.0Hz）、2～4s（0.25～0.50Hz）、4～8s（0.125～0.250Hz）等 3個(gè)周期進(jìn)行計(jì)算。根據(jù)瑞雷面波的特性（圖6），上述3個(gè)周期的波分別對深度0～2、1～4、1～8km的介質(zhì)變化敏感。

圖6 瑞雷面波波速在不同周期的敏感性（據(jù)Liu等（2014））

3.2.1 獲取參考格林函數(shù)REF

獲取參考格林函數(shù)REF有2種方法，即絕對疊加法和相對疊加法（Lecocq et al，2014）。由于水庫蓄水導(dǎo)致的地下介質(zhì)波速變化并不劇烈，因此我們選擇絕對疊加法，即疊加所有的CCF成為REF。圖7為獲得的28對臺(tái)站的參考格林函數(shù)，此處進(jìn)行了反序疊加，旨在展示的方便。

3.2.2 計(jì)算走時(shí)延遲

本文使用移動(dòng)窗口互譜方法計(jì)算相對走時(shí)變化，該方法由 Ratdomopurbo等（1995）提出，最早被Poupinet等（1984）用于研究地震對之間的相對速度變化，其優(yōu)點(diǎn)是在頻率域進(jìn)行計(jì)算，可以明確相關(guān)函數(shù)中相關(guān)信號的帶寬。Brenguier等（2008a、2008b）利用該方法先后研究了法國富爾奈斯火山地下介質(zhì)的相對波速變化及2004年P(guān)arkfield地區(qū)6.0級地震前后沿圣安德烈斯斷層地下介質(zhì)的相對波速變化。Clarke等（2011）詳細(xì)闡述了移動(dòng)窗口互譜方法的原理和步驟，并驗(yàn)證了該方法監(jiān)測地殼波速變化的分辨率和準(zhǔn)確率。

首先，我們將序列CCFs與REF分成許多重疊的窗，去平均、尖滅處理，然后進(jìn)行傅里葉變換得到Fcur（f）和Fref（f）。互譜X（f）被定義為

式中，f為頻率，星號表示復(fù)共軛性。在頻率域，兩序列的相似程度通過能量密度間的相關(guān)系數(shù)C（f）來評估。由式（5）推導(dǎo)得出

圖7 28對臺(tái)站的參考格林函數(shù)（反序疊加）

式中，φ（f）為解纏相位，它與頻率f成線性比例關(guān)系

對于不同的頻率f和相位φ，可通過最小二乘加權(quán)線性回歸獲得m，對應(yīng)的誤差為em，對應(yīng)的權(quán)重為wj

式中，Cj是相關(guān)系數(shù)；Xj為互譜。這種權(quán)重公式的優(yōu)勢在于，能夠在相關(guān)系數(shù)相對為常數(shù)而互譜能量變化的情況下給出較為不同的權(quán)重。根據(jù)式（7）中Δt與m的關(guān)系，將m與em分別除以2π，即可獲得臺(tái)站對間對應(yīng)于不同時(shí)間錯(cuò)動(dòng)的走時(shí)延遲Δt及誤差eΔt。

3.3 計(jì)算平均相對波速變化

假設(shè)相對波速變化Δv/v在空間上是均勻的，那么它與相對走時(shí)變化Δt/t是相反數(shù)，即Δv/v=-Δt/t（Ratdomopurboet al，1995）?？紤]到臺(tái)站間距、信號的能量強(qiáng)弱、CCFs與 REF的相關(guān)程度等因素，本文選擇了互相關(guān)正負(fù)錯(cuò)動(dòng)的8～20s部分（圖5），此部分為面波及尾波。我們將28組臺(tái)站對間不同時(shí)間錯(cuò)動(dòng)的走時(shí)延遲Δt及對應(yīng)的誤差eΔt進(jìn)行平均，得到庫區(qū)地下介質(zhì)平均的走時(shí)延遲及平均誤差對于3個(gè)周期，我們分別通過最小二乘加權(quán)線性回歸獲得為了減小誤差，擬合時(shí)要求直線強(qiáng)制通過原點(diǎn)（Lecocqet al，2014），則

其中，pi為權(quán)重

b對應(yīng)的方差為

最終，根據(jù)式（9），我們獲得了庫區(qū)地下介質(zhì)平均的相對波速變化（圖 8）。

4 結(jié)論與討論

本文的研究時(shí)間段包含了紫坪鋪水庫3次大規(guī)模蓄水和2次泄水的過程，通過噪聲互相關(guān)提取了庫區(qū)28組臺(tái)站對間的經(jīng)驗(yàn)格林函數(shù)CCFs。圖5中以BAJ-GHS臺(tái)站對為例展示了臺(tái)站對間的互相關(guān)結(jié)果，這對臺(tái)站跨越了西南至東北的整個(gè)庫區(qū)（圖3），十分具有代表性。從圖5中可以看出，水位加卸載時(shí)CCFs與REF的互相關(guān)系數(shù)明顯降低，并在時(shí)間上有一定延遲，這定性地表征了水位影響了庫區(qū)地下介質(zhì)波速的變化。本文利用互相關(guān)結(jié)果，通過移動(dòng)窗口互譜方法，分別對 1～2s（0.5～1.0Hz）、2～4s（0.25～0.50Hz）、4～8s（0.125～0.250Hz）等3個(gè)周期，計(jì)算了28組臺(tái)站對間不同時(shí)間錯(cuò)動(dòng)的走時(shí)延遲，3個(gè)周期分別對應(yīng)著0～2、1～4、1～8km的深度。最終，使用面波及尾波部分，通過加權(quán)線性回歸捕捉到了庫區(qū)地下介質(zhì)平均的相對波速變化（圖8）。可以看出，庫區(qū)地下介質(zhì)相對波速變化不超過0.1%。水庫水位的加卸載主要通過壓力和滲透2種作用影響地下介質(zhì)。Biot（1956）認(rèn)為，孔隙流體會(huì)增加壓縮波的速度并降低剪切波的速度。庫區(qū)以碳酸鹽巖為主，具備較好的滲透條件。圖8的3個(gè)周期中，1～2s對應(yīng)的2km內(nèi)的介質(zhì)相對波速變化最大，變化幅度超過了0.05%，這是由于在同樣的水位壓力下，淺層的滲透作用更強(qiáng)烈。

圖8 庫區(qū)地下介質(zhì)平均相對波速變化

4.1 第1次蓄水及泄水

由圖8可見，在2005年9月30日第1次大規(guī)模蓄水前，水位與相對波速變化沒有明顯關(guān)系。蓄水后數(shù)10天內(nèi)，水位由760.36m迅速升高至835.91m，3個(gè)周期對應(yīng)的3個(gè)深度的介質(zhì)相對波速隨之同時(shí)迅速降低，說明此時(shí)庫水壓力作用是導(dǎo)致相對波速變化的主要原因，而滲透作用只影響了2km以內(nèi)的淺層介質(zhì)。之后，水庫開始泄水，對應(yīng)的相對波速也逐漸升高。

4.2 第2次蓄水及泄水

第2次蓄水開始于2006年8月15日，2006年10月14日水位上升到最高值875.18m，前述3個(gè)深度的介質(zhì)相對波速隨之降低，存在明顯的時(shí)間延遲。其中，1～2、2～4s對應(yīng)的相對波速幾乎同時(shí)降到最小值，而4～8s的相對波速在近1個(gè)月后才降到最小值。這說明在承受同樣的庫水壓力下，滲透作用已成為影響相對波速變化的主要因素，且已影響至深度4km左右的斷層。在第2次泄水過程中，3個(gè)深度對應(yīng)的相對波速變化略有延遲，但幾乎同時(shí)升高至最大值，這說明此時(shí)滲透作用已影響至深度8km左右的斷層。

4.3 第3次蓄水

第3次蓄水于2007年12月12日上升到873.39m，3個(gè)周期的相對波速幾乎同時(shí)降到最小值，與水位變化非常一致，互相關(guān)系數(shù)分別達(dá)-0.81，-0.66，-0.76（圖9）。此時(shí)庫水壓力和滲透作用共同影響著介質(zhì)相對波速的變化，水已滲透至深度為8km以下的斷層。

圖9 第3次蓄水時(shí)的水位與相對波速變化的最小二乘擬合曲線

綜上，本文分別對紫坪鋪水庫3次蓄水及2次泄水過程進(jìn)行了分析，討論了不同階段影響相對波速變化的主要因素，追蹤了水在不同深度介質(zhì)中的滲透過程。結(jié)果表明：在紫坪鋪水庫大規(guī)模蓄水前，相對波速變化與水位變化沒有明顯的相關(guān)性；而在之后的水位加卸載過程中，表現(xiàn)出了明顯的負(fù)相關(guān)性，即水位升高，相對波速降低，水位降低，相對波速升高。分析認(rèn)為，第1次蓄水時(shí)，壓力是導(dǎo)致介質(zhì)相對波速變化的主要因素；而后2次蓄水時(shí)，滲透作用成為主要因素。在第2次蓄水至最高點(diǎn)時(shí)，滲透作用影響至深度4km左右的斷層，在第2次泄水至最低點(diǎn)時(shí)，滲透作用已影響至深度8km左右的斷層。

盧顯等（2010）通過研究水庫水位與水庫區(qū)域地震頻度的關(guān)系發(fā)現(xiàn)，在2005年9月、2006年9月和2007年9月這3個(gè)水位快速上升的月份，地震頻次也相應(yīng)增加，且地震的震源深度絕大部分集中在5km處。根據(jù)本文的研究，這3個(gè)月份對應(yīng)著相對波速由0值左右迅速降低的過程。因此對于地下淺層介質(zhì)相對波速連續(xù)變化的研究，有可能作為水庫地震預(yù)測的依據(jù)，具有十分重要的意義。

致謝：本文使用的大量連續(xù)波形及水庫資料均由四川省地震局提供。研究過程中，本文得到了馬延路博士、楊志高博士、周龍泉研究員、杜方研究員、戴仕貴高工、謝蓉華高工和韓進(jìn)高工的支持與幫助，在此一并表示衷心的感謝。同時(shí)，非常感謝審稿專家提出的寶貴意見。