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

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

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

2）中國石家莊050021河北地質大學

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

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

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

2）中國石家莊050021河北地質大學

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

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

0 引言

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

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

1 資料收集

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

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

圖1　相關地磁臺站分布Fig.1 Distribution of related geomagnetic station

表1　臺站參數Table 1 The parameters of stations

表2　2011—2014年每年8月地磁數據Table 2 Geomagnetic data of 2011—2014 in August

表3　2011—2014年每年12月地磁數據Table 3 Geomagnetic data of 2011—2014 in December

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

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

圖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 計算方法

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

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

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

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

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

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

在無偏條件下使估計方差達到最小，經過推導得到權重系數的方程組。

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

克里金法的半變異函數定義為

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

3 資料分析

由于11個地磁臺站海拔高程差異較大，插值前需對地磁觀測數據進行高度改正。具體思路為：選定一個高度作為標準高度，選擇相應經緯度作為計算點得出該地正常場，計算不同高度修正值，若臺站海拔高于標準點，則實測值需要加上改正值；反之，減去改正值。利用surfer軟件，對河北地區2011—2014年地磁資料應用克里金法進行插值計算，插值間隔為0.072×0.072，使x分成100個網格，而y分為82個網格，得到磁場總強度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

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

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

4 結論

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

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

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

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

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

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

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

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

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

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

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

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

張小濤，王莉森，韓麗萍.河北省近年內地磁場變化規律的初步分析［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—），男，助理工程師，主要從事地震監測與異常跟蹤分析工作。E-mail： dongbo0002@126.com