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

火箭軍工程大學907教研室， 西安　710025

時空Kriging算法在區(qū)域地磁場插值中的應(yīng)用及改進

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

火箭軍工程大學907教研室， 西安710025

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

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

1引言

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

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

時空Kriging方法對Kriging進行了時空域擴展，從而充分利用數(shù)據(jù)中的時間域信息來進行插值.目前，時空Kriging方法已經(jīng)在空氣質(zhì)量檢測(De Iaco et al., 2003; Nunes and Soares, 2005)、雨水建模(李莎等，2011; li et al., 2011; Tang, 2002)、風力數(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)甚至用時空Kriging來分析瘧疾門診病人的數(shù)據(jù)，將時空Kriging應(yīng)用到社科數(shù)據(jù)的分析當中，并取得了良好效果.然而，在地磁學領(lǐng)域，時空Kriging方法至今沒有應(yīng)用，使得插值過程中時間域信息被忽視.

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

2時空Kriging算法及其改進

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

2.1普通Kriging

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

(1)

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

(2)

定義變差函數(shù)為

(3)

由于

(4)

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

(5)

寫成矩陣形式為

(6)

其中，

2.2時空Kriging

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

(7)

-Z(s,h)]2,

(8)

則可以通過Kriging方程組

(9)

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

(10)

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

(11)

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

2.3改進時空Kriging算法

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

圖1　Kriging時空擴展演示圖Fig.1　Illustration of expiation of Spatial Temporal Kriging

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

(12)

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

L2范數(shù)的定義為

(13)

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

(14)

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

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

(15)

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

3實驗

3.1數(shù)據(jù)

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

圖2　臺站分布圖Fig.2　Illustration of monitoring stations

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

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

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

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

(16)

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

(17)

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

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

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

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

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

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

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

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

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

3.4地磁圖制備

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

圖3　地磁場數(shù)據(jù)的預(yù)處理(a) (b) 為滿洲里地磁臺MZL的預(yù)處理過程； (c) (d)為銀川地磁臺YCB的預(yù)處理過程. (a) (c)為原始數(shù)據(jù)及擬合的趨勢項； (b) (d)為擬合余項及其正規(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)滿洲里地磁臺MZL和(b)銀川地磁臺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ù)擬合.(a)(c)是時空Kriging的擬合，(b)(d)是改進后時空Kriging的擬合結(jié)果.(d)中觀測值分布平穩(wěn)，能夠較好地集中在擬合曲線周圍.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)做對比試驗，結(jié)果如圖6所示，依次為最近鄰、雙調(diào)和樣條、普通Kriging、時空Kriging、改進時空Kriging插值所得的地磁圖.從結(jié)果可以看出，圖6a中最近鄰方法的分塊效應(yīng)明顯，其他方法更為平滑；圖6b中雙調(diào)和樣條方法結(jié)果出現(xiàn)了沿緯度方向的條帶，而忽視了地磁場在經(jīng)度方向的差異；圖6c中Kriging方法已可以較好地反映地磁場的分布，準確性優(yōu)于前二者；通過圖6c和圖6d的對比可以看出，圖的左上角處，后者紋理更細致、分辨率更高，時空Kriging在處理插值中常見的邊界問題時，效果優(yōu)于普通Kriging；圖6d與圖6e直觀效果差異不大.

4結(jié)果分析

均方差為(Mean Square Error, MSE)

4.1交叉驗證

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

圖6　不同方法插值得到的地磁圖(a) 最近鄰; (b) 雙調(diào)和樣條; (c) Kriging; (d) 時空Kriging; (e) 改進的時空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　每個時間點處的交叉驗證誤差Fig.7　Statistics of cross validation at each time

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

表1　時間域的統(tǒng)計分析

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

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

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

5結(jié)論

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

致謝文中用到了中國地磁臺網(wǎng)等32個地磁臺站的地磁場數(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.

附中文參考文獻

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

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

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

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

10.3969/j.issn.0001-5733.2009.07.016.

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

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

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

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

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

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

徐文耀. 2003. 地磁學. 北京: 地震出版社.

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

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

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

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

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

趙明, 曾小蘋, 林云芳. 1997. 對1989年3月13日磁暴的非線性分析. 地震地磁觀測與研究, 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

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

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

doi:10.6038/cjg20160518 中圖分類號P318

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

陳鼎新， 劉代志， 曾小牛等. 2016. 時空Kriging算法在區(qū)域地磁場插值中的應(yīng)用及改進.地球物理學報，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.