基于速率累加法提取地磁諧波振幅比異常指標(biāo)

李 霞 馮麗麗 趙玉紅 劉 磊 茍智棟

1 中國地震局青海格爾木青藏高原內(nèi)部地球動力學(xué)野外科學(xué)觀測研究站，西寧市柴達(dá)木路242號，810017 2 青海省地震局，西寧市柴達(dá)木路242號，810017

地震孕育的過程是應(yīng)力應(yīng)變在局部地區(qū)(震源體及附近)緩慢積累的過程，而巖石破裂實驗表明，巖石電導(dǎo)率與應(yīng)力應(yīng)變關(guān)系密切[1-4]，因此地震在孕育過程中必然伴隨著巖石電導(dǎo)率的變化。地磁場的內(nèi)場強度既和外場強弱有關(guān)，又和地下物質(zhì)電導(dǎo)率的大小有關(guān)，所以由應(yīng)力累積引起的地下物質(zhì)電導(dǎo)率的變化必然會導(dǎo)致外場和內(nèi)場比值發(fā)生改變，這種改變是可以通過地磁日變化的諧波振幅比反映出來的[5]。馮志生等[6-8]對不同地震震前地磁諧波振幅比異常變化特征進(jìn)行研究發(fā)現(xiàn)，地磁諧波振幅比異常特征與地電阻率趨勢性異常特征具有相似性。

隨著觀測數(shù)據(jù)的累積和地磁諧波振幅比方法的推廣，研究發(fā)現(xiàn)，地磁諧波振幅在地震前確實存在異常[9-14]。但該方法應(yīng)用于日常震情跟蹤工作時存在以下2個問題：1)目前分析主要采用分鐘值采樣數(shù)據(jù)，不能反映淺層介質(zhì)的性質(zhì)變化；2)無定量化指標(biāo)。近年來，速率累加方法在跨斷層水準(zhǔn)觀測資料的前兆異常識別中已經(jīng)得到應(yīng)用，該方法可以使某時段內(nèi)較為頻繁的速率異常變化更加直觀和突出，更易于判定異常和提取異常指標(biāo)[15]。

本文利用諧波振幅比法分析青海地區(qū)磁通門秒采樣觀測資料，利用速率累加方法分析連續(xù)數(shù)據(jù)的逐日變化速率，再根據(jù)需要對分析時段進(jìn)行累加計算，識別數(shù)據(jù)資料的異常變化，進(jìn)一步分析異常變化與強震的關(guān)系，梳理總結(jié)青海地區(qū)磁通門秒值諧波振幅比方法的異常指標(biāo)，并將研究結(jié)果應(yīng)用于青海地區(qū)的震情跟蹤中。

1 主要方法與計算數(shù)據(jù)

1.1 方法與計算步驟

1.1.1 地磁諧波振幅比方法

地磁諧波振幅比YZHx(NS) 和YZHy(EW)的定義為：

(1)

式中，Z(ω)、Hx(ω)、Hy(ω)分別為地磁場的垂直分量、水平NS向分量和WE向分量的振幅譜，可利用地磁三分量的min采樣數(shù)據(jù)計算獲得。對于隨時間發(fā)生周期變化的不均勻場源，在地球介質(zhì)為均勻各向同性平面導(dǎo)體的條件下，有：

(2)

式中，μ為磁導(dǎo)率，σ為電導(dǎo)率，ω為圓頻率，λ為變化磁場的波數(shù)。

(3)

地磁諧波振幅比可以表達(dá)為磁測深視電阻率：

(4)

由式(4)可知，YZHx(NS)和YZHy(EW)與地下介質(zhì)電阻率ρ呈正比。具體計算步驟為：1)計算每天地磁三分量采樣數(shù)據(jù)譜幅度，計算各段各周期諧波振幅比(復(fù)數(shù)余)；2)按照虛實部分別計算各周期的諧波振幅比96段的均值及相位；3)讀入逐頻值文件，按照指定頻帶分別計算諧波振幅比虛實部均值和幅度；4)讀入諧波振幅比逐日值文件，按一定窗長逐日滑動平均，消除噪聲；5)讀入消除噪聲的諧波比逐日值文件，按365 d窗長逐日滑動平均，消除年變。

1.1.2 速率累加方法

速率累加方法是對地磁諧波振幅比計算結(jié)果進(jìn)行二次計算的數(shù)學(xué)方法，其目的是通過二次計算的結(jié)果圖像更加直觀地反映地磁諧波振幅比異常，同時能夠獲得清晰的地磁諧波振幅比異常特征[11]。具體計算步驟為：1)計算地磁諧波振幅比南北向和東西向各周期逐日值的一階差分，獲得各周期一階差分序列，該序列可反映其變化速率；2)對一階差分序列進(jìn)行正負(fù)判斷，大于0取值為1，小于0取值為-1，0取值不變，形成歸一化序列；3)將南北向與東西向各周期歸一化序列對應(yīng)數(shù)值相乘，形成歸一化序列的積序列；4)將各周期積序列進(jìn)行累加計算，獲得各周期速率累加序列，該序列可用于判斷對應(yīng)周期是否存在不同步異常；5)將所有周期累加序列對應(yīng)數(shù)值求和，得到所有周期速率累加和，用于分析不同步異常的起始時間。

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

2015年前青海地區(qū)共有5個地磁觀測點，分別為格爾木、都蘭、德令哈、大武和西寧，“十二五”背景場項目新增設(shè)金銀灘、貴德和玉樹巴塘3個觀測點，8個地磁觀測點共有16套地磁三分量觀測儀器(其中磁通門秒采樣觀測儀14套)，已積累至少5 a的分鐘值采樣及秒采樣連續(xù)觀測資料。分析研究區(qū)5級以上地震與地磁諧波振幅比曲線變化特征的關(guān)系，針對不同方向或不同周期出現(xiàn)的不同步變化進(jìn)行科學(xué)總結(jié)，從時間、空間、強度方面總結(jié)規(guī)律特征。

2 結(jié)果分析

格爾木臺FHDZ-M15[1]諧波振幅比曲線自2016-05起部分速率累加曲線出現(xiàn)轉(zhuǎn)平或連續(xù)波動變化，170 s和260 s的周期曲線首先出現(xiàn)異常變化，自2016-05開始曲線轉(zhuǎn)為平穩(wěn)的波動變化；隨后，290 s的周期曲線在2017年出現(xiàn)轉(zhuǎn)平變化，50 s的周期曲線在2017-10轉(zhuǎn)為平穩(wěn)變化。異常出現(xiàn)24個月后發(fā)生稱多5.3級地震，異常出現(xiàn)27個月后發(fā)生治多5.1級地震，具體見圖1(a)。

圖1 格爾木臺地磁諧波振幅比Fig.1 Geomagnetic harmonic amplitude ratio at Golmud station

格爾木臺GM4[B]諧波振幅比曲線自2016年起部分速率累加曲線出現(xiàn)轉(zhuǎn)平或連續(xù)波動變化，20 s的周期曲線在2016年初出現(xiàn)下降趨勢，同時260 s和290 s的周期曲線在2016年初轉(zhuǎn)為平穩(wěn)的波動變化。隨后，80 s的周期曲線在2017年轉(zhuǎn)平，并呈小幅度下降趨勢；230 s的周期曲線在2017-08轉(zhuǎn)平并呈下降趨勢。異常出現(xiàn)29個月后發(fā)生稱多5.3級地震，異常出現(xiàn)32個月后發(fā)生治多5.1級地震，具體見圖1(b)。

大武臺GM4[8]諧波振幅比曲線自2016-04起部分速率累加曲線出現(xiàn)轉(zhuǎn)平或連續(xù)波動變化，九寨溝7.0級地震前20 s、50 s、80 s、170 s周期曲線均在震前出現(xiàn)轉(zhuǎn)平變化，持續(xù)時間為1.5 a；稱多5.3級地震前20 s、50 s、80 s、110 s、140 s和170 s周期曲線波動變化平穩(wěn)，持續(xù)時間為1 a左右；夏河5.7級地震前20 s、50 s、80 s、110 s、140 s和260 s周期曲線出現(xiàn)轉(zhuǎn)折和下降變化。上述曲線變化的共同特征為周期較小的曲線波動變化更明顯。異常出現(xiàn)16個月后發(fā)生九寨溝7.0級地震，異常出現(xiàn)25個月后發(fā)生稱多5.3級地震，異常出現(xiàn)42個月后發(fā)生夏河5.7級地震，具體見圖2(a)。

圖2 大武臺地磁諧波振幅比Fig.2 Geomagnetic harmonic amplitude ratio at Dawu station

大武臺GM4[b]諧波振幅比曲線自2016年起部分速率累加曲線出現(xiàn)轉(zhuǎn)平或連續(xù)波動變化，九寨溝7.0級地震前20 s、50 s、110 s周期曲線出現(xiàn)顯著轉(zhuǎn)折下降變化，持續(xù)時間為1.5 a左右；80 s、170 s的周期曲線波動變化平穩(wěn)，持續(xù)時間近1 a。稱多5.3級地震前50 s、80 s、110 s、140 s周期曲線變化平穩(wěn)，持續(xù)時間為1 a左右；夏河5.7級地震前20 s、50 s、80 s、110 s、140 s周期曲線出現(xiàn)轉(zhuǎn)折和下降變化，且周期較小曲線波動變化更明顯。異常出現(xiàn)20個月后發(fā)生九寨溝7.0級地震，異常出現(xiàn)29個月后發(fā)生稱多5.3級地震，異常出現(xiàn)46個月后發(fā)生夏河5.7級地震，具體見圖2(b)。

德令哈臺GM4[a]諧波振幅比曲線自2016年起部分速率累加曲線出現(xiàn)轉(zhuǎn)平或連續(xù)波動變化，門源6.4級地震前80 s、140 s、170 s、200 s、260 s周期曲線出現(xiàn)顯著轉(zhuǎn)折下降變化；稱多5.3級地震前20 s、50 s、80 s、110 s、140 s、170 s、200 s、230 s周期曲線變化平穩(wěn)，持續(xù)時間為1 a左右，在此基礎(chǔ)上前50 s、80 s、200 s和230 s周期曲線還出現(xiàn)小幅度的波動變化；甘州5.0級地震前20 s、50 s、80 s、110 s、140 s、200 s、230 s周期曲線波動變化平穩(wěn)，持續(xù)時間為1a，另外290 s周期曲線在震前4個月(2019-05)出現(xiàn)快速下降變化。異常出現(xiàn)29個月后發(fā)生稱多5.3級地震，異常出現(xiàn)45個月后發(fā)生甘州5.0級地震，具體見圖3(a)。

圖3 德令哈臺地磁諧波振幅比Fig.3 Geomagnetic harmonic amplitude ratio at Delhi station

德令哈臺GM4[b]諧波振幅比曲線自2016年起部分速率累加曲線出現(xiàn)轉(zhuǎn)平或連續(xù)波動變化，門源6.4級地震前20 s、140 s、170 s周期曲線出現(xiàn)顯著轉(zhuǎn)折下降變化，震前50 s、80 s周期曲線均處于平穩(wěn)變化階段，地震發(fā)生在曲線轉(zhuǎn)折拐點處；稱多5.3級地震前20 s、110 s、140 s、170 s、230 s、260 s、290 s周期曲線變化平穩(wěn)，持續(xù)時間1.5 a，其中震前20 s周期曲線波動變化顯著；甘州5.0級地震前20 s周期曲線波動變化平穩(wěn)，此外震前110 s、140 s、170 s、230 s、260 s周期曲線震前4個月(2019-05)異常形態(tài)恢復(fù)至正常狀態(tài)。異常出現(xiàn)29個月后發(fā)生稱多5.3級地震，異常出現(xiàn)45個月后發(fā)生甘州5.0級地震，具體見圖3(b)。

都蘭臺GM4[6]諧波振幅比曲線自2016年起部分速率累加曲線出現(xiàn)轉(zhuǎn)平或連續(xù)波動變化，門源6.4級地震前20 s、110 s、170 s周期曲線出現(xiàn)顯著轉(zhuǎn)折下降變化，地震發(fā)生在曲線轉(zhuǎn)折拐點處；稱多5.3級地震前140 s、170 s、200 s、290 s周期曲線變化平穩(wěn)，持續(xù)時間為1.5 a，其中震前170 s與260 s周期曲線波動變化顯著；甘州5.0級地震前20 s、50 s、80 s、110 s、140 s、170 s、200 s、290 s周期曲線波動變化平穩(wěn)，時間持續(xù)1 a以上。異常出現(xiàn)29個月后發(fā)生稱多5.3 級地震，異常出現(xiàn)45個月后發(fā)生甘州5.0級地震，具體見圖4(a)。

圖4 都蘭臺地磁諧波振幅比Fig.4 Geomagnetic harmonic amplitude ratio at Dulan station

都蘭臺GM4[7]諧波振幅比曲線自2016年起部分速率累加曲線出現(xiàn)轉(zhuǎn)平或連續(xù)波動變化，門源6.4級地震前170 s、230 s、260 s、290 s周期曲線出現(xiàn)顯著轉(zhuǎn)折下降變化，地震發(fā)生在曲線趨勢下降過程中；稱多5.3級地震前140 s、170 s、200 s、230 s、290 s周期曲線變化平穩(wěn)，持續(xù)時間為1.5 a，其中震前170 s、200 s與230 s周期曲線波動變化顯著；甘州5.0級地震前110 s、170 s、200 s、230 s、290 s周期曲線波動變化平穩(wěn)，大部分持續(xù)時間1 a以上。異常出現(xiàn)29個月后發(fā)生稱多5.3級地震，異常出現(xiàn)45個月后發(fā)生甘州5.0級地震，具體見圖4(b)。

金銀灘臺GM4[3]諧波振幅比曲線自2016年起部分速率累加曲線出現(xiàn)轉(zhuǎn)平或連續(xù)波動變化，門源6.4級地震前無顯著異常出現(xiàn)；阿左旗5.0級和九寨溝7.0級地震前20 s、80 s、110 s、140 s周期曲線變化平穩(wěn)，持續(xù)時間為1.5 a，其中震前20 s與140 s周期曲線小幅度波動變化顯著；甘州5.0級和夏河5.7級地震前異常形態(tài)不顯著，僅110 s和140 s的周期曲線波動變化平穩(wěn)，持續(xù)時間1 a左右。異常出現(xiàn)18～20個月后發(fā)生阿左旗5.0級和九寨溝7.0級地震，異常出現(xiàn)45～46個月后發(fā)生甘州5.0級和夏河5.7級地震，具體見圖5(a)。

金銀灘臺GM4[4]諧波振幅比曲線自2016年起部分速率累加曲線出現(xiàn)轉(zhuǎn)平或連續(xù)波動變化，門源6.4級地震前僅50 s、260 s周期曲線出現(xiàn)明顯下降趨勢；阿左旗5.0級和九寨溝7.0級地震前50 s、80 s、110 s、170 s、260 s周期曲線變化平穩(wěn)，持續(xù)時間為1.5 a，其中震前50 s與170 s周期曲線在變化初期還伴隨著下降趨勢；甘州5.0級和夏河5.7級地震前20 s、80 s、140 s和260 s周期曲線波動變化平穩(wěn)，持續(xù)時間1.5 a左右，且200 s周期曲線在震前10個月轉(zhuǎn)折變平。異常出現(xiàn)18～20個月后發(fā)生阿左旗5.0級和九寨溝7.0級地震，異常出現(xiàn)45～46個月后發(fā)生甘州5.0級和夏河5.7級地震，具體見圖5(b)。

貴德臺GM4[2]諧波振幅比曲線自2016年起部分速率累加曲線出現(xiàn)轉(zhuǎn)平或連續(xù)波動變化，門源6.4級地震前僅20 s、50 s、80 s、200 s周期曲線出現(xiàn)明顯下降趨勢；阿左旗5.0級和九寨溝7.0級地震前20 s、50 s、80 s、230 s周期曲線變化平穩(wěn)，持續(xù)時間為1.5 a，其中20 s、80 s與230 s周期曲線在變化初期還伴隨著波動變化；甘州5.0級和夏河5.7級地震前20 s、50 s、80 s、110 s和230 s周期曲線變化平穩(wěn)，持續(xù)時間1.5 a左右，且170 s周期曲線在震前10個月轉(zhuǎn)折變平。異常出現(xiàn)18～20個月后發(fā)生阿左旗5.0級和九寨溝7.0級地震，異常出現(xiàn)45～46個月后發(fā)生甘州5.0級和夏河5.7級地震，具體見圖6(a)。

圖6 貴德臺地磁諧波振幅比Fig.6 Geomagnetic harmonic amplitude ratio at Guide station

貴德臺GM4[4]諧波振幅比曲線在門源6.4級、阿左旗5.0級和九寨溝7.0級地震前未出現(xiàn)顯著異常；甘州5.0級和夏河5.7級地震前速率累加曲線50 s、80 s、110 s、170 s和290 s波動變化平穩(wěn)，持續(xù)時間1.5 a左右，且50 s周期曲線在異常出現(xiàn)初期伴隨著快速下降變化。異常出現(xiàn)17～18個月后發(fā)生甘州5.0級和夏河5.7級地震，具體見圖6(b)。

西寧臺FGM[4]諧波振幅比曲線在甘州5.0級和夏河5.7級地震前未出現(xiàn)顯著異常；阿左旗5.0級和九寨溝7.0級地震前速率累加曲線20 s、50 s、80 s、110 s、140 s、170 s自2016-04保持平穩(wěn)變化，其中20 s、50 s和110 s周期曲線還伴隨著小幅下降變化。異常出現(xiàn)14～16個月后發(fā)生阿左旗5.0級和九寨溝7.0級地震，具體見圖7(a)。

圖7 西寧臺和巴塘臺地磁諧波振幅比Fig.7 Geomagnetic harmonic amplitude ratio at Xi’ning and Batang station

玉樹巴塘臺FGM[1]諧波振幅比曲線在雜多6.2級地震前未出現(xiàn)顯著異常；稱多5.3級地震前20 s、80 s、170 s和260 s周期曲線轉(zhuǎn)平，持續(xù)時間1.5 a，期間50 s周期曲線于2018-09出現(xiàn)轉(zhuǎn)折變化；丁青5.1級地震前僅20 s周期曲線波動變化平穩(wěn)，其他周期曲線均保持正常上升趨勢。異常出現(xiàn)17個月后發(fā)生稱多5.3級地震，異常出現(xiàn)37個月后發(fā)生丁青5.1級地震，具體見圖7(b)。

統(tǒng)計地震前9組異常的基本參數(shù)信息可知，除2015年雜多6.2級地震前無異常外，其他地震震前均提取到異常信息，且異常與地磁暴無關(guān)。由于門源地震發(fā)生在觀測時段1 a后，無法記錄到異常開始時間，因此未統(tǒng)計異常持續(xù)時間，其他震例中異常持續(xù)時間最短為14個月，最長為46個月，異常周期數(shù)量較多，詳細(xì)參數(shù)信息見表1。

表1 2015～2020年青海磁通門地磁諧波振幅比異常參數(shù)與地震對應(yīng)統(tǒng)計

本文在總結(jié)異常的過程中發(fā)現(xiàn)如下規(guī)律：1)震中距較遠(yuǎn)的臺站異常曲線數(shù)量較多，異常形態(tài)顯著，如門源6.4級、九寨溝7.0級、稱多5.3級、甘州5.0級及夏河5.7級地震震前表現(xiàn)出的異常特征，具體見圖8。2)異常出現(xiàn)后2 a內(nèi)發(fā)震的頻次較高(圖9(a))，對應(yīng)地震中異常曲線數(shù)量超過1/2(圖9(b))，大部分分布在20～170 s周期內(nèi)(圖9(c))，且震級與震中距成正比，即震級越高震中距越遠(yuǎn)(圖9 (d))。但由于樣本數(shù)量有限，無法給出定量化的描述，且震級與異常曲線數(shù)量、震級與異常持續(xù)時間及震中距與異常持續(xù)時間之間并沒有發(fā)現(xiàn)明顯相關(guān)性，均比較離散。

圖8 異常臺站與震中分布Fig.8 Distribution of anomalous stations and epicenters

圖9 異常參數(shù)統(tǒng)計Fig.9 Statistical table of abnormal parameters

3 指標(biāo)體系

對2015年以來青海地區(qū)14套磁通門觀測資料進(jìn)行地磁諧波振幅比計算，并利用速率累加法提取異常信息，從8組震例信息中總結(jié)異常規(guī)律，確定青海地區(qū)的異常判別指標(biāo)和預(yù)測規(guī)則。

1)判別指標(biāo)。速率累加曲線持續(xù)上升為正常，持續(xù)下降和水平波動為異常，異常時間大于1 a，10條周期曲線中至少5條出現(xiàn)異常，且大部分分布在20～170 s周期范圍內(nèi)為異常成立。

2)預(yù)測規(guī)則。時間：速率累加異常出現(xiàn)后2 a內(nèi)發(fā)震概率較高；地點：震中距與震級成正比，即震中距越大震級越高；強度：M≥5。

4 結(jié) 語

針對當(dāng)前研究中地磁諧波振幅比的成果主要為分鐘值采樣數(shù)據(jù)及無定量化指標(biāo)2個關(guān)鍵問題，重新分析青海地區(qū)磁通門秒采樣觀測資料，并利用速率累加方法分析連續(xù)數(shù)據(jù)的逐日變化速率，再根據(jù)需要對分析時段進(jìn)行累加計算，使原有異常變化更加直觀和突出，異常變化幅度顯著，更易于定量判定和異常指標(biāo)提取。在此基礎(chǔ)上進(jìn)一步分析異常變化與強震的關(guān)系，總結(jié)青海地區(qū)磁通門秒值諧波振幅比方法的異常指標(biāo)，并將研究結(jié)果應(yīng)用于青海地區(qū)的震情跟蹤中。需要強調(diào)的是，本文結(jié)果僅基于5 a數(shù)據(jù)計算得出，因為統(tǒng)計樣本數(shù)量有限，不能完全反映客觀規(guī)律。隨著數(shù)據(jù)資料的進(jìn)一步積累，異常指標(biāo)將會得到完善。