張志宏, 郭安寧, 李夢瑩, 黃明威, 楊牧萍

(1.遼寧省地震局， 沈陽 110034; 2.中國地震局蘭州地震研究所， 蘭州 730000)

經(jīng)過幾十年探索和實踐，中外已發(fā)現(xiàn)很多與地震相關(guān)的前兆現(xiàn)象[1-3]。地磁學(xué)是地震前兆觀測的重要手段，其主要研究內(nèi)容是地球自身及周圍空間磁場的起源、變化、分布規(guī)律及相關(guān)現(xiàn)象。地磁學(xué)通過實驗、觀測和理論研究，認識到了地球磁場的分布及變化規(guī)律。地磁背景場分布與變化的時、空都具有延拓性，且可以用物理意義明確的函數(shù)來表達。從而可以從空間有限個點的觀測值推算其他的點值，用過去的變化預(yù)測未來的趨勢。近年來，地震與地磁場的關(guān)系的研究取得了一些進展。目前中國地磁預(yù)測的地震方法有20多個，主要有：加卸載響應(yīng)比法[4]、地磁逐日比法[5]、日變低點位移法[6]、諧波振幅比法[7]、地磁Z分量日變化空間相關(guān)法[8]和地磁每日一值空間相關(guān)法[9]。這6種不同的分析方法在預(yù)測地震三要素方面都有各自的優(yōu)缺點：日變低點位移法預(yù)測強震發(fā)生的時間有較高的報準率，但是預(yù)測地震發(fā)生的地點區(qū)域有其局限性[6]。地磁Z分量日變化空間相關(guān)法對數(shù)據(jù)的質(zhì)量要求很高；地磁加卸載響應(yīng)比法和地磁逐日比法在預(yù)測地震發(fā)生的時間和震級具有獨特的優(yōu)勢，震中范圍預(yù)測存在不足[10-11]；地磁諧波振幅比法預(yù)測地震的震中大都在異常臺站的周圍，但是異常持續(xù)時間較長，不利于捕捉震前短期信號[12-13]。因此，對一次地震進行綜合預(yù)測成為一個新的探索方向。郭安寧等[14]研究表明多方法分析異?？梢暂^好地把握發(fā)震短臨時間節(jié)點。解滔等[15]利用電磁多方法綜合提取到了較好的2008年汶川MS8.0地震前電磁異常信號。

基于此，現(xiàn)分析2015—2019年東北地區(qū)固定地磁臺連續(xù)觀測波形數(shù)據(jù)，將地磁加卸載響應(yīng)比法和地磁諧波振幅比法方法結(jié)合各自在預(yù)測地震優(yōu)勢，綜合研究2018年5月28日松原MS5.7地震前后數(shù)據(jù)的變化特征，作為地磁多方法預(yù)測地震的震例。

1 震中附近地質(zhì)構(gòu)造背景

如圖1所示為東北地區(qū)橫跨中朝地臺和天山—興安地槽兩大地質(zhì)構(gòu)造單元，以赤峰—開原斷裂為界，北部的天山—興安地槽自西向東可分為大興安嶺華力西褶皺帶、內(nèi)蒙華力西褶皺帶、松遼拗陷和吉黑華力西褶皺帶；南部的中朝地臺自西向東可分為內(nèi)蒙地軸燕山臺褶皺帶、華北斷拗和膠遼臺隆[16]。2018年5月28日吉林松原發(fā)生MS5.7地震，震源深度為10 km，震中位于松遼構(gòu)造盆地[17]。

紅色實心圓為東北地區(qū)公元1100年以來的6級以上地震；藍色三角為分析地磁異常所使用的臺站；黃色五角星為2018年吉林松原MS5.7地震的震中位置圖1 東北地區(qū)地質(zhì)構(gòu)造與地震活動和地磁臺站分布情況Fig.1 Geological structure, seismic activity and distribution of geomagnetic stations in Northeast China

2 計算原理與地磁資料處理

2.1 加卸載響應(yīng)比法計算原理

地震的孕育和發(fā)生是一個復(fù)雜的非線性過程，因此，20世紀90年代初，尹祥礎(chǔ)等[18]提出了定量反映非線性系統(tǒng)不穩(wěn)定程度的參數(shù)——加卸載響應(yīng)比Y，其可表述為

(1)

式(1)中：“+”和“-”分別為應(yīng)力加載和卸載，某區(qū)域在特定時段內(nèi)所有正地震(個數(shù)為N+)能量E的m次冪的和與負地震(個數(shù)為N-)能量的m次冪的和的比值定義為加卸載響應(yīng)比Y。當孕震介質(zhì)處于低應(yīng)力的穩(wěn)態(tài)彈性階段時，Y接近于1。反之，Y隨著失穩(wěn)的過程逐漸增大。自此，加卸載響應(yīng)比理論在各種地球物理參量得到了廣泛應(yīng)用[19]。馮志生等[20]研究發(fā)現(xiàn)地磁垂直分量Z日變相比其他分量震前異常更為顯著。

地磁加卸載響應(yīng)比計算公式為

(2)

式(2)中：RZ(max)為地磁垂直分量日變化幅度極大值；RZ(min)為出現(xiàn)之后的第一個極小值。

2.2 克里金插值法原理

1951年南非工程師Krige D G首先提出了克里金法(Kriging)，相比三角網(wǎng)線內(nèi)插法、最小曲率法、反距離加權(quán)法和多項式擬合等常規(guī)插值法，克里金插值法在保持原場特征的同時更加合理地解釋了地球物理異常[21]。

該方法是運用變異函數(shù)理論對有限區(qū)域內(nèi)的區(qū)域變化的變量無偏最優(yōu)化估計。其優(yōu)點在于考慮待插值點與鄰近觀測數(shù)據(jù)點的位置的同時，對鄰近點之間的空間位置也有所參考，使其誤差接近于零[22]。

設(shè)定一個研究區(qū)D，區(qū)域化變量為Z(x,y)[Z(x,y)∈D]在采樣i(i=1,2,…,n)點處的屬性值為Zi(i=1,2,…,n)，待插值Z0點的插值結(jié)果則可以表述為

(3)

式(3)中：λi(i=1,2,…,m)為待確定權(quán)重系數(shù)；m為已知采樣點。

根據(jù)Kriging原理，Z(x,y)滿足二階平穩(wěn)假設(shè)，即

Cou{Zi,Zi+h}=E[Zi，Zi+h]-m2=C(h)

(4)

根據(jù)無偏性要求可推導(dǎo)得到權(quán)重系數(shù)方程為

(5)

式(5)中：h為滯后距離；Nh為滯后距離數(shù)量；Nk為求取半變異函數(shù)值的距離為h的樣本對數(shù)目。求解權(quán)重系數(shù)可得λ到待插點估計值Z0。

克里金插值法的半變異函數(shù)可表述為

(6)

為了彌補研究區(qū)固定地磁臺網(wǎng)離散的缺點，對加卸載響應(yīng)比值進行密集規(guī)則網(wǎng)格化插值處理，得到直觀的二維等值線空間分布。

2.3 地磁諧波振幅比計算原理

根據(jù)電磁感應(yīng)理論，對于隨時間周期變化的不均勻場源，視地球介質(zhì)為均勻各向同性的平面導(dǎo)體的條件下，地磁諧波振幅比的定義為

(7)

式(7)中：YZHx(NS)和YZHy(EW)分別為南北向和東西向的地磁諧波振幅比；Z(ω)、Hx(ω)、Hy(ω)分別為地磁場的垂直分量、南北向水平分量及東西向水平分量的譜值。

地磁諧波振幅比也可以表達為磁測深視電阻率,即

(8)

式(8)中：ρax和ρay分別為北南向和東西向的深部電阻率；ω為圓頻率;μ0為磁導(dǎo)率;λ為常數(shù)。

由式(8)可知，地磁諧波振幅比能夠反映介質(zhì)深部電阻率的變化。即地磁諧波振幅比算法具有明確的物理意義。

2.4 地磁數(shù)據(jù)處理方法

選取震中距700 km內(nèi)11個地磁固定臺的7套FHD磁力儀、2套GM4型磁通門磁力儀、3套M15磁力儀。地磁加卸載響應(yīng)比異常屬于中短期異常，因此，選取了2018年上半年的地磁Z分量日變化分鐘值采樣數(shù)據(jù)進行48階傅氏擬合濾波后計算日變幅度，得到各個臺站日變幅度逐日變化序列；其次根據(jù)日變幅度逐日變化序列計算各臺站的地磁加卸載響應(yīng)比值，提取所有臺站加卸載響應(yīng)比值，采用Surfer 進行克里金網(wǎng)格化插值處理，得到空間插值等值線數(shù)據(jù)。地磁諧波振幅比計算時間為2016年6月—2019年6月，諧波振幅比具體計算分為5個步驟：①計算每天地磁三分量采樣數(shù)據(jù)譜幅度；②獲得2個方向 5～65 min各周期的地磁諧波振幅比；③按 10 min的頻帶寬度計算各頻帶地磁諧波振幅比的頻帶均值，獲得 10、20、30、40、50、60 min的地磁諧波振幅比；④計算 10、20、30、40、50、60 min的逐日變化序列；⑤采用窗長一年的滑動平均法消除各頻帶地磁諧波比逐日變化序列的年變化。

3 地磁異常分析

3.1 加卸載響應(yīng)比異常分析

地磁場是地球內(nèi)部磁性巖石及分布在地球內(nèi)部和外部電流體系產(chǎn)生的各種磁場成分的疊加。由于磁場起源的不同，各種磁場成分的空間分布和時間變化規(guī)律也大不相同，按照場源的劃分，此次地震場源為外源磁場。外源場起源于地表以上的空間電流系，這些電流體系隨著時間變化較快，變化磁場的畸變具有短時性。因次，地磁加卸載響應(yīng)比值異常表現(xiàn)為短時且跨越式的特征。

如表1所示計算了2008 年以來東北地區(qū)11個固定臺站觀測連續(xù)波形，剔除磁暴引起的日變幅畸變異常數(shù)據(jù)后得到東北地區(qū)加卸載響應(yīng)比均值P，結(jié)果表明東北地區(qū)地磁場P值背景場區(qū)間是[1.53，1.75]。

表1 東北地區(qū)地磁加卸載響應(yīng)比均值P

圖2 東北地區(qū)2018年1—6月地磁加卸載響應(yīng)比時序曲線Fig.2 Time series curve of geomagnetic loading and unloading response ratio in Northeast China from January to June 2018

研究表明，地震對應(yīng)的地磁加卸載異常表現(xiàn)為震前中短期的高值異常，中國大陸地區(qū)異常閾值為3.0對應(yīng)地震的概率較高[23]。如圖2所示為2018 年5月28日MS5.7地震前6個月內(nèi)的加卸載響應(yīng)比時序曲線。圖3(a)～圖3(g)是2018 年5月28日MS5.7地震前6個月內(nèi)加卸載響應(yīng)比高于閾值3.0的7次異常的時空分布圖。圖3(a)中，2018年1月5日距離此次地震震中約137 km的吉林三崗臺加卸載響應(yīng)比超閾值異常，異常值為3.1。經(jīng)調(diào)查異常日為磁靜日，三崗臺周邊環(huán)境及觀測系統(tǒng)正常。圖3(b)中，1月18日內(nèi)蒙錫林浩特及周邊地區(qū)及遼寧鐵嶺地區(qū)加卸載響應(yīng)比超閾值異常，其中遼寧鐵嶺臺為3.1、錫林浩特為3.2。圖3(c)中，2月11日德都地區(qū)加卸載響應(yīng)比存在高值區(qū)，德都臺計算結(jié)果為3.1。圖3(d)中，2月14日遼寧營口、朝陽地區(qū)存在高值區(qū)，其中響應(yīng)比值均為3.0。圖3(e)中，2月20日出現(xiàn)了較大面積的高值異常區(qū)，東北東部大部分地區(qū)處于高值水平，大連、營口及鐵嶺地區(qū)的加卸載響應(yīng)比值處于3.2～4.0的高值水平。圖3(f)中，2月26日黑龍江德都比值為3.0。圖3(g)中，4月24日東北西部地區(qū)出現(xiàn)了大范圍的加卸載超閾值異常，其中距離震源最近的三崗臺為3.8、最遠的錫林浩特比值達到了4.9的高值，通河的比值為3.3，此次異常是2018年5月28日松原MS5.7地震前最后一次異常。

3.2 地磁諧波振幅比異常分析

選取研究區(qū)距離震中最近的三崗和鐵嶺臺的Z、H、D三分量分鐘值數(shù)據(jù)作為計算諧波振幅比。

圖3 2018年5月28日松原MS5.7地震前地磁加卸載響應(yīng)比異常時空演化圖Fig.3 Temporal and spatial evolution of the abnormal geomagnetic load-unload response ratio before the Songyuan MS5.7 earthquake on May 28, 2018

如圖4所示，鐵嶺臺諧波振幅比分布在(0.150，0.425)區(qū)間，YZHx和YZHy兩個方向 2016—2017年 6 個周期變化具有明顯的同步性。2018年上半年 10、20 、30 min 周期開始呈下降趨勢(圖中紅色箭頭)，而 40、50、60 min 正好反向(圖中黑色箭頭)。2018 年 5 月 28 日松原MS5.7 地震發(fā)生在兩個方向10、20 、30 min 周期趨勢下降階段，震后兩個方向6個周期恢復(fù)上升準同步形態(tài)。

如圖5所示，三崗臺計算結(jié)果顯示諧波振幅比分布在(0.046，0.210)區(qū)間內(nèi)，諧波振幅比YZHx和YZHy兩個方向 2016—2017年主要表現(xiàn)為準同步趨勢上升特征。YZHx20、50 min 周期2018年初開始表現(xiàn)為趨勢下降異常(圖5中紅色箭頭)，YZHy50 min同樣表現(xiàn)為趨勢下降異常(圖中黑色箭頭)。2018年5月28日松原MS5.7地震后兩個方向50 min 周期異常恢復(fù)，YZHx20 min 周期異?；謴?fù)后再次出現(xiàn)下降異常。

圖4 2016—2019年鐵嶺臺諧波振幅比Fig.4 Tieling station harmonic amplitude ratio from 2016 to 2019

4 討論

震磁前兆信息的研究主要集中在基本磁場和變化磁場。按照場源的劃分，研究區(qū)固定地磁臺觀測數(shù)據(jù)屬于變化磁場。地震的孕育不僅是震源區(qū)孕震的過程，同時也和區(qū)域應(yīng)力場及多種環(huán)境因素的變化密切相關(guān)。變化磁場是聯(lián)系地震孕育的內(nèi)、外環(huán)境的媒介，其“感應(yīng)-觸發(fā)”作用，對地震預(yù)報有著重要的意義。地磁諧波振幅比和地磁加卸載響應(yīng)比可以理解為“感應(yīng)”和“觸發(fā)”的關(guān)系。

諧波振幅比計算結(jié)果表明，三崗臺和鐵嶺臺諧波振幅比兩個方向的不同周期存在不同步異常，異常特征為下降—轉(zhuǎn)折—恢復(fù)，2018年5月28日MS5.7地震發(fā)生異?；謴?fù)過程階段，震后基本恢復(fù)年變形態(tài)。2018年5月28日松原MS5.7 震中位于NE向扶余—肇東斷裂帶和NW向第二松花江斷裂帶交匯部位，在研究東北地區(qū) NW 向地質(zhì)構(gòu)造時，郯廬斷裂帶北段和兩個分支斷裂帶是研究者關(guān)注的重點。如圖6所示，三崗臺及鐵嶺臺從廣義的地質(zhì)構(gòu)造來看都屬于郯廬斷裂帶北段，此斷裂帶是中國大陸歷史強震活躍區(qū)和現(xiàn)今地震重點危險區(qū)之一，第二松花江斷裂 EW 向的東端與依蘭—伊通斷裂交叉。李君等[24]利用雙差定位方法對此次地震進行了重新定位，并使用CAP方法求解了震源機制解，結(jié)果表明此次地震的發(fā)震斷層為松花江第二斷裂。松原地區(qū)應(yīng)力場受到大尺度板塊構(gòu)造應(yīng)力控制，同時受到區(qū)域構(gòu)造活動的影響[25]。臺站處于高導(dǎo)異常體邊界附近是地磁諧波振幅比異常的機理解釋[26]，由此可推斷，三崗臺和鐵嶺臺可能處于同一高導(dǎo)異常體邊界帶上。

圖5 2016—2019年三崗臺諧波振幅比Fig.5 Sangang station harmonic amplitude ratio from 2016 to 2019

f1為松花江第二斷裂；f2為扶余—肇東斷裂；f3為依蘭—伊通斷裂；f4為密山—敦化斷裂圖6 2018年5月28日松原MS5.7地震地磁諧波振幅比異常臺站分布圖Fig.6 Distribution map of the abnormal stations of the Songyuan MS5.7 earthquake and geomagnetic harmonic amplitude ratio on May 28, 2018

2018年5月28日松原MS5.7地震前地磁加卸載響應(yīng)比出現(xiàn)了7次超閾值異常。其中，1月5日除三崗臺外研究區(qū)空間磁場日變處于正常背景。與其類似的還有2月11日和2月26日。1月18日錫林浩特臺周圍出現(xiàn)約2.5×105km2異常區(qū)，同時鐵嶺臺出現(xiàn)了小范圍異常。2月14日營口臺、朝陽臺周邊存在小范圍異常區(qū)。這5次跨越式的短時異常從此次地震的發(fā)震區(qū)域構(gòu)造角度可靠性不足。而2月20日和4月26日的異常臺站超過了3個，此次地震震中屬于兩次異常的異?？刂茀^(qū)域。因此，2月20日和4月24日可能是“觸發(fā)”階段。

綜上所述，2018年5月28日松原MS5.7地震區(qū)域應(yīng)力場的“感應(yīng)”可認為開始于2015年的下半年，“觸發(fā)”時間節(jié)點分別為2018年的2月和4月末。

5 結(jié)論

通過2018年5月28日松原MS5.7地震前的地磁數(shù)據(jù)分析可以得出以下結(jié)論。

(1)地震前地磁加卸載響應(yīng)比異常是中短期異常，震前會有多次高值跨越短時異常，地震發(fā)生的位置在研究區(qū)地磁加卸載響應(yīng)比值均值2倍以上區(qū)間。

(2)地震諧波振幅比異常在震前1年內(nèi)出現(xiàn)，異常特征為不同方向同周期不同步異常，其形態(tài)表現(xiàn)為下降—轉(zhuǎn)折—恢復(fù)，地震發(fā)生在恢復(fù)過程中。

(3)磁加卸載響應(yīng)比和地磁諧波振幅比兩種分析方法各有優(yōu)勢，前者是臨震短期信號，而后者可能反映的是區(qū)域電性結(jié)構(gòu)趨勢異常。結(jié)合兩種異常分析方法把握地震的震中和發(fā)震時間有重要的意義。