張永奇,韓美濤,韓曉飛,鄭增記,翟宏光

(1.陜西省地震局，陜西 西安 710068；2.中國地震局地質(zhì)研究所 地震動力學國家重點實驗室，北京 100029)

Kriging法在區(qū)域重力場插值中的適用性研究

張永奇1,2,韓美濤1,韓曉飛1,鄭增記1,翟宏光1

(1.陜西省地震局，陜西 西安 710068；2.中國地震局地質(zhì)研究所 地震動力學國家重點實驗室，北京 100029)

由于重力觀測網(wǎng)點位分布不均勻，密度不夠，在精細刻畫區(qū)域重力場變化時受到很大限制，因此需采用合理的插值算法。文中介紹Kriging插值算法和變異函數(shù)理論模型的基礎上，結合陜西地區(qū)2011—2012年離散重力變化數(shù)據(jù)，利用最小二乘法得到各個理論變異函數(shù)模型的擬合參數(shù)，并將實驗變異函數(shù)模型用于Kriging插值算法，同時綜合考慮交叉驗證精度統(tǒng)計、插值精度統(tǒng)計結果，以此來研究Kriging插值算法的適用性，研究結果表明：基于最小二乘法獲得的球形模型、指數(shù)模型擬合參數(shù)精度高，高斯模型、對數(shù)模型精度最差；基于球形模型、指數(shù)模型的實驗變異函數(shù)用于Kriging插值算法得到的插值結果精度最高、圖像平滑自然、異常區(qū)明顯，是一種適合陜西地區(qū)重力離散數(shù)據(jù)進行插值計算的有效算法。

Kriging法；變異函數(shù)模型；陜西地區(qū)；重力場；插值

重力資料在地震長期預測預報中占有重要地位，合理的重力網(wǎng)布設是獲取真實可靠資料的基礎。但由于重力網(wǎng)布設易受各種因素的限制，點位分布往往呈現(xiàn)不均勻、密度不夠等特點，這對于精細刻畫某個區(qū)域重力場變化來說，顯得力不從心。因此，針對離散觀測值采用合理的插值算法進行插值計算就成為一個重要的解決途徑。目前，重力場插值應用較多的算法有最小二乘配置法、樣條函數(shù)法、Kriging法等[1-3]。Kriging法是一種依據(jù)已知點進行空間插值的地質(zhì)學統(tǒng)計方法，該方法是由南非礦產(chǎn)地理學家P.G.Krige最先提出，該法綜合考慮空間點的結構性和隨機性，估值具有最優(yōu)、線性、無偏等特點，在地學領域得到廣泛的應用，并取得很好的效果[4-9]。目前，基于Kriging法進行數(shù)據(jù)插值的研究主要集中在幾種插值算法的比較，研究表明Kriging法具有明顯優(yōu)越性[3,6,9-10]；基于不同變異函數(shù)的Kriging法進行重力場插值的研究相對較少。變異函數(shù)是Kriging法進行插值的核心，常用的變異函數(shù)理論模型有球面模型、高斯模型、指數(shù)模型等[12-14]，變異函數(shù)模型選擇是否恰當，直接影響插值效果的好壞，甚至影響計算結果的可靠性，因此對不同變異函數(shù)模型的適用性進行研究，有利于提高Kriging法插值的精度和可靠性?；谏鲜龇治觯疚脑诮榻BKriging法及變異函數(shù)模型的基礎上，以陜西區(qū)域離散重力變化數(shù)據(jù)作為研究對象，對基于不同變異函數(shù)的Kriging法插值結果進行精度統(tǒng)計，以此評價變異函數(shù)模型的適用性，最終得到適合該區(qū)域重力場插值的可靠算法。

1 數(shù)據(jù)處理理論與模型

1.1 普通Kriging法

假定Z(x)服從二階平穩(wěn)假設，其在空間位置x1,x2,…,xn有n個觀測值Z(xi),i=1,2,…,n,則在未知點x0處，其觀測值Z(x0)可由已知的n個觀測值Z(xi)估計。

(1)

式中:λi為待定的權重系數(shù)，為達到線性無偏估計，使估計方差最小，權系數(shù)λi可由Kriging方程組決定。

(2)

式中:μ是拉格朗日乘數(shù);γij表示第i,j個觀測點間向量所對應的變異函數(shù)值;γi0表示第i個觀測點到估值點間向量對應的變異函數(shù)值，這些量可通過實驗變異函數(shù)擬合理論模型計算得到。在式(2)約束條件下，求出λj，再回代式(1)，即可根據(jù)已知點觀測數(shù)據(jù)內(nèi)插出待定點的觀測估值[14]。

1.2 變異函數(shù)

(3)

式中:γ*(h)為理論變異函數(shù)值γ(h)的估計值;N為被向量h相隔的實驗數(shù)據(jù)對的數(shù)目。

按式(3)計算出空間相距h距離的變異函數(shù)值γ*(h)后，再選擇合適的變異函數(shù)理論模型進行擬合，進而求出適合研究區(qū)域的變異函數(shù)理論模型。

1.3 變異函數(shù)理論模型

變異曲線圖表示的是一定滯后距h的變異函數(shù)值γ*(h)與該h的對應圖，如圖1所示。圖中C0稱為塊金效應，表示h很小時兩點x和x+h之間空間信息的變化；a稱為變程，其大小反映了研究對象中某一區(qū)域變化量的變化程度。CS稱為總基臺值，反映了某區(qū)域化變量在研究范圍內(nèi)變異的強度，它是最大滯后距的可遷移性變異函數(shù)的極限值，當h→∞時，γ(∞)=C(0)=Var[Z(x)]=CS，即當h→∞時，變異函數(shù)值近于先驗方差C(0)，當無塊金效應C0時，CS=C，當有塊金效應時，CS=C+C0；而C稱為基臺值，它是先驗方差與塊金效應(常數(shù))之變差，C=CS-C[14]。

圖1 變異曲線圖

變異函數(shù)理論模型分為有基臺值和無基臺值兩大類。其中有基臺類模型包括球形模型、高斯模型、指數(shù)模型、純塊金模型、線性模型(有基臺)，有基臺模型一般認為區(qū)域變量滿足二階平穩(wěn)假設；無基臺類模型包括冪函數(shù)模型、對數(shù)模型、線性模型(無基臺)等，一般認為只滿足本征假設[14]。具體模型如下：

1)球面模型:

(4)

2)高斯模型:

(5)

3)指數(shù)模型:

(6)

4)純塊金效應模型:

(7)

該模型對純隨機變量才適用。

5)有基臺的線性模型:

(8)

6)冪模型:

γ(h)=hθ,0<θ<2.

(9)

當θ=1，即為無基臺線性模型，γ(h)=ah,h>0，a為常數(shù)。

7)對數(shù)模型:

γ(h)=log(h).

(10)

由于h→0時，log(h)→-∞，這與變異函數(shù)性質(zhì)不符合，因此在區(qū)域化變量為離散點時不適用。

利用式(3)求出各γ(h)，按最小二乘擬合法求出式(4)～式(10)中的參數(shù)C0，C，a，則可確定出樣本空間的實驗變異函數(shù)模型,也就是適合研究區(qū)域的變異函數(shù)模型[11-14]。

2 實例及結果分析

2.1 研究區(qū)概況

本文選取陜西地區(qū)2011～2012年離散重力變化數(shù)據(jù)作為研究對象，重力網(wǎng)范圍為32.42°～36.25°N，106.15°～110.56°E，總共包括221個點，其中關中地區(qū)測點密度最大，平均點間距約10～15 km；陜南地區(qū)測點密度次之，平均點間距約20～30 km；陜北地區(qū)點位分布最為稀疏，平均點間距達到40～50 km。陜西地區(qū)地形地貌差異較大，陜北地區(qū)以黃土塬為主，海拔大致為800～1 300 m；關中地區(qū)為斷陷盆地，平均海拔約520 m；陜南地區(qū)為秦巴山地，海拔基本上為1 500～2 000 m。為了同時顧及點位分布不均以及地形差異的影響，本文選擇Kriging法對離散觀測值進行插值，其目的既考慮觀測值的空間性也考慮觀測點間的相關性。通過對不同變異函數(shù)模型的參數(shù)進行擬合，獲得不同的實驗變異函數(shù)，并將其用于Kriging插值算法，通過精度對比分析，得到可靠性和精度都比較滿意的Kriging插值算法。

2.2 研究結果與分析

首先利用研究區(qū)分布的211個離散點數(shù)據(jù)，建立適合本研究區(qū)的實驗變異函數(shù)。根據(jù)離散點空間分布特征及其相關性，獲得實驗變異函數(shù)參數(shù)，具體參數(shù)見表1。由表1知，變異函數(shù)最大滯后距為1.9°，滯后組數(shù)劃分為25組，滯后寬度為0.076°，變異函數(shù)高值、滯后方向、滯后方向公差分別為644、60°、90°。其次求取變異函數(shù)模型參數(shù)，選擇變異函數(shù)理論模型(球形模型、指數(shù)模型、冪函數(shù)模型、塊金效應模型等)，并給定初始參數(shù)值，利用最小二乘擬合法對各模型參數(shù)進行迭代處理，設置迭代次數(shù)為50，經(jīng)過計算獲得各組合變異函數(shù)模型(塊金效應+其它模型)參數(shù)并繪制擬合圖，如表2、圖2所示。表2中參數(shù)就是實驗變異函數(shù)模型所需的參數(shù)C0，C，a，將參數(shù)代入模型(4)～模型(10)，就可以獲得適合研究區(qū)的變異函數(shù)模型，根據(jù)模型(1)～模型(3)就可以求取研究區(qū)任一點的重力異常值。從圖2直觀可知，球形模型、指數(shù)模型、線性模型、冪函數(shù)模型擬合效果較好，高斯模型、對數(shù)模型擬合效果較差，研究表明對數(shù)模型不適合進行離散點擬合[15]。

表1 實驗變異函數(shù)參數(shù)

表2 不同變異函數(shù)模型擬合參數(shù)

備注：線性模型斜率slope=225;冪函數(shù)模型因變量θ=0.692

圖2 不同模型的變異函數(shù)擬合圖

為了進一步比較不同變異函數(shù)模型的適用性，采用交叉驗證的方法進行分析。將221個離散點數(shù)據(jù)全部作為檢核點，計算擬合點值與離散數(shù)據(jù)點值的殘差精度以及擬合點自身精度，精度統(tǒng)計如表3所示。分析表3發(fā)現(xiàn)，從擬合點的均方根來分析，精度從高到低順序為：高斯模型>指數(shù)模型>球形模型=冪函數(shù)模型>線性模型>對數(shù)模型；但從擬合點值與離散數(shù)據(jù)點值的殘差均方根來分析，精度從高到低順序為：球形模型=冪函數(shù)模型>線性模型>指數(shù)模型>高斯模型>對數(shù)模型；從標準差來分析，也可以得到和均方根大致相同的結論，因此綜合4個精度指標考慮，球形模型、冪函數(shù)模型擬合精度較高，更加適合研究區(qū)離散數(shù)據(jù)Kriging法插值，而高斯模型、對數(shù)模型擬合精度最差。

表3 不同變異函數(shù)模型的交差驗證精度統(tǒng)計

利用上述各個實驗變異函數(shù)模型進行Kriging法插值，總共獲得8 700個節(jié)點的插值數(shù)據(jù)，緯度方向插值節(jié)點為100行，經(jīng)度方向插值節(jié)點為87列。根據(jù)插值結果繪制插值等值線圖，超出插值邊界最大值的數(shù)據(jù)采用白化值，如圖3所示，圖3中(a)～(f)分別代表球形模型、高斯模型、指數(shù)模型、線性模型、冪函數(shù)模型、對數(shù)模型的插值結果。同時統(tǒng)計基于不同變異函數(shù)模型對應的Kriging法插值精度，具體統(tǒng)計指標有插值花費時間、插值最大值、最小值、均值、標準差，具體見表4。從圖3直觀看出，(a)、(c)～(e)結果較好，很好地體現(xiàn)了重力場異常變化的細節(jié)特征，是對離散觀測值的有效插值和擬合，而(b)、(f)插值效果最差，(b)反映的是研究區(qū)的趨勢變化，其中的細節(jié)體現(xiàn)不完整，而(f)完全不適合離散點插值，對重力場異常變化反映不明顯，而且量值遠遠大于采用其它模型進行插值的結果。從細節(jié)上分析，基于球形模型和指數(shù)模型的插值結果，“牛眼”較少，曲線比較光滑，由南向北過渡比較自然；而基于線性模型和冪函數(shù)模型的插值結果反映細節(jié)比較清晰，形成很多局部異常，因此基于這兩種模型得到的插值結果最好是進行高通濾波，得到更加真實的異常變化圖。

由表4知，幾種模型花費的時間差別不大，對數(shù)模型用時最多，精度最差。從插值的最大、最小值來看都沒有超過原始離散點的最大值72.7、最小值-52.7，表明插值數(shù)據(jù)可靠性較高，只進行內(nèi)插計算沒有外推計算。從標準差來看高斯模型精度最高，但這是一種“虛假”的高精度，并不適合研究區(qū)的重力異常插值；其余5種模型中球形模型和指數(shù)模型精度較高，與圖3結果表現(xiàn)一致，是真實反映插值效果的指標；冪函數(shù)模型、線性模型精度稍差，應該是局部區(qū)域插值精度差導致的。

表4 不同變異函數(shù)模型的插值精度統(tǒng)計

圖3 基于不同變異函數(shù)模型的kriging重力場插值圖(ugal)

3 結 論

本文針對陜西地區(qū)2011～2012年211個離散重力變化值，基于不同變異函數(shù)模型構建Kriging插值法進行模型適用性研究，得出如下結論：

1)對于本次試驗數(shù)據(jù)，最大滯后距為1.9°，滯后組數(shù)分為25組，滯后方向為60°；經(jīng)反復測試變異函數(shù)模型選擇純塊金模型和其它模型的組合，擬合效果較好，單獨模型擬合效果不佳；且利用最小二乘擬合迭代計算，迭代次數(shù)設置為50次，基本上可以獲得最優(yōu)解。

2)利用最小二乘擬合獲得的各種模型的參數(shù)差異較大，同時發(fā)現(xiàn)數(shù)據(jù)具有一定程度的各向異性，量值皆為2，各向異性方向變化范圍為15.61°～16.33°。變異函數(shù)模型參數(shù)精度統(tǒng)計使用交叉驗證方法，交叉驗證時采用全部離散點數(shù)據(jù)，其優(yōu)點是避免由于選擇部分離散點而帶入人為的統(tǒng)計偏差。交叉驗證結果顯示球形模型、冪函數(shù)模型精度最高，均方根和標準誤差均最小，說明這兩種模型建立參數(shù)選取合理，適合用于陜西地區(qū)離散重力變化網(wǎng)格化插值。

3)利用基于不同變異函數(shù)模型的Kriging法進行插值時，考慮插值曲線的光滑性，整個研究區(qū)域劃分為8 700個節(jié)點，插值得到各個節(jié)點的重力變化值。為了減小邊緣效應的影響，超出插值最大值的區(qū)域采用白化值。由圖3、表4可以得出：采用高斯模型、對數(shù)模型的Kriging法進行插值效果最差，不適合陜西地區(qū)離散重力變化插值；插值效果最好的是采用球形模型、指數(shù)模型，精度最高、誤差最小、圖像顯示曲線光滑自然，且沒有出現(xiàn)插值外推現(xiàn)象。

[1] 祝意青，胡斌，李輝，等.用雙三次樣條函數(shù)模擬青藏高原東北緣重力場動態(tài)圖像[J].中國地震，2005，21(2)：165-171.

[2] 祝意青，陳斌，張希，等.景泰5.9級地震前后的重力變化研究[J].中國地震，2001，17(4)：356-363.

[3] 李振海，汪海洪.重力數(shù)據(jù)網(wǎng)格化方法比較[J].大地測量與地球動力學，2010，30(1)：140-144.

[4] 陳小斌.中國陸地現(xiàn)今水平形變狀況及其驅(qū)動機制[J].中國科學(D輯：地球科學)，2007，37(8)：1056-1064.

[5] 朱守彪，蔡永恩，石耀霖,等.青藏高原及鄰區(qū)現(xiàn)今地應變率場的計算及其結果的地球動力學意義[J].地球物理學報，2005，48(5)：1053-1061.

[6] 曾懷恩，黃聲享.基于Kriging方法的空間數(shù)據(jù)插值研究[J].測繪工程，2007，16(5)：5-8.

[7] 劉曉剛，吳曉平，王寶軍，等.衛(wèi)星重力梯度數(shù)據(jù)的網(wǎng)格化方法研究[J].大地測量與地球動力學，2010，30(6)：60-65.

[8] 劉曉霞，江在森，武艷強，等.Kriging方法在GPS速度場網(wǎng)格化和應變率場計算中的適用性[J].武漢大學學報(信息科學版)，2014，39(4)：457-461.

[9] 劉兆平，楊進，武煒，等.地球物理數(shù)據(jù)網(wǎng)格化方法的選取[J].物探與化探，2010，34(1)：93-97.

[10] 姚道榮，鐘波，汪海洪，等.最小二乘配置與普通Kriging法的比較[J].大地測量與地球動力學，2008，28(3)：77-82.

[11] 熊俊楠，馬洪濱.變異函數(shù)的自動擬合研究[J].測繪信息與工程，2008，33(1)：27-29.

[12] 李章林，王平，李冬梅，等.實驗變差函數(shù)計算方法的研究與運用[J].國土資源信息化，2008(2)：10-14.

[13] 郭泉河，李秀海.不同變異函數(shù)的泛Kriging法的GPS高程擬合結果[J].黑龍江工程學院學報(自然科學版)，2011， 25(4)：26-28.

[14] 沈云中，陶本藻.實用測量數(shù)據(jù)處理方法[M].北京：測繪出版社，2012：194-214.

[15] 侯景儒，尹鎮(zhèn)南，李維明，等.實用地質(zhì)統(tǒng)計學[M].北京：地質(zhì)出版社，1998：31-68.

Research on the applicability of Kriging methodin regional gravity field interpolation

ZHANG Yongqi1,2, HAN Meitao1, HAN Xiaofei1, ZHENG Zengji1，ZHAI Hongguang1

(1. Shaanxi Provincial Earthquake Administration, Xi’an 710068,China;2. State Key Laboratory of Earthquake Dynamics，Institute of Geology, China Earthquake Administration，Beijing 100029，China)

Because the gravity observation network distribution is not uniform and the point density is not enough, it is very limited to describe the variation of regional gravity field, which needs to use the reasonable interpolation algorithm. Based on the detailed introduction of the Kriging interpolation algorithm and the variation function theory, this paper unifies the Shaanxi area 2011～2012 discrete variation of gravity data, obtains the fitting parameters of theoretical variation function model by least square method, and uses the experimental variation model for Kriging interpolation algorithm. Considering the cross validation accuracy statistics, the interpolation precision results, the application to study Kriging interpolation algorithm, the result shows that: the fitting parameter accuracy of spherical model and exponential model are high based on the least square method, and the accuracy of Gauss model and logarithmic model is the worst; the experimental variation function based on the exponential model and spherical model used in the Kriging interpolation algorithm, can get the highest accuracy of interpolation, meanwhile the image is smooth and the abnormal area is obvious, so it is a suitable interpolation calculation algorithm for gravity discrete data of Shaanxi area.

Kriging method；variation function model；Shaanxi area；gravity field；interpolation

2016-10-28

2016年度中國地震局“三結合”課題(162704)；陜西省起航與創(chuàng)新基金資助項目(201514)

張永奇(1985-)，男，工程師,博士研究生.

著錄：張永奇,韓美濤,韓曉飛，等.Kriging法在區(qū)域重力場插值中的適用性研究[J].測繪工程，2018，27(2)：1-6.

10.19349/j.cnki.issn1006-7949.2018.02.001

P207

A

1006-7949(2018)02-0001-06

張德福]