曹穎 黃江培 付虹
云南省地震局,昆明市北辰大道148號 650224
目前,對水庫誘發地震的形成有多種認識和解釋,但都一致認為,水庫所在區域地質和地球物理環境、庫區從地表到深部的巖石組成和物質結構狀態、區域斷裂構造分布、產狀、力學性質、發育程度、現今活動性等,都是發生水庫誘發地震的重要因素(任金衛,1993)。近年來,隨著數字地震觀測技術的發展,深部地球物理探測研究是獲得區域深部構造環境和孕震環境的重要手段。層析成像方法不僅能用于一般區域,對于水庫地區也同樣適用。利用層析成像方法探測地下速度結構,并結合地震震源分布揭示地下介質的特點,為研究水庫誘發地震的形成提供了新的途徑。王亮等(2015)利用速度結構和小震精定位結果得出紫坪鋪水庫水的滲透作用對汶川地震的發生沒有直接影響的結論;鐘羽云等(2010)利用速度結構和震源參數結果研究認為,由于水庫蓄水后庫水下滲,地震大多發生在低速異常區內,水庫誘發地震之初的幾年中震源深度有一個逐漸變大的過程;李強等(2009)對三峽水庫壩址及鄰區中上地殼P波速度結構研究發現,庫水滲透作用對地殼淺層速度結構產生影響。
tomoDD雙差地震層析成像方法由Zhang等(2003)首先提出,該方法避免了復雜速度結構引起的射線路徑異常對定位結果的影響,同時,利用觀測走時、走時殘差進行地震定位和速度結構聯合反演,不僅能夠獲得較常規地震層析成像更為精確的震源密集區附近的速度結構,還能夠獲得地震發生的絕對位置,其地震定位的精度與雙差地震定位的精度相當。在國外雙差層析成像方法被廣泛應用于地震區精細速度結構研究,并已獲得了從區域尺度到幾百米尺度的高分辨率成像和地震定位結果(Shelly et al,2006;Zhang et al,2004;Okada et al,2005;Pei et al,2010),在我國也有很多應用(于湘偉等,2010;李海鷗等,2011;鄧文澤等,2014;王小娜等,2015;呂子強等,2016)。
云南小灣水庫是滇西瀾滄江中下游河段梯級電站的龍頭水庫,水庫設計最大壩高292m,設計正常蓄水位高程1242m,水庫總庫容約149×108m3,水庫區跨越大理、臨滄、保山3個地區,由西支干流瀾滄江和東支支流黑惠江組成,庫長分別為178、124km,水庫靠近庫壩區的最高海拔約1160m,靠近庫尾及庫區正北的黑惠江流域最高海拔約1210m,總體呈東南低、西北高的地形特點,水庫于2008年12月16日開始蓄水。小灣水庫位于構造活動較強烈的地區,在唐古拉-昌都-蘭坪-思茅褶皺系和岡底斯-念青唐古拉褶皺系的接合部位(毛玉平等,2004),庫區出露的地層主要是變質巖和沉積巖,局部亦有巖漿巖發育(任金衛,1993)。庫區及其周邊主要斷裂有怒江斷裂、保山-施甸斷裂、瀾滄江斷裂、蘭坪-云龍斷裂、維西-喬后斷裂、紅河斷裂、程海-賓川斷裂、無量山斷裂、南汀河斷裂、昌寧斷裂、柯街斷裂等(圖1(b))。歷史上該區域地震多發,距壩址僅75km處有7級強震的活動記錄。對于小灣水庫的研究已有很多,鄔成棟等(2010)利用小灣電站水庫誘發地震監測臺網記錄的數字波形資料,使用遺傳算法反演了小灣水庫近場328個中小地震的地震矩、應力降和拐角頻率等震源參數;李永莉等(2012)計算了小灣水庫蓄水前后的水庫地震波速比變化;姜金鐘等(2016)利用地震精定位方法分析了小灣水庫的地震活動性;柯乃琛等(2016)使用 Velest程序計算得到小灣庫區蓄水前的最小一維速度結構,然后以該最小一維速度結構作為初始模型,使用Simul2000程序計算了小灣水庫蓄水前后庫區上地殼介質三維速度結構和震源位置。在上述研究的基礎上,本文使用雙差地震層析成像方法(tomoDD)聯合絕對到時、相對到時計算小灣水庫庫區及周邊區域水庫蓄水后不同時間段內的地震重定位結果和三維P波速度結構,分析由蓄水所導致的該區域地下介質波速的變化,以期深化對小灣水庫地震發震構造、發震機理等的認識。

圖1 研究區域活動斷裂、地震分布、網格劃分和所用臺站分布
為監測小灣水庫及其周邊的地震活動,2005年6月1日建成并運行小灣水庫地震臺網12個一期臺站,臺站包圍并均勻分布于庫壩區,2009年9月1日開始運行二期臺站4個,主要分布于庫首區,2010年4月小灣一、二期臺站合并運行。2012年8月1日漫灣電站水庫地震臺網建成,并與小灣水庫臺網聯合運行。但從2014年開始水庫臺網運行率下降,基本無法記錄到庫壩區的地震,所以在本文中加入了距研究區內地震約250km范圍內的云南區域地震臺網的49個臺站,包括30個固定臺站、3個騰沖火山臺網臺站、3個下關小孔徑臺網臺站以及13個流動臺站(圖1(a))。本文選取的研究范圍為小灣水庫庫區及其周邊地區(24.5°~25.5°N,99.0°~100.5°E),數據為至少有 6個臺站記錄到的地震事件,時間段為小灣水庫蓄水后的 2008年12月16日~2016年12月31日,經過篩選共挑選出 11249個ML≥-0.5地震(圖1(b)),再去除走時曲線中離散較大的震相,最后得到77306個P波絕對到時資料(圖2)。

圖2 所選取的走時-震中距曲線
由小灣水庫水位資料(圖3)可知,小灣水庫開始蓄水后水位快速上升,至2009年7月,水位由蓄水前的996m上升至1100m,2010年10月水位上升至1210m,由圖3可見,2次水位快速上升期間庫區內部分區域的地震活動頻次明顯增多。隨后,水位以年為周期在1200~1250m間變化(姜金鐘等,2016)。根據水位變化,選取2008年12月16日~2011年6月30日及2011年7月1日~2016年12月31日2個時間段來分析研究區地下速度結構,其中,第1個時間段是蓄水后水位快速上升至平穩變化的時間段,在這個時間段內地震活動明顯增多,第2個時間段是每年水位周期變化的時間段。由柯乃琛(2016)的研究可知,從小灣水庫開始蓄水至2011年6月地震頻次大量增加,此后一直保持在較高的水平。第1個時間段內符合選取條件的地震事件有4474個,并產生31921條絕對P波到時,最大震級地震為2010年6月1日施甸ML4.8地震;第2個時間段內有6775個地震,產生45385條絕對P波到時,最大震級地震是2015年10月30日昌寧MS5.1地震,其震中距水庫大壩僅77km,距瀾滄江僅11km。在此基礎上進行地震對匹配,選擇每次地震與2個地震之間的最大距離為20km,每個地震最多可與30個地震組成地震對,最終構建了相關到時數據。其中,2008年12月16日~2011年6月30日有270189條P波相對到時,2011年7月1日~2016年12月31日有381176條P波相對到時。

圖3 小灣水庫水位變化與地震月頻次
反演中選取的坐標原點為(25°N,99.75°E)。在劃分網格時,較小的網格分辨率較高,但若網格太小又會影響地震層析成像圖像的質量,所以,根據所選地震事件和臺站的分布情況,在反演之前進行了大量的分辨率測試,以尋求最佳的網格分布,我們分別測試了5km×5km、10km×10km、15km×15km、20km×20km的網格間隔。雖然2個時間段使用的臺站不一樣,但差別并不大,且2個時間段的地震分布密集區基本相同,均為黑惠江和昌寧-施甸一帶,2個時間段的橫向、縱向分辨率均可達到10km,研究區東側和北側外圍地震分布較稀疏,故使用 20km×20km的間隔(圖1(b))。2個時間段垂直向網格相同,均為 0、3、5、10、15、20、25、30km。選擇一維模型作為初始模型,初始一維速度模型為使用Kissling方法得到的小灣水庫最小一維速度模型(表1)(柯乃琛等,2016),2個時間段均使用該初始模型。空間任意點的速度利用線性插值求得(Thurber,1983)。
由于阻尼參數和平滑權重的大小對反演結果的穩定性有較大影響,因此,對不同平滑權重和阻尼參數進行了權衡分析,構建不同平滑權重、阻尼參數的解的方差與數據方差之間的均衡曲線,圖4給出了2個時間段的均衡曲線,2個時間段的平滑權重和阻尼參數的搜索范圍一致,分別為0~600、10~2000。圖4中所示的解的方差和數據方差僅有相對意義,解的方差包含了事件定位和速度模型參數的影響。由圖4可見,對于不同的平滑權重,最佳阻尼參數約為750;由阻尼參數為750時的一系列平滑權重的數據方差與速度模型方差的均衡曲線得到的最優平滑權重約為150。

表1 一維P波速度和波速比模型表
2個時間段的數據經過11次迭代后,第1個時間段2008年12月16日~2011年6月30日的數據在最后一次反演后拾取的到時殘差(RMSCT)從0.36s降到0.15s;第2個時間段2011年7月1日~2016年12月30日的從0.28s降到0.06s。圖5為2個時間段定位模型的殘差變化圖,以0.1s為間隔統計所有地震的均方根殘差。由圖5可見,初始均方根殘差分布較為分散,最后均方根殘差則緊密地分布于0s附近。由于初始定位結果采用簡單的一維速度模型并且用于定位的數據類型較少,僅使用了絕對到時數據,因此地震定位精度較低;而使用地震波絕對、相對到時數據聯合反演震源參數以及速度結構會使得反演后的地震定位及速度模型精度都有了顯著的改進。

圖4 平滑權重和阻尼參數均衡曲線

圖5 蓄水后2個時間段的殘差變化
從研究區域(圖1(b))2個時間段地震重定位結果來看(圖6),研究區內主要有 A、B、C、D、E等地震聚集區域,A區為黑惠江,B區為小灣水庫回水瀾滄江段,C區為瀾滄江昌寧段,D區為施甸一帶,E區為距瀾滄江僅11km的2015年10月30日昌寧MS5.1地震余震序列分布區。第1個時間段的地震主要分布在 A、B、C、D區,由姜金鐘等(2016)及柯乃琛(2016)的研究結果可知,A、B、C三個區域蓄水前地震活動頻次不高,蓄水后頻次開始增加,并隨水位漲落而起伏。D區地震大多為2010年6月1日施甸ML4.8地震余震序列,且D區的地震一直以來都較多,與水庫蓄水前后水位變化間的關系不明顯。由圖6(c)可見,A、B、C三個區域大部分地震的震源深度不超過10km;D區地震震源深度為0~30km,與姜金鐘等(2016)的結果相符。

圖6 重定位后地震分布
第2個時間段的地震分布與第1個時間段大致相同,分別是A、B、C、D、E區,D區地震一直較多,其中包含了2012年9月11日施甸MS4.7地震余震序列,E區主要為距瀾滄江僅11km的2015年10月30日昌寧MS5.1地震序列,該地震為水庫蓄水后發生的震級最大的地震。由圖6(b)可見,A、B、C區地震的震源深度大多小于10km,與第1個時間段差不多,但A區發生的2016年1月4日ML4.0地震的震源深度約為15km。C區昌寧MS5.1地震的余震序列震源深度大多小于10km。D區施甸MS4.7地震余震序列的震源深度為0~30km。2個時間段內小灣水庫庫壩區地震活動均未有明顯增強,表明小灣水庫大壩位置地質構造較為穩定。
由姜金鐘等(2016)的研究結果可知,A、B、C三個區域地震的震源深度平均值均小于10km,D區地震的震源深度平均值大于10km,小于15km,姜金鐘等(2016)將此結果與天然構造地震的深度結果(張國民等,2002;楊智嫻等,2004;Yang et al,2005)進行對比,得出 A、B、C三個區域存在明顯的水庫觸發地震,而D區域的地震則屬于構造地震。
利用棋盤測試方法評價成像結果的分辨率。在初始一維速度模型上分別增加±6%的正負相間速度擾動,從而獲得初始三維速度模型,根據實際地震數據正演計算合成走時數據集,使用合成走時數據集進行雙差地震層析成像,通過對反演結果中每一個節點的速度與理論值進行對比,了解成像結果的分辨率。對2個時間段都進行了棋盤測試,由于20、25、30km深度處的反演結果基本為空白或失真,所以分辨率測試只選取了0、3、5、10、15km剖面。圖7為2個時間段的棋盤測試結果。由圖7可見,由于不同深度處射線的分布不同,因此,0km深度處恢復不好;3km深度處的中心區域恢復較好;5、10、15km深度處恢復都較好。
選擇0、3、5、10、15km深度處水平層面的速度分布進行討論,對于2.1節所討論的地震集中分布的區域將著重分析。由0km深度水平層面P波速度結構(圖8)可知,2個時間段研究區的速度異常與地表地形起伏之間有密切關系,高速區基本沿著瀾滄江分布,而瀾滄江兩岸均為高山。
由圖8可見,2個時間段3km深度處地震均集中分布在 A、B、D區,D區第1個時間段地震大多為2010年6月1日施甸ML4.8地震的余震序列,第2個時間段地震分布較為離散,其中包括2012年9月11日施甸MS4.7地震的余震序列。所有這些地震集中分布的區域均為高速區,第1個時間段該高速區的 P波速度約為4.5km/s,小灣水庫庫區的 P波速度為3.8~4.0km/s;第2個時間段地震分布集中區域的P波速度為4.2~4.5km/s,水庫庫區P波速度為 4.2~4.5km/s。
2個時間段5km深度處地震仍然主要分布于A、B、D區,D區2個時間段的地震仍然主要是2次施甸地震的余震序列,第2個時間段E區地震主要為距瀾滄江僅11km的2015年10月30日昌寧MS5.1地震的余震序列。A區、B區和庫壩區第1個時間段P波速度為5.8~6.0km/s,第2個時間段 A區和庫壩區為5.6~5.8km/s,B區為5.9~6.0km/s;C區第1個時間段速度為6.0~6.1km/s,第 2個時間段為 5.8~6.0km/s;D區第 1個時間段速度為 5.9~6.0km/s,第 2個時間段為 5.8~5.9km/s。

圖7 不同深度處的P波速度棋盤測試結果
2個時間段10km深度處A區基本為低速區,P波速度均為5.8~6.0km/s;庫壩區及B區為高速區,B區第1個時間段P波速度為6.2~6.4km/s,第2個時間段為6.1~6.4km/s,庫壩區第1個時間段P波速度約為6.6km/s,第2個時間段為6.2~6.4km/s;C區第1個時間段 P波速度約為5.8km/s,第2時間段約為5.7km/s;D區第1個時間段 P波速度為6.1~6.3km/s,第2個時間段為5.8~6.2km/s;E區2個時間段P波速度均約為6.0km/s。

圖8 不同深度的P波速度、P波速度差和地震震中分布
至15km深度處,D區2個時間段地震分布仍相對較多,而0~10km深度地震分布集中的A、B、C區在15km深度處卻很少有地震發生。該深度處第1個時間段A區P波速度為5.9~6.0km/s,第 2個時間段為 6.0~6.2km/s;B區第 1個時間段 P波速度約為5.8km/s,第 2個時間段約為6.4km/s;庫壩區2個時間段P波速度為6.1~6.3km/s;C區第1個時間段P波速度約為5.8km/s,第2個時間段約為6.2km/s;D區第1個時間段 P波速度為6.0~6.2km/s,第2個時間段為 6.2~6.4km/s;E區第 1個時間段 P波速度約為5.6km/s,第2個時間段為5.6~5.8km/s。
為了對比2個時間段研究區內P波速度的變化,將第2個時間段的速度減去第1個時間段的速度得到圖8(c)所示結果。由圖8(c)可見,A區在0、3、5、10km深度處 P波速度均為下降,本文認為這可能是由于黑惠江是一條年輕的支流,蓄水后庫水沿著該區域巖石破碎形成的通道向四周進行滲透,使得周邊介質空隙含水率增加,強烈地改變了介質的性質,造成了此處P波速度呈現低速異常變化。對于B區,0~10km深度內,在地震分布的區域,P波速度隨深度增加而降低,至5km深度處,速度降低幅度很小,至10km深度處,基本為P波速度上升,可能是由于該區域地下巖體主要為二長花崗巖,庫水不易滲透,最多只能滲透至5km深度。C區在0、3、5、10km深度處P波速度均為下降,C區位于瀾滄江斷裂與柯街斷裂的交界處,由柯乃琛(2016)對小灣水庫蓄水水位與地震活動性間的相關性研究可知,蓄水回水位沿著瀾滄江超過了C區后,C區的地震大量增加,說明該區域P波速度降低的原因是庫水沿瀾滄江斷裂與柯街斷裂的巖石破碎帶向下滲透從而改變了地下介質的性質。D區在0、3、5、10km深度處P波速度均為下降,但是D區遠離瀾滄江,蓄水前后地震活動都較活躍,故本文認為該區域構造活動較為活躍,與水庫蓄水無關。E區在0、3、10km深度處P波速度降低,5km深度處地震主要分布于速度降ΔvP上升的區域,E區位于柯街斷裂上,E區發生的地震與蓄水回水無關,初步判定是構造地震。水庫庫壩區的P波速度變化沒有規律,表明了該區域地質構造的穩定性。
圖9為 A-A′、B-B′、C-C′、D-D′、E-E′剖面的 P波速度結構分布。由圖9可見,2個時間段內地表低速層基本不變,A-A′剖面上的小灣水庫回水瀾滄江段的地震集中區域在第1個時間段震源深度基本在5km以內,第2個時間段的地震震源深度大于5km,小于10km,P波速度未有太大的變化。B-B′和C-C′剖面均橫穿過黑惠江上的地震集中區域,2個剖面上地震聚集區的震源深度在第1個時間段內均小于5km,而在第2個時間段內大于 5km,小于10km,在地震聚集區的P波速度稍微變小。D-D′剖面上紅色箭頭所標注地即為瀾滄江保山段上地震聚集區域,該地震聚集區的震源深度在2個時間段內均很淺,大多并未超過10km,第2個時間段震源聚集區域低速區變厚。E-E′剖面上紅色箭頭所標注地為2015年10月30日昌寧MS5.1地震余震區,在第1個時間段內該區域并沒有地震,余震深度小于10km,P波速度基本無變化,第1個時間段內E-E′剖面西端的地震為2012年9月11日施甸MS4.7地震的一部分余震,該區域波速相對較低。
本文認為由不同地震聚集區在2個時間段的震源深度及P波速度變化可得到一些認識,即黑惠江上的地震集中區(A區)和小灣水庫回水瀾滄江段的地震集中區(B區),隨著蓄水時間的延長,庫水沿著破碎帶向下滲透,地震震源深度隨之增加。瀾滄江保山段上地震聚集區(C區)的震源深度并未隨蓄水時間延長而增加,2個時間段的震源深度分布相差不大,表明該區域的地下通道和巖性均更利于庫水滲透,蓄水回水至該段后,庫水滲透至較深處,而在10km深度以下地質構造較穩定,故庫水難以向下滲透。

圖9 P波速度沿垂直剖面深度分布
已有研究認為,只有當水體和庫區的地質構造條件、水文地質條件、巖性條件和應力條件等因素適當地結合后才有可能觸發較大的地震(丁原章等,1989;李安然等,1992;Simpson et al,1988)。水庫水體的加載和滲透主要對庫底巖石及斷層產生3種作用(陳翰林等,2009),即彈性效應、孔隙壓變化及斷層弱化等,其中,水庫蓄水載荷會導致斷層面上正應力和剪應力增大,正應力增大會使得深部斷層增強,而剪應力增加,斷層變弱還是變強則取決于斷層走向與區域應力場方向之間的關系。水體載荷的影響深度非常淺,載荷和滲透引起的孔隙壓變化則首先增加淺部的孔隙壓,繼而使得孔隙壓向深部擴散,從而改變介質的強度和斷層的摩擦阻力。所以,結合地震的震源深度和速度變化可初步判斷不同的地震聚集區在小灣水庫蓄水后地震增加的原因。在本文中,基于已有的認識可知:
(1)A、B區的地震在第1個時間段內震源深度小于5km,在第2個時間段內大于5km,但并未超過10km,第2個時間段的P波速度相比第1個時間段下降,說明這2個區域的地震與水庫蓄水有關,是由于水體滲透導致孔隙壓變化,并隨著時間的變化孔隙壓變化朝著更深的部位擴散,從而導致介質變化,P波速度降低。
(2)C區的地震震源深度小于10km,大多在5km以內,并未隨時間而向深部發展,第2個時間段P波速度相比第1個時間段下降,在蓄水回水位到達該區域后地震開始增多,說明C區地震的發生是由于庫水的滲透作用,該區域地下10km深度以上滲透通道有利于庫水的滲透,而10km深度以下地質構造較為穩定。
(3)小灣水庫蓄水后,在水庫大壩附近并沒有明顯的地震增多現象,2個時間段內的P波速度也無規律性的變化,表明了小灣水庫大壩位置地質構造的穩定性。
(4)D區歷來地震多發,由震源深度以及速度結構變化可知,震源深度為0~30km,P波速度變化無規律,說明蓄水后施甸一帶發生的地震,如2010年6月1日施甸 ML4.8、2012年9月11日施甸MS4.7等地震均屬于構造地震,與水庫蓄水關系不大。
(5)發生在E區的2015年10月30日昌寧MS5.1地震及其余震震中位于柯街斷裂與瀾滄江斷裂的交匯處,震源深度為0~10km,大多在5km以內,由主震震源機制解結果①云南省地震局監測中心,2014,大震應急數據產品產出判斷,NE向柯街斷裂為發震斷裂,具正斷性質。在第1個時間段內該區域基本無地震,2個時間段內該區域的P波速度上升。初步判定該地震是與小灣水庫蓄水有關的構造地震,但需進一步對該震群單獨作計算分析才能得到更準確的結果。
致謝:中國科技大學的張海江教授為本研究提供了tomoDD程序,兩位審稿專家為本文提出了寶貴的修改意見,在此一并表示感謝。