基于泛克里格方法的位場(chǎng)擴(kuò)邊處理

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

(1.中南大學(xué) 地球科學(xué)與信息物理學(xué)院，長(zhǎng)沙 410083；2.中南大學(xué) 有色金屬成礦預(yù)測(cè)與地質(zhì)環(huán)境監(jiān)測(cè)教育部重點(diǎn)實(shí)驗(yàn)室，長(zhǎng)沙 410083；3.湖南省永吉高速公路建設(shè)開(kāi)發(fā)有限公司，吉首 416000)

0 引言

一些常見(jiàn)的位場(chǎng)數(shù)據(jù)處理與轉(zhuǎn)換可以在空間域或者頻率域內(nèi)進(jìn)行[1-2](位場(chǎng)導(dǎo)數(shù)[3]、位場(chǎng)延拓[4]等)。

由于在頻率域的處理相對(duì)于空間域的處理更為簡(jiǎn)潔[5]，所以目前的位場(chǎng)數(shù)據(jù)處理與轉(zhuǎn)換更多的是在頻率域內(nèi)進(jìn)行的。

在頻率域內(nèi)對(duì)位場(chǎng)數(shù)據(jù)進(jìn)行快速傅立葉變換時(shí)要求數(shù)據(jù)為規(guī)則網(wǎng)格數(shù)據(jù)，為了減弱邊界效應(yīng)以及Gibbs效應(yīng)[6]，要求網(wǎng)格數(shù)據(jù)在x、y上的點(diǎn)數(shù)均為2的冪次方，并且擴(kuò)邊后的邊界值均等于同一數(shù)值[7]。所以當(dāng)原始的位場(chǎng)網(wǎng)格數(shù)據(jù)不滿足上述要求時(shí)，就必須采用適當(dāng)?shù)臄U(kuò)邊算法對(duì)原始數(shù)據(jù)進(jìn)行擴(kuò)邊處理以滿足頻率域快速傅立葉變換的要求。目前比較常用的擴(kuò)邊算法是余弦擴(kuò)邊，該方法利用余弦公式對(duì)數(shù)據(jù)區(qū)域進(jìn)行擴(kuò)展，擴(kuò)展后的邊界值均為零，從而滿足了頻率域快速傅立葉變換的要求。

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

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

1.1 非平穩(wěn)隨機(jī)變量的漂移與波動(dòng)

在普通克里金方法中，區(qū)域變量Z(x)的數(shù)學(xué)期望是假設(shè)不變的，而在泛克里格方法中隨機(jī)變量Z(x)在研究區(qū)域內(nèi)的數(shù)學(xué)期望E[Z(x)]是變化的，即：

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

(1)

非平穩(wěn)隨機(jī)變量包括兩個(gè)部分：①漂移；②波動(dòng)。其中漂移是指隨機(jī)變量在點(diǎn)x處的數(shù)學(xué)期望，即式(1)中的m(x)，而波動(dòng)R(x)是指點(diǎn)x處隨機(jī)變量與漂移的差。漂移與波動(dòng)的關(guān)系為式(2)。

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

(2)

從實(shí)際意義上來(lái)說(shuō)，漂移可以理解為變量在研究區(qū)域的某種明確的變化趨勢(shì)，而波動(dòng)則是隨機(jī)變量在m(x)附近的隨機(jī)擺動(dòng)，波動(dòng)的數(shù)學(xué)期望為零，即：

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

(3)

漂移一般用一次或二次多項(xiàng)式表示，即：

(4)

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

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

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

(5)

1.2 泛克里格方程組

區(qū)域變量Z(x)在待插值點(diǎn)x0處的真實(shí)值Z(x0)是未知的，因此需要用x0影響范圍內(nèi)的測(cè)量值的加權(quán)平均來(lái)估計(jì)待插值點(diǎn)的值：

(6)

而待插值點(diǎn)x0處的真實(shí)值Z(x0)與估計(jì)值Z*(x0)之差應(yīng)該滿足無(wú)偏條件，即：

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

(7)

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

(8)

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

(9)

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

文獻(xiàn)[9]推導(dǎo)了用變差函數(shù)表示的估計(jì)方差表達(dá)式：

(10)

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

(11)

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

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

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

表1 理論模型的參數(shù)

(12)

2 理論模型實(shí)驗(yàn)

筆者設(shè)計(jì)了圖1所示的理論模型對(duì)本文提出的方法效果進(jìn)行驗(yàn)證，為了更加明顯地突出邊界處理的效果，特意在數(shù)據(jù)區(qū)域邊界放置了四個(gè)長(zhǎng)方形異常體2、3、4、5，模型參數(shù)如表1所示。

用文獻(xiàn)[15]介紹的式(13)、式(14)計(jì)算了模型的重力場(chǎng)及其垂向一階導(dǎo)數(shù)理論值，網(wǎng)格密度為101×101。

(13)

(14)

分別用余弦法、泛克里格法、徑向基函數(shù)法、自然鄰點(diǎn)插值法、線性插值三角網(wǎng)法對(duì)重力場(chǎng)數(shù)據(jù)進(jìn)行擴(kuò)邊，擴(kuò)邊后的邊界值統(tǒng)一設(shè)定為零，然后利用離散的網(wǎng)格數(shù)據(jù)基于ISVD算法[13]換算出其垂向一階導(dǎo)數(shù)。圖2、圖3是基于各種擴(kuò)邊方法通過(guò)頻率域處理?yè)Q算垂向一階導(dǎo)數(shù)與理論值的對(duì)比及誤差絕對(duì)值對(duì)比。表2統(tǒng)計(jì)了基于各種擴(kuò)邊方法換算出的垂向一階導(dǎo)數(shù)與理論值的誤差。

圖2 基于各種擴(kuò)邊方法通過(guò)頻率域處理?yè)Q算垂向一階導(dǎo)數(shù)與理論值的對(duì)比Fig.2 The compare between transformed vertical first derivative and the theoretical value based on different edge extending methods

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

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

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

3 實(shí)際應(yīng)用

我們選取了Gooper于2008年在SEG網(wǎng)站上提供的南非Witwatersrand Basin某區(qū)域的開(kāi)源重力數(shù)據(jù)，利用泛克里格擴(kuò)邊方法對(duì)原始數(shù)據(jù)做了擴(kuò)邊處理，然后用ISVD算法換算了其一、二、三階垂向?qū)?shù)，如圖3所示。由圖3可以看出，一階垂向?qū)?shù)的灰度圖相對(duì)原始重力異常而言能夠觀察到非常明顯的地質(zhì)體主體框架，右上角的環(huán)狀地質(zhì)體邊界清晰可辨，其左側(cè)和下方的地質(zhì)結(jié)構(gòu)體相對(duì)原始重力異常圖有明顯的突顯；而二階和三階垂向?qū)?shù)則突出了更多的地質(zhì)細(xì)節(jié)，地質(zhì)體與周?chē)鷩鷰r區(qū)分明顯，能夠觀測(cè)到明顯的褶皺和斷層的分布和走向。

圖4 實(shí)測(cè)重力異常數(shù)據(jù)及其一、二、三階垂向?qū)?shù)Fig.4 The measured gravity anomaly data and its first, second and third order vertical derivative

4 總結(jié)

鑒于傳統(tǒng)的余弦位場(chǎng)擴(kuò)邊方法不能保證擴(kuò)邊后位場(chǎng)數(shù)據(jù)的連續(xù)且平滑的要求，筆者旨在從目前常用的地球物理數(shù)據(jù)網(wǎng)格化方法中尋找一種合適的擴(kuò)邊方法，通過(guò)參考多篇相關(guān)文獻(xiàn)，篩除了效果不佳的方法，從中挑出了效果較好的泛克里格法、徑向基函數(shù)法、自然鄰點(diǎn)插值法、線性插值三角網(wǎng)法這幾種方法與余弦擴(kuò)邊法對(duì)位場(chǎng)數(shù)據(jù)做了擴(kuò)邊處理，通過(guò)換算其垂向一階導(dǎo)數(shù)對(duì)比了各自的精度。結(jié)果表明，泛克里格擴(kuò)邊方法的精度明顯高于余弦擴(kuò)邊法，而與其他擴(kuò)邊方法相比，其精度也是最高的。通過(guò)實(shí)測(cè)數(shù)據(jù)的應(yīng)用結(jié)果表明，基于泛克里格擴(kuò)邊的位場(chǎng)數(shù)據(jù)的高階導(dǎo)數(shù)邊值處理是比較成功的。