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

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

1）中國(guó)石家莊050021河北省地震局

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

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

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

1）中國(guó)石家莊050021河北省地震局

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

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

地磁場(chǎng)；克里金插值法；等值線

0 引言

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

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

1 資料收集

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

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

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

表1　臺(tái)站參數(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個(gè)地磁臺(tái)站相距較遠(yuǎn)，地磁觀測(cè)數(shù)據(jù)相差較大，為了直觀對(duì)比各臺(tái)數(shù)據(jù)變化，同軸繪制日均值曲線圖。將2011—2014年每年8月和12月F、Z日均值減去該月的月均值，得到圖2、圖3。

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

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

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

2 計(jì)算方法

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

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

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

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

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

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

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

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

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

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

3 資料分析

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

圖4　2014年8月和12月地磁場(chǎng)總強(qiáng)度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月地磁場(chǎng)垂直分量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變化趨勢(shì)較平穩(wěn)，由北向南呈減小趨勢(shì)。

由圖5可以看出：地磁場(chǎng)垂直分量Z的變化與F基本一致。Z值等值線與緯度基本平行，且絕對(duì)值隨緯度的增加而增加。在（114°—116°E，37°—39°N）地磁場(chǎng)垂直分量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值變化趨勢(shì)較平穩(wěn)，由北向南呈減小趨勢(shì)。

4 結(jié)論

根據(jù)上述資料處理結(jié)果可知，河北地區(qū)近期地磁場(chǎng)總強(qiáng)度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值隨時(shí)間增加而增大，且呈逐年增大趨勢(shì)，年增幅分別約為35 nT和70 nT，未發(fā)現(xiàn)趨勢(shì)性變化；③F值、Z值分布與該區(qū)內(nèi)臺(tái)站分布狀況有關(guān)。地磁臺(tái)站比較集中在河北區(qū)的東北部和西南部，其等值線分布比較密集，而地磁臺(tái)站分布較少的西北部和東南部等值線分布相對(duì)較離散；④本文在分析河北地區(qū)地磁場(chǎng)變化特征時(shí)，只對(duì)該區(qū)地磁場(chǎng)的日均值和月均值進(jìn)行了相關(guān)研究，在今后工作中將進(jìn)一步考慮地磁場(chǎng)的日變化、季節(jié)變化和年變化等特征。

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

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

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

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

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

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

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

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

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

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

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

張小濤，王莉森，韓麗萍.河北省近年內(nèi)地磁場(chǎng)變化規(guī)律的初步分析［J］.地震地磁觀測(cè)與研究，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)測(cè)與異常跟蹤分析工作。E-mail： dongbo0002@126.com