巖石圈磁場數(shù)據(jù)網(wǎng)格化方法研究

張忠龍 倪 喆 趙育飛

1）中國貴州 562400 黔西南民族職業(yè)技術(shù)學(xué)院

2）中國昆明 650224 云南省地震局

0 引言

地震發(fā)生受到構(gòu)造運(yùn)動控制（中國科學(xué)院地質(zhì)研究所，1997），構(gòu)造運(yùn)動伴隨著巖石圈磁場的變化，巖石圈磁場及其變化是地球深部物理過程的重要信息來源之一（丁鑒海等，1994）。地磁總強(qiáng)度數(shù)據(jù)通過日變通化、長期變改正處理后剝離得到巖石圈磁場，巖石圈磁場位場分離后得到地殼深部、淺表和居里面的能量加卸載信息，直接反映了地球內(nèi)部各深度的物理過程和構(gòu)造活動能量，比較真實可靠。但獲取的地磁總強(qiáng)度數(shù)據(jù)是離散分布的點位，而實際采集的數(shù)據(jù)點位有限，難以直觀反映巖石圈磁場變化，進(jìn)而不利于判定巖石圈磁場異常區(qū)，需要對巖石圈磁場數(shù)據(jù)進(jìn)行網(wǎng)格化處理。網(wǎng)格化數(shù)學(xué)模型的好壞不但影響到網(wǎng)格化巖石圈磁場數(shù)據(jù)的精度和可信程度，而且會進(jìn)一步影響解釋處理圖件的質(zhì)量、效果和可靠性。常用網(wǎng)格化數(shù)學(xué)模型有最小曲率、克里格、改進(jìn)Shepard、反距離加權(quán)和徑向基函數(shù)等（李振海等，2010）。文中將以上5 種數(shù)學(xué)模型應(yīng)用于小江斷裂帶巖石圈磁場數(shù)據(jù)網(wǎng)格化，并采用實驗方法對其精度進(jìn)行比較分析。

1 數(shù)學(xué)模型基本原理

1.1 最小曲率

最小曲率法用公式表示為

1.2 克里格

克里格模型根據(jù)所給數(shù)據(jù)趨勢，構(gòu)造一個只與點間空間距離相關(guān)的半方差函數(shù)，從而求取各擬合點的權(quán)系數(shù)，據(jù)此進(jìn)行數(shù)值內(nèi)插（刁鑫鵬等，2012）。

對于區(qū)域化變量Z(xi,yi)，假設(shè)在n個位置取樣，有Z(x1,y1)，Z(x2,y2)，…，Z(xn,yn)，則點(xp,yp)處的估計量為

式中，λi為待定克里格權(quán)系數(shù)，為滿足無偏、方差最小的條件，需

式中，μ為拉格朗日乘子，γ(xi,xj)為變異函數(shù)。求出λi便可求得未采樣點值。

理論變異函數(shù)模型較多，常用線性模型、指數(shù)模型、高斯模型和球面模型等（Pan，1997）。

1.3 改進(jìn)Shepard

Shepard 是一種改進(jìn)的加權(quán)平均法，采用1964 年提出的局部逼近模型，以R為半徑的擬合區(qū)分為2 個環(huán)帶，并分別定義權(quán)函數(shù)（管澤霖等，1994）。

式中，ri為已知點與待定點之間的平面距離。

擬合數(shù)學(xué)模型為

1.4 反距離加權(quán)

反距離加權(quán)的基本思想是將權(quán)函數(shù)看作距離倒數(shù)的二次方（鄭曉月等，2008），則

式中，Z(x,y)表示坐標(biāo)點(x,y)上的插值結(jié)果；Zi為第i個點的已知值；di為第i點與插值點之間的距離。

1.5 徑向基函數(shù)

徑向基函數(shù)是多個數(shù)據(jù)網(wǎng)格化方法的組合，其基函數(shù)由單個變量的函數(shù)構(gòu)成（陳歡歡等，2007），一個點(x,y)的基函數(shù)形式往往是hi(x,y)=h(di)，其中di表示由點(x,y)到第i個數(shù)據(jù)點的距離。基函數(shù)類型有如下幾種。

（1）倒轉(zhuǎn)復(fù)二次函數(shù)，公式如下

（2）復(fù)對數(shù)，公式如下

（3）復(fù)二次函數(shù)，公式如下

（4）自然三次樣條函數(shù)，公式如下

（5）薄板樣條法函數(shù)，公式如下

2 數(shù)學(xué)模型精度評定與實驗

2.1 精度評定

數(shù)學(xué)模型精度的檢驗方法主要有交叉驗證和驗證方法（楊元德等，2013），即預(yù)留一個或多個數(shù)據(jù)樣本點，對預(yù)留數(shù)據(jù)點進(jìn)行預(yù)測，計算所有估計值與實際值的誤差，以此來判斷估值方法的優(yōu)劣。交叉驗證的原理是：使用全部數(shù)據(jù)評價自相關(guān)模型，逐一刪除每個數(shù)據(jù)點，并預(yù)測被刪除點的值。驗證方法的原理是：刪除部分?jǐn)?shù)據(jù)（稱之為驗證數(shù)據(jù)），使用剩余數(shù)據(jù)研究趨勢及預(yù)測的自相關(guān)模型。

綜上，選取采樣點計算的均方根預(yù)測誤差（RMSPE，Root Mean Square Predictive Error）來評定各種插值方法的精度，其值越小說明插值結(jié)果越接近真實值，結(jié)果越可靠；使用插值數(shù)據(jù)殘差均方根（RMS）來評價網(wǎng)格數(shù)據(jù)與插值數(shù)據(jù)的一致性。

均方根預(yù)測誤差的計算式為

式中，k表示用于估計網(wǎng)格數(shù)據(jù)精度的樣本數(shù)，zj表示第j（1 ≤j≤k）個樣點插值后的值，表示第j個樣點實測值。

插值數(shù)據(jù)殘差均方根計算式為

式中，n表示全部用于插值的樣本數(shù)據(jù)點數(shù)，zi表示第i個點插值后的值，表示第i點的實測值。

2.2 實驗數(shù)據(jù)

小江斷裂地磁總強(qiáng)度加密區(qū)緯度跨度約23°—28°，經(jīng)度跨度約101°—104°，測點間距約25 km，共布設(shè)了216 個測點。選取該測區(qū)2017 年和2018 年各2 期巖石圈磁場數(shù)據(jù)作為本文的實驗數(shù)據(jù)。測點分布和實驗數(shù)據(jù)信息見圖1 和表1。

圖1 小江斷裂總強(qiáng)度加密區(qū)測點分布Fig.1 Distribution of measuring points in the total intensity of Xiaojiang fault

表1 小江斷裂總強(qiáng)度加密區(qū)數(shù)據(jù)信息Table 1 Data information on the total intensity of the Xiaojiang fault

對表1 中每一期地磁總強(qiáng)度數(shù)據(jù)進(jìn)行預(yù)處理，剔除變遷測點和孤立異常點數(shù)據(jù)。對預(yù)處理后的總強(qiáng)度數(shù)據(jù)進(jìn)行日變通化改正和長期變改正，用長期變改正得到的磁場數(shù)據(jù)減去國際地磁參考場模型數(shù)據(jù)得到巖石圈磁場數(shù)據(jù)。由于每期數(shù)據(jù)獨立測量，實驗數(shù)據(jù)不具有相關(guān)性。選用均方根預(yù)測誤差（RMSPE）和插值數(shù)據(jù)殘差均方根（RMS）等評定指標(biāo)，對5 種數(shù)學(xué)模型的網(wǎng)格化數(shù)據(jù)結(jié)果進(jìn)行精度評定，評定結(jié)果見表2—表5。結(jié)果表明：①克里格與反距離加權(quán)插值方法的數(shù)據(jù)網(wǎng)格化精度優(yōu)于最小曲率、徑向基函數(shù)和改進(jìn)Shepard 等3 種插值方法；②克里格與反距離加權(quán)插值方法的數(shù)據(jù)網(wǎng)格化精度相當(dāng)，2 種插值方法均可用于巖石圈磁場數(shù)據(jù)的網(wǎng)格化。

表2 2017 年3 月數(shù)據(jù)不同數(shù)學(xué)模型精度評定結(jié)果（單位：nT）Table 2 Data accuracy of different mathematical model assessment results in March 2017（Unit:nT）

表3 2017 年8 月數(shù)據(jù)不同數(shù)學(xué)模型精度評定結(jié)果（單位：nT）Table 3 Data accuracy of different mathematical model assessment results in August 2017（Unit:nT）

表4 2018 年3 月數(shù)據(jù)不同數(shù)學(xué)模型精度評定結(jié)果（單位：nT）Table 4 Data accuracy of different mathematical model assessment results in March 2018（Unit:nT）

表5 2018 年8 月數(shù)據(jù)不同數(shù)學(xué)模型精度評定結(jié)果（單位：nT）Table 5 Data accuracy of different mathematical model assessment results in August 2018（Unit:nT）

進(jìn)一步從網(wǎng)格化圖形質(zhì)量、數(shù)據(jù)點位空間相關(guān)性等結(jié)果來比較克里格插值與反距離加權(quán)插值法的優(yōu)劣。如圖2 所示，反距離加權(quán)插值主要考慮數(shù)據(jù)的光滑性，而數(shù)據(jù)點位空間相關(guān)性不足，網(wǎng)格化結(jié)果失真；相反，克里格插值網(wǎng)格化圖形平滑，能更好反映數(shù)據(jù)點的空間相關(guān)性，巖石圈磁場變化的分布和大小特征明顯，克里格插值優(yōu)于反距離加權(quán)插值法。

圖2 2018 年3 月小江斷裂巖石圈磁場數(shù)據(jù)克里格插值和反距離加權(quán)插值結(jié)果(a)克里格插值；(b)反距離加權(quán)插植Fig.2 The meshing results of the Kriging model and inverse distance weighted method for the magnetic field data of the Xiaojiang fault in March 2018

3 結(jié)論

通過比較小江斷裂地磁總強(qiáng)度加密區(qū)巖石圈磁場數(shù)據(jù)的格網(wǎng)化結(jié)果可以發(fā)現(xiàn)，克里格與反距離加權(quán)插值的數(shù)據(jù)網(wǎng)格化精度優(yōu)于最小曲率、徑向基函數(shù)和改進(jìn)Shepard 等插值方法；克里格插值方法的網(wǎng)格化數(shù)據(jù)精度略優(yōu)于反距離加權(quán)插值方法。進(jìn)一步從網(wǎng)格化圖形質(zhì)量、數(shù)據(jù)點位空間相關(guān)性等結(jié)果來比較克里格插值與反距離加權(quán)插值法的優(yōu)劣。結(jié)果表明克里格插值兼顧數(shù)據(jù)平滑性和各實測點與待估點之間的空間位置關(guān)系，并較好利用已有實測數(shù)據(jù)空間分布的結(jié)構(gòu)特征，避免了系統(tǒng)誤差，巖石圈磁場變化的分布和大小特征明顯，得出克里格插值是巖石圈磁場數(shù)據(jù)網(wǎng)格化的最優(yōu)數(shù)學(xué)模型的結(jié)論。