改進的地磁低點位移法與原方法對比分析

馬 亮

1 甘肅省地震局，蘭州市東崗西路450號，730000

地磁低點位移法是一種使用地磁日變數據預報地震的方法。陳紹明[1]在研究1966-03邢臺6.8級、7.2級地震余震的前兆時發現地磁異常現象，并于1972年提出“日變低點位移法”。丁鑒海等[2-3]在此基礎上通過研究全國地磁臺站日變低點位移的空間分布特點，于1975年提出“低點位移平面圖”方法，并以此開展地震預報實踐工作。1986年陳紹明[4]在研究華北地區地磁日變曲線的反相位現象時提出“日變反向”的概念。丁鑒海[3]研究發現，低點時間均出現在臺站當地時間12時左右，于1988年提出低點時間的“經度效應”，并在1994年用數學歸納法推理低點位移法的判斷標準，從而解決該方法有經驗、無指標的缺點。2009年王亞麗等[5]通過建立地磁低點時間的期望值與地理經度的線性回歸方程指出，低點時間標準差與緯度間具有正相關關系。2017年許康生等[6]在研究岷縣-漳縣M6.6地震前地磁低點位移現象時發現，單個臺站的異常幅度與該臺站的震源距有關。2018年郭增建等[7]提出地磁低點位移以磁暴倍九法作為補充預測的觀點。賈昕曄等[8]使用“理論低點時間”與“實測低點時間”來區分地方時與世界時，并繪制低點時間垂直分量等值線圖。2019馬亮[9]提出隨經度變化的地磁臺站當地時的計算公式。戴勇等[10]使用克里格插值法求出岷縣-漳縣M6.6地震前的低點位移事件的突變界線，這是首次脫離手繪方式求取突變界線。本文將使用數學語言及MATLAB語言，對地磁低點位移法與其搭載軟件進行改進，通過建立突變界線與臺網地磁低點時間之間的幾何解析關系來實現自動成圖。

1 數據選取與問題提出

中國地震臺網中心調用全國124個地磁臺站的數據來監測地磁低點位移現象，圖1為臺站空間分布情況。東五區和東六區臺站密度小，經向間距大；臺灣省、上海市與港澳地區未建設地磁臺站。該臺網中觀測儀包含FHD-2質子矢量磁力儀、FHDZ-M15自動化臺站系統、GM4磁通門磁力儀、CTM-DI磁通門磁力儀、MINGEO磁通門磁力儀、Mag-01H磁通門地磁經緯儀、G856質子(旋進)磁力儀和GSM-19FD overhauser磁力儀，本文選用的數據均為該臺網的地磁低點時間的垂直分量。由于地磁變化常受到空間磁暴的影響，本文對中國科學院國家空間科學中心空間環境預報中心(http://www.sepc.ac.cn/Kp3HPred_chn.php)提供的磁暴信息進行篩除，以排除磁暴干擾。

采用等距離圓錐投影以確保經向距離不變形圖1 地磁臺站空間分布Fig.1 Spatial distribution of geomagnetic stations

王亞麗等[5]在研究低點時間的經度效應時發現，垂直分量的低點時間期望值與經度之間具有很強的負相關關系，這反映出低點時間具有顯著的地方時依賴性。由于每個測點的低點時間約為測點地方時12:00，經度相差15°的兩個地磁臺站的世界時剛好相差1 h(圖1)。如果兩個臺站的經向間距超過一個時區，突變界線的走向會出現較大誤差；如果兩個臺站的經向間距超過兩個時區，在未出現地磁異常的情況下也可能存在地磁低點位移現象。為了消除經度效應帶來的誤差，需采用標準公式將世界時統一轉化為地方時，再劃分突變界線。

中國地震臺網中心一般通過地震前兆信息處理與分析軟件EIS2000調用地磁數據并繪制低點時間分布圖。EIS2000中地磁低點時間采用格林威治時間，手繪突變界線并進行平滑處理，但該方法會產生兩種誤差。第一種誤差由經度效應引起，同一時區內最西側與最東側的格林威治時間會自動產生1 h偏差(圖1)，每個臺站格林威治時間的時差隨經度而不是隨時區序數而變化，尤其是突變界線沿南北向延展時，格林威治時間對突變界限的走向會產生較大影響。第二種誤差為手繪突變界線產生的隨機誤差，同時在平滑處理過程中又會產生二次誤差。第一種誤差較為簡單，可用臺站當地時間代替格林威治時間予以消除；第二種誤差處理難度較大，本文將采用新思路進行處理。

2 Voronoi剖分的應用

從集合的觀點來看，突變界線上任意一點到相鄰兩個臺站的距離相等，因此突變界線一定沿這兩個臺站的垂直平分線延伸。為獲取垂直平分線的坐標，本文引入Voronoi剖分技術，用Voronoi剖分線代替突變界線的可能走向線。Voronoi圖(以下簡稱V圖)與Delaunay三角剖分圖互為偶圖。以某一散點為核心的V圖內部部分被稱為該點的值守范圍，V圖中相鄰兩個散點之間的垂直平分線稱為Dirichlet鑲嵌線[11]，圖2為甘肅-青海地磁臺網的Delaunay三角剖分網和Dirichlet鑲嵌線。按照地磁臺站低點時間的大小對每個臺站的值守范圍進行涂色，如果Dirichlet鑲嵌線兩側的色差超過2 h，地磁低點位移的突變界線便會沿該條Dirichlet鑲嵌線延展。

圖2 甘肅-青海地磁臺網的Delaunay三角網和Voronoi圖Fig.2 Delaunay triangulation and Voronoi diagram of the Gansu-Qinghai geomagnetic station network

3 新時間格式與原時間格式的計算結果對比

從地磁日變曲線可以直觀地看出，地磁低點時間大多位于當地時間12:00前后，且服從正態分布[9]。為消除采用格林威治時間對突變界線的走向造成的誤差，本節地磁低點時間統一采用當地時間，與格林威治時間的轉化關系為：

Tl=T0+N

(1)

式中，Tl、T0分別為臺站當地時間與格林威治時間，N為臺站所在時區的序數。若將Tl表示為T0與臺站經度的函數，則式(1)可寫為[9]：

(2)

式中，lon為地磁臺站的經度。為避免低點時間隨經度變化而對突變界線的走向產生影響，本節的地磁低點時間統一采用式(2)進行求解，并非采用當地時區時間。由于T0已知，因此只需知道臺站經度便可計算地磁低點時間的當地時間，進而用集合思想計算出突變界線的坐標數組。

為評價使用式(2)的新算法相對于原算法的優勢，本文以2017年九寨溝7.0級地震發生前07-20出現的地磁低點位移現象為例，分析時間格式的改進對地磁低點位移突變界線走向的影響。從圖1中可以看出，全國地磁臺站的空間分布跨越5個時區，如果采用格林威治時間計算突變界線，東五區和東七區臺站的地磁低點時間在正常情況下將會存在2 h偏差，為了直觀地表現臺站間的偏差，通過最近鄰插值法將全國124個采用格林威治時間的地磁臺的低點時間繪制成空間等值線分布圖(圖3)。

由于世界時會隨經度而變化，在地磁數據無異常時，臺站經向跨度較大的地區低點時間可能會相差2 h以上，從而產生“假異常”現象。圖3中新疆5個地磁臺外邊緣的突變界線就屬于“假異常”突變界線，“假異常”是由世界時的固有性質而產生，與場地條件無關，采用本文方法中式(2)將世界時轉化為當地時后即可消除。

圖3 采用格林威治時間的突變界線Fig.3 Abrupt boundaries based on GMT

各臺站地質構造、地下電性結構存在差異，臺站資料差異對突變界線的影響已體現在原始數據中。因此只要使用相同的數據，無論采用原始手繪方法還是本文方法，在去除手繪的隨機誤差后，突變界線走向一致。手繪的原始突變界線為折線，但經過后期平滑處理后可轉化為平滑曲線。

新疆、西藏、青海、甘肅地區地磁臺站稀疏，各臺站經向跨度大，按照格林威治時間計算的突變界線的誤差會被放大，為體現與式(2)地方時的差異，本節按式(2)將圖3中中國西部臺站的原有低點時間轉化為地方時，并求出兩者的偏差值(表1)。

表1 西部地磁臺網的格林威治時間與地方時的低點時間對比Tab.1 Comparison of low point-time between GMT and local time of western geomagnetic network

從表1可以看出，獅泉河臺與喀什臺的低點時間與背景值之差超過60 min，各臺站格林威治地磁低點時間差異非常大，且經度跨度越大，地磁低點時間的差異越大。將格林威治時間轉化為式(2)地方時后，各臺站的低點時間約為12:00，與文獻[3]研究結果一致，表明新方法可有效消除“經度效應”。但文獻[3]未定義地方時的計算公式，該研究中“地方時”概念可能源于臺站的時區時間，因此本文建議將各臺站按式(2)計算的地方時的地磁低點時間長期觀測的期望定義為地磁低點時間的背景值，以消除原有方法的“經度效應”。

4 突變界線的定位

現虛擬一個由8個地磁臺站組成的地磁臺網，臺站緯、經度分別用矩陣lat和lon表示，低點時間用矩陣data表示，突變界線在MATLAB中的坐標計算代碼如下：

lat=[y1,y2,y3,y4,y5,y6,y7,y8];

lon=[x1,x2,x3,x4,x5,x6,x7,x8];

data=[a1,a2,a3,a4,a5,a6,a7,a8];

[X,Y]=meshgrid(linspace(y1,y8,n),

linspace(x1,x8,n));

Z=griddata(lat,lon,data,X,Y, ‘nearest’);

[C,h,CF]=contourf(Y,X,Z, 2.5);

其中，命令linspace(y1,y8,n)用來產生y1到y8之間等距離的n個點。[X,Y]=meshgrid(X,Y)(X=[x1,x2,x3],Y=[y1,y2,y3])函數用來產生網格矩陣，X為網格矩陣的橫坐標范圍，Y是網格矩陣的縱坐標范圍；返回網格矩陣的坐標[X,Y]。Z=griddata(lat,lon,data,X,Y, ‘n’)為散點插值函數，用于在散點[lat,lon]之間進行插值[11]，data為該點的色值，[X,Y]表示插值的網格(密度遠大于[lat,lon])，n為插值類型；返回插值結果矩陣Z。函數[C,h,CF]=contourf(Y,X,Z,m)可用來繪制二維等值線圖，Z為平面內網格點[Y,X]處的色值，m為梯度線的個數，添加色彩后會產生m+1種顏色。同時返回與命令contourc中相同的等高線矩陣C，C也可被命令clabel使用；返回包含patch圖形對象的句柄向量h；返回填充色值的矩陣CF。

突變界線可通過設置參數m，將等值線階數降為1階來實現。將contourf命令中的m值設為[11]：

m=n-1

(3)

(4)

式中，m為等值線條數，n為梯度數，d為總梯度差，L為步長。例如在某次地震事件中，低點時間等值線條數為m，所涂顏色共有n種，網內地磁低點時間的最早值與最晚值之差為d，兩個區域之間的低點時間相差步長為L，且L≡2。例如網內低點時間的最早值與最晚值分別為08:00和15:00，此時d=7，將d和L≡2代入式(3)和式(4)中求出m值為2.5。執行contourf命令可輸出以m為參數的地磁低點時間的等值線，突變界線為其中一條或一段。

式(3)中m并不是突變界線的條數，而是確定等值線的參數。在地磁臺網中，所有臺站的低點時間的相互差異程度決定m值的范圍，當0

在2015年皮山6.5級地震前06-21地磁低點位移事件中，m值為6，等值線有5條，突變界線有1條。在2016年門源6.4級地震前2015-12-23甘肅、青海地區出現的地磁低點位移現象中，m值為2，等值線有2條，突變界線有1條。在2016年甘肅金塔4.7級地震前02-28地磁低點位移事件中，m值為3，等值線有3條，突變界線有2條，本節案例參見文獻[11]。

5 數據分析

中國地震臺網中心于2017-07-20在華北、長江中下游地區出現的地磁低點位移現象中發現，低點時間被分為內外兩個區域，按原數據的格林威治時間統計的內區域平均低點時間為6.73時，外區域平均低點時間為3.84時，兩者相差2.89 h；將格林威治時間轉化為式(2)地方時后，內區域的平均低點時間為14.39時(表2)。外區域的平均低點時間為11.01時，與內區域相差3.38 h。桃源臺儀器出現故障，因此用邵陽臺地磁數據代替。在地磁低點位移現象發生后的第19天，即北京時間2017-08-08 21:19，四川省阿壩州九寨溝縣發生MS7.0地震，比文獻[3]中統計的極可能發震區間的最小值小4 d，震中距突變界線約185 km。2017-08-08的K指數為2，未受到外空間磁暴干擾，地磁低點數據真實可靠，因此認為，07-20出現的地磁低點位移現象與此次地震具有較好的相關性，圖4為中國地震臺網中心采用人工繪制的突變界線。

表2 內區域全部臺站的低點時間Tab.2 Low-point time of all stations in inner region

圖4 2017-07-20出現的地磁低點位移現象(采用當地時間)Fig.4 The geomagnetic low-point displacement phenomenon on July 20, 2017(use local time)

突變界線由乳山、濟南、紅山、定襄、太原、臨汾、萬州、石柱、涪陵、重慶、桃源、九峰、高郵、海安臺共同組成，外區域西端存在榆林、漢中、成都、美姑、貴陽、邵陽6個控制性臺站，對突變界線的走向起決定性作用。為了直觀地表現兩種形式的低點時間在空間分布上的差異以及§4中自動定位技術提取與人工畫線的差異，采用式(2)地方時結合§5方法求出全國124個臺站的低點時間等值線圖(圖5)。

圖5 采用當地時間繪制的突變界線Fig.5 Abrupt boundaries based on local time

由圖5可見，最近鄰插值結果中共有4條梯度線，第1條與人工繪制的突變界線形狀相似，第2條在襄陽市、南陽市周圍，第3條在山東省無棣大山地震臺周圍，第4條在新疆維吾爾自治區喀什市周圍。由于單個地磁臺站低點時間的偏移未形成大區域，因此認為，圖5中第2～4條梯度線不是地磁低點位移突變界線，只有第1條梯度線符合突變界線的定義。

圖5為使用式(2)地方時的結果，喀什臺與都蘭臺的相對偏差減小為185 min，烏什臺、紅淺臺、烏魯木齊臺與都蘭臺的低點時間色值相同。圖3為使用格林威治時間的結果，喀什臺和都蘭臺低點時間色值不同，相對偏差為274 min，烏什臺、紅淺臺、烏魯木齊臺和都蘭臺低點時間色值不同。上述分析表明，新方法可減小具有經向跨度的臺站之間低點時間的偏差，可過濾掉新疆地區的“假異常”突變界線。

對比圖4與圖5可知，前者是平滑曲線，后者為折線。圖4中東北地區的第2條突變界線在人工繪制時未被發現，而圖5中已繪制東北地區的突變界線。在華北、華中地區，兩個圖的突變界線走向大體一致，在河北北部略有差異。圖5中的突變界線在內區域存在兩個孤島，由于單個地磁臺站低點時間的偏移未形成大區域，因此本文未把孤島外輪廓作為地磁低點位移突變界線。

6 結 語

本文提出的方法在不改變地磁低點位移法基本原理的情況下可有效排除干擾，減少誤差。九寨溝地震前的地磁低點位移現象分析表明，本文提出的突變界線定位法與原有算法相比具有以下優點：

1)原有方法會產生“經度效應”， 新方法可很好地消除。

2)原有方法的突變界線可以反映突變界線的大致走向，但無法獲取突變界線內區域的孤島；而新方法可以精確地繪制突變界線內區域中的孤島。

3)原有地磁低點時間的數據格式無法統計與地磁低點背景時間的偏差，并且使用概率統計法識別異常會受來自外空場的干擾所影響，從而較難識別真正的異常；而新方法中的低點時間采用式(2)地方時，可以更精確地識別異常和來自外空場的干擾，過濾 “假異常”突變界線。

4)如果將線條看作點的集合，且假定原方法手繪的突變界線足夠精確，則原方法的突變界線為新方法的突變界線的非空真子集，后者包含前者，人工劃線時容易遺漏一些隱藏的突變界線，而新方法不會遺漏隱藏的突變界線。