時(shí)空Kriging算法在區(qū)域地磁場(chǎng)插值中的應(yīng)用及改進(jìn)

陳鼎新， 劉代志， 曾小牛， 孟亮， 楊曉君

火箭軍工程大學(xué)907教研室， 西安　710025

時(shí)空Kriging算法在區(qū)域地磁場(chǎng)插值中的應(yīng)用及改進(jìn)

陳鼎新， 劉代志*， 曾小牛， 孟亮， 楊曉君

火箭軍工程大學(xué)907教研室， 西安710025

摘要本文以區(qū)域地磁場(chǎng)為應(yīng)用背景，研究了時(shí)間域信息對(duì)地磁場(chǎng)插值效果的影響.用時(shí)空Kriging算法進(jìn)行插值時(shí)，在時(shí)空Kriging變差函數(shù)的擬合過(guò)程中，將角度信息引入向量距離，定義了新的向量距離形式，提出了改進(jìn)的時(shí)空Kriging算法，并對(duì)權(quán)重的選取進(jìn)行討論.對(duì)經(jīng)度87.2°E—126.6°E，緯度19.0°N—49.6°N范圍內(nèi)的32個(gè)臺(tái)站時(shí)均值地磁場(chǎng)數(shù)據(jù)進(jìn)行插值實(shí)驗(yàn).在添加了時(shí)間信息之后，插值結(jié)果在兼顧經(jīng)度方向和緯度方向變化的同時(shí)，邊值問(wèn)題也得以緩解，明顯優(yōu)于傳統(tǒng)方法.交叉驗(yàn)證結(jié)果表明，隨著時(shí)間域信息和向量角度的加入，Kriging算法的性能得到改進(jìn)，地磁場(chǎng)的插值精度依次提高.對(duì)比實(shí)驗(yàn)中，改進(jìn)Kriging算法的精度最高、性能最穩(wěn)定，MAE集中在1 nT左右，MSE集中在[0,5]的范圍內(nèi).

關(guān)鍵詞地磁場(chǎng)； Kriging； 時(shí)空Kriging； 變差函數(shù)； 向量角度

1引言

地磁場(chǎng)由95%的主磁場(chǎng)和5%的變化場(chǎng)組成，其中，主磁場(chǎng)由液態(tài)地核的磁流體發(fā)電機(jī)效應(yīng)和地殼的磁性物質(zhì)產(chǎn)生，變化場(chǎng)受太陽(yáng)活動(dòng)、地球自轉(zhuǎn)、公轉(zhuǎn)等因素影響，變化過(guò)程具有時(shí)間的連續(xù)性(徐文耀，2003).地磁場(chǎng)研究在空間天氣預(yù)報(bào)、地磁導(dǎo)航等領(lǐng)域有著重要應(yīng)用(徐文耀，2009；袁楊輝，2012).目前，描述全球尺度地磁場(chǎng)的模型有GMM、WMM和CM4等，但這些模型的精度都不高(例如，CM4的時(shí)間精度最高為6小時(shí))，在局地分析、地磁實(shí)時(shí)導(dǎo)航等領(lǐng)域不能滿足要求(劉代志等，2005；易世華等，2013；牛超等，2014).因此，如何在局部地區(qū)獲得有效的地磁場(chǎng)分析數(shù)據(jù)，成為關(guān)鍵技術(shù)之一(趙明等，1997；郭才發(fā)等，2011).

空間相關(guān)的插值方法有Kriging(Goldenberg, 2006；張曉明和趙剡，2009； Huang et al., 2012)、多項(xiàng)式插值(孫涵等，2010)、最近鄰等，已較多地應(yīng)用于地磁場(chǎng)的建模分析中.然而，地磁臺(tái)站信號(hào)記錄的地磁信號(hào)，為時(shí)間序列，地磁場(chǎng)數(shù)據(jù)在時(shí)間上具有連續(xù)性(李夕海等，2009；曾小牛，2011；曾小牛，2014).以上方法只關(guān)注臺(tái)站空間位置對(duì)插值的影響，而忽略了時(shí)間域的信息.

時(shí)空Kriging方法對(duì)Kriging進(jìn)行了時(shí)空域擴(kuò)展，從而充分利用數(shù)據(jù)中的時(shí)間域信息來(lái)進(jìn)行插值.目前，時(shí)空Kriging方法已經(jīng)在空氣質(zhì)量檢測(cè)(De Iaco et al., 2003; Nunes and Soares, 2005)、雨水建模(李莎等，2011; li et al., 2011; Tang, 2002)、風(fēng)力數(shù)據(jù)插值(Gneiting, 2002)以及地下水建模(Snepvangers et al., 2003; Jost et al., 2005; Cichota et al., 2006; Lark et al., 2006)等領(lǐng)域得到了廣泛應(yīng)用.Gething等(2007)甚至用時(shí)空Kriging來(lái)分析瘧疾門(mén)診病人的數(shù)據(jù)，將時(shí)空Kriging應(yīng)用到社科數(shù)據(jù)的分析當(dāng)中，并取得了良好效果.然而，在地磁學(xué)領(lǐng)域，時(shí)空Kriging方法至今沒(méi)有應(yīng)用，使得插值過(guò)程中時(shí)間域信息被忽視.

筆者曾研究過(guò)用時(shí)空Kriging方法來(lái)分析地磁場(chǎng)數(shù)據(jù)所涉及的精度問(wèn)題、插值過(guò)程中數(shù)據(jù)庫(kù)容量的選擇等.本文重點(diǎn)討論從空間維度向時(shí)間-空間維度擴(kuò)展過(guò)程中，Kriging方法自身發(fā)生的變化.在時(shí)空變差函數(shù)的擬合過(guò)程中，將角度信息添加到向量距離的定義當(dāng)中，提出了改進(jìn)的時(shí)空Kriging算法，并對(duì)向量距離中角度信息的權(quán)重選擇進(jìn)行討論.

2時(shí)空Kriging算法及其改進(jìn)

時(shí)空Kriging的核心思想，是由Kriging方法在空間域插值的基礎(chǔ)上，向時(shí)間-空間域進(jìn)行擴(kuò)展.

2.1普通Kriging

Kriging 方法是以南非礦業(yè)工程師D.G.Krige(克里金)名字命名的一項(xiàng)實(shí)用空間估計(jì)技術(shù)，是地質(zhì)統(tǒng)計(jì)學(xué)的重要組成部分(Cressie, 1988; Bajat et al., 2013).假設(shè)在待估計(jì)點(diǎn)x的臨域內(nèi)共有n個(gè)實(shí)測(cè)點(diǎn)，即x1,x2,…xn，其樣本值為Z(xi).普通Kriging的插值公式為(Cressie, 1990; Du and Yang, 2012)

(1)

其中的權(quán)重系數(shù)λi通過(guò)Kriging方程組解算，公式為

(2)

定義變差函數(shù)為

(3)

由于

(4)

式(2)可以轉(zhuǎn)化為

(5)

寫(xiě)成矩陣形式為

(6)

其中，

2.2時(shí)空Kriging

時(shí)空Kriging是在Kriging算法的基礎(chǔ)上，引入時(shí)間域的連續(xù)性與相關(guān)性(Rouhani and Myers, 1990).用A=(si,tj)示時(shí)空域中某一點(diǎn)的坐標(biāo)，則該點(diǎn)處特征量的值可以表示為鄰域內(nèi)所有點(diǎn)的加權(quán)和，即

(7)

-Z(s,h)]2,

(8)

則可以通過(guò)Kriging方程組

(9)

來(lái)求解λ.現(xiàn)實(shí)中，由于時(shí)間和空間中變量的量綱不同，時(shí)空聯(lián)合域的距離很難統(tǒng)一表示，直接求變差函數(shù)的時(shí)空聯(lián)合分布是困難的(王建民等，2014).然而，通過(guò)一些典型的模型(Journel, 1986;Porcuetal, 2008)，例如Product-Sum(李莎等，2011；宋瑩，2013)、Product-Integration(DeIacoetal，2002；李莎等，2012)及Cressie-Huang模型(CressieandHuang，1999)等，可以用時(shí)間域條件分布函數(shù)和空間域條件分布函數(shù)

(10)

來(lái)表示時(shí)空域的變差函數(shù)(CressieandHuang，1999)，公式為

(11)

將式(11)帶入到式(9)中，即可解得時(shí)空域Kriging方程組.

2.3改進(jìn)時(shí)空Kriging算法

隨著定義范圍的不同，式(3)中Z(x)與式(10)中Z(x)的含義也發(fā)生了變化.如圖1所示，式(3)中Z(x)表示空間中點(diǎn)(x)處地磁場(chǎng)的值，而式(10)中Z(x)則表示，點(diǎn)(x)處各時(shí)間采樣點(diǎn)所組成的序列.

圖1　Kriging時(shí)空擴(kuò)展演示圖Fig.1　Illustration of expiation of Spatial Temporal Kriging

換言之，Z(x)表示對(duì)應(yīng)的空間切片，即Z(x)=(Z(x,t1),Z(x,t2),…,Z(x,tn))，Z(x)亦表示對(duì)應(yīng)的時(shí)間切片，即Z(t)=(Z(x1,t),Z(x2,t),…,Z(xm,t)).Z(x)和Z(t)實(shí)質(zhì)上是由切片上各點(diǎn)的值組成的向量，因此，式(10)也可表示為

(12)

其中，d(Z(xi)-Z(xi+hs))和d(Z(tj)-Z(tj+ht))表示向量間的距離.當(dāng)取L2范數(shù)來(lái)表示d(Z(xi)-Z(xi+hs))和d(Z(tj)-Z(tj+ht))時(shí)，則得到通常的時(shí)空Kriging方法；在此定義下，再取為某一常數(shù)，則得到空間域的普通Kriging方法.然而，在時(shí)間-空間域中，L2范數(shù)不一定是最好的方案.

L2范數(shù)的定義為

(13)

D(Z1,Z2)=p‖Z1-Z2‖2+(1-p)

(14)

式(14)等號(hào)右邊第一項(xiàng)為L(zhǎng)2范數(shù)定義的向量距離，第二項(xiàng)為向量之間的角度，通過(guò)設(shè)置參數(shù)p來(lái)控制二者的權(quán)重.p=1時(shí)，向量距離中的角度項(xiàng)被去掉，只保留了幅度部分，此時(shí)的向量距離與向量中值濾波器(VectorMedianFilter，VMF)中的距離定義相同(Astolaetal., 1990)；p=0時(shí)，空間距離被去掉，只考慮角度差異，此時(shí)的向量距離與向量角度濾波器(VectorDirectionalFilter，BVDF)中的距離定義相同(Trahaniasetal., 1996).

將式(14)代入式(12)，得到新向量距離定義下的條件變差函數(shù)為

(15)

進(jìn)而構(gòu)造出時(shí)空域聯(lián)合變差函數(shù)γst(hs,ht).

3實(shí)驗(yàn)

3.1數(shù)據(jù)

實(shí)驗(yàn)采用的數(shù)據(jù)，是中國(guó)地磁臺(tái)網(wǎng)及其他幾個(gè)單獨(dú)臺(tái)站(共32個(gè)觀測(cè)臺(tái)站)從Jan. 1st2004到Sep.30th2004的地磁場(chǎng)日均值序列.臺(tái)站位置如圖2所示，五角星分別表示32個(gè)觀測(cè)臺(tái)站的地理位置，覆蓋了經(jīng)度87.2°E—126.6°E，緯度19.0°N—49.6°N的范圍.實(shí)驗(yàn)過(guò)程中，取式(14)中p=0.3，其理由將在4.2節(jié)中作詳細(xì)討論和說(shuō)明.

圖2　臺(tái)站分布圖Fig.2　Illustration of monitoring stations

3.2數(shù)據(jù)預(yù)處理

地磁場(chǎng)的記錄數(shù)據(jù)中，含有一些季節(jié)性的積累誤差，在插值之前，應(yīng)該首先對(duì)這些成分進(jìn)行擬合、消除.Kriging與時(shí)空Kriging方法都需要滿足二次平穩(wěn)假設(shè)，因此，首先對(duì)每一個(gè)臺(tái)站的時(shí)間序列進(jìn)行正規(guī)化，以使數(shù)據(jù)獨(dú)立不相關(guān).

預(yù)處理分為兩步：

(1) 消除數(shù)據(jù)中的趨勢(shì)項(xiàng).各臺(tái)站的積累誤差都不相同，對(duì)各臺(tái)站分別進(jìn)行擬合，得到趨勢(shì)項(xiàng)，再進(jìn)行消除，從而得到余項(xiàng)εi(t)為

(16)

(2) 用各自的均值、方差對(duì)余項(xiàng)εi(t)進(jìn)行正規(guī)化，公式為

(17)

預(yù)處理的結(jié)果如圖3所示，分別以滿洲里地磁臺(tái)(117.4°E,49.6°N)和銀川地磁臺(tái)(106.8°E,38.3°N)數(shù)據(jù)為例進(jìn)行說(shuō)明.圖3a和圖3c中，數(shù)據(jù)的趨勢(shì)項(xiàng)可以通過(guò)3次多項(xiàng)式進(jìn)行擬合.對(duì)擬合的余項(xiàng)進(jìn)行正規(guī)化，結(jié)果如圖3b和圖3d.

預(yù)處理前后數(shù)據(jù)的自相關(guān)函數(shù)如圖4所示.自相關(guān)函數(shù)是平穩(wěn)的重要判據(jù)之一(李莎等，2012).在預(yù)處理之后，自相關(guān)函數(shù)在非零點(diǎn)處能夠迅速地減小，且速度遠(yuǎn)大于原始數(shù)據(jù).由此可見(jiàn)，預(yù)處理后數(shù)據(jù)的平穩(wěn)性增加，使得數(shù)據(jù)能更好地滿足時(shí)空Kriging的二次平穩(wěn)條件.

3.3構(gòu)建時(shí)空變差函數(shù)

3.3.1條件變差函數(shù)的擬合

時(shí)空變差函數(shù)的構(gòu)建，是整個(gè)工作的核心.變差函數(shù)類型的選擇以及擬合的精度，都會(huì)影響到時(shí)空Kriging插值的精度.根據(jù)式(10)和式(15)分別擬合兩種距離定義下的條件變差函數(shù).通常用到的變差函數(shù)有sphere, index, Gauss, power等，根據(jù)觀測(cè)值的分布情況選擇合適的函數(shù)進(jìn)行擬合(De Iaco et al., 2003；李莎等，2012).根據(jù)采樣數(shù)據(jù)的形狀，試驗(yàn)中空間變差函數(shù)均采用sphere函數(shù)，時(shí)間變差函數(shù)采用Gauss函數(shù)進(jìn)行擬合，如圖5所示.從圖中對(duì)比可以看出，圖5d中觀測(cè)值比圖5c中更加規(guī)則，使得擬合時(shí)估計(jì)值能夠更好地跟蹤到觀測(cè)值，擬合更準(zhǔn)確.

3.3.2時(shí)空變差函數(shù)的求取

本文采用Product-Sum模型，利用條件分布函數(shù)來(lái)構(gòu)建時(shí)空域變差函數(shù).其模型公式為

γst(hs,ht)=(k1Ct(0)+k2)γs(hs)

+(k2Cs(0)+k3)γt(ht)-k3γs(hs)γt(ht),其中，

3.4地磁圖制備

在時(shí)空Kriging方法中，取數(shù)據(jù)庫(kù)長(zhǎng)度為20天.當(dāng)計(jì)算某一個(gè)時(shí)間點(diǎn)t處的值時(shí)，時(shí)間區(qū)間[t-20,t-1]內(nèi)的所有觀測(cè)值都納入數(shù)據(jù)庫(kù).在87.2°E—126.6°E,19.0°N—49.6°N的范圍，取間隔1°進(jìn)行插值.若地理坐標(biāo)為A(Lon1,Lat1)和B(Lon2,Lat2)，兩點(diǎn)間地理距離為

圖3　地磁場(chǎng)數(shù)據(jù)的預(yù)處理(a) (b) 為滿洲里地磁臺(tái)MZL的預(yù)處理過(guò)程； (c) (d)為銀川地磁臺(tái)YCB的預(yù)處理過(guò)程. (a) (c)為原始數(shù)據(jù)及擬合的趨勢(shì)項(xiàng)； (b) (d)為擬合余項(xiàng)及其正規(guī)化結(jié)果.Fig.3　Preprocessing of geomagnetic field data(a) (b) illustrate the preprocessing of data set in monitoring station MZL, while (c) (d) illustrate the station YCB. (a) (c) are the original data and its fitted value, and (b) (d) are residual and its standardized data.

圖4　(a)滿洲里地磁臺(tái)MZL和(b)銀川地磁臺(tái)YCB數(shù)據(jù)的自相函數(shù)圖Fig.4　Autocorrelogram of data sets in monitoring stations (a) MZL and (b) YCB

圖5　變差函數(shù)擬合結(jié)果(a)(b)為空間變差函數(shù)擬合，(c)(d)為時(shí)間變差函數(shù)擬合.(a)(c)是時(shí)空Kriging的擬合，(b)(d)是改進(jìn)后時(shí)空Kriging的擬合結(jié)果.(d)中觀測(cè)值分布平穩(wěn)，能夠較好地集中在擬合曲線周?chē)?Fig.5　Illustration of fitted variogram functions(a)(b) depict the fitting of space conditional variogram. (c)(d) illustrate the fitting of time conditional variogram. (a)(c) are the results of Spatial Temporal Kriging while (b)(d) is of the improved one. In (d), distribution of sampled data is more stable, and concentrated near the fitted curve.

distance=R×arccos(sin(N1)sin(N2)

+cos(N1)cos(N2)cos(E2-E1)),

用不同方法插值，并用最近鄰法(Barber et al., 1996)和雙調(diào)和樣條插值(Sandwell, 1987)做對(duì)比試驗(yàn)，結(jié)果如圖6所示，依次為最近鄰、雙調(diào)和樣條、普通Kriging、時(shí)空Kriging、改進(jìn)時(shí)空Kriging插值所得的地磁圖.從結(jié)果可以看出，圖6a中最近鄰方法的分塊效應(yīng)明顯，其他方法更為平滑；圖6b中雙調(diào)和樣條方法結(jié)果出現(xiàn)了沿緯度方向的條帶，而忽視了地磁場(chǎng)在經(jīng)度方向的差異；圖6c中Kriging方法已可以較好地反映地磁場(chǎng)的分布，準(zhǔn)確性優(yōu)于前二者；通過(guò)圖6c和圖6d的對(duì)比可以看出，圖的左上角處，后者紋理更細(xì)致、分辨率更高，時(shí)空Kriging在處理插值中常見(jiàn)的邊界問(wèn)題時(shí)，效果優(yōu)于普通Kriging；圖6d與圖6e直觀效果差異不大.

4結(jié)果分析

均方差為(Mean Square Error, MSE)

4.1交叉驗(yàn)證

時(shí)間序列長(zhǎng)度為274天，計(jì)算時(shí)取某一天t的前20天[t-20.t-1]實(shí)測(cè)值作為時(shí)空Kriging的數(shù)據(jù)庫(kù).對(duì)[21,274]區(qū)間內(nèi)的每一天，分別用上述5種方法進(jìn)行插值的交叉驗(yàn)證，結(jié)果如圖7所示.不論MAE或MSE，最近鄰、雙調(diào)和樣條和Kriging的誤差都大于時(shí)空Kriging和改進(jìn)的時(shí)空Kriging.而改進(jìn)的時(shí)空Kriging比時(shí)空Kriging結(jié)果更穩(wěn)定，誤差變化幅度小，MAE集中在1左右，MSE集中在[0,5]的范圍內(nèi).

圖6　不同方法插值得到的地磁圖(a) 最近鄰; (b) 雙調(diào)和樣條; (c) Kriging; (d) 時(shí)空Kriging; (e) 改進(jìn)的時(shí)空Kriging.Fig.6　Interpolation image of different methods(a) Nearest; (b) Double harmonic spline; (c) Kriging; (d) Spatial Temporal Kriging; (e) Improved Spatial Temporal Kriging.

圖7　每個(gè)時(shí)間點(diǎn)處的交叉驗(yàn)證誤差Fig.7　Statistics of cross validation at each time

對(duì)圖7中各個(gè)時(shí)刻的MAE和MSE，沿時(shí)間軸對(duì)誤差進(jìn)行統(tǒng)計(jì)分析，其結(jié)果如表1所示.在引入了時(shí)間域的信息之后，時(shí)空Kriging比之前的空間域傳統(tǒng)插值方法，誤差減小.但時(shí)空Kriging誤差的方差比普通Kriging略有增加，反映出誤差的波動(dòng)較大.這種現(xiàn)象在經(jīng)過(guò)改進(jìn)之后得到解決，改進(jìn)后時(shí)空Kriging誤差的各項(xiàng)統(tǒng)計(jì)量均達(dá)到最小，插值的精度、穩(wěn)定性明顯提高.

表1　時(shí)間域的統(tǒng)計(jì)分析

圖8　權(quán)重p對(duì)時(shí)空Kriging的影響Fig.8　Choose of weight p in Improved Spatial Temporal Kriging

4.2改進(jìn)時(shí)空Kriging的權(quán)值選取

式(14)中，權(quán)重p決定了角度信息在向量距離中的比重，進(jìn)而影響到變差函數(shù)擬合的精度.取間隔0.05，在0≤p≤1的范圍內(nèi)考察權(quán)重對(duì)擬合、插值的影響，結(jié)果如圖8所示.當(dāng)p=0時(shí)，向量距離中只有角度信息，幅值信息的缺乏導(dǎo)致大的誤差.之后，誤差隨p的增大而減??；p=0.3時(shí)，插值誤差達(dá)到最小，此時(shí),Mean of MAEs=0.976177666， Mean of MSEs=1.961507859；之后，誤差隨p增大而增大，p=1時(shí)，角度信息缺失，此時(shí)向量距離退化為歐式距離，插值結(jié)果與時(shí)空Kriging相同.在p=0.3處，幅值信息與角度信息達(dá)到平衡，使得時(shí)間切片間和空間切片間的差異在向量距離中有合適的表達(dá)，從而提高變差函數(shù)擬合與時(shí)空Kriging插值的精度.因此，在本文的實(shí)驗(yàn)中，均采用p=0.3.

5結(jié)論

對(duì)Kriging方法進(jìn)行時(shí)間-空間域的擴(kuò)展，其核心是對(duì)其變差函數(shù)進(jìn)行時(shí)間-空間域的擴(kuò)展.在加入了時(shí)間域信息之后，時(shí)空Kriging方法的交叉驗(yàn)證誤差減小，插值精度提高；在Kriging變差函數(shù)的時(shí)空域擴(kuò)展過(guò)程中，將角度信息引入向量距離，從而提出了改進(jìn)的時(shí)空Kriging算法.實(shí)驗(yàn)結(jié)果表明，新的向量距離能夠更有效地表征數(shù)據(jù)，改進(jìn)的時(shí)空Kriging方法插值精度進(jìn)一步得到提高.

致謝文中用到了中國(guó)地磁臺(tái)網(wǎng)等32個(gè)地磁臺(tái)站的地磁場(chǎng)數(shù)據(jù)，在此一并表示感謝.

References

Astola J, Haavisto J, Neuvo Y. 1990. Vector median filters. Proceedings of the IEEE, 78(4): 678-689.

Cichota R, Hurtado A L B, de Jong van Lier Q. 2006. Spatio-temporal variability of soil water tension in a tropical soil in Brazil.Geoderma, 133(3-4): 231-243.

Cressie N. 1988. Spatial prediction and ordinary kriging.MathematicalGeology, 20(4): 405-421.

Cressie N. 1990. The origins of kriging.MathematicalGeology, 22(3): 239-252.

Cressie N, Huang H C. 1999. Classes of nonseparable, spatio-temporal stationary covariance functions.JournaloftheAmericanStatisticalAssociation, 94(448): 1330-1339.Davis B M, Borgman L E. 1982. A note on the asymptotic distribution of the sample variogram.MathematicalGeology, 14(2): 189-193.

De Iaco S, Myers D E, Posa D. 2002. Nonseparable space-time covariance models: some parametric families.MathematicalGeology, 34(1): 23-42.

De Iaco S, Myers D E, Posa D. 2003. The linear coregionalization model and the product-sum space-time variogram.MathematicalGeology, 35(1): 25-38. Du R Q, Yang J. 2012. Application of ordinary kriging method in data processing of magnetic survey.∥Proceedings of the 2012 7th International Conference on Computer Science & Education. Melbourne, Australia: IEEE, 771-774.

Gething P W, Akinson P M, Noor A M, et al. 2007. A local space-time kriging approach applied to a national outpatient malaria data set.Computers&Geosciences, 33(10): 1337-1350.

Gneiting T. 2002. Nonseparable, stationary covariance functions for space-time data.JournaloftheAmericanStatisticalAssociation, 97(458): 590-600. Goldenberg F. 2006. Geomagnetic navigation beyond the magnetic compass.∥ Proceedings of the 2006 IEEE/ION Position, Location, and Navigation Symposium. Coronado, CA: IEEE, 684-694.

Guo C F, Zhang L J, Cai H. 2011. Analysis of geomagnetic navigation accuracy under magnetic storms.Chin.J.SpaceSci. (in Chinese), 31(3): 372-377.

Huang H, Liu Q, Wang J, et al. 2012. Kriging-based study on the visualization of magnetic method data.∥ Proceedings of the 2012 International Conference on Image Analysis and Signal Processing. Hangzhou, China: IEEE, 1-6.

Jost G, Heuvelink G B M, Papritz A. 2005. Analysing the space-time distribution of soil water storage of a forest ecosystem using spatio-temporal kriging.Geoderma, 128(3-4): 258-273.

Journel A G. 1986. Geostatistics: Models and tools for the earth sciences.MathematicalGeology, 18(1): 119-140.

Knotters M, Brus D J, Oude Voshaar J H. 1995. A comparison of kriging, co-kriging and kriging combined with regression for spatial interpolation of horizon depth with censored observations.Geoderma, 67(3-4): 227-246. Lark R M, Bellamy P H, Rawlins B G. 2006. Spatio-temporal variability of some metal concentrations in the soil of eastern England, and implications for soil monitoring.Geoderma, 133(3-4): 363-379.

Laurenceau J, Sagaut P. 2008. Building efficient response surfaces of aerodynamic functions with kriging and cokriging.AIAAJournal, 46(2): 498-507.

Li S, Shu H, Dong L. 2011. Research and realization of Kriging interpolation based on spatial-temporal variogram.ComputerEngineeringandApplications(in Chinese), 47(23): 25-26.

Li S, Shu H, Xu Z. 2012. Interpolation of temperature based on spatial-temporal kriging.GeomaticsandInformationScienceofWuhanUniversity(in Chinese), 37(2): 237-241.

Li X H, Liu G, Liu D Z, et al. 2009. Discrimination of nuclear explosions and earthquakes using the nearest support vector feature line fusion classification algorithm.ChineseJ.Geophys. (in Chinese), 52(7): 1816-1824, doi: 10.3969/j.issn.0001-5733.2009.07.016.Liu D Z, Wang R M, Li X H, et al. 2005. Onset point identification of single-channel seismic signal based on wavelet packet and the AR model.ChineseJ.Geophys. (in Chinese), 48(5): 1098-1102.Niu C, Liu D Z, Li X H, et al. 2014. Feature analysis of geomagnetic variation time series based on complexity theory.ChineseJ.Geophys. (in Chinese), 57(8): 2713-2723, doi: 10.6038/cjg20140829. Nunes C, Soares A. 2005. Geostatistical space-time simulation model for air quality prediction.Environmetrics, 16(4): 393-404.

Porcu E, Mateu J, Saura F. 2008. New classes of covariance and spectral density functions for spatio-temporal modelling.Stoch.Environ.Res.RiskAssess., 22(S1): 65-79.

Rouhani S, Myers D E. 1990. Problems in space-time kriging of geohydrological data.MathematicalGeology, 22(5): 611-623.

Sandwell D T. 1987. Biharmonic spline interpolation of GEOS-3 and SEASAT altimeter data.GeophysicalResearchLetters, 14(2): 139-142.

Snepvangers J J J C, Heuvelink G B M, Huisman J A. 2003. Soil water content interpolation using spatio-temporal kriging with external drift.Geoderma, 112(3-4): 253-271.

Song Y. 2013. Construction and application of the spatio-temporal variability function of groundwater regionalized variable[Master thesis] (in Chinese). Beijing: Beijing Forestry University.

Song Y, Zhang Z, Zhang H X. 2013a. The simulation parameter interpolation of spatial and temporal variability of groundwater.JournalofAridLandResourcesandEnvironment(in Chinese), 27(8): 120-124.Song Y, Zhang Z, Zhang H X, et al. 2013b. Evaluation of the spatial and temporal variability of environmental chloride ion in groundwater by simulation.ChineseJournalofUndergroundSpaceandEngineering(in Chinese), 9(2): 457-461.

Sun H, Feng Y, An Z Z, et al. 2010. Analysis of variation in geomagnetic field of Chinese mainland based on comprehensive model CM4.ActaPhysicaSinica(in Chinese), 59(12): 8941-8953.

Tang Y B. 2002. Comparison of semivariogram models for Kriging monthly rainfall in eastern China.JournalofZhejiangUniversityScience, 3(5): 584-590.Trahanias P E, Karakos D G, Venetsanopoulos A N. 1996. Directional processing of color images: theory and experimental results.IEEETransactionsonImageProcessing, 5(6): 868-880.

Wang J M, Zhang J, Deng Z B, et al. 2014. Slope deformation analyses with space-time Kriging interpolation method.JournalofChinaCoalSociety(in Chinese), 39(5): 874-879.

Xu A P, Hu L, Shu H. 2011. Extension and implementation from spatial-only to spatiotemporal Kriging interpolation.JournalofComputerApplications(in Chinese), 31(1): 273-276.

Xu W Y. 2003. Geomagnetism (in Chinese). Beijing: Seismological Press.

Xu W Y. 2009. Physics of Electromagnetic Phenomena of the Earth (in Chinese). Hefei: University of Science and Technology of China Press.

Yi S H, Liu D Z, He Y L, et al. 2013. Modeling and forecasting of the variable geomagnetic field by support vector machine.ChineseJ.Geophys. (in Chinese), 56(1): 127-135, doi: 10.6038/cjg20130113.

Yuan Y H. 2012. The research of the geomagnetic variation field in the geomagnetic navigation[Ph. D thesis] (in Chinese). Wuhan: Huazhong University of Science and Technology.Zeng X N, Li X H, Liu D Z, et al. 2011. Regularization analysis of integral iteration method and the choice of its optimal step-length.ChineseJ.Geophys. (in Chinese), 54(11): 2943-2950, doi: 10.3969/j.issn.0001-5733.2011.11.024. Zeng X N, Liu D Z, Li X H, et al. 2014. An Improved iterative Wiener filter for downward continuation of potential fields.ChineseJ.Geophys. (in Chinese), 57(6): 1958-1967, doi: 10.6038/cjg20140626.

Zhang Q Y, Wu J T. 2015. Image super-resolution using windowed ordinary Kriging interpolation.OpticsCommunications, 336: 140-145.

Zhang X M, Zhao Y. 2009. Local geomagnetic field mapping based on kriging interpolation method.ElectronicMeasurementTechnology(in Chinese), 32(4): 122-125.

Zhao M, Zeng X P, Lin Y F. 1997. Nonlinear analysis of magnetic storm in March 13, 1989.ObservationandResearchonSeismologicalandGeomagnetism(in Chinese), 18(3): 56-60.

附中文參考文獻(xiàn)

郭才發(fā), 張力軍, 蔡洪. 2011. 磁暴期間的地磁導(dǎo)航精度分析. 空間科學(xué)學(xué)報(bào), 31(3): 372-377.

李莎, 舒紅, 董林. 2011. 基于時(shí)空變異函數(shù)的Kriging插值及實(shí)現(xiàn). 計(jì)算機(jī)工程與應(yīng)用. 47(23): 25-26.

李莎, 舒紅, 徐正全. 2012. 利用時(shí)空Kriging進(jìn)行氣溫插值研究. 武漢大學(xué)學(xué)報(bào). 信息科學(xué)版. 37(2): 237-241.

李夕海, 劉剛, 劉代志等. 2009. 基于最近鄰支撐向量特征線融合算法的核爆地震識(shí)別. 地球物理學(xué)報(bào), 52(7): 1816-1824, doi:

10.3969/j.issn.0001-5733.2009.07.016.

劉代志, 王仁明, 李夕海等. 2005. 基于小波包分解及AR模型的單通道地震波信號(hào)初至點(diǎn)檢測(cè). 地球物理學(xué)報(bào), 48(5): 1098-1102.

牛超, 劉代志, 李夕海等. 2014. 基于復(fù)雜度理論的地磁變化場(chǎng)時(shí)間序列特征分析. 地球物理學(xué)報(bào), 57(8): 2713-2723, doi: 10.6038/cjg20140829.

宋瑩. 2013. 地下水區(qū)域化變量時(shí)空變異函數(shù)構(gòu)建及其應(yīng)用研究[碩士論文]. 北京: 北京林業(yè)大學(xué).

孫涵, 馮彥, 安振昌等. 2010. 利用地磁場(chǎng)綜合模型CM4分析中國(guó)大陸地區(qū)地磁場(chǎng)變化. 物理學(xué)報(bào), 59(12): 8941-8953.

王建民, 張錦, 鄧增兵等. 2014. 時(shí)空Kriging插值在邊坡變形監(jiān)測(cè)中的應(yīng)用. 煤炭學(xué)報(bào), 39(5): 874-879.

徐愛(ài)萍, 胡力, 舒紅. 2011. 空間克里金插值的時(shí)空擴(kuò)展與實(shí)現(xiàn). 計(jì)算機(jī)應(yīng)用, 31(1): 273-276.

徐文耀. 2003. 地磁學(xué). 北京: 地震出版社.

徐文耀. 2009. 地球電磁現(xiàn)象物理學(xué). 合肥: 中國(guó)科學(xué)技術(shù)大學(xué)出版社.

易世華, 劉代志, 何元磊等. 2013. 變化地磁場(chǎng)預(yù)測(cè)的支持向量機(jī)建模. 地球物理學(xué)報(bào), 56(1): 127-135, doi: 10.6038/cjg20130113.袁楊輝. 2012. 地磁導(dǎo)航中地磁變化場(chǎng)的研究[博士論文]. 武漢: 華中科技大學(xué).

曾小牛, 李夕海, 劉代志等. 2011. 積分迭代法的正則性分析及其最優(yōu)步長(zhǎng)的選擇. 地球物理學(xué)報(bào), 54(11): 2943-2950, doi: 10.3969/j.issn.0001-5733.2011.11.024.

曾小牛, 劉代志, 李夕海等. 2014. 位場(chǎng)向下延拓的改進(jìn)迭代維納濾波法. 地球物理學(xué)報(bào), 57(6): 1958-1967, doi: 10.6038/cjg20140626.

張曉明, 趙剡. 2009. 基于克里金插值的局部地磁圖的構(gòu)建. 電子測(cè)量技術(shù), 32(4): 122-125.

趙明, 曾小蘋(píng), 林云芳. 1997. 對(duì)1989年3月13日磁暴的非線性分析. 地震地磁觀測(cè)與研究, 18(3): 56-60.

(本文編輯劉少華)

Application and improvement of spatial temporal Kriging in geomagnetic field interpolation

CHEN Ding-Xin, LIU Dai-Zhi*, ZENG Xiao-Niu, MENG Liang， YANG Xiao-Jun

Staffroom907,PLARocketForceEngineeringUniversity,Xi′an710025,China

AbstractThis paper utilized Spatial Temporal Kriging method to improve the interpolation of regional geomagnetic field, by taking time domain information into consideration. Direction information was added to the variogram function of Spatial Temporal Kriging. A new kind of vector distance was defined to improve the method, meanwhile, selection of the weight in vector distance was discussed. Interpolation of geomagnetic field data from 32 monitoring stations in the region of longitude 87.2°E—126.6°E, latitude 19.0°N—49.6°N illustrated that, the results of methods which made use of information in time domain were much better than the traditional ones. Changes of geomagnetic field along longitude direction and latitude direction were both expressed adequately, while the boundary effect decreased. Illustration of the cross validation of geomagnetic field data indicated that, with the addition of time domain information and vector direction, the performance of Kriging method was improved in sequence. The improved Spatial Temporal Kriging method had the most accurate and stable performance in contrast test, with MAE concentrated at 1nT, and MSE concentrated in the interval of [0,5].

KeywordsGeomagnetic field; Kriging; Spatial Temporal Kriging; Variogram; Vector direction

基金項(xiàng)目國(guó)家自然科學(xué)基金(41374154, 61304240)及中國(guó)博士后科學(xué)基金(2014M552589)聯(lián)合資助.

作者簡(jiǎn)介陳鼎新，男，1986年生，在讀博士研究生，主要從事地磁信號(hào)處理研究.E-mail:chendx12@mails.tsinghua.edu.cn *通訊作者劉代志，男，1960年生，博士生導(dǎo)師，主要從事國(guó)家安全地球物理研究.E-mail:daizhiliu@163.com

doi:10.6038/cjg20160518 中圖分類號(hào)P318

收稿日期2015-07-30, 2016-03-31收修定稿

陳鼎新， 劉代志， 曾小牛等. 2016. 時(shí)空Kriging算法在區(qū)域地磁場(chǎng)插值中的應(yīng)用及改進(jìn).地球物理學(xué)報(bào)，59(5):1743-1752,doi:10.6038/cjg20160518.

Chen D X, Liu D Z, Zeng X N, et al. 2016. Application and improvement of spatial temporal Kriging in geomagnetic field interpolation.ChineseJ.Geophys. (in Chinese),59(5):1743-1752,doi:10.6038/cjg20160518.