劉亞靜，李勝男，王誠(chéng)聰

(1.華北理工大學(xué) 礦業(yè)工程學(xué)院，河北 唐山 063210;2.中國(guó)科學(xué)院 東北地理與農(nóng)業(yè)生態(tài)研究所，吉林 長(zhǎng)春 130000)

0 引 言

局部多項(xiàng)式插值法是一種局部加權(quán)最小二乘擬合，根據(jù)有限個(gè)已知采樣點(diǎn)數(shù)據(jù)，采用多個(gè)多項(xiàng)式來(lái)進(jìn)行曲面的擬合[1]，每一個(gè)插值點(diǎn)的預(yù)測(cè)值都對(duì)應(yīng)一個(gè)多項(xiàng)式，利用局部多項(xiàng)式插值法得到的曲面，能夠較好地描述短程變異[2]。在測(cè)量過(guò)程中，由于受到諸多因素的影響，會(huì)使測(cè)量數(shù)據(jù)中含有粗差，使得觀測(cè)值將不服從正態(tài)分布，若直接使用這些數(shù)據(jù)進(jìn)行曲面插值，其結(jié)果很難與真實(shí)情況相吻合[3]?？共罟烙?jì)解決了這一問(wèn)題，即使用等價(jià)權(quán)函數(shù)來(lái)代替一般權(quán)函數(shù)，對(duì)含粗差的數(shù)據(jù)進(jìn)行降權(quán)處理，則可以有效提高模型的抗差能力。

在現(xiàn)有研究中，吳富梅等利用水準(zhǔn)網(wǎng)數(shù)據(jù)，對(duì)選權(quán)迭代法的幾種權(quán)函數(shù)在不同顯著性水平下的抗差性進(jìn)行了比較，認(rèn)為當(dāng)數(shù)據(jù)不穩(wěn)定時(shí)，IGGⅢ和丹麥法的抗差效果略優(yōu)于其他權(quán)函數(shù)[4]。徐波、楊勇喜等通過(guò)坐標(biāo)參數(shù)轉(zhuǎn)換實(shí)例對(duì)選權(quán)迭代法進(jìn)行分析，發(fā)現(xiàn)通過(guò)IGG函數(shù)進(jìn)行定權(quán)可有效降低含粗差觀測(cè)值的權(quán)值，在進(jìn)行三維坐標(biāo)轉(zhuǎn)換時(shí)能獲取更高精度的轉(zhuǎn)換參數(shù)[5-6]。孫同賀等將具有穩(wěn)健初值的選權(quán)迭代法融入到DEM粗差探測(cè)中，證實(shí)穩(wěn)健初值的選權(quán)迭代法具有很強(qiáng)的穩(wěn)健性和粗差探測(cè)能力，可以實(shí)現(xiàn)對(duì)粗差的探測(cè)和剔除[7]。HOLLAND等提出選權(quán)迭代法，可通過(guò)反復(fù)迭代，逐步減小粗差觀測(cè)值的權(quán)重，從而進(jìn)行粗差定位和剔除[8]。穩(wěn)健估計(jì)選權(quán)迭代法具有很好的抗差能力[9]，但其相關(guān)研究大多是針對(duì)水準(zhǔn)網(wǎng)測(cè)量和坐標(biāo)轉(zhuǎn)換方面的研究，對(duì)于多項(xiàng)式插值采樣點(diǎn)的定權(quán)鮮有闡述。

一般認(rèn)為DEM誤差來(lái)源于源數(shù)據(jù)的采樣誤差和DEM插值過(guò)程中的插值誤差[10]。在進(jìn)行多項(xiàng)式插值時(shí)，常用的采樣點(diǎn)定權(quán)函數(shù)包括反距離加權(quán)法和局部距離比法，但對(duì)于含粗差的觀測(cè)數(shù)據(jù)，這2種傳統(tǒng)的方法會(huì)將粗差當(dāng)作一般數(shù)據(jù)處理，導(dǎo)致在構(gòu)建曲面模型時(shí)粗差引入，使得曲面發(fā)生歪曲[11-14]。因此文中將穩(wěn)健估計(jì)選權(quán)迭代法與局部多項(xiàng)式插值進(jìn)行結(jié)合，分別采用傳統(tǒng)的反距離加權(quán)法和穩(wěn)健估計(jì)選權(quán)迭代法中的Huber權(quán)函數(shù)法、Andrews權(quán)函數(shù)法、IGG方案對(duì)局部多項(xiàng)式插值的采樣點(diǎn)進(jìn)行定權(quán)，通過(guò)均方根誤差(root mean squared error，RMSE)對(duì)各方法的擬合精度進(jìn)行比較分析，并且對(duì)比分析采用不同定權(quán)方法得到的插值曲面和預(yù)測(cè)高程相似度，探討所引進(jìn)的3種定權(quán)方法中，在數(shù)據(jù)含有粗差時(shí)采用不同定權(quán)方法降低含粗差觀測(cè)值的權(quán)[15-17]，對(duì)于含粗差數(shù)據(jù)插值擬合效果最優(yōu)的方法，以提高其插值精度，對(duì)局部多項(xiàng)式插值的定權(quán)模型改進(jìn)，為局部多項(xiàng)式插值曲面的建立提供保證。

1 局部插值理論模型與精度評(píng)定

1.1 動(dòng)態(tài)球半徑選點(diǎn)法

動(dòng)態(tài)球半徑選點(diǎn)法是以內(nèi)插點(diǎn)為球心，R為半徑，選取各內(nèi)插點(diǎn)周圍可用來(lái)計(jì)算的采樣點(diǎn)[18-19](圖1)。半徑R可通過(guò)構(gòu)造函數(shù)來(lái)確定

圖1 動(dòng)態(tài)球半徑選點(diǎn)法原理示意Fig.1 Schematic diagram of dynamic sphere radius point selection method

(1)

式中N為采樣點(diǎn)總數(shù)，個(gè)；A為研究區(qū)域面積，m2；n為規(guī)定界限球內(nèi)的采樣點(diǎn)數(shù),個(gè)。

1.2 二次曲面模型

高精度曲面的繪制必須保證每一個(gè)數(shù)據(jù)點(diǎn)在權(quán)值其有效范圍內(nèi)平滑過(guò)渡[11]，采用二階可導(dǎo)的二次曲面來(lái)建立局部曲面模型，其二次曲面模型表達(dá)式如下

z=f(x,y)+ε

(2)

f(x,y)=a0+a1x+a2y+a3xy+a4x2+a5y2

(3)

式中z為插值點(diǎn)高程值，m;ε為擬合殘差;f(x,y)為二次曲面點(diǎn)的內(nèi)插值,m，ai(i=0,1,2,…,5)為多項(xiàng)式系數(shù)；x和y分別為平面采樣點(diǎn)的橫縱坐標(biāo)。

通過(guò)最小二乘原理，求解多項(xiàng)式系數(shù)，矩陣表達(dá)形式如下

Z=B×A+ε

(4)

A=(BTPB)-1(BTPε)

(5)

矩陣表示為

1.3 精度評(píng)定

選用均方根誤差(RMSE)進(jìn)行擬合優(yōu)度的評(píng)價(jià)，當(dāng)預(yù)測(cè)值與真實(shí)值相同RMSE=0，誤差越大，該值越大。均方根誤差計(jì)算公式為

(6)

2 采樣點(diǎn)權(quán)函數(shù)

采樣點(diǎn)權(quán)函數(shù)可以用來(lái)反映采樣點(diǎn)與內(nèi)插點(diǎn)的相關(guān)程度。對(duì)區(qū)域內(nèi)的采樣點(diǎn)賦權(quán)，權(quán)重(Pi)的確定與該點(diǎn)到內(nèi)插點(diǎn)的距離(di)有關(guān)，距離越小，該點(diǎn)對(duì)內(nèi)插點(diǎn)的影響就越大，權(quán)重也越大?？刹捎梅淳嚯x加權(quán)法、Huber權(quán)函數(shù)、Andrews權(quán)函數(shù)和IGG方案4種方法確定采樣點(diǎn)的權(quán)重。

多具有不規(guī)則狀晶形，與石英、云母等脈石礦物關(guān)系密切，少部分與黃鐵礦連體，白鎢礦中含鎢約為60.81%，粒度主要為中細(xì)粒，占75.14%，部分為微細(xì)粒，占12.02%，少量為粗粒，約占12.82%。礦石中的白鎢礦解離度較好，-0.074 mm含量占65%的磨礦細(xì)度下，92.73%的白鎢礦已單體解離，6.49%的白鎢礦與脈石礦物共生。因此，白鎢礦主要分為兩類，一類多數(shù)為單體解離，但粒度大小差別較大。另一類粒度較細(xì)，與脈石礦物與石英、長(zhǎng)石等脈石礦物共生或呈港灣狀接觸，或包裹于脈石礦物顆粒中。

2.1 反距離加權(quán)法

反距離加權(quán)法是常用的采樣點(diǎn)定權(quán)法。該方法依據(jù)相近相似的原理[20]，采樣點(diǎn)權(quán)重的大小會(huì)隨著其與插值點(diǎn)之間距離的增大而減小，距離插值點(diǎn)越近，采樣點(diǎn)的權(quán)重就越大[14]。

(7)

(8)

式中di為內(nèi)插點(diǎn)到采樣點(diǎn)的距離，m；X，Y，Z分別為內(nèi)插點(diǎn)的坐標(biāo)；xi,yi,zi分別為采樣點(diǎn)的坐標(biāo)。

2.2 穩(wěn)健估計(jì)權(quán)函數(shù)

采用穩(wěn)健估計(jì)權(quán)函數(shù)進(jìn)行采樣點(diǎn)的定權(quán)，是以內(nèi)插點(diǎn)與采樣點(diǎn)之間的距離為自變量，內(nèi)插點(diǎn)的值為因變量[18]，通過(guò)采樣點(diǎn)權(quán)函數(shù)計(jì)算其權(quán)重。權(quán)函數(shù)分別選用Huber權(quán)函數(shù)、Andrews權(quán)函數(shù)和IGG方案，其中Huber函數(shù)包含正常域和可疑域，只對(duì)可疑域內(nèi)的觀測(cè)值進(jìn)行降權(quán)處理，不能將有害信息剔除掉；Andrews函數(shù)包含可疑域和淘汰域，將所有的觀測(cè)值都進(jìn)行降權(quán)處理；而IGG方案同時(shí)擁有正常域、可疑域和淘汰域，充分利用所有觀測(cè)數(shù)據(jù)[4]。

2.2.1 Huber權(quán)函數(shù)

Huber權(quán)函數(shù)確定的權(quán)因子為

(9)

式中ω為采樣點(diǎn)權(quán)因子；u為采樣點(diǎn)殘差。

2.2.2 Andrews權(quán)函數(shù)

(10)

式中ω為采樣點(diǎn)權(quán)因子；u為采樣點(diǎn)殘差。

2.2.3 IGG方案

IGG方案是由周文江教授最初提出的[21]，該方案將數(shù)據(jù)劃分為正常、可用和有害數(shù)據(jù)，同時(shí)，將權(quán)分區(qū)為保權(quán)區(qū)、降權(quán)區(qū)及淘汰區(qū)，充分考慮測(cè)量數(shù)據(jù)的實(shí)際情況，很好地將不同權(quán)重函數(shù)和容差標(biāo)準(zhǔn)用于不同的測(cè)量數(shù)據(jù)。由IGG方案確定的權(quán)因子[16]為

(11)

式中vi為第i個(gè)測(cè)量的殘差；σ0為標(biāo)準(zhǔn)偏差；k，r為閾值的調(diào)制因子，其中k=1.0～1.5，r=2.5～3.0，該方案k=1.5，r=3.0。

由式(11)可知，IGG方案將權(quán)重確定分為3部分，采用不間斷抗差標(biāo)準(zhǔn)[22]：當(dāng)殘差不大于kσ0時(shí)，采用最小二乘法，其權(quán)重仍為1；當(dāng)殘差介于kσ0與rσ0時(shí)，則采用殘差反比例減少法，降低對(duì)其權(quán)重，以減弱對(duì)參數(shù)識(shí)別的影響；當(dāng)殘差大于rσ0時(shí)，其權(quán)重為0，即淘汰此測(cè)量數(shù)據(jù)，充分體現(xiàn)IGG方案的抗差能力。

3 案例分析

3.1 數(shù)據(jù)來(lái)源

實(shí)驗(yàn)采用《ArcGIS地理信息系統(tǒng)空間分析實(shí)驗(yàn)教程》中的實(shí)驗(yàn)數(shù)據(jù)。該組數(shù)據(jù)共有170個(gè)采樣點(diǎn)，分布均勻，屬性字段分別是橫坐標(biāo)x，縱坐標(biāo)y，以及高程z。部分采樣點(diǎn)數(shù)據(jù)見(jiàn)表1。

表1 部分采樣點(diǎn)數(shù)據(jù)Table 1 Partial sampling points data

為驗(yàn)證Huber函數(shù)法、Andrews函數(shù)法及IGG方案對(duì)含有粗差數(shù)據(jù)的抗差性強(qiáng)度區(qū)別，在實(shí)驗(yàn)數(shù)據(jù)中隨機(jī)加入3組粗差。運(yùn)用MATLAB軟件進(jìn)行曲面插值和擬合，并在三維空間中顯示出來(lái)(圖2)，同時(shí)將穩(wěn)健估計(jì)采樣點(diǎn)定權(quán)3種方法預(yù)測(cè)得到的高程數(shù)據(jù)進(jìn)行對(duì)比分析，通過(guò)散點(diǎn)圖表示(圖3)。

圖2 曲面擬合結(jié)果對(duì)比Fig.2 Comparison of curved surface fitting results

圖3 插值結(jié)果相似度對(duì)比Fig.3 Similarity comparison of interpolation results

3.2 插值曲面對(duì)比分析

對(duì)比分析圖2中的擬合曲面效果，可以看出選權(quán)迭代法對(duì)含有粗差的數(shù)據(jù)進(jìn)行定權(quán)的插值結(jié)果與無(wú)粗差數(shù)據(jù)時(shí)結(jié)果很接近，效果較好，圖中標(biāo)注的部分為插值曲面中差異較明顯的區(qū)域。

3.3 預(yù)測(cè)高程相似度對(duì)比分析

將采用反距離加權(quán)法對(duì)無(wú)粗差數(shù)據(jù)計(jì)算出的高程，與采用穩(wěn)健估計(jì)權(quán)函數(shù)對(duì)于含粗差數(shù)據(jù)計(jì)算出的高程，通過(guò)散點(diǎn)圖(圖3)來(lái)反應(yīng)其結(jié)果的相似程度?？v軸代表數(shù)據(jù)無(wú)粗差時(shí)采用反距離權(quán)重法定權(quán)得到的高程值，橫軸表示采用選權(quán)迭代法計(jì)算出的高程值。

圖3中3幅圖像數(shù)據(jù)點(diǎn)的分布均呈現(xiàn)出y=x的趨勢(shì)，但由IGG方案定權(quán)后進(jìn)行插值擬合的效果更好，說(shuō)明當(dāng)數(shù)據(jù)存在粗差時(shí)，3種選權(quán)迭代法中，采用IGG方案對(duì)采樣點(diǎn)進(jìn)行定權(quán)預(yù)測(cè)的內(nèi)插點(diǎn)高程值，與數(shù)據(jù)無(wú)粗差時(shí)常用的反距離權(quán)重法進(jìn)行定權(quán)計(jì)算得到的高程結(jié)果最為接近。

3.4 擬合優(yōu)度對(duì)比

當(dāng)測(cè)量數(shù)據(jù)含有粗差時(shí)，通過(guò)對(duì)比4種采樣點(diǎn)定權(quán)方法得到的均方根誤差(RMSE)評(píng)價(jià)其曲面擬合優(yōu)度(見(jiàn)表2)，RMSE越大，則說(shuō)明預(yù)測(cè)值與真實(shí)值的誤差越大[23-26]，即抗差性越弱。由表2可知，當(dāng)高程數(shù)據(jù)含有粗差，采用傳統(tǒng)的反距離加權(quán)法進(jìn)行插值時(shí)，RMSE值為1.433 3，相比之下，Huber函數(shù)法和Andrews函數(shù)法插值得到的擬合優(yōu)度更高，分別是0.742 8和0.683 2，表明穩(wěn)健估計(jì)選權(quán)迭代法具有一定的抗差能力。當(dāng)采用IGG方案對(duì)采樣點(diǎn)進(jìn)行定權(quán)時(shí)，RMSE最小為0.531 8，可以看出該方法的抗差能力最強(qiáng)，因此可以推斷當(dāng)觀測(cè)數(shù)據(jù)存在誤差時(shí)，所用的3種方法中，IGG方案對(duì)采樣點(diǎn)定權(quán)所得到的高程值與真實(shí)值最為接近，擬合度最高，插值效果最好。

表2 含粗差時(shí)擬合優(yōu)度對(duì)比Table 2 Analysis of goodness of fit with gross error

4 結(jié) 論

1)在局部多項(xiàng)式插值中，當(dāng)數(shù)據(jù)存在粗差時(shí)，利用傳統(tǒng)的反距離權(quán)重法定權(quán)會(huì)將粗差當(dāng)作一般數(shù)據(jù)進(jìn)行處理，因此在構(gòu)建曲面模型時(shí)會(huì)因粗差的存在而使曲面發(fā)生扭曲，與真實(shí)情況下的插值結(jié)果差別較大。運(yùn)用選權(quán)迭代法將穩(wěn)健估計(jì)思想運(yùn)用到空間確定性插值中，在一定程度上能使插值結(jié)果不受粗差的影響。

2)通過(guò)實(shí)例得出，在對(duì)含粗差數(shù)據(jù)進(jìn)行多項(xiàng)式插值時(shí)，采用穩(wěn)健估計(jì)選權(quán)迭代中的IGG方案對(duì)采樣點(diǎn)定權(quán)，經(jīng)過(guò)多次迭代，可使含有粗差的觀測(cè)值權(quán)降為零或接近零，減弱粗差對(duì)插值結(jié)果的影響，且IGG方案對(duì)采樣點(diǎn)定權(quán)的方法優(yōu)于傳統(tǒng)的反距離權(quán)重法與選權(quán)迭代的Huber函數(shù)法和Andrews函數(shù)法。

3)將抗差估計(jì)的思想融入到局部多項(xiàng)式插值模型中，設(shè)計(jì)出對(duì)粗差具有抵抗能力的插值模型，對(duì)于提高局部多項(xiàng)式插值的準(zhǔn)確性具有很強(qiáng)的現(xiàn)實(shí)意義，但由于IGG方案中k,r的2個(gè)閾值調(diào)制因子取值的不同，會(huì)導(dǎo)致插值結(jié)果有所不同，因此對(duì)該方案插值模型中的參數(shù)取值仍需進(jìn)行優(yōu)化。