基于泛克里格方法的位場擴邊處理

李 盼， 戴前偉，2， 呂宏安

(1.中南大學 地球科學與信息物理學院，長沙 410083；2.中南大學 有色金屬成礦預測與地質環境監測教育部重點實驗室，長沙 410083；3.湖南省永吉高速公路建設開發有限公司，吉首 416000)

0 引言

一些常見的位場數據處理與轉換可以在空間域或者頻率域內進行[1-2](位場導數[3]、位場延拓[4]等)。

由于在頻率域的處理相對于空間域的處理更為簡潔[5]，所以目前的位場數據處理與轉換更多的是在頻率域內進行的。

在頻率域內對位場數據進行快速傅立葉變換時要求數據為規則網格數據，為了減弱邊界效應以及Gibbs效應[6]，要求網格數據在x、y上的點數均為2的冪次方，并且擴邊后的邊界值均等于同一數值[7]。所以當原始的位場網格數據不滿足上述要求時，就必須采用適當的擴邊算法對原始數據進行擴邊處理以滿足頻率域快速傅立葉變換的要求。目前比較常用的擴邊算法是余弦擴邊，該方法利用余弦公式對數據區域進行擴展，擴展后的邊界值均為零，從而滿足了頻率域快速傅立葉變換的要求。

雖然余弦擴邊方法有著簡單易實現的特點，但是作為一種優秀的擴邊算法，應該使得擴展后的數據區域能夠與原始數據區域平滑、連續地銜接，并且擴展后的數據能盡可能地接近真實位場值，很明顯余弦擴邊并不能保證擴邊區域與原始區域的平滑連接，更不能保證擴展區域的數據接近真實的位場值，因此,尋找一種更加合理的擴邊方法是有必要的。擴邊算法的實質是利用空間插值算法估計出未知區域的位場值，前人具體研究過各種插值方法在地球物理數據網格化中的應用[8-10]，其中克里格方法是公認的效果最好的空間插值方法之一，作為普通克里格的改進，泛克里格方法考慮到了非平穩隨機變量的漂移與波動，充分利用了數據點的空間相關性，能夠盡可能準確地對未知點進行估計。其他常見的空間插值方法有反距離加權法、最小曲率法、徑向基函數法、自然鄰點插值法、線性插值三角網法等方法。其中反距離加權方法容易出現“牛眼”現象[11]，不宜用來做擴邊處理；而最小曲率法雖然能很好地保證數據的圓滑，但是在數據稀疏區域插值結果可能會出現嚴重的“震蕩”現象[12]，嚴重影響擴邊效果，需要對其穩定性做進一步改進。筆者設計了理論模型對泛克里格擴邊方法的效果進行了驗證，同時為了說明該方法的效果，還將其與余弦擴邊方法以及其他幾種效果較好的徑向基函數法、自然鄰點插值法、線性插值三角網法的擴邊效果進行了比較，對原始數據進行擴邊，然后利用涉及到頻率域處理的ISVD算法[13]換算位場垂向一階導數，通過比較換算結果與理論值的誤差來說明各個擴邊方法的優劣性。

1 泛克里格插值方法原理[14]

1.1 非平穩隨機變量的漂移與波動

在普通克里金方法中，區域變量Z(x)的數學期望是假設不變的，而在泛克里格方法中隨機變量Z(x)在研究區域內的數學期望E[Z(x)]是變化的，即：

E[Z(x)]=m(x)

(1)

非平穩隨機變量包括兩個部分：①漂移；②波動。其中漂移是指隨機變量在點x處的數學期望，即式(1)中的m(x)，而波動R(x)是指點x處隨機變量與漂移的差。漂移與波動的關系為式(2)。

Z(x)=m(x)+R(x)

(2)

從實際意義上來說，漂移可以理解為變量在研究區域的某種明確的變化趨勢，而波動則是隨機變量在m(x)附近的隨機擺動，波動的數學期望為零，即：

E[R(x)]=E[Z(x)]-E[m(x)]=0

(3)

漂移一般用一次或二次多項式表示，即：

(4)

式中：fl(x)為已知函數；μl為未知系數。

在曲面插值中，漂移通常采用式(3)表示。

m(x,y)=μ0+u1x+u2y+u3xy+u4x2+u5y2

(5)

1.2 泛克里格方程組

區域變量Z(x)在待插值點x0處的真實值Z(x0)是未知的，因此需要用x0影響范圍內的測量值的加權平均來估計待插值點的值：

(6)

而待插值點x0處的真實值Z(x0)與估計值Z*(x0)之差應該滿足無偏條件，即：

E[Z*(x0)-Z(x0)]≡0

(7)

結合式(6)、式(7)及式(1)、式(4)可得：

(8)

如果要使式(8)對任何一組系數都成立，則必須滿足：

(9)

式(9)即為用k+1個方程表示的無偏條件。

文獻[9]推導了用變差函數表示的估計方差表達式：

(10)

采用拉格朗日乘數法，將無偏條件式(9)代入式(10)得目標函數式(11)。

(11)

將式(11)對λi及μl求一階偏導，令其為“0”，

便可以得到泛克里格方程組(12)。

圖1 理論模型的平面圖及其重力場值的灰度圖Fig.1 The model plane position and its gravity anomaly gray image

表1 理論模型的參數

(12)

2 理論模型實驗

筆者設計了圖1所示的理論模型對本文提出的方法效果進行驗證，為了更加明顯地突出邊界處理的效果，特意在數據區域邊界放置了四個長方形異常體2、3、4、5，模型參數如表1所示。

用文獻[15]介紹的式(13)、式(14)計算了模型的重力場及其垂向一階導數理論值，網格密度為101×101。

(13)

(14)

分別用余弦法、泛克里格法、徑向基函數法、自然鄰點插值法、線性插值三角網法對重力場數據進行擴邊，擴邊后的邊界值統一設定為零，然后利用離散的網格數據基于ISVD算法[13]換算出其垂向一階導數。圖2、圖3是基于各種擴邊方法通過頻率域處理換算垂向一階導數與理論值的對比及誤差絕對值對比。表2統計了基于各種擴邊方法換算出的垂向一階導數與理論值的誤差。

圖2 基于各種擴邊方法通過頻率域處理換算垂向一階導數與理論值的對比Fig.2 The compare between transformed vertical first derivative and the theoretical value based on different edge extending methods

圖3 基于各種擴邊方法換算出的垂向一階導數與理論值的誤差絕對值Fig.3 The absolute value of the error of transformed vertical first derivative and the theoretical value based on different edge extending methods

圖2是基于各種擴邊方法對原始數據進行了擴邊之后，利用ISVD算法換算出的位場垂向一階導數，并基于式(14)計算出的理論值作為對比。通過對比可以看出來基于離散值換算出的結果與理論值是有一定誤差的，誤差主要出現于曲線的“峰頂”和“谷底”處，基于各種擴邊方法換算的結果都非常接近，主要區別在于數據邊界位置，即圖2兩個紅色橢圓所圈出的區域，在其他地方，這五種擴邊方法換算結果曲線都比較吻合，而在邊界出則出現了比較明顯的分離，其中青色虛線所代表的余弦擴邊法出現了比較明顯的偏差，其余四種方法比較接近。圖3給出了各種方法的換算結果與理論值的誤差絕對值的曲線圖，由圖3可以看出，在靠中間的區域，所有方法的誤差曲線都非常接近，只存在細微的差距，主要區別在邊界處，余弦法的表現是比較差的，出現了明顯的誤差陡增的現象，這主要是余弦法在擴邊時，對于原始數據與擴邊數據的銜接處處理得不夠光滑，使得連接處出現“溝壑”現象，使得數據嚴重失真，這必然造成換算結果的誤差增大。而其他四種方法相對來說表現比較接近，在對應圖2橢圓所圈出的邊界處只出現了較小幅度的誤差增加，整體來說，沒有出現明顯的誤差陡增，比較平穩。仔細對比的話，泛克里格方法對應的誤差曲線處于最下方，也就是說誤差最小。表2計算了各種擴邊方法換算結果與理論值的最大誤差與均方根誤差，由表2可以看出，泛克里格方法無論是最大誤差還是均方根誤差都是最小的，而余弦法則是表現最差的。

表2 基于各種擴邊方法換算出的垂向一階導數與理論值的誤差Tab.2 The error of transformed vertical first derivative and the theoretical value based on different edge extending methods

3 實際應用

我們選取了Gooper于2008年在SEG網站上提供的南非Witwatersrand Basin某區域的開源重力數據，利用泛克里格擴邊方法對原始數據做了擴邊處理，然后用ISVD算法換算了其一、二、三階垂向導數，如圖3所示。由圖3可以看出，一階垂向導數的灰度圖相對原始重力異常而言能夠觀察到非常明顯的地質體主體框架，右上角的環狀地質體邊界清晰可辨，其左側和下方的地質結構體相對原始重力異常圖有明顯的突顯；而二階和三階垂向導數則突出了更多的地質細節，地質體與周圍圍巖區分明顯，能夠觀測到明顯的褶皺和斷層的分布和走向。

圖4 實測重力異常數據及其一、二、三階垂向導數Fig.4 The measured gravity anomaly data and its first, second and third order vertical derivative

4 總結

鑒于傳統的余弦位場擴邊方法不能保證擴邊后位場數據的連續且平滑的要求，筆者旨在從目前常用的地球物理數據網格化方法中尋找一種合適的擴邊方法，通過參考多篇相關文獻，篩除了效果不佳的方法，從中挑出了效果較好的泛克里格法、徑向基函數法、自然鄰點插值法、線性插值三角網法這幾種方法與余弦擴邊法對位場數據做了擴邊處理，通過換算其垂向一階導數對比了各自的精度。結果表明，泛克里格擴邊方法的精度明顯高于余弦擴邊法，而與其他擴邊方法相比，其精度也是最高的。通過實測數據的應用結果表明，基于泛克里格擴邊的位場數據的高階導數邊值處理是比較成功的。