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

（四川省地震局，四川 成都 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)造活動(dòng)強(qiáng)烈，地震地質(zhì)背景復(fù)雜，是高烈度地震危險(xiǎn)區(qū)，穿過水庫庫區(qū)的主要活動(dòng)斷裂包括大渡河斷裂、磨西斷裂，其中磨西斷裂歷史上1786年曾發(fā)生73/4級(jí)地震，震中距大壩僅54 km。大崗山水庫具有高壩大庫，且位于地震危險(xiǎn)性極高的磨西8級(jí)地震潛在震源區(qū)內(nèi)的特點(diǎn)，對(duì)大崗山水庫地震活動(dòng)監(jiān)測(cè)和相關(guān)研究工作十分必要。

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

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

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

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

對(duì)兩事件i和j到同一臺(tái)站k的方程作差，有：

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

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

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

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

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

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

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

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

圖2 重定位前射線路徑

2 結(jié)果及分析

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

從2014年2月1日至2017年12月31日，研究區(qū)內(nèi)臺(tá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級(jí)、2016年3月24日4.0級(jí)、2017年1月4日4.0級(jí)。大崗山水庫分兩個(gè)階段蓄水，在第二階段蓄水位達(dá)到正常蓄水位后，水位保持相對(duì)恒定。在2015年4月、2016年3-4月、2017年3月出現(xiàn)了地震活動(dòng)增強(qiáng)的現(xiàn)象，如圖3、圖4所示。

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

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

2.2 反演參數(shù)選擇

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

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

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

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

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

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

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

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

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

第一階段為蓄水前，此階段參與反演地震數(shù)2 766條，重定位后獲得2 406次地震的精確定位結(jié)果，均方根殘差從0.415 1 s降低到0.1 s，X、Y、Z三個(gè)分向震源位置估算誤差降低到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ū)相當(dāng)，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)系不大，推測(cè)存在小型隱伏構(gòu)造。

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

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

第二階段為開始蓄水至正常蓄水位，此階段參與反演地震2 700個(gè)重定位后獲得2 497次地震的精確定位結(jié)果，均方根殘差從0.384 4 s降低到0.08 s，X、Y、Z三個(gè)分向震源位置估算誤差降低到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個(gè)重定位后獲得4 678次地震的精確定位結(jié)果，均方根殘差從0.374 s降低到0.071 1 s，X、Y、Z三個(gè)分向震源位置估算誤差降低到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ū)地震分布與前兩個(gè)階段相似，呈帶狀分布，如圖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個(gè)重定位后獲得2 128次地震的精確定位結(jié)果，均方根殘差從0.365 s降低到0.065 9 s，X、Y、Z三個(gè)分向震源位置估算誤差降低到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)

針對(duì)蓄水后B區(qū)在蓄水后地震活動(dòng)增強(qiáng)的現(xiàn)象，我們做了Y=5 km、X=-5 km兩個(gè)剖面，如圖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 四個(gè)階段庫區(qū)東西向P波速度和S波速度剖面

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

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

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

3 討論與結(jié)論

綜上所述，本文利用水庫地震臺(tái)網(wǎng)2014年2月1日至2017年12月31日在研究區(qū)（101.9°～102.5°E，29.2°～29.9°N）內(nèi)記錄到的ML≥0級(jí)地震事件12 764次，采用雙差層析成像的方法分別對(duì)蓄水前、蓄水過程、蓄水后等四個(gè)階段進(jìn)行反演，獲得了水庫蓄水區(qū)及其臨區(qū)三維速度結(jié)構(gòu)和精確的地震定位結(jié)果。棋盤測(cè)試結(jié)果顯示，各個(gè)階段磨西斷裂西側(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)存在明顯的橫向不均勻性，各個(gè)階段A區(qū)P波、S波速度均高于B區(qū)，結(jié)合地質(zhì)資料，B區(qū)分布有3條斷裂，且大渡河斷裂穿過大渡河，河水可能沿?cái)嗔褲B透，造成該區(qū)域P波、S波速度較低。B區(qū)各個(gè)階段在0～5 km的淺層P波、S波速度均為未發(fā)生較大變化，這與常見的水庫蓄水后速度結(jié)構(gòu)發(fā)生較大變化不同[5-6，10-12]，分析原因，因大渡河斷裂穿過大渡河，在蓄水前，河水可能就已經(jīng)滲透到地下，孔隙水已達(dá)到飽和，而在蓄水后，淺層的水含量未發(fā)生較大變化，因此蓄水前后速度結(jié)構(gòu)未發(fā)生較大變化。

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