董 博紀春玲張環(huán)曦李 鳳周安聘牛淑瑜趙雨晨劉 靜劉 檀朱音杰

1）中國石家莊050021河北省地震局

2）中國石家莊050021河北地質(zhì)大學

采用克里金插值法分析河北地區(qū)地磁場變化特征

董 博1），2）紀春玲1）張環(huán)曦1）李 鳳1）周安聘1）牛淑瑜1）趙雨晨1）劉 靜1）劉 檀1）朱音杰1）

1）中國石家莊050021河北省地震局

2）中國石家莊050021河北地質(zhì)大學

通過克里金插值方法，對河北地區(qū)2011—2014年地磁場中總強度F值及垂直分量Z值數(shù)據(jù)進行計算，分析該區(qū)近期地磁場變化特征。結(jié)果表明：在地震活動較為平靜狀態(tài)下，河北地區(qū)地磁場變化與地理緯度改變具有密切關(guān)系，與地理經(jīng)度變化關(guān)系不大；F值與Z值等值線分布與地磁臺站分布有關(guān)，且隨時間增加而增大。

地磁場；克里金插值法；等值線

0 引言

眾所周知，地球磁場隨時間緩慢變化，對這種地球基本磁場長期變化的研究，是地球物理學的重要課題之一。地磁場由不同磁場疊加，按其性質(zhì)，可分為來源于地球內(nèi)部并構(gòu)成地磁場主要成分的穩(wěn)定磁場和起源于地球外部的變化磁場。同時，地磁場有各類長期與短期變化，只有研究并充分認識其規(guī)律，才能正確提取異常場。因此，分析、研究某一地區(qū)地磁場變化特征，可以進一步認識該區(qū)地磁場活動規(guī)律，排除非震干擾信息，捕獲真正有用的地震前兆信息（張小濤等，2008）。

由于地磁場的分布和變化具有空間上的相關(guān)性和時間上的延拓性，可以根據(jù)時空分布特征，采用物理意義明確的數(shù)學方式來描述（薛革等，1995）。本文主要利用精度較高的克里金插值方法，對2011—2014年河北及周邊地區(qū)質(zhì)量較好且連續(xù)率較高的11個定點地磁臺站觀測資料進行分析。河北省及周邊地區(qū)地磁臺站分布較合理，地磁觀測數(shù)據(jù)可以反映該區(qū)域地磁場時空變化。通過對河北地區(qū)地磁場進行時空分布特征的分析計算，得到該地區(qū)近期地磁場時空分布演變規(guī)律及特征，不僅對了解地磁場變化形態(tài)具有意義，且有助于識別地磁異常信號，為提取更有價值的地震地磁異常信息做好準備工作（張小濤等，2008）。

1 資料收集

河北省現(xiàn)有10個地磁觀測臺站，因黃壁莊臺、文安臺、順平臺近年才開展數(shù)字化觀測，觀測資料大量缺記，故選用其他7個連續(xù)性較好、數(shù)據(jù)質(zhì)量可靠的地磁觀測資料作為研究對象。為了有效減小邊界效應(yīng)，并考慮地震臺站間距及分布，增加河北周邊地區(qū)4個地震臺（遼寧朝陽臺、河南?？h臺、山西大同臺、山東濟南臺）地磁觀測資料參與研究。各臺站參數(shù)見表1，地震臺分布見圖1。

大量震例表明，地磁垂直分量Z震前異常比其他分量明顯多（解用明，1996），因此本文選取地磁場總強度F及垂直分量Z研究河北地磁場近期變化，并考慮地震活動較為平靜時期及最新資料兩個因素。統(tǒng)計可知，每年8月地磁暴相對較多，地磁場擾動相對較大，而12月地磁場相對平穩(wěn)，且2011—2014年每年8月和12月全球地震，尤其大震發(fā)生次數(shù)相對較少，故選取該時段8月和12月地磁F值、Z值月均值進行分析，見表2、表3。

圖1　相關(guān)地磁臺站分布Fig.1 Distribution of related geomagnetic station

表1　臺站參數(shù)Table 1 The parameters of stations

表2　2011—2014年每年8月地磁數(shù)據(jù)Table 2 Geomagnetic data of 2011—2014 in August

表3　2011—2014年每年12月地磁數(shù)據(jù)Table 3 Geomagnetic data of 2011—2014 in December

因所選11個地磁臺站相距較遠，地磁觀測數(shù)據(jù)相差較大，為了直觀對比各臺數(shù)據(jù)變化，同軸繪制日均值曲線圖。將2011—2014年每年8月和12月F、Z日均值減去該月的月均值，得到圖2、圖3。

由圖2—圖3可以看出：①2011—2014年，地磁擾動較大的每年8月，日均值曲線變化幅度較大，最大幅度約300 nT；地磁場相對平靜的12月，日均值曲線變化幅度相對較小，變化范圍基本在50 nT以內(nèi)；②承德臺、昌黎臺地磁觀測數(shù)據(jù)日變化幅度相對較小，其他臺站相對較大；各地磁臺觀測數(shù)據(jù)變化幅度不同，長期變化趨勢一致，能較客觀記錄并反映河北地區(qū)地磁場近期變化特征。

圖2　2011—2014年每年8月各臺F、Z日均值Fig.2 The daily mean value of F and Z of 2011—2014 in August

圖3　2011—2014年每年12月各臺F、Z日均值Fig.3 The daily mean value of F and Z of 2011—2014 in December

2 計算方法

地磁數(shù)據(jù)的觀測和采樣往往在不規(guī)則且離散的測網(wǎng)（或測線）進行。為合理并正確解釋地球物理異常，在保持原有場特征前提下，有必要進行均勻且較密集規(guī)則網(wǎng)格的數(shù)據(jù)插值，形成直觀的二維等值線，為進一步異常分析提供依據(jù)（劉強等，2013）。

常用格網(wǎng)化插值方法包括三角網(wǎng)線性內(nèi)插方法、徑向基函數(shù)法、最小曲率法、克里金法、反距離加權(quán)法和多項式擬合法等（唐澤圣，1999；王建等，2004；王廣君等，2005）。反距離加權(quán)法缺點是，在格網(wǎng)區(qū)域內(nèi)會產(chǎn)生圍繞觀測點的“牛眼”，給磁法數(shù)據(jù)解釋帶來不便（李富等，2009）。三角網(wǎng)線性內(nèi)插方法僅用三角形的3個頂點線性內(nèi)插出三角形內(nèi)所有點的值，觀測點分布稀疏情況下很難保證插值精度，一般適合于小比例尺的地層模型和斷層（張琚等，2006）。徑向基函數(shù)法、最小曲率法和多項式擬合法均屬于平滑插值，一般不適合具有局部異常的磁法數(shù)據(jù)插值。當空間連續(xù)性變化的屬性不規(guī)則時，克里金插值法能夠在保持原有場特征前提下，更合理地對地球物理異常做出正確的解釋（陳歡歡等，2007；江玉樂等，2007）。

局部地磁場中的地球主磁場、地磁異常場及擾動磁場與克里金法的3種成分假設(shè)具有一一對應(yīng)關(guān)系：地球主磁場是局部地磁場中與恒定均值或趨勢有關(guān)的結(jié)構(gòu)性成分；地磁異常場為與空間變化有關(guān)的隨機變量，而且地磁異常場滿足二階平穩(wěn)假設(shè)；擾動磁場則可以視為與空間無關(guān)的隨機磁場噪聲。因此，克里金法能較為準確地反映局部地磁場的分布空間相關(guān)特性規(guī)律，可用于局部磁場空間插值（劉強等，2013）。

克里金法（Kriging）是由南非金礦工程師Krige D G 于1951年提出的一種空間插值方法，法國地質(zhì)學家Dr Mathhemn在1962年加以完善，并以克里金命名，最早應(yīng)用于礦產(chǎn)資源儲量估算。該方法建立在變異函數(shù)理論分析基礎(chǔ)上，對有限區(qū)域內(nèi)的區(qū)域化變量取值進行無偏最優(yōu)估計（best linear unbiased estimator，BLUE）?？死锝鸱ㄊ蔷€性的，因為其估計值是根據(jù)已有資料的加權(quán)線性結(jié)合獲得的。與其他估計方法相比，克里金插值法的平均殘差或誤差接近于零，使誤差方差最小是其顯著特點。該方法與傳統(tǒng)插值方法的不同之處在于，估計原觀測樣本數(shù)值時，不僅考慮待插值點與鄰近觀測數(shù)據(jù)點的空間位置，還考慮各鄰近點之間的空間位置關(guān)系，而且根據(jù)已有觀測數(shù)據(jù)空間分布的特點，使其估計比傳統(tǒng)方法更加貼近實際。

克里金法的基本原理為：設(shè)研究區(qū)為A，點承載的區(qū)域化變量為Z（x，y）［Z（x，y）∈A］在采樣點i（i=1，2，…，n）處的屬性值為Zi（i=1，2，…，n），待插值點Z0的屬性值插值結(jié)果是m個已知采樣點屬性值的加權(quán)和，即

式中：λi（i=1，2，…，m）為待確定權(quán)重系數(shù)。根據(jù)克里金原理，Z（x，y）在研究區(qū)內(nèi)滿足二階平穩(wěn)假設(shè)，即：①Z（x，y）的數(shù)學期望存在且等于常數(shù)，即E［Z（x，y）］=m；②Z（x，y）的協(xié)方差函數(shù)存在且相等（即只依賴于空間滯后距離h，而與位置無關(guān)），即

在無偏條件下使估計方差達到最小，經(jīng)過推導(dǎo)得到權(quán)重系數(shù)的方程組。

式中：μ為拉格朗日常數(shù)；γ（hi-hj）為樣本點間的半變異函數(shù)值；γ（h0-hi）為已知點與待插值點間的半變異函數(shù)值。求出權(quán)重系數(shù)λi（i=1，2，…，m），即可求得待插值點屬性值的估計值Z0。

克里金法的半變異函數(shù)定義為

式中：h為上文提到的滯后距離；Nh是滯后距離數(shù)量；Nk是用來計算半變異函數(shù)值的距離為h的樣本對數(shù)目；γ（h）表示半變異函數(shù)值只與h有關(guān)。由于距離為h的樣本點對通常有很多個，通過公式（4）得到的半變異函數(shù)值與h是一對多的關(guān)系，據(jù)此擬合半變異函數(shù)的實際形式比較困難。通常做法是，選取一定步長來對點分組，分別計算每組的平均半變異值，通過選取合適的擬合模型擬合參數(shù)，得出半變異函數(shù)的確切函數(shù)形式，帶入公式（3），求解待插值點的估計值。

3 資料分析

由于11個地磁臺站海拔高程差異較大，插值前需對地磁觀測數(shù)據(jù)進行高度改正。具體思路為：選定一個高度作為標準高度，選擇相應(yīng)經(jīng)緯度作為計算點得出該地正常場，計算不同高度修正值，若臺站海拔高于標準點，則實測值需要加上改正值；反之，減去改正值。利用surfer軟件，對河北地區(qū)2011—2014年地磁資料應(yīng)用克里金法進行插值計算，插值間隔為0.072×0.072，使x分成100個網(wǎng)格，而y分為82個網(wǎng)格，得到磁場總強度F和垂直分量Z的二維等值線和三維等值線。因2011—2014年地磁場總強度F和垂直分量Z的變化趨勢基本一致，在此只給出2014年F、Z等值線圖，見圖4、圖5。

圖4　2014年8月和12月地磁場總強度F二維、三維等值線（a） 2014年8月F值二維等值線；（b） 2014年12月F值二維等值線；（c） 2014年8月F值三維等值線；（d） 2014年12月F值三維等值線Fig.4 The 2D and 3D contour map of geomagnetic field total intensity F of 2014 in August and December

圖5　2014年8月和12月地磁場垂直分量Z值二維、三維等值線（a） 2014年8月Z值二維等值線；（b） 2014年12月Z值二維等值線；（c） 2014年8月Z值三維等值線；（d） 2014年12月Z值三維等值線Fig.5 The 2D and 3D contour map of geomagnetic field vertical component Z of 2014 in August and December

為避免邊界效應(yīng)，只分析圖形內(nèi)部等值線變化。由圖4可以看出：F值隨緯度增加而增加；在（114°—116°E，37°—39°N）F變化率較大，（118°—120°E，36°—41 °N）變化率較??；（117°E，36.5°N）的F值較周圍地區(qū)偏低，為該區(qū)極小值；（116.9°E，41°N）F值較周圍地區(qū)偏高，為該區(qū)極大值；2011—2014年F變化趨勢較平穩(wěn)，由北向南呈減小趨勢。

由圖5可以看出：地磁場垂直分量Z的變化與F基本一致。Z值等值線與緯度基本平行，且絕對值隨緯度的增加而增加。在（114°—116°E，37°—39°N）地磁場垂直分量Z變化率較大，（118°—120°E，36°—39°N）變化率較??；（117°E，36.5°N）的Z值較周圍地區(qū)偏低，為該區(qū)極小值；（116.9°E，41°N）的Z值較周圍地區(qū)偏高，為該區(qū)極大值；2011—2014年Z值變化趨勢較平穩(wěn)，由北向南呈減小趨勢。

4 結(jié)論

根據(jù)上述資料處理結(jié)果可知，河北地區(qū)近期地磁場總強度F和垂直分量Z的變化具有抑下特征：①F值、Z值變化與緯度密切相關(guān)，且隨緯度增加而增加，與經(jīng)度關(guān)系較??；二者極大值均出現(xiàn)在（116.9°E，41°N）附近，而極小值出現(xiàn)在（117°E，36.5°N）附近，極大值與極小值經(jīng)緯度位置經(jīng)度差異不大，緯度存在較大差異；②F值、Z值隨時間增加而增大，且呈逐年增大趨勢，年增幅分別約為35 nT和70 nT，未發(fā)現(xiàn)趨勢性變化；③F值、Z值分布與該區(qū)內(nèi)臺站分布狀況有關(guān)。地磁臺站比較集中在河北區(qū)的東北部和西南部，其等值線分布比較密集，而地磁臺站分布較少的西北部和東南部等值線分布相對較離散；④本文在分析河北地區(qū)地磁場變化特征時，只對該區(qū)地磁場的日均值和月均值進行了相關(guān)研究，在今后工作中將進一步考慮地磁場的日變化、季節(jié)變化和年變化等特征。

上述對地磁場特征的討論主要總結(jié)了地震活動較為平靜時期的特點，在此基礎(chǔ)上如果出現(xiàn)明顯的地磁場變化，則應(yīng)考慮地震發(fā)生的可能性。由于研究區(qū)域內(nèi)地磁臺站分布較少等原因，分析精度和深度都還不夠，有待于今后更進一步工作。

陳歡歡，李星，丁文秀.Surfer 8.0等值線繪制中的十二種插值方法［J］.工程地球物理學報，2007，4（1）：52-57.

江玉樂，李才明，張朝霞.B樣條分層離散插值在磁測異常處理中的應(yīng)用［J］.西南石油大學學報，2007，29（4）：12-15.

李富，王永華.10種插值方法在物探數(shù)據(jù)處理中的對比——以電法和磁法資料中的應(yīng)用為例［J］.四川地質(zhì)學報，2009，29（4）：474-476.

劉強，陶鈞，劉旭，等.基于克里金插值的磁法數(shù)據(jù)格網(wǎng)化研究［J］.河南科學，2013，31（7）：1 039-1 044.

唐澤圣.三維數(shù)據(jù)場可視化［M］.北京：清華大學出版社，1999.

王廣君，房建成.一種星圖識別的星體圖像高精度內(nèi)插算法［J］.北京航空航天大學學報，2005，31（5）：566-569.

王建，白世彪，陳曄.地理信息制圖［M］.北京：中國地圖出版社，2004.

薛革，陶九慶，張繼紅.山東地區(qū)地磁背景場的變化特征［J］.高原地震，1995，7（4）：50-54.

解用明.河北省近年地磁Z分量長期變特征［J］.山西地震，1996，2：47-50.

張琚，李星.利用IDL進行地學數(shù)據(jù)處理的多種插值法［J］.工程地球物理學報，2006，3（1）：49-53.

張小濤，王莉森，韓麗萍.河北省近年內(nèi)地磁場變化規(guī)律的初步分析［J］.地震地磁觀測與研究，2008，29（6）：17-21.

Abstract

By using the method of Kriging，the paper analyzed and calculated the total intensity F value and vertical component Z value from the year 2011 to 2014 in Hebei Province，and the authors finally got the variable features of geomagnetic field in Hebei area.The results showed that the variation of the geomagnetic background field closely related to the change of geographic latitude and no relationship with the change of geographic longitude；the distribution of isogram of the F and Z values was relation with the distribution of the geomagnetic stations，the values of F and Z were larger with the time.

The analysis of the variable features of geomagnetic field in Hebei area by using the method of Kriging

Dong Bo1），2），Ji Chunling1），Zhang Huanxi1），Li Feng1），Zhou Anpin1），Niu Shuyu1），Zhao Yuchen1），Liu Jing1），Liu Tan1）and Zhu Yinjie1）
1） Earthquake Administration of Hebei Province，Shijiazhuang 050021，China
2） Hebei GEO University，Shijiazhuang 050021，China

geomagnetic field，Kriging，isogram

10.3969/j.issn.1003-3246.2016.03.025

董博（1986—），男，助理工程師，主要從事地震監(jiān)測與異常跟蹤分析工作。E-mail： dongbo0002@126.com