大崗山水庫庫區(qū)蓄水前后三維速度成像研究

莊園旭，丁 勇，杜 瑤，馮 靜，康 宏，王余偉

（四川省地震局，四川 成都 610041）

大崗山水電站是大渡河干流的大型水電工程之一，庫區(qū)位于青藏高原東南緣，地屬川西南高原區(qū)。庫區(qū)內(nèi)山頂海拔一般為2 000 m以上，庫區(qū)右岸最高峰貢嘎山海拔為7 556 m，左岸馬鞍山海拔為4 021 m，壩區(qū)枯水期河水位952～957 m，水面寬40～80 m。電站正常蓄水位1 130 m，最大壩高210 m，總庫容7.42億m3，是典型的高山峽谷型高壩大庫容水庫。大崗山水庫于2014年12月導(dǎo)流洞下閘開始第一階段蓄水，2015年5月導(dǎo)流底孔下閘開始第二階段蓄水。

庫區(qū)位于鮮水河斷裂帶、龍門山斷裂帶和安寧河斷裂帶交匯部位，構(gòu)造活動強烈，地震地質(zhì)背景復(fù)雜，是高烈度地震危險區(qū)，穿過水庫庫區(qū)的主要活動斷裂包括大渡河斷裂、磨西斷裂，其中磨西斷裂歷史上1786年曾發(fā)生73/4級地震，震中距大壩僅54 km。大崗山水庫具有高壩大庫，且位于地震危險性極高的磨西8級地震潛在震源區(qū)內(nèi)的特點，對大崗山水庫地震活動監(jiān)測和相關(guān)研究工作十分必要。

1 方法與數(shù)據(jù)選取

根據(jù)射線理論，地震i到地震臺k的體波到時T可以表示：

式中，τi為地震的發(fā)震時刻，u為慢度場，ds為射線的路徑微元。

對于近震層析成像，震源坐標（x1，x2，x3）、發(fā)震時刻、射線路徑、慢度矢量都是未知的。地震到時和地震定位的關(guān)系是非線性的，為了方便計算，采用Taylor級數(shù)展開的方法使其線性化，同時用三維網(wǎng)格速度模型，就可以把上式方程寫成一個線性的等式：

對兩事件i和j到同一臺站k的方程作差，有：

假設(shè)這兩個事件彼此靠近，可以認為事件到同一臺站的射線路徑幾乎相同并且速度結(jié)構(gòu)是已知的，則等式（3）可以簡化為：

方程（5）既可以用震相目錄的絕對走時，也可以用互相關(guān)數(shù)據(jù)的相對走時差。

雙差層析成像方法[1]是在雙差定位方法[2]和傳統(tǒng)走時層析成像的基礎(chǔ)上提出的一種走時層析成像方法。與傳統(tǒng)的層析成像方法相比，該方法利用絕對到時以及相對到時實現(xiàn)了三維速度結(jié)構(gòu)的反演和地震重定位，因為震源位置和速度結(jié)構(gòu)存在耦合效應(yīng)，并且真實的速度結(jié)構(gòu)未知，利用絕對走時數(shù)據(jù)確定震源區(qū)域以外區(qū)域的速度結(jié)構(gòu)，利用到時差數(shù)據(jù)確定震源區(qū)的精細速度結(jié)構(gòu)，雙差層析成像可以通過適當調(diào)整速度結(jié)構(gòu)來提高地震定位的精度，地震定位的相對精確可以降低由于射線路徑的不同造成的速度結(jié)構(gòu)反演的誤差，從而提高速度結(jié)構(gòu)反演的精度。因此相比于傳統(tǒng)的層析成像方法（僅采用絕對到時）雙差層析成像方法可以得到更精確的速度結(jié)構(gòu)和地震定位結(jié)果。自該方法提出以來，被廣泛應(yīng)用在區(qū)域成像中。鄒振軒等[3]利用浙江數(shù)字地震臺網(wǎng)大量地震記錄，通過地震波雙差走時層析成像，研究浙江省及周邊地區(qū)上地殼精細速度結(jié)構(gòu)和局部各向異性，并探討地殼速度結(jié)構(gòu)與活動斷裂及地質(zhì)構(gòu)造的關(guān)系。鄧文澤等[4]利用川西流動地震臺陣、汶川地震震后應(yīng)急臺網(wǎng)記錄到的P波到時資料，對2008年5月至2008年10月期間發(fā)生的汶川地震余震序列應(yīng)用雙差層析成像方法進行了地震震源和三維P波速度結(jié)構(gòu)的聯(lián)合反演，確認了汶川地震的發(fā)震構(gòu)造特征。Dixit等[5]利用雙差層析成像的方法，獲得了印度Koyna-Warna水庫區(qū)域高分辨率三維速度結(jié)構(gòu)，并討論了庫水滲透、地震分布和速度異常之間的關(guān)系。羅佳宏[6]利用三峽加密觀測臺網(wǎng)記錄的2009年3月到2010年12月高精度地震數(shù)據(jù)，使用雙差層析成像方法，聯(lián)合反演了三峽庫區(qū)上地殼的震源位置參數(shù)和三維上地殼P波和S波速度結(jié)構(gòu)，討論了速度結(jié)構(gòu)、微震分布和斷裂的關(guān)系。

大崗山水庫地震臺網(wǎng)于2013年3月建成投測。大崗山水庫地震臺網(wǎng)由8個測震臺站組成，布局孔徑NS 20 km×EW 40 km，臺站沿庫區(qū)均勻展布，并全部包圍了可能誘發(fā)地震的重點監(jiān)視區(qū)段，區(qū)內(nèi)監(jiān)測震級下限為ML0.5級。大崗山水庫臺網(wǎng)在臺站地震觀測系統(tǒng)上全部采用反饋式短周期地震計和24位數(shù)據(jù)采集器，觀測頻帶2s～40Hz。同時，大渡河流域還分布有瀑布溝水庫臺網(wǎng)、瀘定水庫臺網(wǎng)、黃金坪水庫臺網(wǎng)等，如圖1所示，黑色小方框為研究區(qū)（101.9°～102.5°E，29.2°～29.9°N）。此外，四川區(qū)域數(shù)字地震臺網(wǎng)在距離大崗山大壩100 km范圍內(nèi)布設(shè)有6個寬頻帶地震儀，記錄了大量觀測資料均可作為本研究的有效補充。

圖1 研究區(qū)域位置及地震臺站分布

本文利用大崗山水庫數(shù)字地震臺網(wǎng)及其他臺網(wǎng)2014年2月1日至2017年12月31日在研究區(qū)內(nèi)記錄到的ML≥0級地震事件12 764次，采用雙差層析成像的方法進行反演。在速度結(jié)構(gòu)反演之前，先要在研究區(qū)建立直角坐標系并設(shè)置反演網(wǎng)格點，研究區(qū)坐標系X軸與Y軸分別沿EW向和NS向，沒有進行旋轉(zhuǎn)，坐標中心點為大崗山水庫大壩（102.2191°E、29.4497°N）。

由于研究區(qū)地震分布和臺站分布不均勻，根據(jù)初始數(shù)據(jù)顯示的研究區(qū)地震分布特點和射線覆蓋路徑，如圖2所示，本文采用非均勻網(wǎng)格點設(shè)置，X軸網(wǎng)格點為-40 km、-35 km、-30 km、-27 km、-24 km、-22 km、-20 km、-18 km、-16 km、-13 km、-10 km、-5 km、0 km、5 km、10 km、20 km、30 km、40 km，共18個節(jié)點；Y軸網(wǎng)格點為-40 km、-30 km、-20 km、-10 km、-5 km、0 km、5 km、10 km、13 km、16 km、18 km、20 km、22 km、24 km、27 km、30 km、35 km、40 km、50 km、60 km，共20個節(jié)點，Z軸網(wǎng)格節(jié)點為0 km、2 km、5 km、10 km、18 km、28 km、40 km。網(wǎng)格點原點為大崗山水庫大壩位置。一維速度模型采用趙珠等[7]給出的四川西部地區(qū)P波速度模型，見表1，波速比Vp/Vs取1.73。

表1 研究區(qū)域地殼P波速度模型

圖2 重定位前射線路徑

2 結(jié)果及分析

2.1 大崗山庫區(qū)地震活動性分析

從2014年2月1日至2017年12月31日，研究區(qū)內(nèi)臺網(wǎng)共記錄到12 764次ML≥0的地震事件，其中ML0～0.9地震9 317次，ML1.0～1.9地震3 097次，ML2.0～2.9地震315次，ML3.0～3.9地震32次，ML≥4.0地震3次分別為2016年3月18日4.2級、2016年3月24日4.0級、2017年1月4日4.0級。大崗山水庫分兩個階段蓄水，在第二階段蓄水位達到正常蓄水位后，水位保持相對恒定。在2015年4月、2016年3-4月、2017年3月出現(xiàn)了地震活動增強的現(xiàn)象，如圖3、圖4所示。

圖3 大崗山庫區(qū)蓄水位與地震日頻次關(guān)系圖

圖4 大崗山庫區(qū)N-T圖

2.2 反演參數(shù)選擇

本文分為四個階段對庫區(qū)三維速度結(jié)構(gòu)和地震分布進行精確定位。第一階段為蓄水前2014年2月1日至2014年12月31日，第二階段為蓄水過程2015年1月1日至2015年12月31日，第三階段為蓄水后2016年1月1日至2016年12月31日，第四階段為2017年1月1日至2017年12月31日。

四個階段反演均采用相同的初始速度模型和相同的網(wǎng)格節(jié)點，反演中共設(shè)置了6組迭代，每組迭代2次，迭代的截斷值設(shè)為9倍標準差，并同時反演速度結(jié)構(gòu)和地震位置。另外，由于數(shù)據(jù)量較大，反演解算時使用帶阻尼LSQR算法，該算法中平滑系數(shù)smooth和阻尼系數(shù)damp約束著地震位置和慢度的變化量，對收斂速度和結(jié)果平滑程度影響較大。為此，本文采用L-curve方法進行測試以找到合適的系數(shù)值，通過歸一化模型與走時殘差關(guān)系曲線（圖5）分析顯示，當固定平滑系數(shù)為100模型較平滑，阻尼系數(shù)分別為80、80、140、60時走時殘差也相對較小，因此在反演成像中采用這兩個值，在合成分辨率測試中smooth取10，damp取值和反演過程取值相同。

圖5 阻尼系數(shù)選擇曲線

2.3 棋盤測試和反演結(jié)果

為確保成像結(jié)果的可靠性，在實際數(shù)據(jù)反演之前，采用棋盤格測試[8]的方法對解的質(zhì)量進行分析，本文棋盤測試中網(wǎng)格設(shè)置與反演時網(wǎng)格設(shè)置相同，速度擾動取正常值的±5%，圖中黑色表示正速度異常，白色表示負速度異常，棋盤測試結(jié)果如圖6、圖10、圖14、圖18所示，由于地震分布集中在鮮水河斷裂西側(cè)并且臺站沿大渡河兩岸分布，造成無地震和臺站區(qū)域的檢測板測試恢復(fù)較差，本文中檢測板測試結(jié)果恢復(fù)較好的范圍為X軸-30～10 km，Y軸-10～30 km區(qū)域，在深度上，2～10 km深度范圍具有較好的模型分辨率，本文只使用恢復(fù)較好區(qū)域成像結(jié)果進行解釋。

圖6 第一階段不同深度V p（左）V s（右）棋盤測試

圖10 第二階段不同深度V p（左）V s（右）棋盤測試

圖14 第三階段不同深度V p（左）V s（右）棋盤測試

圖18 第四階段不同深度V p（左）V s（右）棋盤測試

第一階段為蓄水前，此階段參與反演地震數(shù)2 766條，重定位后獲得2 406次地震的精確定位結(jié)果，均方根殘差從0.415 1 s降低到0.1 s，X、Y、Z三個分向震源位置估算誤差降低到336 m、229.5 m、459 m。不同深度震源位置和速度結(jié)構(gòu)切片如圖7所示，從圖中可以看出，大崗山庫區(qū)的速度結(jié)構(gòu)存在明顯的橫向不均勻性。地震較為集中的X=-30 km到X=-10 km，Y=5 km到Y(jié)=25 km的磨西斷裂西側(cè)區(qū)域（以下稱為A區(qū)），在5 km以上的淺層區(qū)域P波速度和S波速度均表現(xiàn)為高速異常，5～10 km深度范圍高Vp范圍減小，低Vp范圍增加。X=-10 km到X=0 km，Y=0 km到Y(jié)=10 km的水庫蓄水區(qū)（以下稱為B區(qū)），5 km以上P波、S波速度都較低，5～10 km范圍表現(xiàn)出P波速度與A區(qū)相當，S波依然表現(xiàn)為低速。

圖7 第一階段不同深度V p（左）V s（右）速度結(jié)構(gòu)

第一階段地震主要分布在A區(qū)，呈叢集狀分布，定位結(jié)果與杜瑤等[9]使用雙差定位法得到的結(jié)果基本一致。震源深度在10 km以內(nèi)，呈叢集狀分布，且位于高Vp、高Vs區(qū)域，如圖8、圖9所示。從地震分布上看，與磨西斷裂關(guān)系不大，推測存在小型隱伏構(gòu)造。

圖8 大崗山庫區(qū)第一階東西向P波速度和S波速度剖面

圖9 大崗山庫區(qū)第一階南北向P波速度和S波速度剖面

第二階段為開始蓄水至正常蓄水位，此階段參與反演地震2 700個重定位后獲得2 497次地震的精確定位結(jié)果，均方根殘差從0.384 4 s降低到0.08 s，X、Y、Z三個分向震源位置估算誤差降低到227.7 m、206.4 m、333.1 m。不同深度震源位置和速度結(jié)構(gòu)切片如圖11所示，A區(qū)在5 km以上的淺層區(qū)域P波速度和S波速度均表現(xiàn)為高速異常，5～10 km深度范圍高Vp范圍減小，低Vp范圍增加。B區(qū)5 km以上P波、S波速度都較低，5～10 km范圍P波、S波均表現(xiàn)為低速。A區(qū)地震分布與第一階段相似，呈帶狀分布，震源深度主要在7 km以內(nèi)，且處于P波、S波高速區(qū)，根據(jù)圖12、圖13震源深度特征依然能看到與第一階段類似的隱伏斷裂特征。B區(qū)出現(xiàn)明顯的地震叢集，震源深度在2～7 km范圍內(nèi)，如圖22、圖23所示。

圖11 第二階段不同深度V p（左）V s（右）速度結(jié)構(gòu)

圖12 大崗山庫區(qū)第二階段東西向P波速度和S波速度剖面

圖13 大崗山庫區(qū)第二階段南北向P波速度和S波速度剖面

第三階段參與反演地震4 980個重定位后獲得4 678次地震的精確定位結(jié)果，均方根殘差從0.374 s降低到0.071 1 s，X、Y、Z三個分向震源位置估算誤差降低到144.9 m、136.2 m、201.2 m。不同深度震源位置和速度結(jié)構(gòu)切片如圖15所示，速度分布與第二階段類似，A區(qū)在5 km以上的淺層區(qū)域P波速度和S波速度均表現(xiàn)為高速異常，5～10 km深度范圍高Vp范圍減小，低Vp范圍增加。B區(qū)5 km以上P波、S波速度都較低，5～10 km范圍P波、S波均表現(xiàn)為低速。A區(qū)地震分布與前兩個階段相似，呈帶狀分布，如圖16、圖17所示。此階段B區(qū)地震數(shù)目增加，叢集性更明顯，震源深度主要分布在2～8 km范圍內(nèi)，如圖22、圖23所示。

圖15 第三階段不同深度V p（左）V s（右）速度結(jié)構(gòu)

圖16 大崗山庫區(qū)第三階段東西向P波速度和S波速度剖面

圖17 大崗山庫區(qū)第三階段南北向P波速度和S波速度剖面

第四階段參與反演地震2 311個重定位后獲得2 128次地震的精確定位結(jié)果，均方根殘差從0.365 s降低到0.065 9 s，X、Y、Z三個分向震源位置估算誤差降低到208.3 m、207.6 m、282.0 m。不同深度震源位置和速度結(jié)構(gòu)切片如圖19所示，A區(qū)、B區(qū)速度分布和地震分布與第二、三階段相似，在此不再贅述。

圖19 第四階段不同深度V p（左）V s（右）速度結(jié)構(gòu)

針對蓄水后B區(qū)在蓄水后地震活動增強的現(xiàn)象，我們做了Y=5 km、X=-5 km兩個剖面，如圖22、圖23所示，第一階段地震分布散亂，從第二階段開始，B區(qū)出現(xiàn)了地震叢集，震源深度主要集中在2～8 km，B區(qū)在0～5 km范圍P波速度和S波速度均沒有明顯變化，而在5～10 km出現(xiàn)的S波高速異常，是反演過程中地震叢集過于集中該區(qū)域射線密度大而其他區(qū)域射線分布少造成的。

圖22 四個階段庫區(qū)東西向P波速度和S波速度剖面

圖23 四個階段庫區(qū)南北向P波速度和S波速度剖面

圖20 大崗山庫區(qū)第四階段東西向P波速度和S波速度剖面

圖21 大崗山庫區(qū)第四階段南北向P波速度和S波速度剖面

3 討論與結(jié)論

綜上所述，本文利用水庫地震臺網(wǎng)2014年2月1日至2017年12月31日在研究區(qū)（101.9°～102.5°E，29.2°～29.9°N）內(nèi)記錄到的ML≥0級地震事件12 764次，采用雙差層析成像的方法分別對蓄水前、蓄水過程、蓄水后等四個階段進行反演，獲得了水庫蓄水區(qū)及其臨區(qū)三維速度結(jié)構(gòu)和精確的地震定位結(jié)果。棋盤測試結(jié)果顯示，各個階段磨西斷裂西側(cè)的A區(qū)（X=-30 km到X=-10 km，Y=5 km到Y(jié)=25 km）和大崗山水庫蓄水去的B區(qū)（X=-10 km到X=0 km，Y=0 km到Y(jié)=10 km）均有較高分辨率。

從蓄水各階段成像結(jié)果來看，研究區(qū)速度結(jié)構(gòu)存在明顯的橫向不均勻性，各個階段A區(qū)P波、S波速度均高于B區(qū)，結(jié)合地質(zhì)資料，B區(qū)分布有3條斷裂，且大渡河斷裂穿過大渡河，河水可能沿斷裂滲透，造成該區(qū)域P波、S波速度較低。B區(qū)各個階段在0～5 km的淺層P波、S波速度均為未發(fā)生較大變化，這與常見的水庫蓄水后速度結(jié)構(gòu)發(fā)生較大變化不同[5-6，10-12]，分析原因，因大渡河斷裂穿過大渡河，在蓄水前，河水可能就已經(jīng)滲透到地下，孔隙水已達到飽和，而在蓄水后，淺層的水含量未發(fā)生較大變化，因此蓄水前后速度結(jié)構(gòu)未發(fā)生較大變化。

從蓄水各階段地震精確定位結(jié)果來看，A區(qū)地震分布集中且呈帶狀分布，震中分布方向與磨西斷裂走向一致，從各階段震源深度分布可以清晰分辨出A區(qū)存在一條走向與磨西斷裂一致傾向WE的隱伏斷裂，該區(qū)域地震主要受本底構(gòu)造活動影響。B區(qū)蓄水前地震稀少，蓄水后即出現(xiàn)了地震活動增強的現(xiàn)象，震中分布呈叢集狀，震源深度主要在2～8 km范圍且各階段震源深度未發(fā)生明顯變化，P波、S波速度也未發(fā)生明顯變化，因此B區(qū)出現(xiàn)的地震叢集可能是由于庫水荷載引起局部應(yīng)力變化造成的[13]。