巖石圈磁場數據網格化方法研究

張忠龍 倪 喆 趙育飛

1）中國貴州 562400 黔西南民族職業技術學院

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

0 引言

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

1 數學模型基本原理

1.1 最小曲率

最小曲率法用公式表示為

1.2 克里格

克里格模型根據所給數據趨勢，構造一個只與點間空間距離相關的半方差函數，從而求取各擬合點的權系數，據此進行數值內插（刁鑫鵬等，2012）。

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

式中，λi為待定克里格權系數，為滿足無偏、方差最小的條件，需

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

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

1.3 改進Shepard

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

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

擬合數學模型為

1.4 反距離加權

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

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

1.5 徑向基函數

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

（1）倒轉復二次函數，公式如下

（2）復對數，公式如下

（3）復二次函數，公式如下

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

（5）薄板樣條法函數，公式如下

2 數學模型精度評定與實驗

2.1 精度評定

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

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

均方根預測誤差的計算式為

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

插值數據殘差均方根計算式為

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

2.2 實驗數據

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

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

表1 小江斷裂總強度加密區數據信息Table 1 Data information on the total intensity of the Xiaojiang fault

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

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

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

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

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

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

圖2 2018 年3 月小江斷裂巖石圈磁場數據克里格插值和反距離加權插值結果(a)克里格插值；(b)反距離加權插植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 結論

通過比較小江斷裂地磁總強度加密區巖石圈磁場數據的格網化結果可以發現，克里格與反距離加權插值的數據網格化精度優于最小曲率、徑向基函數和改進Shepard 等插值方法；克里格插值方法的網格化數據精度略優于反距離加權插值方法。進一步從網格化圖形質量、數據點位空間相關性等結果來比較克里格插值與反距離加權插值法的優劣。結果表明克里格插值兼顧數據平滑性和各實測點與待估點之間的空間位置關系，并較好利用已有實測數據空間分布的結構特征，避免了系統誤差，巖石圈磁場變化的分布和大小特征明顯，得出克里格插值是巖石圈磁場數據網格化的最優數學模型的結論。