顧及方向遮蔽性的反距離權(quán)重插值法

李正泉，吳堯祥

浙江省氣候中心，浙江 杭州310017

1 引 言

空間插值常用于將離散點(diǎn)數(shù)據(jù)轉(zhuǎn)換為連續(xù)曲面，以便更為直觀地考察數(shù)據(jù)要素的空間分布模式。在等值線繪圖中，離散點(diǎn)空間插值是繪制等值線的關(guān)鍵基礎(chǔ)，插值方法選取不僅影響要素值的空間分布形式，而且還影響著等值線繪制的視覺效果［1－5］?，F(xiàn)今多被采用的空間插值方法有反距離權(quán)重法、克里金法、樣條法、三角剖分法、多項(xiàng)式法、趨勢面法等。反距離權(quán)重算法簡便、普適性強(qiáng)，但僅考慮插值點(diǎn)與樣本點(diǎn)的空間距離，沒有包含方向相關(guān)［6－7］；克里金以空間統(tǒng)計學(xué)作為理論基礎(chǔ)，使用變異函數(shù)度量樣點(diǎn)間的空間相關(guān)性，可對插值誤差進(jìn)行逐點(diǎn)理論估計，但其變異函數(shù)需人為經(jīng)驗(yàn)選定，且當(dāng)變異函數(shù)為多個組合時，其計算量劇增［8－9］；樣條法更為適合高密度的樣點(diǎn)內(nèi)插，而對于稀疏有限的樣點(diǎn)該方法多不適合［10－11］；多項(xiàng)式和趨勢面的空間插值更多依賴于插值要素固有的空間分布趨勢，若樣點(diǎn)布局遵從插值要素的空間變化，則可得到理想的插值效果，反之則難以保證插值的合理性［12－13］；基于三角剖分的插值方法可獲得高精度的插值曲面，在地表特征表達(dá)和地形等高線的繪制中經(jīng)常使用［14－16］，但該方法也存在類似于克里金的插值缺陷，即算法復(fù)雜計算量大，且其插值曲面的等值線棱角過于明顯。簡而言之，各種空間插值方式均有優(yōu)缺，并不存在普適最優(yōu)的插值方法。

在眾多空間插值方法中，反距離權(quán)重因其插值原理易于理解、算法簡單便于實(shí)現(xiàn)，且插值后能夠保留原樣點(diǎn)真值，常被作為離散點(diǎn)空間分析的傳統(tǒng)方法之一。但反距離權(quán)重未考慮插值點(diǎn)與樣本點(diǎn)之間的方向關(guān)系，當(dāng)樣點(diǎn)多集中于插值點(diǎn)某個方向時，該方向的樣點(diǎn)總權(quán)重就會大大增加，而其他方向上的樣點(diǎn)因權(quán)重過小而常被忽視，空間插值缺少方向均衡性，其插值合理性亦受到質(zhì)疑［17－18］。另一方面，反距離權(quán)重的等值線繪圖通常有一弊病，即圖中常會有孤立圓區(qū)存在，且易在樣本極值點(diǎn)周圍產(chǎn)生同心圓區(qū)（“牛眼”現(xiàn)象）［19－22］，這嚴(yán)重影響了等值線圖的審美效果。在大多情況下，為追求繪圖美觀不得不摒棄反距離權(quán)重而選用其他插值方式，不惜以數(shù)據(jù)失真來換取整幅等值線圖的協(xié)調(diào)美觀。本文旨在對反距離權(quán)重插值算法進(jìn)行優(yōu)化改進(jìn)，減少或消除反距離權(quán)重等值線繪圖中不合理的孤立圓和同心圓現(xiàn)象的存在。通過增加一個可反映插值樣點(diǎn)方位的調(diào)和權(quán)重系數(shù)K，克服反距離權(quán)重插值中僅考慮插值點(diǎn)與樣本點(diǎn)之間距離而未考慮它們之間方位的插值缺陷，由此提出了調(diào)和反距離權(quán)重（adjusted inverse distance weighted，AIDW）空間插值方式。

2 調(diào)和反距離權(quán)重插值

2.1 反距離權(quán)重算法

反距離權(quán)重（inverse distance weighted，IDW）算法于1968年由Shepard提出，1985年Watson等將其應(yīng)用于空間插值的等值線繪制［23－24］，繼而IDW算法被廣泛應(yīng)用于各行業(yè)領(lǐng)域的空間分析與制圖［25］。IDW算法公式通常被表示為

式中，Z是插值點(diǎn)估計值；Zi是第i個樣本點(diǎn)觀測值；di是插值點(diǎn)與第i個樣本點(diǎn)之間的歐氏距離；n是用于估算插值點(diǎn)值的樣點(diǎn)數(shù)；p是冪指數(shù)。在空間分析和等值線繪圖軟件（如ArcGIS、Surfer等）的計算程序中，p值通常被默認(rèn)為2，此時IDW被稱為距離平方反比法。文獻(xiàn)［26］表明隨著p值的增大，IDW的插值結(jié)果越具有平滑效果。

從IDW算法（式（1））可看出，插值點(diǎn)的估算值僅與插值樣本距插值點(diǎn)的遠(yuǎn)近相關(guān)，并未考慮樣本點(diǎn)的方位性（即樣本點(diǎn)被表示為各向同性）。IDW插值的基本假設(shè)是樣點(diǎn)在插值區(qū)呈均勻分布，當(dāng)樣點(diǎn)在各方向較均勻分布時，該插值算法十分可靠。然而在眾多情況下，樣點(diǎn)在各向分布并非均勻，甚至?xí)霈F(xiàn)樣點(diǎn)集中于某一方向的現(xiàn)象，違背了基本假設(shè)，其插值合理性就難被保證。針對IDW這一插值局限，本文提出了調(diào)和反距離權(quán)重（AIDW）插值算法。

2.2 調(diào)和反距離權(quán)重算法

2.2.1 基本假設(shè)

AIDW算法與IDW基本類似，關(guān)鍵在于增加了可反映插值點(diǎn)與樣本點(diǎn)方位關(guān)系的調(diào)和權(quán)重系數(shù)K，其基本假設(shè)是：距插值點(diǎn)近的樣本點(diǎn)，對其后方的樣本點(diǎn)有遮蔽效應(yīng)，當(dāng)兩樣本點(diǎn)與插值點(diǎn)的連線夾角α＜360°／n（n為插值搜索鄰域內(nèi)的樣點(diǎn)個數(shù)）時，遮蔽效應(yīng)存在，當(dāng)α≥360°／n時，遮蔽效應(yīng)消失。在AIDW插值過程中，受遮蔽影響的樣本點(diǎn)，其插值權(quán)重將被削弱，削弱的程度取決于該樣點(diǎn)K值的大?。▍⒁?.2.2節(jié)）。

按上述假設(shè)：圖1（a）所示的5個樣點(diǎn)在方向上均勻地分布在插值點(diǎn)（中心點(diǎn)）周圍，任意兩樣點(diǎn)與插值點(diǎn)的連線夾角均大于或等于72°（即α≥360°／5），即認(rèn)為該5個樣點(diǎn)間相互不存在遮蔽效應(yīng)；在圖1（c）中，任意兩樣點(diǎn)與插值點(diǎn)的連線夾角均小于72°，即認(rèn)為距插值點(diǎn)的近樣點(diǎn)，對其后的樣點(diǎn)均具有遮蔽效應(yīng)；在大多情況下，樣點(diǎn)在插值點(diǎn)周圍的分布應(yīng)類似圖1（b），既不像圖1（a）均勻分布，也不像圖1（c）集中分布。圖1（b）中Z1、Z3對任一樣點(diǎn)均無遮蔽，Z2對Z4、Z5有遮蔽，Z4對Z5也有遮蔽。

圖1 插值搜索鄰域中樣點(diǎn)散布示例（5個樣點(diǎn)）Fig.1 Examples of sample points scattering in searching neighborhood（five points）

2.2.2 算法表達(dá)式

將IDW傳統(tǒng)的算法思想與本文的基本假設(shè)結(jié)合，提出了AIDW算法，該算法可被表達(dá)為

式中，Z、Zi、di、n及p的定義與式（1）相同；i為插值搜索鄰域內(nèi)的樣點(diǎn)距插值點(diǎn)遠(yuǎn)近的排序序號（排列順序?yàn)橛山吝h(yuǎn)）；ki是序號i樣點(diǎn)的方位調(diào)和權(quán)重系數(shù)，它表示序號i樣點(diǎn)受其他樣點(diǎn)遮蔽的綜合效應(yīng)。因序號為1的樣點(diǎn)距插值點(diǎn)最近，任何樣點(diǎn)都不會對其有遮蔽影響，故設(shè)定k1＝1。sinpθij是序號i樣點(diǎn)受序號j樣點(diǎn)遮蔽程度的計算式

式中，1≤j≤i－1，θij是序號i、j兩樣點(diǎn)連線與其中線（過插值點(diǎn)）的夾角（銳角或直角）；aij是序號i、j兩樣點(diǎn)與插值點(diǎn)的連線夾角（圖2）。依據(jù)AIDW的基本假設(shè)，當(dāng)αij≥360°／n時，序號j樣點(diǎn)對序號i樣點(diǎn)無遮蔽影響，故規(guī)定此時的sinpθij值為1。

2.2.3 算式sinpθ對樣點(diǎn)遮蔽效應(yīng)的表達(dá)

圖2示例解釋了樣點(diǎn)Zi所受遮蔽效應(yīng)與其位置的關(guān)系以及sinpθ對這種關(guān)系的表達(dá)能力。在Zi以固定a角逐漸移近插值點(diǎn)Zo的過程中（圖2（a）），θ角和sinpθ值在不斷增大，Zj對Zi的遮蔽效應(yīng)逐步降低；當(dāng)Zi移動到距插值點(diǎn)的距離等于Zj距插值點(diǎn)距離時（即兩樣點(diǎn)已無先后，Zj、Zi與插值點(diǎn)構(gòu)成了等腰三角形），θ角為90°，sinpθ＝1，此時表明Zj已不會對Zi產(chǎn)生遮蔽影響。另一種情況，當(dāng)Zi保持與插值點(diǎn)的距離不變，但逐漸向Zj背后靠近時（圖2（b）），顯然樣點(diǎn)Zj對Zi的遮蔽效應(yīng)在逐步增大，這一過程中θ角和sinpθ值也在不斷變小；當(dāng)Zi移動到Zj與插值點(diǎn)連線的正后方，此時θ角變?yōu)?°，sinpθ＝0，表明Zj已完全遮蔽了Zi對插值點(diǎn)的影響。

圖2 sinpθ隨樣點(diǎn)位置的變化示例Fig.2 Examples of sinpθvalue changing with sample point moving

若將樣點(diǎn)Zi的變化位置看作是其他樣點(diǎn)的分布，則圖2解釋了插值樣點(diǎn)在兩種特例分布下θ和sinpθ的極端變化情況。即兩樣點(diǎn)與插值點(diǎn)的距離趨近于等同，θ→90°，sinpθ→1；某樣點(diǎn)位置趨向于另一樣點(diǎn)（與插值點(diǎn)連線）正后方，θ→0°，sinpθ→0。通常情況下，有遮蔽效應(yīng)的兩樣點(diǎn)，它們的θ角多介于0°～90°，sinpθ值多介于0～1。

2.3 AIDW的特點(diǎn)分析

由于AIDW綜合考慮了樣點(diǎn)距離和樣點(diǎn)方位對插值點(diǎn)的共同影響，因而其插值過程（圖3）和插值結(jié)果（圖4）較IDW更趨于合理。從圖3（a）中可看出，插值樣點(diǎn)在各向上的分布，較難滿足樣點(diǎn)均勻布局這一IDW插值基本假設(shè)，此種情況下IDW的插值可靠性難以保證。AIDW插值考慮了樣點(diǎn)間的方位影響，通過計算各樣點(diǎn)間的遮蔽效應(yīng)，削減受遮蔽樣點(diǎn)的插值權(quán)重。按AIDW算法，在圖3（b）中因Z1對Z6、Z3對Z7和Z8、Z4對Z7有遮蔽影響，這些受遮蔽樣點(diǎn)的插值權(quán)重被削減，Z10、Z11、Z12分別被Z4、Z3、Z7完全遮蔽，它們的插值權(quán)重降至為0。依照式（2）和式（3），最終插值點(diǎn)估算值的計算式為

式中，Z為插值點(diǎn)（中心點(diǎn)）估算值；Z1—Z9為樣點(diǎn)觀測值；d1—d9為樣點(diǎn)與插值點(diǎn)的歐氏距離；p是冪指數(shù)；θ角如圖3（b）所示。由此看來，從某種意義上可以說AIDW算法可使插值樣點(diǎn)的分布在各向上更趨于均勻。

圖3 IDW和AIDW兩種插值的樣點(diǎn)處理方式Fig.3 Treating mode of sample points of the IDW and AIDW interpolations

為了便于理解AIDW和IDW兩種插值在等值線繪制中的各自表現(xiàn)，本文用一個簡單特例對其加以表述。在圖4中紫色樣點(diǎn)與插值點(diǎn)（黑色點(diǎn)）呈等距離分布，距離為1；淺橙色樣點(diǎn)與插值點(diǎn)也呈等距離分布，距離為2，且淺橙色樣點(diǎn)位于插值點(diǎn)與紫色樣點(diǎn)連線的正后方，各樣點(diǎn)的觀測值如圖4所示。按IDW算法（p＝1），插值點(diǎn)的估算值為50；按AIDW（p＝1）算法，插值點(diǎn)的估算值為20（淺橙色樣點(diǎn)完全被遮蔽，僅需考慮紫色樣點(diǎn)插值權(quán)重）。從樣點(diǎn)數(shù)據(jù)值的分布可看出（圖4，等值線為30和40），很大程度上AIDW要比IDW的插值結(jié)果更趨于合理。IDW算法使低值點(diǎn)中出現(xiàn)了高值，其等值線繪制也就增大了“牛眼”的出現(xiàn)幾率；而AIDW插值結(jié)果與臨近數(shù)據(jù)值接近，很好地避免了“牛眼”現(xiàn)象（圖4（b））。有關(guān)AIDW和IDW的插值誤差及其等值線繪制的實(shí)例對比分析，參見第4節(jié)。

3 AIDW插值程序設(shè)計

AIDW的插值程序可分為插值前準(zhǔn)備和插值計算兩個過程。插值前準(zhǔn)備主要是用于搜索合適的插值樣點(diǎn)，并為下一步的插值計算提供di和fij值；插值計算過程主要是求算反映樣點(diǎn)遮蔽程度的sinpθij值，并結(jié)合di、Zi值求算插值點(diǎn)的Z值（圖5）。插值搜索鄰域的大小以格點(diǎn)數(shù)（k×k）表示，m是搜索鄰域內(nèi)的樣點(diǎn)數(shù)，n是插值所需的樣點(diǎn)數(shù)，d、f分別為樣點(diǎn)與插值點(diǎn)的歐氏距離和兩樣點(diǎn)間的歐氏距離，i、j、u、v均為插值樣點(diǎn)的序號（i＝1，2，…，n；有1≤j≤i－1、u＜j、i＜v≤n），aij、θij、Zi與式（4）中的定義相同。

sinpθij的求算涉及到兩次運(yùn)用三角余弦定理，一次是判別條件cosαij的求算，另一次是cosθij的求算。為減少程序運(yùn)算過程，若插值搜索鄰域固定樣點(diǎn)數(shù)的設(shè)定大于3個（即n＞3），則兩樣點(diǎn)是否存在遮蔽的判別條件，程序直接使用cosαij與cos（360°／n）對比，而不使用αij角度與360°／n對比，以便去除對αij角度的求解。另一方面，當(dāng)出現(xiàn)某樣點(diǎn)i被另一樣點(diǎn)j完全遮蔽（即cosαij＝1）時，該樣點(diǎn)i與其他樣點(diǎn)的遮蔽效應(yīng)，程序采取了判別賦值方式。因已有sinpθij＝0，故i樣點(diǎn)與其他樣點(diǎn)的遮蔽關(guān)系可不再考慮，方便起見可直接賦值sinpθui＝1和sinpθiv＝1（u＜j、i＜v≤n），這種判別賦值方式亦省去一些不必要的sinpθ求算。

圖5 AIDW的插值程序設(shè)計流程（n＞3）Fig.5 Flow chart of the AIDW interpolation program（n＞3）

4 AIDW空間插值實(shí)例分析

本文對浙江省69個氣象臺站的某年降水量觀測數(shù)據(jù)，分別利用IDW和AIDW插值算法進(jìn)行空間插值，進(jìn)而繪制降水量等值線圖。插值實(shí)例中，IDW和AIDW算法的p取值為2，插值搜索鄰域的固定樣點(diǎn)數(shù)為5，插值范圍為浙江全省，插值空間分辨率即格點(diǎn)大小設(shè)置為1km×1km。IDW插值及等值線繪制由ArcGIS軟件的ArcMap 9.3完成，AIDW插值和等值線繪制由文中的插值程序及文獻(xiàn)［27—28］中的等值線繪制方式完成。

浙江省年降水量IDW和AIDW兩種空間插值的誤差分析見表1。AIDW插值的平均誤差（ME）、平均絕對誤差（MAE）、平均相對誤差（MRE）以及均方根誤差（RMSE），均小于IDW插值方式，尤其是ME和RMSE這兩種誤差更是明顯低于IDW方式，這表明了AIDW相比于IDW，其獲得的空間插值結(jié)果更趨于合理。

表1 IDW和AIDW空間插值交叉驗(yàn)證的誤差分析Tab.1 Errors analysis of cross validation of the IDW and AIDW interpolation

從圖6可看出（等值線基值為1000mm雨量，間隔為50mm雨量）：IDW和AIDW兩者的等值線圖總體分布趨勢基本一致（圖6（a）），但在局部區(qū)域卻存在顯著差異。區(qū)域放大圖顯示：在浙西開化、浙西北安吉和長興、浙東南溫嶺、樂清、溫州、瑞安及平陽等區(qū)，IDW繪圖中出現(xiàn)的同心圓或孤立圓現(xiàn)象（圖6（b）），在 AIDW 繪圖中（圖6（c））得以消除。

圖4 IDW和AIDW插值方式的等值線繪制特例Fig.4 Typical examples of contour plotting for the IDW and AIDW interpolation methods

圖6 IDW和AIDW插值方式的等值線繪制實(shí)例對比Fig.6 Comparison of contour plotting cases of the IDW and AIDW interpolation methods

在浙西開化區(qū)，隨插值點(diǎn)距開化站（高雨量站）的距離變遠(yuǎn)，IDW插值的開化站權(quán)重會逐漸減小，而常山、江山、衢州和淳安等低雨量站的權(quán)重相應(yīng)增加，故越遠(yuǎn)離開化站，其降水插值結(jié)果越小，從而在開化站周圍形成了同心圓等值線（圖6（b））。在 AIDW 插值過程中，由于考慮了站點(diǎn)方位遮蔽效應(yīng)，當(dāng)對開化站左側(cè)區(qū)域插值時，開化站對周圍其他站產(chǎn)生了不同程度的遮蔽，且其他站之間也存在遮蔽，由此常山、江山、衢州、淳安等低雨量站的權(quán)重被大幅削減，故開化站左側(cè)區(qū)的插值未被低雨量站越拉越低，從而避免了同心圓等值線出現(xiàn)（圖6（c））。與開化情況相類似的還有浙西北的安吉和浙東南的溫嶺。在浙西北的長興區(qū)，IDW的等值線圖中有一個孤立圓（圖6（b）），這是因?yàn)楹莺偷虑鍍蓚€高雨量站的存在，使得長興站周圍的插值，隨距長興站越遠(yuǎn)，其插值結(jié)果越高，由此等值線繪制也不得不形成一個孤立圓區(qū)；而在AIDW插值過程中，當(dāng)對長興站左側(cè)區(qū)插值時，湖州和德清兩個高雨量站因受長興站的遮蔽以及湖州、安吉兩站對德清站的遮蔽，湖州和德清兩站的插值權(quán)重（尤其是德清站）明顯被減小，因此長興站左側(cè)區(qū)的插值結(jié)果未被明顯拉高，等值線繪制時此區(qū)也就與其他區(qū)域合并在了一起，并未出現(xiàn)孤立圓現(xiàn)象（圖6（c））。與長興情況相類似的還有浙西南樂清、溫州、瑞安及平陽等區(qū)。

值得注意的是：AIDW插值算法并未消除IDW繪圖中浙西龍游的孤立圓以及浙西北德清和杭州的同心圓現(xiàn)象（如圖6（b）和圖6（c））。仔細(xì)查看可發(fā)現(xiàn)，這些站周圍的各站點(diǎn)幾乎呈均勻布局（上下左右各1個），并不像開化和長興那樣集中于一側(cè)，這也意味著AIDW插值算法即可保留IDW插值思想，又可減小樣點(diǎn)分布不均IDW的插值缺陷。簡而言之，AIDW相對于IDW，其插值結(jié)果更趨于合理，等值線繪制也更接近于人工作業(yè)，如開化、安吉、長興、溫嶺、樂清、溫州、瑞安及平陽等區(qū)的等值線走向判識。

5 結(jié) 論

AIDW算法依據(jù)樣點(diǎn)與插值點(diǎn)的連線夾角a判別樣點(diǎn)間是否有遮蔽效應(yīng)，采用兩樣點(diǎn)連線與其中線（過插值點(diǎn)）夾角的sinpθ值表征樣點(diǎn)受遮蔽程度，以方位權(quán)重系數(shù)調(diào)和距離權(quán)重，插值過程考慮樣點(diǎn)距離和樣點(diǎn)方位共同影響，其空間插值結(jié)果更趨于合理?；贏IDW空間插值曲面的等值線繪制，可有效減少IDW等值線繪制中孤立圓和同心圓的出現(xiàn)幾率，使等值線繪圖更具科學(xué)美觀性。

因AIDW算法涉及cosα和cosθ值求算，故其空間插值運(yùn)算速度較IDW插值略慢。在浙江省69氣象站點(diǎn)的年降水量空間插值實(shí)例中，AIDW和IDW插值運(yùn)算時間均不足1s，但當(dāng)插值的空間分辨率由1km×1km提高到500m×500m，插值搜索鄰域固定樣點(diǎn)數(shù)設(shè)為15個時，浙江全省的AIDW插值運(yùn)算時間為11s，而IDW插值運(yùn)算時間為6s。隨插值格點(diǎn)總數(shù)以及插值搜索鄰域固定樣點(diǎn)數(shù)的增加，AIDW插值運(yùn)算速度會略慢于IDW插值。

［1］ LONGLEY P A，GOODCHILD M F，MAGUIRE D J，et al.Geographical Information Systems：Principles， Techniques，Management and Applications［M］.New York：Wiley，1999：481－492.

［2］ HENGL T A.Practical Guide to Geostatistical Mapping of Environmental Variables［M］.Ispra：Official Publications of the European Communities，2007：11－24.

［3］ LIU Jianjun，CHEN Jun，WANG Donghua，et al.The Study of Description and Application of Contour Adjacency Relationship［J］.Acta Geodaetica et Cartographica Sinica，2004，33（2）：174－178.（劉建軍，陳軍，王東華，等.等高線鄰接關(guān)系的表達(dá)及應(yīng)用研究［J］.測繪學(xué)報，2004，33（2）：174－178.）

［4］ YOU Shucheng，YAN Tailai.A Study on Artificial Neural Network Based Surface Interpolation［J］.Acta Geodaetica et Cartographica Sinica，2000，29（2）：30－34.（尤淑撐，嚴(yán)泰來.基于人工神經(jīng)網(wǎng)絡(luò)面插值的方法研究［J］.測繪學(xué)報，2000，29（2）：30－34.）

［5］ SHI Wenjiao，YUE Tianxiang，SHI Xiaoli.Research Progress in Soil Property Interpolators and Their Accuracy［J］.Journal of Natural Resources，2012，27（1）：163－175.（史文嬌，岳天祥，石曉麗.土壤連續(xù)屬性空間插值方法及其精度的研究進(jìn)展［J］.自然資源學(xué)報，2012，27（1）：163－175.）

［6］ LU G Y，WONG D W.An Adaptive Inverse Distance Weighting Spatial Interpolation Technique［J］.Computer and Geoscience，2008，34：1044－1055.

［7］ DENG Xiaobin.Comparison between Two Space Interpolation Methods Based on ArcGIS［J］.Geospatial Information，2008，6（6）：85－87.（鄧曉斌.基于 ArcGIS兩種空間插值方法的比較［J］.地理空間信息，2008，6（6）：85－87.）

［8］ KNOTTERS M，BRUS D J，VOSHAAR J H O.A Comparison of Kriging，Co－Kriging and Kriging Combined with Regression for Spatial Interpolation of Horizon Depth with Censored Observations［J］. Geoderma，1995，67（3）：227－246.

［9］ NIU Wenjie.Comparison of Theory and Application between Thin Plate Spline and Universal Kriging［J］.Journal of Engineering Graphics，2010，4：123－129.（牛文杰.薄板樣條法和泛克里金法在理論和應(yīng)用方面的比較［J］.工程圖學(xué)學(xué)報，2010，4：123－129.）

［10］ WANG Jingxin.Research on Spline Interpolation Poisedness and Interpolation Approach Problems［D］.Dalian：Dalian University of Technology，2004.（王晶昕.樣條插值適定性與插值逼近問題研究［D］.大連：大連理工大學(xué)，2004.）

［11］ TANG Bo，TONG Ling，KANG Shaozhong.Effects of Spatial Station Density and Interpolation Methods on Accuracy of Reference Crop Evapotranspiration［J］.Transactions of the Chinese Society of Agricultural Engineering，2013，29（13）：60－66.（湯博，佟玲，康紹忠.站點(diǎn)密度及插值方法對ET0空間插值精度的影響［J］.農(nóng)業(yè)工程學(xué)報，2013，29（13）：60－66.）

［12］ LI Xin，CHENG Guodong，LU Ling.Comparison of Spatial Interpolation Methods［J］.Advance in Earth Sciences，2000，15（3）：260－265.（李新，程國棟，盧玲.空間內(nèi)插方法比較［J］.地球科學(xué)進(jìn)展，2000，15（3）：260－265.）

［13］ LIU Jinsong，CHEN Hui，YANG Binyun.Comparison of Interpolation Methods on Annual Mean Precipitation in Hebei Province［J］.Acta Ecologica Sinica，2009，29（7）：3493－3500.（劉勁松，陳輝，楊彬云.河北省年均降水量插值方法比較［J］.生態(tài)學(xué)報，2009，29（7）：3493－3500.）

［14］ ZHANG Meihua，LIANG Wenkang.An Improved Method for Contouring on Network of Triangles［J］.Journal of Computer－Aided Design and Computer Graphics，1997，9（3）：213－217.（張梅華，梁文康.一個三角形網(wǎng)格上等值線圖的繪制算法［J］.計算機(jī)輔助設(shè)計與圖形學(xué)學(xué)報，1997，9（3）：213－217.）

［15］ HU Jinxing，MA Zhaoting，WU Huanping，et al.Massive Data Delaunay Triangulation Based on Grid Partition Method［J］.Acta Geodaetica et Cartographica Sinica，2004，33（2）：163－167.（胡金星，馬照亭，吳煥萍，等.基于格網(wǎng)劃分的海量數(shù)據(jù)Delaunay三角剖分［J］.測繪學(xué)報，2004，33（2）：163－167.）

［16］ ZHANG Yao，F(xiàn)AN Hong，HUANG Wang.The Method of Generating Contour Tree Based on Contour Delaunay Triangulation［J］. Acta Geodaetica et Cartographica Sinica，2012，43（3）：461－467.（張堯，樊紅，黃旺.基于Delaunay三角網(wǎng)的等高線樹生成方法［J］.測繪學(xué)報，2012，43（3）：461－467.）

［17］ DALE Z，CLAIRE P，AMY R，et al.An Experimental Comparison of Ordinary and Universal Kriging and Inverse Distance Weighting［J］.Mathematical Geology，1999，31（4）：375－390.

［18］ AZPURUA M，RAMOS K D.A Comparison of Spatial Interpolation Methods for Estimation of Average Electromagnetic Field Magnitude［J］.Progress In Electromagnetics Research M，2010，14：135－145.

［19］ CAROL A G，RICHARD B F，GARY W H，et al.Comparison of Kriging and Inverse Distance Methods for Mapping Soil Parameters［J］.Soil Science Society of A－merica Journal，1996，60（4）：1237－1247.

［20］ ZHU Huiyi，LIU Shulin，JIA Shaofeng.Problems of the Spatial Interpolation of Physical Geographical Elements［J］.Geographical Research，2004，23（4）：425－432.（朱會義，劉述林，賈紹鳳.自然地理要素空間插值的幾個問題［J］.地理研究，2004，23（4）：425－432.）

［21］ FENG Zhiming，YANG Yanzhao，DING Xiaoqiang，et al.Optimization of the Spatial Interpolation Methods for Climate Resources［J］.Geographical Research，2004，23（3）：357－364.（封志明，楊艷昭，丁曉強(qiáng)，等.氣象要素空間插值方法優(yōu)化［J］.地理研究，2004，23（3）：357－364.）

［22］ NUSRET D，DUG S.Applying the Inverse Distance Weighting and Kriging Methods of the Spatial Interpolation on the Mapping the Annual Precipitation in Bosnis and Rzegovina［C］∥Proceedings of International Congress on Environmental Modelling and Software.Leipzig：［s.n.］，2012.

［23］ SHEPARD D.A Two－dimensional Interpolation Function for Irregularly－Spaced Data［C］∥ Proceedings of the 1968 ACM National Conference.New York：ACM Press，1968：517－524.

［24］ WATSON D F，PHILIP G M.A Refinement of Inverse Distance Weighted Interpolation［J］. Geoprocessing，1985，2：315－327.

［25］ ZENG Hongwei，LI Lijuan，ZHANG Yongxuan，et al.Study on Spatial Interpolation of Precipitation with Large Scale Samples：A Case Study on 2009’s Precipitation of China［J］.Progress in Geography，2011，30（7）：811－818.（曾紅偉，李麗娟，張永萱，等.大樣本降水空間插值研究：以2009年中國年降水為例［J］.地理科學(xué)進(jìn)展，2011，30（7）：811－818.）

［26］ HUSAR R B，F(xiàn)ALKE S R.Uncertainty in the Spatial Interpolation of PM10Monitoring Data in Southern California［J／OL］.［1996－04－05］.http：∥capita.wustl.edu／CAPITA／capitareports／CAInterp／CAINTERP.doc.

［27］ XIE Shunping，TIAN Desen.To Improve Tracing Isopleth by Route Grid Method［J］.Acta Geodaetica et Cartographica Sinica，1995，24（1）：52－56.（謝順平，田德森.完善等值線追蹤的路徑柵格法［J］.測繪學(xué)報，1995，24（1）：52－56.）

［28］ WANG Yu，ZHU Changqing，SHI Wenzhong.Application of B Spline and Smoothing Spline on Interpolating the DEM Based on Rectangular Grid［J］.Acta Geodaetica et Cartographica Sinica，2000，29（3）：240－244.（王昱，朱長青，史文中.B樣條與磨光樣條在基于矩形格網(wǎng)的DEM內(nèi)插中 的應(yīng)用 ［J］.測繪學(xué)報，2000，29（3）：240－244.）