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

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

0 引言

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

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

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

1 資料選取與計算

1.1 資料選取

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

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

1.2 數(shù)據(jù)計算

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

(1)

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

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

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

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

2 計算結(jié)果與極化特征分析

2.1 時序特征分析

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

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

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

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

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

2.2 逐日相關(guān)分析

對尼泊爾地震震中距≤1500km地磁臺站的極化值進行了計算分析， 如圖 7 所示。由于多數(shù)臺站從2015年開始觀測， 故統(tǒng)一從2015年1月開始分析， 全部使用GM4儀器數(shù)據(jù)的計算結(jié)果。從各臺極化曲線中可以看出， 在尼泊爾地震前98d開始出現(xiàn)極化值超閾值的變化， 多數(shù)臺站異常持續(xù)約2d， 且部分臺站在2月13日也出現(xiàn)了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年的滑動平均極化結(jié)果

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

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

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

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

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

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

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

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

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

3 結(jié)論與討論

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

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

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

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

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

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