大點(diǎn)線距數(shù)據(jù)網(wǎng)格化方法研究
——以達(dá)來(lái)廟地區(qū)高精度磁測(cè)數(shù)據(jù)為例

，,2，

(1.中國(guó)地質(zhì)大學(xué)(武漢)資源學(xué)院，湖北 武漢 430074； 2.中國(guó)地質(zhì)大學(xué)地質(zhì)過(guò)程與礦產(chǎn)資源國(guó)家重點(diǎn)實(shí)驗(yàn)室，湖北 武漢 430074)

0 引 言

大點(diǎn)線距數(shù)據(jù)主要是指測(cè)點(diǎn)距、測(cè)線距較大的地球物理數(shù)據(jù)，這些數(shù)據(jù)往往都是大比例尺數(shù)據(jù)。而這些地球物理數(shù)據(jù)在解釋之前，需要經(jīng)過(guò)數(shù)據(jù)網(wǎng)格化處理。不同網(wǎng)格化方法的使用對(duì)數(shù)據(jù)結(jié)果會(huì)產(chǎn)生影響，即空間內(nèi)插模型不同，插值結(jié)果不同。常用的數(shù)據(jù)網(wǎng)格化方法有線性插值法(Watson,1992)、三角網(wǎng)法(劉天佑,2004)、最小曲率法(Briggs,1974)、反距離加權(quán)法(唐澤圣,1999)、克里格法(王仁鐸等,1988)、張力樣條法(唐澤圣,1999;王穎等,2011)等。針對(duì)這些方法，王穎等(2011)利用航空磁測(cè)數(shù)據(jù)網(wǎng)格化測(cè)試認(rèn)為張力樣條法適用于航測(cè)數(shù)據(jù)；劉兆平等(2010)認(rèn)為數(shù)據(jù)量小于1 000個(gè)的情況下克里格法和最小曲率法最可行；而葛志廣等(2010)在利用高精度磁法數(shù)據(jù)網(wǎng)格化測(cè)試后發(fā)現(xiàn)采用反距離加權(quán)插值法、最近鄰點(diǎn)法、徑向基函數(shù)法進(jìn)行網(wǎng)格化,灰色區(qū)域負(fù)異常呈“塊”狀或孤立異常點(diǎn)分布,效果不理想。在大點(diǎn)、線間距的數(shù)值插值方面，前人利用了高精度磁法數(shù)據(jù)或航測(cè)數(shù)據(jù)，使用諸多方法測(cè)試，而在這些數(shù)據(jù)網(wǎng)格化中，等值線往往會(huì)呈現(xiàn)沿測(cè)線方向的“串珠”結(jié)構(gòu)。就這個(gè)問(wèn)題，目前研究不多，謝汝寬等(2013)在航磁數(shù)據(jù)網(wǎng)格化中提出了雙方向測(cè)線型數(shù)據(jù)的網(wǎng)格化方法，取得了較好的效果。

在前人研究的基礎(chǔ)上，使用Oasis montaj平臺(tái)提供的三角網(wǎng)法、距離倒數(shù)加權(quán)法、雙向線法、克里格法對(duì)達(dá)來(lái)廟地區(qū)高精度數(shù)據(jù)進(jìn)行網(wǎng)格化，結(jié)合當(dāng)?shù)氐刭|(zhì)特征，通過(guò)參數(shù)設(shè)置和方法調(diào)整，交叉驗(yàn)證后總結(jié)出該平臺(tái)下適用于該類數(shù)據(jù)網(wǎng)格化的方法，并解決數(shù)據(jù)沿測(cè)線方向呈現(xiàn)的“串珠”結(jié)構(gòu)的問(wèn)題。

所用數(shù)據(jù)為達(dá)來(lái)廟地區(qū)1∶5萬(wàn)高精度磁測(cè)數(shù)據(jù)，研究區(qū)域出露的地層有上石炭—下二疊統(tǒng)寶力高廟組：為灰-灰褐色安山巖、灰紫色英安質(zhì)-流紋質(zhì)火山熔巖及火山碎屑巖；石炭紀(jì)似斑狀黑云母斑巖，第四系全新統(tǒng)主要為角礫層，巖屑層，砂土層。出露的侵入巖為石炭紀(jì)中粗粒黑云母二長(zhǎng)花崗巖、細(xì)粒二長(zhǎng)花崗巖、中粗粒黑云母花崗巖及中細(xì)粒黑云母花崗巖；侏羅紀(jì)灰紫紅、灰肉紅色中細(xì)粒花崗閃長(zhǎng)巖。脈巖類包括中酸性、基性脈巖均有出露。區(qū)內(nèi)與礦床成因關(guān)系密切的侵入巖主要為斑巖類。

區(qū)域上構(gòu)造主要呈北東、北北東向展布，其中巖體和地層受構(gòu)造控制非常明顯，也多呈北東、北北東向產(chǎn)出。其次較發(fā)育的是北西向構(gòu)造，再次為近東西向、近南北向與北北東向構(gòu)造。斷裂構(gòu)造一般數(shù)十千米，小者數(shù)千米—幾百米，主要形成于古生界—中生界。測(cè)線沿北西向布置。測(cè)點(diǎn)距離40 m，測(cè)線間距100 m。

目前，針對(duì)數(shù)據(jù)插值檢驗(yàn)的方法主要是交叉驗(yàn)證法(高艷芳等,2012)。該方法首先需要預(yù)留若干數(shù)據(jù)樣本點(diǎn)，然后對(duì)這些數(shù)據(jù)點(diǎn)作出預(yù)測(cè)，計(jì)算所有估計(jì)值與實(shí)際測(cè)量值的誤差，從而判斷插值方法的優(yōu)劣。交叉驗(yàn)證即通過(guò)對(duì)已知的N個(gè)點(diǎn)的觀測(cè)數(shù)據(jù)，計(jì)算和分析網(wǎng)格化后每個(gè)觀測(cè)點(diǎn)上數(shù)據(jù)的殘差來(lái)進(jìn)行數(shù)據(jù)網(wǎng)格化質(zhì)量評(píng)估。數(shù)據(jù)殘差即網(wǎng)格化后的評(píng)估值Zm,i與實(shí)際觀測(cè)值Ze,i之差。評(píng)估一個(gè)好的插值模型，一般選用可以反映利用樣本點(diǎn)數(shù)據(jù)的估計(jì)靈敏度和極值效應(yīng)的均方根誤差(RMSE)，其值要最小或較小，從而來(lái)確定插值模型的好壞。

(1)

能總體反映估計(jì)誤差的大小的平均相對(duì)誤差(MRE)：

(2)

可以估算可能的誤差范圍的平均絕對(duì)誤差(MAE)：

(3)

1 三角網(wǎng)法

(4)

利用Oasis montaj軟件對(duì)數(shù)據(jù)進(jìn)行網(wǎng)格化處理。在進(jìn)行三角網(wǎng)網(wǎng)格化中需要設(shè)置不同的重復(fù)次數(shù)(因?yàn)門IN程序在任何(X,Y)位置能處理的數(shù)據(jù)點(diǎn)都不超過(guò)1個(gè),所以需要對(duì)重復(fù)Z值進(jìn)行選擇，選項(xiàng)是所有值求和或者使用平均值)和插值方法(本次選線性插值)。利用實(shí)驗(yàn)數(shù)據(jù)，在該軟件下進(jìn)行測(cè)試，結(jié)果如圖1。

通過(guò)對(duì)數(shù)據(jù)網(wǎng)格化后并繪制三角網(wǎng)，三角網(wǎng)呈規(guī)則分布。驗(yàn)證結(jié)果如表1。

根據(jù)判斷原則，從表中信息可以得出：(1) 三角網(wǎng)法在很大程度上尊重了原始數(shù)據(jù)，插值誤差都接近0,應(yīng)用效果較好；(2) 三角網(wǎng)圖中顯示了該方法沿測(cè)線方向進(jìn)行插值，沒(méi)有發(fā)現(xiàn)跳三角的現(xiàn)象；(3) 從垂向一階導(dǎo)數(shù)圖發(fā)現(xiàn)，異常呈線性，三角網(wǎng)法網(wǎng)格化后能反映出異常特征。

圖1 三角網(wǎng)法試驗(yàn)結(jié)果

表1 三角網(wǎng)法交叉驗(yàn)證結(jié)果

2 距離倒數(shù)加權(quán)法

距離倒數(shù)插值方法是加權(quán)移動(dòng)平均方法的一種(張瑞新等,1991)。假設(shè)未知點(diǎn)x0處屬性值是在局部鄰域內(nèi)所有數(shù)據(jù)點(diǎn)的距離加權(quán)平均值，加權(quán)移動(dòng)平均方法的計(jì)算如式(5)：

(5)

(6)

式中，權(quán)重系數(shù)由函數(shù)φ(d(x,xi))計(jì)算，要求當(dāng)d→0時(shí)φ(d)→1，一般取倒數(shù)或負(fù)指數(shù)形式；d-r、e-d、e-d2,其中,φ(d(x,xi))最常見(jiàn)的形式是距離倒數(shù)加權(quán)函數(shù)，形式如式(7)：

(7)

式(7)中，xj為未知點(diǎn)，xi為已知數(shù)據(jù)點(diǎn)。

在Oasis montaj軟件中進(jìn)行距離倒數(shù)加權(quán)法需要進(jìn)行下列參數(shù)設(shè)置：加權(quán)冪、加權(quán)坡度、搜索半徑。

需要說(shuō)明的是：距離倒數(shù)加權(quán)函數(shù)根據(jù)離開(kāi)定義搜索半徑的距離來(lái)指定平均權(quán)重：1/(距離冪+1/坡度)，其中距離的單位是網(wǎng)格圖單元的單位。在該模塊里，系統(tǒng)所使用的冪的默認(rèn)值是2，坡度的默認(rèn)值是1，這樣將生成一個(gè)鐘形權(quán)重函數(shù)。這里權(quán)重在0距離是有限的，所以這里需要坡度>0。坡度降低會(huì)使函數(shù)變平，使原網(wǎng)格圖單元的點(diǎn)權(quán)重增加，圖像變得平滑。

如此便可以通過(guò)對(duì)這些參數(shù)的調(diào)整來(lái)對(duì)數(shù)據(jù)進(jìn)行更好的網(wǎng)格化處理。試驗(yàn)結(jié)果如圖2。

圖2 距離倒數(shù)加權(quán)法試驗(yàn)結(jié)果

根據(jù)上面不同參數(shù)的設(shè)置可以看出：在冪指數(shù)為2、坡度為1時(shí)，搜索半徑不同，產(chǎn)生的效果不同；同樣，當(dāng)冪指數(shù)不同時(shí)，產(chǎn)生的效果也不同。通過(guò)設(shè)置參數(shù)，冪指數(shù)為2、坡度為1時(shí)，當(dāng)搜索半徑為數(shù)據(jù)網(wǎng)格單元的4倍、5倍、6倍時(shí)數(shù)據(jù)效果相近；同樣，在冪指數(shù)為2、坡度為0.5時(shí)，數(shù)據(jù)效果相近。對(duì)其效果進(jìn)行驗(yàn)證(表2)。

表2 距離倒數(shù)加權(quán)交叉驗(yàn)證結(jié)果

通過(guò)判斷原則，分析表2中的信息得出：(1) 冪指數(shù)為2、坡度系數(shù)為1時(shí)，搜索半徑小的情況下，數(shù)據(jù)網(wǎng)格化較好；冪指數(shù)為2，坡度系數(shù)為0.5時(shí)，情況相同。在該方法中搜索半徑影響了數(shù)據(jù)網(wǎng)格化質(zhì)量；(2) 在高精度磁測(cè)數(shù)據(jù)中，可取網(wǎng)格單元的4～5倍來(lái)控制數(shù)據(jù)質(zhì)量。而冪指數(shù)和坡度是控制數(shù)據(jù)網(wǎng)格化效果是否圓滑的參數(shù)。

選冪指數(shù)為2、坡度系數(shù)為1、搜索半徑為88時(shí)，作垂向一階導(dǎo)數(shù)圖(圖3)。

圖3 距離倒數(shù)加權(quán)結(jié)果

從圖3中發(fā)現(xiàn)：使用距離導(dǎo)數(shù)加權(quán)法雖然能在網(wǎng)格化時(shí)突出一些北東向特征，但在作垂向一階導(dǎo)數(shù)后，圖像沿測(cè)線方向呈帶狀分布，應(yīng)用效果并不理想。

3 雙向線法

雙向線網(wǎng)格化又叫做測(cè)線網(wǎng)格化方法(謝汝寬等,2013)，是Oasis montaj軟件中常用的數(shù)據(jù)網(wǎng)格化方法之一，其主要用于對(duì)平行測(cè)線或者近似平行測(cè)線的數(shù)據(jù)網(wǎng)格化。雙向線法中所指的雙方向是對(duì)測(cè)線數(shù)據(jù)先進(jìn)行某一方向進(jìn)行插值，再沿垂直于前一插值方向進(jìn)行插值。在Oasis montaj中，雙向線網(wǎng)格化需要2個(gè)步驟來(lái)完成：(1) 每條線沿原始測(cè)線插值，并使用觀測(cè)數(shù)據(jù)在每個(gè)需要的網(wǎng)格線與測(cè)線交叉點(diǎn)處產(chǎn)生數(shù)據(jù)值；(2) 每條測(cè)線的交叉點(diǎn)再沿測(cè)線方向插值，產(chǎn)生網(wǎng)格點(diǎn)上的值。在這個(gè)過(guò)程中，該軟件中提供的插值方法有線性、最小曲率、Akima、最鄰近法。此外，在二次插值前還可以通過(guò)網(wǎng)格圖的正確方向來(lái)強(qiáng)調(diào)數(shù)據(jù)中的地質(zhì)趨勢(shì)，即調(diào)整網(wǎng)格圖方向。

Akima插值方法是在2個(gè)實(shí)測(cè)點(diǎn)之間進(jìn)行插值，但同時(shí)還需要與這2個(gè)點(diǎn)相近的4個(gè)實(shí)測(cè)點(diǎn)上的觀測(cè)值(Akima,1970)。設(shè)置過(guò)實(shí)測(cè)點(diǎn)Pi(xi)和Pi+1(xi+1)的三次多項(xiàng)式為：

y=p0+p1(x-xi)+p2(x-xi)2+p3(x-xi)3

(8)

式(8)中，p0=y1；p1=t1；p2=[3(y2-y1)/(x2-x1)-2t1-t2]/(x2-x1)；p3=[t1+t2-2(y2-y1)/(x2-x1)]/(x2-x1)2；y1、y2分別為點(diǎn)Pi和Pi+1上的函數(shù)值；t1、t2分別為點(diǎn)Pi和Pi+1上的斜率，而每個(gè)點(diǎn)的斜率t由附近5個(gè)點(diǎn)確定。

根據(jù)研究區(qū)的數(shù)據(jù)情況，數(shù)據(jù)測(cè)線呈北西向，地質(zhì)構(gòu)造為北東向展布。因此，在網(wǎng)格化時(shí)對(duì)數(shù)據(jù)參數(shù)需要進(jìn)行調(diào)整，考慮設(shè)置數(shù)據(jù)方向?yàn)?5°。使用Akima和線性插值對(duì)數(shù)據(jù)進(jìn)行網(wǎng)格化，并作垂向一階導(dǎo)數(shù)圖(圖4)。

通過(guò)試驗(yàn)發(fā)現(xiàn)，該方法適用于測(cè)線型數(shù)據(jù)，且測(cè)線數(shù)據(jù)網(wǎng)格化后效果較好，突出了測(cè)線垂直方向上的異常情況。分別選取雙向Akima插值和雙向線性插值法進(jìn)行驗(yàn)證(表3)。

表3 雙向線法交叉驗(yàn)證結(jié)果

根據(jù)判斷原則，結(jié)合表中信息可知：(1) 2種插值方法的使用，對(duì)數(shù)據(jù)結(jié)果影響不大，但雙向Akima方法較優(yōu)于線性法；(2) 該方法的使用受測(cè)線數(shù)據(jù)方向影響，根據(jù)原理，該方法更適用于垂向測(cè)線數(shù)據(jù)插值，所以在使用中要注意測(cè)線方向的調(diào)整；(3) 從垂向?qū)?shù)圖發(fā)現(xiàn)，2種方法效果相當(dāng)，都能較好地反映北東向異常的特征。

圖4 雙向線法試驗(yàn)結(jié)果

4 克里格法

(9)

式(9)中，h為滯后距，N(h)為數(shù)據(jù)在[f(xi),f(xi+h)]的對(duì)數(shù)。

利用克里格法一般是先對(duì)已知數(shù)據(jù)做實(shí)驗(yàn)變差函數(shù)，再用理論變差函數(shù)與實(shí)驗(yàn)變差函數(shù)進(jìn)行擬合。幾種常用理論變差函數(shù)所用的模型如下。

(1) 球狀：

(變程a′=a) (10)

(2) 指數(shù)模型：

(變程a′=3a) (11)

(3) 高斯：

在Oasis montaj軟件中,使用克里格法網(wǎng)格化同樣是通過(guò)調(diào)整實(shí)驗(yàn)變差函數(shù)，再使用理論變差函數(shù)模型進(jìn)行套合。利用普通克里格法將磁測(cè)數(shù)據(jù)進(jìn)行網(wǎng)格化，數(shù)據(jù)測(cè)線沿北西向展布，進(jìn)行實(shí)驗(yàn)變差函數(shù)結(jié)構(gòu)分析發(fā)現(xiàn)，該結(jié)構(gòu)為球狀模型，因此再利用球狀模型理論變差函數(shù)進(jìn)行套合?？紤]到數(shù)據(jù)各向異性的特征，構(gòu)造走向?yàn)楸睎|向，所以設(shè)置該方向上的權(quán)重，優(yōu)先對(duì)北東向進(jìn)行插值網(wǎng)格化。參數(shù)如下：塊金值為0，拱高為130，變程為750，方向?yàn)?5，方向權(quán)重為4 (各向異性)。

圖5 變差函數(shù)圖

經(jīng)過(guò)變差函數(shù)調(diào)整和其他參數(shù)設(shè)置并對(duì)結(jié)果作垂向一階導(dǎo)數(shù)(圖6)。

圖6 克里格法試驗(yàn)結(jié)果

利用克里格法發(fā)現(xiàn)，通過(guò)調(diào)整變差函數(shù)和方向，消除了沿測(cè)線方向上的部分“串珠”結(jié)構(gòu)，而從導(dǎo)數(shù)圖像也能看出其北東向異常特征(表4)。

表4 克里格法交叉驗(yàn)證結(jié)果

經(jīng)過(guò)驗(yàn)證，克里格法誤差也比較小，能突出方向上的特征。

5 結(jié) 論

利用內(nèi)蒙古達(dá)來(lái)廟地區(qū)高精度磁測(cè)數(shù)據(jù)進(jìn)行網(wǎng)格化測(cè)試，通過(guò)利用三角網(wǎng)法、距離倒數(shù)加權(quán)法、雙向線法、克里格法對(duì)高精度磁測(cè)數(shù)據(jù)的網(wǎng)格效果分析得出以下幾點(diǎn)結(jié)論。

(1) 三角網(wǎng)法最能反映原始數(shù)據(jù)情況，在高精度磁法數(shù)據(jù)使用線性三角網(wǎng)法能快速計(jì)算，得到網(wǎng)格化結(jié)構(gòu)。但會(huì)產(chǎn)生一些沿測(cè)線方向的“串珠”結(jié)構(gòu)，所以采用線性的方法插值效果較好，尤其是對(duì)于大點(diǎn)線類型的數(shù)據(jù)，線性三角網(wǎng)法效果較好。

(2) 距離倒數(shù)加權(quán)法受到冪指數(shù)和搜索半徑的影響較大。大點(diǎn)線距數(shù)據(jù)由于數(shù)據(jù)本身測(cè)點(diǎn)測(cè)線距離比較大，所以要調(diào)整合適的冪指數(shù)和搜索半徑，距離越大，搜索半徑也要隨之提高，這樣才能更好地避免數(shù)據(jù)呈“串珠”狀。經(jīng)過(guò)測(cè)試后發(fā)現(xiàn)，根據(jù)達(dá)來(lái)廟地區(qū)地質(zhì)構(gòu)造概況，該區(qū)所使用的高精度磁測(cè)數(shù)據(jù)，如果采用距離倒數(shù)加權(quán)法，應(yīng)該使用網(wǎng)格圖單元大小4倍的搜索半徑，從而取得相對(duì)合理的結(jié)果。對(duì)于大點(diǎn)線數(shù)據(jù)網(wǎng)格化，該方法并不作為最佳方法來(lái)選擇。

(3) 雙向線法是一種基于測(cè)線的網(wǎng)格化方法，在本次研究高精度數(shù)據(jù)網(wǎng)格化中，因數(shù)據(jù)是測(cè)線數(shù)據(jù)，選用該方法能得出較好的結(jié)果。而要利用雙向線法必須選擇合適的雙向插值方法，在大點(diǎn)線距數(shù)據(jù)中，選用線性插值、Akima插值效果較好。此外，測(cè)線方向的控制對(duì)該方法應(yīng)用效果有很大影響。由于構(gòu)造方向呈北東向，所以在方法使用時(shí)應(yīng)考慮測(cè)線方向性問(wèn)題。對(duì)于大點(diǎn)線數(shù)據(jù)，如果能控制好測(cè)線方向和選擇合適的插值方法，該方法較為適用。

(4) 克里格法的使用必須考慮到變差函數(shù)。在本次研究高精度磁測(cè)數(shù)據(jù)中，通過(guò)調(diào)整合適的變差函數(shù)來(lái)得到較好的插值模型，能較好地使用克里格法。在大點(diǎn)線距數(shù)據(jù)網(wǎng)格化的使用中，套合變差函數(shù)模型必須考慮方向性，對(duì)于各向異性的數(shù)據(jù)，使用方向權(quán)重能較好地體現(xiàn)出構(gòu)造方向上的異常特征。

(5) 大點(diǎn)線數(shù)據(jù)的網(wǎng)格化，要使數(shù)據(jù)最符合原始異常可以選擇線性三角網(wǎng)法，要得到最優(yōu)網(wǎng)格效果可以選擇使用雙向線法，而克里格法可以突出一些方向性，距離倒數(shù)加權(quán)法可以不常用。

郭良輝,孟小紅,郭志宏,等.2005.地球物理不規(guī)則分布數(shù)據(jù)的空間網(wǎng)格化法[J].物探與化探,29(5):438-442.

葛志廣,宋俊杰.2010.高精度磁法數(shù)據(jù)網(wǎng)格化方法的選取[J].工程地球物理學(xué)報(bào),7(2):169-172.

高艷芳,陳實(shí),馮斌.2012.交叉驗(yàn)證在離散數(shù)據(jù)網(wǎng)格化時(shí)的應(yīng)用[J].物探化探計(jì)算技術(shù),34(5):619-621.

劉天佑.2004.應(yīng)用地球物理[M].武漢:中國(guó)地質(zhì)大學(xué)出版社.

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

蘇姝,林愛(ài)文,劉慶華.2004.普通Kriging法在空間內(nèi)插中的運(yùn)用[J].江南大學(xué)學(xué)報(bào)：自然科學(xué)版,3(1):18-21.

唐澤圣.1999.三維數(shù)據(jù)場(chǎng)可視化[M].北京:清華大學(xué)出版社.

王仁鐸,胡光道.1988.線性地質(zhì)統(tǒng)計(jì)學(xué)[M].北京:地質(zhì)出版社.

王穎,祝民強(qiáng),喬康寧.2011.航測(cè)數(shù)據(jù)處理中的空間插值方法比較[J].測(cè)繪與空間地理信息，34(1):234-237,241.

謝汝寬,王平,郭華,等.2013.考慮航磁水平梯度變化的ΔT網(wǎng)格化方法研究[J].地球物理學(xué)報(bào),56(2):660-666.

張瑞新,于汝綬.1991.距離加權(quán)法的完善及其在層狀礦床儲(chǔ)量計(jì)算中的應(yīng)用前景[J].中國(guó)礦業(yè)大學(xué)學(xué)報(bào),20(4):45-53.

AKIMA H.1970.A new method of interpolation and smooth curve fitting based on local procedures[J].Journal of the ACM,17(4):589-602.

BRIGGS C I.1974.Machine contouring using minimum curvature[J].Geophysics,39(1):39-48.

WATSON F D.1992.Contouring:A Guide to the Analysis and Display of Spatial Data[M].Oxford,UK：Pergamon Press.