幾種插值方法在微DEM構(gòu)建中的應(yīng)用

宋向陽,吳發(fā)啟

（西北農(nóng)林科技大學(xué)資源環(huán)境學(xué)院,陜西楊陵 712100）

DEM是地形分析的主要數(shù)據(jù),其精度對(duì)地形分析結(jié)果的精度具有明顯影響[1]。目前,利用離散高程點(diǎn)數(shù)據(jù)進(jìn)行空間插值來構(gòu)建DEM是生成地形分析數(shù)據(jù)的主要途徑[2],也是ArcGIS軟件地形分析的重要功能之一[3]。一直以來,空間插值方法對(duì)DEM精度的影響都得到研究人員的廣泛關(guān)注。早在1988年,Hutchinson提出利用高程信息生成DEM的Hutchinson算法[4],并在其基礎(chǔ)上開發(fā)出專業(yè)化軟件ANUDEM[5];鄧小波等[6]研究了不同插值方法對(duì)氣候要素分析的影響,認(rèn)為改進(jìn)插值方法和優(yōu)化插值參數(shù)能夠改善插值效果。另外,湯國安等[7]的研究表明,不同的插值方法生成的DEM具有明顯的差異,而這種差異對(duì)地形分析結(jié)果的影響比較顯著;孟慶香等[8]研究基于GIS的黃土高原降水量和年均溫空間插值,結(jié)果表明不同插值方法構(gòu)建的降水量和年均溫分布圖差異較大;李新等[9]通過比較研究不同空間插值方法差異,結(jié)果表明要得到理想的空間插值效果,必須針對(duì)研究區(qū)的實(shí)際情況,在對(duì)實(shí)測(cè)樣點(diǎn)數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析和反復(fù)試驗(yàn)比較的基礎(chǔ)上選擇最佳的插值方法。前人的這些研究成果主要是針對(duì)宏觀尺度的DEM構(gòu)建及分析,而關(guān)于微觀尺度DEM構(gòu)建方法的研究還幾乎是一個(gè)空白。

眾所周知,黃土高原是我國乃至世界上最為嚴(yán)重的水土流失地區(qū)之一。在該區(qū)坡耕地的土壤侵蝕量占到總侵蝕量的50%～60%[10]。雖該地區(qū)造成土壤侵蝕的因素眾多,但人為土地管理方式是影響土壤侵蝕的主要原因之一。這些管理往往在耕地表面留下土粒級(jí)和土塊級(jí)的高低起伏,直接影響到產(chǎn)流、匯流和泥沙運(yùn)動(dòng)。因此,這種微地形特征(或稱為地表糙度)的研究,對(duì)揭示水土流失成因有重要作用。鑒于此,本文在模擬該區(qū)比較常見的三種耕作措施,即等高耕作、人工掏挖和人工鋤耕,并以直線坡面為對(duì)照措施,利用非接觸式激光掃描儀測(cè)量這幾種措施下的地表高程的基礎(chǔ)上,結(jié)合GIS軟件中數(shù)據(jù)分析工具,以均方根預(yù)測(cè)誤差、誤差平均值和誤差圖等指標(biāo)對(duì)測(cè)定的地表高程數(shù)據(jù)進(jìn)行探索性分析,以便確定構(gòu)建微DEM的最佳插值方法,服務(wù)于黃土坡耕地土壤侵蝕機(jī)理的研究。

1 研究方法

本研究借用了以下幾種前人的研究方法來建立微DEM。

1.1 反距離加權(quán)插值法

反距離加權(quán)插值法(Inverse Distance Weighting,簡稱IDW)[11]是利用鄰近已知點(diǎn)的數(shù)值進(jìn)行加權(quán)運(yùn)算,所需的權(quán)重根據(jù)距離的冪來確定,離插值點(diǎn)越近的樣本點(diǎn)賦予的權(quán)重越大。權(quán)重顯著影響內(nèi)插結(jié)果,其選擇標(biāo)準(zhǔn)是最小平均絕對(duì)誤差[12]。Husar[13]等人的研究結(jié)果表明,距離的冪越高,內(nèi)插結(jié)果越具有平滑的效果。該方法的插值計(jì)算是一個(gè)均分過程,要求采樣點(diǎn)均勻分布,并且密集程度足以滿足在分析中反映局部表面變化。

1.2 局部多項(xiàng)式插值法

局部多項(xiàng)式插值法(Local Polynomial Interpolation,簡稱LPI)[14]是一種局部加權(quán)最小二乘擬合法,根據(jù)有限的采樣數(shù)據(jù),采用多個(gè)多項(xiàng)式來擬合表面,它引入了“距離權(quán)”的概念,對(duì)于未知樣點(diǎn)的計(jì)算,考慮在局部范圍內(nèi)全部已知樣點(diǎn)對(duì)其貢獻(xiàn),距未知樣點(diǎn)近的點(diǎn)權(quán)重較大,距未知點(diǎn)遠(yuǎn)的點(diǎn)權(quán)重較小,每一個(gè)未知樣點(diǎn)的預(yù)測(cè)值都對(duì)應(yīng)一個(gè)多項(xiàng)式,每個(gè)多項(xiàng)式都處在特定重疊的鄰近區(qū)域內(nèi),通過最小二乘法求解鄰域內(nèi)多項(xiàng)式組成的方程組來得到擬合表面。

1.3 徑向基函數(shù)插值法

徑向基函數(shù)插值法(Radial Basis Function,簡稱RBF)[15]屬于人工神經(jīng)網(wǎng)絡(luò)方法的一種,該方法是根據(jù)有限采樣數(shù)據(jù),選擇合適的徑向基函數(shù)生成一個(gè)具有最小曲率,且到各樣點(diǎn)的Z值的距離最小的曲面。該方法所擬合的表面經(jīng)過所有樣點(diǎn)數(shù)據(jù),并可以計(jì)算出高于或低于樣點(diǎn)Z值的預(yù)測(cè)值,適用于采樣點(diǎn)數(shù)據(jù)集大、表面變化平緩的情況,當(dāng)局部變異性大,且無法確定樣點(diǎn)數(shù)據(jù)的準(zhǔn)確性,或樣點(diǎn)數(shù)據(jù)具有很大不確定性時(shí),該方法并不適用。

1.4 克里格插值法

克里格插值法(Kriging)[16]是用協(xié)方差函數(shù)和變異函數(shù)來確定高程變量隨空間距離而變化的規(guī)律,以距離為自變量的變異函數(shù),計(jì)算相鄰高程值關(guān)系權(quán)值,在有限區(qū)域內(nèi)對(duì)區(qū)域化變量進(jìn)行無偏最優(yōu)估計(jì)的一種方法。從數(shù)學(xué)的角度來講,克里格插值法包括析取克里格(Disjunctive Kriging)、普通克里格(Ordinary Kriging)和概率克里格(Probability Kriging)等。對(duì)不同的樣本數(shù)據(jù),按照其數(shù)據(jù)特征的不同應(yīng)該選擇不同的插值方法才能取得較好的插值效果,如普通克里格法前提條件是樣本數(shù)據(jù)符合正態(tài)分布;若樣本不服從正態(tài)分布時(shí),一般選用析取克里格;當(dāng)同一事物的兩種屬性存在相關(guān)關(guān)系,且一種屬性不易獲取時(shí),可選用協(xié)同克里格方法,借助另一屬性實(shí)現(xiàn)該屬性的空間內(nèi)插。

2 試驗(yàn)數(shù)據(jù)的采集與處理

2.1 試驗(yàn)設(shè)計(jì)

試驗(yàn)在中國科學(xué)院黃土高原土壤侵蝕與旱地農(nóng)業(yè)國家重點(diǎn)實(shí)驗(yàn)室降雨大廳進(jìn)行。試驗(yàn)小區(qū)規(guī)格為1.0 m×1.0 m×0.5 m(以下簡稱:小區(qū)),每個(gè)措施設(shè)置3個(gè)重復(fù)。小區(qū)地表處理⑴土樣準(zhǔn)備:將過篩(篩孔0.5 cm)的土樣分層填裝在小區(qū)中,并將土壤容重控制在1.20～1.30 g/cm3,含水量控制在10%左右;⑵耕作措施的布設(shè):在填好土的侵蝕槽中人工構(gòu)建不同的耕作措施,包括:①等高耕作:選用橫坡耕作方式,壟高7～10 cm,壟距為30 cm;②人工掏挖:采用镢頭掏挖地表,深度5～8 cm,間距20～25 cm;③人工鋤耕:沿地表以傳統(tǒng)方式鋤耕,深度4～5 cm;④CK:平整坡面,不做任何處理(見圖1)。

2.2 數(shù)據(jù)采集

本研究數(shù)據(jù)利用非接觸式激光掃描儀對(duì)微地表進(jìn)行掃描后獲得,每個(gè)處理下數(shù)據(jù)為46行,40列,共1 840個(gè)高程值,行列方向上空間間隔均為0.02 m,其存儲(chǔ)格式如表1。

2.3 數(shù)據(jù)處理

數(shù)據(jù)處理在ArcGIS 9.3中進(jìn)行,使用探索性數(shù)據(jù)分析工具Geostatistical Analyst-Explore Data所提供的一系列圖形工具探測(cè)數(shù)據(jù)分布、全局和局部異常值(過大值或過小值)、尋求全局的變化趨勢(shì)、研究空間自相關(guān),得到不同耕作措施下地表高程值統(tǒng)計(jì)特征,如表2。

圖1 耕作措施布設(shè)

表1 高程采樣數(shù)據(jù)

表2 不同耕作措施下地表相對(duì)高程統(tǒng)計(jì)特征表

從表2中可以看出,標(biāo)準(zhǔn)差分別為CK(0.004 8)＜等高耕作(0.012 4)＜人工鋤耕(0.012 7)＜人工掏挖(0.018 3),表明采樣點(diǎn)數(shù)據(jù)具有較高的精度。同時(shí)反映出規(guī)則采樣對(duì)地表微地形起伏狀況變化不明顯的CK表現(xiàn)能力最強(qiáng),其次是等高耕作,而人工鋤耕和人工掏挖相對(duì)較弱;地表微地形起伏狀況數(shù)據(jù)表明CK的起伏度(0.023)＜人工鋤耕(0.063)＜等高耕作(0.075)＜人工掏挖(0.094),同時(shí)最大高程和最低高程數(shù)據(jù)都出現(xiàn)在人工掏挖措施下的采樣數(shù)據(jù)中,分別為0.270和0.176,這表明人工掏挖措施微地形局部變異較大;空間變異系數(shù)的變化分別為等高耕作(0.99)＞人工掏挖(0.97)＞CK(0.87)＞人工鋤耕(0.72),表明采樣數(shù)據(jù)均具有較強(qiáng)的空間相關(guān)性;峰度和偏態(tài)數(shù)據(jù)均不為0,這表明采樣數(shù)據(jù)均不符合正態(tài)分布,但等高耕作措施下的地表微地形數(shù)據(jù)更接近正態(tài)分布,其偏態(tài)和峰度分別為0.003 2和1.811。

結(jié)果表明,由于偏態(tài)和峰度以及直方圖、QQplot圖等表明采樣點(diǎn)數(shù)據(jù)不符合正態(tài)分布。所以,不宜采用采用普通克里格、簡單克里格和泛克里格插值法。但是,高程采樣點(diǎn)數(shù)據(jù)足以反映地表微地形局部表面變化且具有較強(qiáng)的空間相關(guān)性,達(dá)到IDW插值法、局部多項(xiàng)式插值法、徑向基函數(shù)插值方法和析取克里格插值法對(duì)空間數(shù)據(jù)的要求,因此,理論上可以采用這些方法進(jìn)行插值分析。

3 結(jié)果與分析

3.1 鄰域搜索半徑的設(shè)置

對(duì)于高程這種區(qū)域化變量來說,鄰域搜索半徑的大小不僅決定著鄰域采樣數(shù)據(jù)對(duì)預(yù)測(cè)點(diǎn)的權(quán)重,而且影響插值精度和插值表面的光滑度。因此,根據(jù)采樣點(diǎn)數(shù)據(jù)特征,設(shè)置鄰域搜索合理的鄰域搜索半徑是空間插值的基礎(chǔ),也是插值精度的保證。

對(duì)于等高耕作,由于采樣數(shù)據(jù)橫向相鄰點(diǎn)高程相似度較高,區(qū)域扇形設(shè)置應(yīng)更加扁平且平行于等高線方向,扇區(qū)橢圓長半軸(Major semiaxis)為 0.075,短半軸(Minor semiaxis)為0.03,角度(Angle)為90,步長0.02,搜索半徑內(nèi)使用最小點(diǎn)數(shù)為橫向相鄰的2個(gè)點(diǎn),最大為橫向4個(gè)。

對(duì)于人工掏挖、人工鋤耕和直線坡面,由于采樣數(shù)據(jù)橫向與縱向相鄰點(diǎn)高程相似度差異不大,區(qū)域扇形設(shè)置應(yīng)接近圓形。扇區(qū)橢圓長半軸為0.04,短半軸為0.04,角度為0,步長為0.02,搜索半徑內(nèi)使用最小點(diǎn)數(shù)為中心4個(gè)點(diǎn),最大點(diǎn)數(shù)為12個(gè)。

3.2 不同插值方法比較

通過上面的分析,在了解了數(shù)據(jù)特征的基礎(chǔ)上,本文擬選擇反距離加權(quán)插值法、局部多項(xiàng)式插值法、析取克里格插值法和徑向基函數(shù)插值法,分別對(duì)4種微地形數(shù)據(jù)進(jìn)行插值分析。通過多次計(jì)算,比較誤差平均值、均方根預(yù)測(cè)誤差、標(biāo)準(zhǔn)誤差預(yù)測(cè)圖以及使用軟件中的Compare功能等,得出針對(duì)不同耕作措施下4種空間插值方法的最優(yōu)參數(shù)和最佳插值結(jié)果,如表3。

表3 空間插值方法比較分析

由表3可以看出,不同插值方法對(duì)于不同地表微DEM構(gòu)建的差異較大。對(duì)于等高耕作來說,球面模型的析取克里格插值法優(yōu)于其他三種,其絕對(duì)誤差平均值為DK(0.19)＜LPI(0.32)＜IDW(1.20)＜RBF(1.34),同時(shí)其均方根預(yù)測(cè)誤差也最小(2.26)。這表明析取克里格插值方法對(duì)于具有空間相關(guān)性但不符合正態(tài)分布的高程數(shù)據(jù)是最優(yōu)插值方法;對(duì)于人工掏挖來說,Power=2時(shí)的局部多項(xiàng)式插值法明顯優(yōu)于其他方法,其絕對(duì)誤差平均值為 LPI(0.07)＜RBF(0.13)＜DK(0.25)＜IDW(0.63),均方根預(yù)測(cè)誤差為LPI(3.31)＜RBF(3.60)＜DK(6.02)＜IDW(6.21)。這表明地表微地形存在局部短程變異時(shí),應(yīng)該優(yōu)先考慮局部多項(xiàng)式插值方法;對(duì)于人工鋤耕來說,基于規(guī)則樣條函數(shù)的徑向基函數(shù)插值方法是最優(yōu)插值方法。其絕對(duì)誤差平均值為 RBF(0.05)＜LPI(0.10)＜IDW(0.28)＜DK(1.43),均方根預(yù)測(cè)誤差為RBF(4.46)＜DK(5.76)＜LPI(6.49)＜IDW(7.24)。這表明基于規(guī)則樣條函數(shù)的徑向基函數(shù)插值方法更適合于起伏狀況變化不大的微地形表面;對(duì)于CK來說,各種插值方法精度差異不大,其絕對(duì)誤差平均值為RBF(0.02)＜DK(0.03)＜LPI(0.07)＜IDW(0.10),均方根預(yù)測(cè)誤差為LPI(1.45)=IDW(1.45)＜DK(1.47)＜RBF(1.49)。這表明當(dāng)高程采樣點(diǎn)數(shù)據(jù)足以表現(xiàn)地表微地形局部變化時(shí),對(duì)于微地形起伏狀況無明顯變化的區(qū)域,插值方法對(duì)插值表面精度影響不大。

圖2顯示了利用不同插值方法對(duì)不同地表微地形插值計(jì)算的效果。從圖中可以看出不同插值方法在微DEM構(gòu)建過程中精度差異。對(duì)于等高耕作和人工掏挖來說,局部多項(xiàng)式插值法和析取克里格插值法效果相近且都較好,而基于規(guī)則樣條函數(shù)的徑向基函數(shù)插值法和反距離加權(quán)插值法效果較差,這和表3中結(jié)果相同。這表明局部多項(xiàng)式插值法和析取克里格插值法對(duì)于微地形變化具有一定方向性以及存在局部短程變異的地表能夠取得較高的精度和更加光滑的表面;對(duì)于人工鋤耕來說,基于規(guī)則樣條函數(shù)的徑向基函數(shù)插值方法效果最佳,其次是局部多項(xiàng)式插值方法和析取克里格插值方法,反距離加權(quán)插值方法效果最差。這一結(jié)果與表3中反映的精度差異一致,這表明地表微地形起伏變化較為平緩區(qū)域應(yīng)該優(yōu)先考慮基于規(guī)則樣條函數(shù)的徑向基函數(shù)插值方法;對(duì)于CK來說,析取克里格插值表面最為光滑,局部多項(xiàng)式插值方法次之,徑向基函數(shù)插值方法和反距離加權(quán)插值方法較差。雖然表3中結(jié)果表明各插值方法對(duì)CK的差異不大,但是通過圖2表明插值精度較高的插值方法取得的微地形表面不一定光滑。因此,在插值過程中還要在插值精度較高的方法中選擇效果最佳的插值方法。

圖2 利用不同插值方法對(duì)不同地表微地形插值計(jì)算的效果

3.3 最優(yōu)插值方法的選擇

通過對(duì)比表3中不同耕作措施下各插值方法的誤差平均值、均方根預(yù)測(cè)誤差,來定量判斷各插值方法的準(zhǔn)確性。同時(shí),對(duì)比圖2中插值表面的光滑程度來綜合考慮不同耕作措施條件下構(gòu)建微DEM的最佳插值方法。

4 結(jié)論

優(yōu)選空間插值方法和最優(yōu)參數(shù),能夠更加客觀、準(zhǔn)確地反映不同耕作措施下地表微地形起伏狀況的空間分布特點(diǎn)和變化趨勢(shì),便于構(gòu)建微地形高精度DEM。因此,本文對(duì)ArcGIS 9.3中4種常用空間插值方法進(jìn)行比較研究,通過對(duì)誤差平均值和均方根預(yù)測(cè)誤差等指標(biāo)的比較,得到不同耕作措施下構(gòu)建厘米級(jí)微DEM的最佳插值方法。

結(jié)果表明,反距離加權(quán)插值方法在4種耕作措施下誤差平均值和均方根預(yù)測(cè)誤差皆比較大,因而不適合進(jìn)行微地形的構(gòu)建。這主要是由于該方法僅考慮距離作為區(qū)域變量相關(guān)性的度量,而忽略了樣點(diǎn)的空間分布格局;局部多項(xiàng)式插值方法在人工掏挖措施下兩指標(biāo)都比較小,這主要是由于該方法能夠較好地保留地表微地形細(xì)節(jié)變化,鄰域重疊保持了插值表面的連續(xù)性,更適合于存在局部短程變異的數(shù)據(jù)集;析取克里格插值方法在等高耕作措施下兩指標(biāo)都比較小,這主要是因?yàn)樵摲椒紤]了已知樣本點(diǎn)的空間分布及與未知樣點(diǎn)的空間方位關(guān)系,更適合于微地形起伏狀況具有一定方向性變化的數(shù)據(jù)集;基于徑向基函數(shù)的規(guī)則樣條函數(shù)插值法在人工掏挖措施下兩指標(biāo)都比較小,其主要原因是規(guī)則樣條函數(shù)更適合于表面變化平緩的微地形的擬合;以上幾種插值方法在直線坡面耕作措施下,兩指標(biāo)都最小,這主要是由于采樣數(shù)據(jù)足夠精確情況下,各種插值方法表現(xiàn)差異不大所致。

因此,等高耕作措施微地表采用基于球面模型的析取克里格插值法精度較高;人工掏挖措施微地表采用局部多項(xiàng)式插值法效果較好;人工鋤耕措施微地表采用規(guī)則樣條函數(shù)的徑向基函數(shù)插值方法精度較好;對(duì)于直線坡面上述幾種插值方法的精度差異不大。

[1]李志林,朱慶.數(shù)字高程模型[M].2版.武漢:武漢大學(xué)出版社,2007.

[2]Wackernagel H.Multivariate Geostatistics:An Introduction with application[M].3rd Edition.Berlin:Springer Verlag,2003:416-429.

[3]王勁峰,李連發(fā),葛詠,等.地理信息空間分析的理論體系探討[J].地理學(xué)報(bào),2000,55(1):92-98.

[4]Hutchinson M F.Anew procedure for girding elevation and strearm line data with automatic removal of spurious pits[J].Journal or Hydrology,1989,106(3/4):211-232.

[5]Hutchinson M F.ANUDEM version 5.1 User Guide[C]//Centre for Resource and Environmental Studies.The Australian National University,Canberra,2004:1-22.

[6]鄧小波,羅宇翔,于飛,等.西南復(fù)雜山地農(nóng)業(yè)氣候要素空間插值方法比較[J].中國農(nóng)業(yè)氣象,2008,29(4):458-462.

[7]楊昕,湯國安,劉學(xué)軍,等.數(shù)字地形分析的理論方法與應(yīng)用[J].地理學(xué)報(bào),2009,64(9):1058-1070.

[8]孟慶香,劉國彬,楊勤科.基于GIS的黃土高原氣象要素空間插值方法[J].水土保持研究,2010,17(1):10-14.

[9]李新,陳國棟,盧玲.空間內(nèi)插方法比較[J].地球科學(xué)研究,2000,15(3):260-264.

[10]Tang Keli,Zhen Fenli,Ca Xuan.Soil erosion on the sloping farmland in the loess plateau of China[C]//Proceedings of the Fourth International Symposium on River Sedimentation.Beijing:China Ocean Press,1989.

[11]常文淵,戴新剛,陳洪武.地質(zhì)統(tǒng)計(jì)學(xué)在氣象要素場(chǎng)插值的實(shí)例研究[J].地球物理學(xué)報(bào),2004,47(6):982-990.

[12]方書敏,錢正堂,李遠(yuǎn)平.甘肅省降水的空間內(nèi)插方法比較[J].干旱區(qū)資源與環(huán)境,2005,19(3):47-50.

[13]Husar R B,Falke S R.Uncertainty in the spatial interpolation of PM 10 monitoring data in Southern California[EB/OL].http://capita.wustl.edu/CA PITA/CapitaReports/CaInterp/CaINTERP.HTM L,1997-03-03/1999-10-25.

[14]Fan Jianqing,Irè ne Gijbels.Local polynomial modelling and its applications[M].Chapman&Hall/CRC,2003:159-216.

[15]周志剛,陳麗紅.RBF神經(jīng)網(wǎng)絡(luò)及其在數(shù)值計(jì)算中的應(yīng)用[J].復(fù)雜系統(tǒng)與復(fù)雜性科學(xué),2006,2(3):64-68.

[16]曾懷恩,黃聲享.基于Kriging方法的空間數(shù)據(jù)插值研究[M].測(cè)繪工程,2007,16(5):5-13.