2015年尼泊爾MW7.8地震前磁擾動極化異常特征分析

管貽亮 董曉娜 尹玉振 馮麗麗 殷海濤*

1)山東省地震局， 濟南 250014 2)青海省地震局， 西寧 810001

0 引言

地震的孕育和發生往往伴隨著產生電磁輻射， 目前認為超低頻(ULF)信號中可能包含較多的震磁異常信息。國外學者在1988年Spitak地震(Molchanovetal.， 1992； Kopytenkoetal.， 1993)、 1989年Loma Prieta地震(Fraser-Smithetal.， 1990)和1993年Guam地震(Hayakawaetal.， 1996)等多個地震震前發現了電磁輻射異常。伴隨著孕震過程， 能量不斷積累， 產生的感應磁場信號也不斷增強。但相對于強大的背景場而言， 變化磁場較為微弱， 且ULF信號從震源區到地表不斷衰減， 增加了電磁輻射異常信號的提取難度。研究發現， 來自巖石圈的電磁信號主要在垂直分量上體現， 而外源場信號主要體現在水平分量上(Molchanovetal.， 1995a， b)， 極化法利用二者的比值進行震磁異常信息的提取。該方法最早由Hayakawa等(1996)提出， 其發現1993年關島MS8.0 地震前關島臺的數據出現了頻帶范圍為0.002～0.005Hz、 幅值約為0.1nT的輻射異常。Hattori等(2002)發現在2次Kagoshima地震前存在明顯的極化值增高現象； Molchanov等(2003)基于俄羅斯Kamchatka半島建立的綜合觀測臺分析了M4～6.2震群前短時間(2～6d)的極化值波動情況， 并認為極化值與空間磁場活動無相關性； Hobara等(2004)對0.0017～2Hz頻段的極化值進行了精細化分析， 在Izu和Guam 2個地震震前提取到了0.01～0.022Hz和0.05～0.1Hz超正常值約2倍的極化異常。近年來， 國內學者先后對新疆(馮志生等， 2010； Lietal.， 2011； 張敏等， 2020)、 青海(廖曉峰等， 2019； 李霞等， 2021)和云南(李琪等， 2015； 樊文杰等， 2018)等地區的資料進行了極化分析， 在震例積累方面成果頗豐。

雖然國內外學者對極化特征及其映震效果展開了一系列探索并取得了不錯的效果， 但目前對該方法依舊存在一定的質疑。爭議之一集中于極化異常與空間磁場變化存在一定的關系， 而與地震活動關聯不大。Thomas等(2009)和Masci(2011)對地磁水平和垂直分量與地磁K指數(以0～9表征所在區域的地磁場擾動程度，K值越大表明磁場空間活動越強)的相關性進行了研究， 認為極化值增高與外部空間磁場活動有關， 而與地震活動無關。爭議之二是受限于數據采樣率、 臺站數量等， 相關研究多是基于單一臺站的時序分析， 由于計算參數不統一、 異常提取標準不一致和分析數據時間短等， 導致異常信度較低。一系列的爭議也限制了該算法的應用和發展， 而從能量積累、 信號衰減、 空間干擾等角度分析， 小地震引起的電磁輻射異常可被觀測到或提取出的可能性較小， 因此對大地震前的極化異常分析更有可信度。本研究在盡可能保證數據信度(相同儀器型號、 統一計算參數、 多臺站、 長時間尺度等)的前提下， 試圖從時間和空間2個維度對極化異常特征進行研究， 采取定量的方法深度分析極化異常與地磁活動、 地震活動的關系。

2015年4月25日14時11分在尼泊爾(28.15°N， 84.71°E)發生MW7.8地震， 震源深度8.2km(USGS)， 震源機制結果顯示該地震為低角度逆沖型地震。該地震的震中位于喜馬拉雅造山帶， 這一地帶長約2500km， 呈EW向展布， 是印度板塊和歐亞板塊碰撞和擠壓的產物， 長期處于構造強烈活動的狀態， 也是目前陸內地震多發區。本次地震引發了多次中強度余震， 且均為淺源地震， 對尼泊爾地區、 中國西藏地區以及周邊國家造成了重大的人員傷亡和財產損失。然而研究表明， 該地震所釋放的能量遠低于空區理論預測的釋放能量(劉靜等， 2015)， 說明該地區的地震危險性依舊很高， 需要利用多種手段和方法對該地區開展地震活動性研究。諸多學者從構造活動(楊曉平等， 2016)、 發震背景(鄧起東等， 2014； 滕吉文等， 2017)以及應力變化(占偉等， 2015； 萬永革等， 2017； 王東振等， 2017； 張迎峰等， 2017)等方面對該地震進行了分析， 但對于地震前后的空間磁場變化情況研究較少， 尤其是缺乏深層次的特征分析。本文以震中附近臺站的極化異常特征和同步性轉折為切入點， 對地震孕育發生過程中的震磁關系進行研究， 基于中國地磁臺站的極化數據對ULF信號的異常輻射范圍進行空間特征分析， 從地質構造和動力學演化方向對震磁極化異常演化過程進行了探索性研究。

1 資料選取與計算

1.1 資料選取

近年來， 隨著地磁臺網規模不斷擴大， 地震監測能力得到有效提升的同時也產出了大量的秒值數據， 為提取ULF信息提供了很好的數據基礎。為了研究尼泊爾地震前構造活動可能引起的區域磁場波動， 本文對中國部分地磁臺站數據進行了極化計算分析。其中， 震中距拉薩地磁臺628km， 距獅泉河臺652km， 且2個臺站均有GM4和M15 2種型號的儀器， 數據質量較好， 故重點對其進行了時序異常特征的提取分析。在對中國地磁臺站數據進行計算時， 為了減少單一臺站數據畸變對空間結果的影響， 基于數據有效性進行了數據篩選和質量提升。臺站數據連續率 <20% 的資料不參與計算， 對突跳數據(>5倍均方差)進行插值處理。經分析， 2014—2015年中國大陸地區地磁臺站的數據中符合計算標準的臺站共有48個， 分布如圖 1 所示。從圖中可以看出， 中西部臺站相對較多， 尤其是川滇地區較為密集， 東北和東南地區地磁臺站相對較少。

圖 1 參與計算的地磁觀測臺站空間分布Fig. 1 The spatial distribution of geomagnetic observation stations used in the calculation.

1.2 數據計算

極化法是在譜分析的基礎上， 將磁場垂直分量和水平分量的頻譜振幅相比(Hayakawaetal.， 1996)， 具體定義為

(1)

其中，Yzh為極化值，Z(ω)為地磁垂直分量的譜，H(ω)為地磁水平分量的譜， 由NS向和EW向的譜值計算得出。

圖 2 拉薩臺2015—2017年極化各階段的計算結果Fig. 2 The polarization calculation results of each stage of Lasa station from 2015 to 2017.a 地磁K指數日值曲線； b 垂直分量頻譜曲線； c水平分量頻譜曲線； d 極化值日值曲線與富氏擬合曲線； e 5日滑動極化值與5日滑動K指數對比曲線

計算得到的極化日值與K指數的相關系數為-0.38， 呈較明顯的負相關， 說明極化高值并不是由空間磁場活動所引起的。為便于提取高值異常， 以富氏擬合值+2倍均方差為閾值對低比值進行剔除， 然后再次計算極化日均值。從長時間尺度來看， 極化日均值存在明顯的年變特征， 當數據時間較短時這種年變趨勢會影響對異常的判斷甚至產生假異常， 結果如圖2d所示。作者認為這也可能是Hayakawa等(1996)在得出 “震前極化值趨勢性增大”結論時受到爭議的原因之一， 遂對其做富氏擬合去除年變， 去除年變后計算5日滑動平均值以濾除短期噪聲。5日滑動極化值與5日滑動K指數的相關系數為-0.46， 呈比較明顯的負相關， 如圖2e所示。基于對以上計算過程的認知和計算結果的逐步分析， 認為極化法能夠抑制外源場信號， 突出與震源有關的局部磁異常信息， 可進行震磁異常信息提取， 詳細計算流程如圖 3 所示。

圖 3 極化值計算流程Fig. 3 The flow chart of polarization calculation.

2 計算結果與極化特征分析

2.1 時序特征分析

分別計算了距離震中最近的拉薩和獅泉河2個臺站GM4和M15儀器的極化值， 如圖 4 和圖 5 所示。從極化時間序列中可以看出： 距離震中最近的拉薩臺在2015年1月17日(震前98d)出現超閾值(2倍均方差， 下同)異常， 異常持續3d； 獅泉河臺在1月18日出現超閾值現象。查看2個臺站極化異常時段的原始數據和觀測日志發現， 數據沒有受到明顯的外部干擾。極化異常期間2個臺站的K指數處于明顯的低值區， 說明該時段外空磁場變化比較平靜， 可排除外部干擾和空間磁場活動的影響， 此次同步性異常可能是尼泊爾震前區域磁場活動的真實反映。

圖 4 拉薩臺的極化結果與地磁K指數曲線對比圖Fig. 4 The comparison chart of polarization results and K index of Lhasa station.a GM4儀器計算的極化值； b M15儀器計算的極化值； c 拉薩臺的地磁K指數

圖 5 獅泉河臺的極化結果與地磁K指數曲線對比圖Fig. 5 The comparison chart of polarization results and K index of Shiquanhe station.a GM4儀器計算的極化值； b M15儀器計算的極化值； c 獅泉河臺的地磁K指數

根據滑動平均結果能更加清晰地對數據變化趨勢進行分析， 如圖 6 所示。2個臺站的極化值在2015年1月18日均出現極化超閾值現象， 其中拉薩臺的異常幅度和異常持續時間均大于獅泉河臺。2個臺站的極化值在2015年2月13日前后也出現了極化值增大的現象， 部分儀器的極化值超過了閾值， 但整體幅度明顯低于18日。分析認為， 此次上升可能是18日磁場變化的延續， 極化異常是成組出現的， 反映了震前磁場一段時間內的變化。對比相同臺站不同儀器的極化結果發現， M15儀器的極化背景值和異常幅度均高于GM4儀器， 其異常靈敏度更優。因此在進行多臺分析， 尤其是定量分析時， 使用同型號儀器的數據更有說服力。

圖 6 5日滑動平均極化結果Fig. 6 The polarization results of 5-day moving average.a 拉薩臺GM4儀器； b 拉薩臺M15儀器； c 獅泉河臺GM4儀器； d 獅泉河臺M15儀器

2.2 逐日相關分析

對尼泊爾地震震中距≤1500km地磁臺站的極化值進行了計算分析， 如圖 7 所示。由于多數臺站從2015年開始觀測， 故統一從2015年1月開始分析， 全部使用GM4儀器數據的計算結果。從各臺極化曲線中可以看出， 在尼泊爾地震前98d開始出現極化值超閾值的變化， 多數臺站異常持續約2d， 且部分臺站在2月13日也出現了1次同步性的高值異常， 異常幅度略低于18日， 與拉薩臺和獅泉河臺的變化基本一致。

圖 7 尼泊爾MW7.8地震前震中附近臺站的極化異常時序圖Fig. 7 The timing chart of polarization abnormality of nearby stations before the Nepal MW7.8 earthquake.道孚臺(a)、 肅北臺(b)、 格爾木臺(c)、 都蘭臺(d)、 喀什臺(e)和且末臺(f)2015年的滑動平均極化結果

多臺極化值曲線具有同步性變化， 反映了磁場變化的空間輻射性。這種時間尺度上的一致性是空間磁場變化的正常體現， 而在幅度上的差異性可能隱藏了一定的異常信息。以拉薩臺和獅泉河臺為例， 分別計算了震前和震后98d(異常距離發震時間)時間序列的相關系數， 震前數據的相關系數是0.88， 震后為0.72， 可以看出2個臺站震后數據變化的同步性明顯低于震前數據的同步性。以此為出發點， 為了更加直觀地顯示極化值相關性隨時間的變化關系， 計算了上述兩兩臺站間的逐日滑動相關系數。具體計算過程為： 以一定時間長度作為窗長(N)計算得到第1天的相關系數， 窗長不變向后逐日滑動計算每天的系數， 即用第M-N到M天的數據計算第M天的相關系數。本研究中選取N=60為基準時間， 盡可能在保證總體數據穩定性的同時突出樣本的差異性。為了保證異常期間的數據穩定， 從2014年10月開始計算， 相關系數的起始時間為2015年12月1日。由于道孚臺和都蘭臺自2015年開始觀測， 極化異常期間無法計算相關系數， 故不作分析。對拉薩臺、 獅泉河臺、 肅北臺、 格爾木臺、 喀什臺和且末臺6個臺站的結果進行分析， 如圖 8 所示， 每個子圖分別為該臺站與其余5個臺站的逐日滑動相關系數曲線。 為了便于分析， 對每個子圖中的5條曲線取均值形成1條曲線， 如圖 9 所示。

圖 8 臺站間逐日滑動相關曲線Fig. 8 The daily correlation curve between the individual stations.拉薩臺(a)、 獅泉河臺(b)、 肅北臺(c)、 格爾木臺(d)、 喀什臺(e)和且末臺(f)分別與其余5個臺站的逐日滑動相關曲線

圖 9 臺站間逐日滑動相關均值曲線Fig. 9 The daily mean correlation curve between the individual stations.拉薩臺(a)、 獅泉河臺(b)、 肅北臺(c)、 格爾木臺(d)、 喀什臺(e)和且末臺(f)分別與其余5個臺站的逐日滑動相關均值曲線

從相關曲線中可以看出， 多數臺站在極化異常后出現逐日相關增強的現象， 其中拉薩臺、 喀什臺、 格爾木臺和且末臺4個臺站變化明顯。這種發生在地震孕育過程中的多臺間相關性增強的現象， 可能與地震能量積累過程中所引起的空間磁場增強有關。其后在尼泊爾地震發生前幾天， 相關性出現了急劇的下降—上升轉折， 拉薩臺、 喀什臺和且末臺的變化尤為明顯， 從曲線上看地震發生在轉折上升期間， 這種高度的一致性轉折可能代表了地震應力釋放過程中的磁場變化。

分析各臺站的地磁K指數和原始數據， 在極化異常到地震發生這段時間磁場活動均比較平靜， 只有2015年3月17日最大K指數為7， 各臺原始數據均呈平穩變化， 故相關性出現轉折的原因可排除磁場波動和數據畸變的影響。分析認為震前可能存在外場磁信號對區域磁場觀測產生了影響， 該信號強度較大， 影響了區域內正常的磁場變化， 使不同臺站出現了較強的同步性變化， 即震前的高同步性可能是與震源相關的外場信號施加造成的， 是地震孕震過程中磁場變化的真實反映。鑒于逐日相關計算是一個逐日累積計算的過程， 每天的數據都影響整體的結果， 將其用作異常判據需要滿足一定的條件， 即數據不可存在大范圍缺數， 此外需要排除高壓直流、 磁暴等確定的外部干擾。

2.3 空間特征及動力演化分析

本文計算了中國范圍內符合要求的48個地磁臺站的極化值相對變化量， 選取震中附近25個臺站的數據進行空間特征分析。對2015年1月18日的異常進行了時間和尺度上的分析， 發現該組異常主要從17日開始在部分臺站出現， 持續至19日， 具體如圖 10 所示， 其中紅色的零值線為空間異常閾值線。17日有10個臺站(40%)的極化值超過閾值， 伴隨著地震的孕育、 震磁信號的傳播， 異常區域隨著時間的推移不斷擴大， 18日異常空間范圍和異常幅度均增大， 18個臺站(72%)超過閾值。到19日異常幅度逐漸降低， 空間范圍有所縮小， 15個臺站(60%)超閾值， 20日異常基本消失。異常變化過程中空間范圍和異常幅度均呈現增長—減少的變化趨勢， 但整個異常空間分布始終圍繞震中。該異常變化在時間上的推移和空間上的一致性是區域磁場變化的真實反映， 存在比較明顯的外部作用。而引起該變化的外部場具有很強的動力學特征， 可能與構造活動和地殼應力變化有關。大范圍的異常與尼泊爾地震的強度也是適應的， 下文將嘗試從區域構造活動和地震的孕育過程角度對極化異常的時空演化特征進行分析。

圖 10 極化逐日空間異常分布圖Fig. 10 The spatial distribution map of polarization abnormality.a 2015年1月17日極化異常空間分布圖； b 2015年1月18日極化異常空間分布圖； c 2015年1月19日極化異常空間分布圖

從異常空間分布圖中可以看出， 空間異常主要呈SE-NW向展布， 異常分布區域與喜馬拉雅構造帶主峰北側青藏高原構造分布基本一致， 如圖 1 所示。震前磁場異常可能反映了區域構造的動力學演化過程。對地震孕育過程所引起的電磁信號而言， 無論用何種震磁效應解釋， 磁場信號的產生和傳播必然與地質結構及構造活動有著密不可分的關系。尼泊爾地震發生前后， 多位學者對這一區域的地質結構背景(楊曉平等， 2016)、 動力演化特征(鄧起東等， 2014； 滕吉文等， 2017)等進行了研究， 通過重力、 GPS等手段對應力場變化進行了分析(占偉等， 2015； 王東振等， 2017)， 證實了震前印度板塊與歐亞板塊之間存在持續擠壓變形， 這次大地震是活動帶上長期以來深部物質與能量強烈交換、 運動導致應力積累和釋放的產物。大地電磁測深結果也揭示青藏高原中下地殼普遍具有高導特征， 研究認為可能是地殼中巖石部分熔融和含鹽流體共同作用的結果(金勝等， 2010)。綜合前人在該區域的研究成果， 結合前面對震前極化增強的認識， 分析認為本次極化的空間異常與尼泊爾地震前的構造活動有關。板塊運動引起地殼應力變化從而激發區域磁場變化， 而該地區存在大量的熔融物質等流體為磁場的孕育和傳播提供了基礎和通道。反之， 大范圍的磁擾動極化異常也揭示了該地區震前的構造動力學演化過程， 為研究震前構造活動提供了一種地球物理方案。不可否認的是， 臺站的分布密度、 臺基巖性等客觀因素會影響極化異常的空間分布， 尤其是對定量分析造成了一定的限制。雖然能夠在震后對異常演化過程進行推演， 但基于極化空間異常的地震預測依舊存在諸多困難， 如本研究中并沒有發現震中距與極化異常幅度間的量化關系， 后續需要加大對異常信息的再挖掘力度。

3 結論與討論

(1)通過對極化算法的研究可以確認， 該方法能夠抑制外源場信號， 突出與震源有關的局部磁異常信息。本文的研究工作驗證了磁暴等外空磁場變化不會引起極化值的高值異常， 可用于提取震磁異常信息。但所需數據的時長至少應在1a以上， 以保證結果的可靠性。

(2)對尼泊爾MW7.8地震震中1500km范圍內的臺站極化值進行計算分析， 發現震前98d出現了同步性高值異常， 鑒于這種區域同步性變化提出了極化逐日相關分析。臺站間逐日滑動相關分析結果顯示， 尼泊爾地震前存在明顯的極化相關增強現象， 而地震發生在同步性轉折階段， 分析認為該同步性變化可能是區域磁場在震源場外加作用下引起的。

(3)基于極化值相對變化量進行的空間異常分析顯示， 尼泊爾地震前中國有超過20%的臺站出現持續3d以上的異常， 且異常空間分布具有明顯的區域性， 可為后續震前磁擾動異常的判別提供一些啟示。

(4)研究認為極化異常的產生和分布與地質構造的動力學演化密切相關， 地質活動產生的應力積累是引起磁場變化的主要原因， 而區域構造結構則為磁場的發展和傳播提供了通道。因此， 將地球物理場變化與地質特征相結合進行異常分析， 可更加清晰地解釋地震的孕育和發展過程， 提高極化異常的可信度。

(5)震磁機理研究和震前異常信息的提取依舊是震磁研究的難點。本研究中沒有發現震中距與極化異常幅度間的量化關系， 分析認為與站點構造環境、 孕震背景、 空間磁場變化復雜性和震例較少等多方面因素有關。極化方法的應用還需要更多的震例積累和理論分析， 以獲取更多的證據， 提高可靠性。

致謝馮志生研究員、 朱培育、 李霞、 賀曼秋、 廖曉峰、 何暢、 梵文杰、 艾薩·伊斯馬伊力、 袁文秀等在方法研究和數據計算方面提供了幫助和指導； 審稿專家為本文提出了寶貴的修改意見。在此一并表示感謝!