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

安艷茹 張曉東

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

2）中國地震臺網中心，北京市西城區三里河南橫街5號 100045

0 引言

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

1 紫坪鋪水庫及地震臺網概況

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

圖1 紫坪鋪水庫水位變化

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

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

圖2 庫壩區地下水類型（據王云基（2001））

圖3 紫坪鋪水庫臺網臺站分布

表1 紫坪鋪水庫臺網臺站及儀器參數表

2 觀測資料

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

圖4 各臺數據連續情況

3 理論基礎與數據處理

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

本文使用MSNoise軟件（Lecocq et al，2014）處理數據。該軟件已應用于奧克蘭火山地區（Kasper，2014）及地震前后（Berkeley Seismological Laboratory，2014）的分析研究中，并得到了穩定可靠的結果。

3.1 噪聲互相關提取經驗格林函數CCFs

3.1.1 單臺數據預處理

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

3.1.2 互相關計算提取經驗格林函數

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

式中，x為點的位置；t為時間；i為虛數單位；an為模態激發函數；un和ωn為地球的本征方程和本征頻率。發散場的一個重要性質就是模態振幅是不相關的隨機變量

式中，δnm為克羅內克函數；F為頻譜能量密度。因式（2）的交叉項在取平均時消失，所以場中x、y點的相互關系簡化為

式中，τ為時間錯動。比較發現，式（3）與x、y點間真實的格林函數相比，僅差一個振幅因子F。因此，根據足夠長時間內得到的場和場的關系，可從發散場中提取2點的格林函數，這正是噪聲成像的理論基礎。之后，該方法在海洋聲學（Roux et al，2004）和地震學（Shapiro et al，2004、2005）等領域得到了迅速的發展。

數學上，頻率域與時間域的互相關計算是等價的。若2個序列x（t）、y（t）的傅里葉變換分別是X（f）、Y（f），則頻率域的互相關公式為

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

本文將1天的數據分成48段×30min（50%重疊），在頻率域對每個臺站對的每段數據做互相關運算，并疊加形成每天的經驗格林函數CCF。在測量波速時，會因為與REF相關度較差的CCF而產生強烈的變化；因此我們分別以10、30、50天作為滑動窗，疊加獲得 CCFs，以提高相關度，獲得更為可靠的波速變化。圖5中，以BAJ-GHS臺站對為例，展示了互相關提取的經驗格林函數以50天為滑動窗疊加后的能量干涉圖，并繪出了不同長度的滑動窗獲得的CCFs與REF的互相關系數曲線。由圖5可見，10天滑動窗的互相關系數曲線抖動較明顯，30、50天滑動窗的互相關系數曲線較為平滑，能夠突出顯著的變化。

圖5 BAJ-GHS臺站對50天疊加的CCFs的能量干涉圖（a）、紫坪鋪水庫水位變化（b）、CCFs與REF互相關系數的變化（c）

3.2 計算走時延遲Δt

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

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

3.2.1 獲取參考格林函數REF

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

3.2.2 計算走時延遲

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

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

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

圖7 28對臺站的參考格林函數（反序疊加）

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

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

式中，Cj是相關系數；Xj為互譜。這種權重公式的優勢在于，能夠在相關系數相對為常數而互譜能量變化的情況下給出較為不同的權重。根據式（7）中Δt與m的關系，將m與em分別除以2π，即可獲得臺站對間對應于不同時間錯動的走時延遲Δt及誤差eΔt。

3.3 計算平均相對波速變化

假設相對波速變化Δv/v在空間上是均勻的，那么它與相對走時變化Δt/t是相反數，即Δv/v=-Δt/t（Ratdomopurboet al，1995）。考慮到臺站間距、信號的能量強弱、CCFs與 REF的相關程度等因素，本文選擇了互相關正負錯動的8～20s部分（圖5），此部分為面波及尾波。我們將28組臺站對間不同時間錯動的走時延遲Δt及對應的誤差eΔt進行平均，得到庫區地下介質平均的走時延遲及平均誤差對于3個周期，我們分別通過最小二乘加權線性回歸獲得為了減小誤差，擬合時要求直線強制通過原點（Lecocqet al，2014），則

其中，pi為權重

b對應的方差為

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

4 結論與討論

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

圖8 庫區地下介質平均相對波速變化

4.1 第1次蓄水及泄水

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

4.2 第2次蓄水及泄水

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

4.3 第3次蓄水

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

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

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

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

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