插值條件下格網(wǎng)DEM坡度計(jì)算模型的噪聲誤差分析

盧華興，劉學(xué)軍，王永君，田梓文

1.東南大學(xué) 交通學(xué)院，江蘇 南京 210018；2.南京師范大學(xué) 地理科學(xué)學(xué)院，江蘇 南京 210097；3.國(guó)家海洋局 第一海洋研究所，山東 青島 266100

1 引 言

坡度（slope）是表達(dá)地形曲面結(jié)構(gòu)和形態(tài)的基本參數(shù)之一，也是地表過(guò)程模擬、土壤侵蝕模型、土地利用分類等環(huán)境模型的基本變量，在地貌、地理、土壤、資源與環(huán)境等領(lǐng)域有著廣泛的應(yīng)用［1］。目前在基于格網(wǎng)DEM數(shù)據(jù)的數(shù)字地形分析中，格網(wǎng)DEM的坡度計(jì)算誤差主要來(lái)源于DEM數(shù)據(jù)誤差、DEM 數(shù)據(jù)結(jié)構(gòu)和分析算法［2］，國(guó)內(nèi)外學(xué)者從不同的角度對(duì)這一問(wèn)題進(jìn)行了探討和分析，如DEM分辨率對(duì)坡度計(jì)算精度的影響［3－4］、DEM 誤差 在 地 形 分 析 中 的 傳 遞 分 析［5－6］、不同的地形分析算法的精度評(píng)價(jià)［7－8］等，然而這些研究的不足之處在于并沒(méi)有考慮DEM誤差的空間自相關(guān)特性。事實(shí)上，DEM作為地形曲面的離散化表達(dá)，由于插值技術(shù)的存在，DEM誤差在空間上是相關(guān)的，而DEM誤差的空間自相關(guān)性已經(jīng)被證明存在并極大影響DEM及數(shù)字地形分析的不確定性，如文獻(xiàn)［9］給出了坡度中誤差、DEM中誤差及DEM誤差相關(guān)系數(shù)之間的計(jì)算公式，并從試驗(yàn)上證明當(dāng)DEM誤差相互獨(dú)立時(shí)，坡度誤差可達(dá)33%，遠(yuǎn)大于通過(guò)計(jì)算值和實(shí)測(cè)值對(duì)坡度誤差的估計(jì)；文獻(xiàn)［10］以隨機(jī)模擬的方式，通過(guò)改變DEM誤差的空間自相關(guān)大小，發(fā)現(xiàn)等高線密度（與面積百分比）隨著莫蘭指數(shù)I增大而降低，而河網(wǎng)密度（百分比）卻在I＝0.4附近達(dá)到了最?。晃墨I(xiàn)［11］采用指數(shù)自相關(guān)函數(shù)模型，從試驗(yàn)角度證實(shí)了DEM誤差空間自相關(guān)性對(duì)坡度精度的影響；文獻(xiàn)［2］采用隨機(jī)模擬并結(jié)合試驗(yàn)驗(yàn)證的方式，通過(guò)選擇兩種空間自相關(guān)函數(shù)模型（線形模型、指數(shù)模型）論述在兩種空間自相關(guān)函數(shù)模型下的窗口分析的協(xié)方差矩陣坡度誤差的影響。然而，以上模擬或試驗(yàn)方法在數(shù)字地形分析精度評(píng)價(jià)方面多數(shù)采用經(jīng)驗(yàn)空間自相關(guān)模型，經(jīng)驗(yàn)?zāi)Ｐ痛嬖谌缦氯毕荩孩倏臻g自相關(guān)模型選擇及其參數(shù)的確定上具有盲目性，無(wú)法忠實(shí)地描述DEM誤差空間自相關(guān)性規(guī)律；②空間自相關(guān)模型在整個(gè)數(shù)字地形分析精度評(píng)價(jià)過(guò)程中一成不變，無(wú)法描述DEM格網(wǎng)的誤差空間分異特征。

本文首先從DEM統(tǒng)一插值模型入手，推導(dǎo)插值條件下的格網(wǎng)DEM鄰近格網(wǎng)之間噪聲誤差的局部自相關(guān)性模型，從理論上描述了格網(wǎng)DEM噪聲誤差對(duì)坡度精度影響，最后從試驗(yàn)角度具體分析3種典型插值條件下的坡度計(jì)算的局部空間自相關(guān)性對(duì)數(shù)字地形分析精度的影響，并評(píng)價(jià)了不同地形分析算法對(duì)DEM噪聲誤差的抗差能力。

2 格網(wǎng)DEM插值的噪聲誤差分析

DEM插值是指由采樣點(diǎn)的高程信息按照特定的空間自相關(guān)性規(guī)則去估計(jì)待測(cè)點(diǎn)高程的數(shù)學(xué)方法，其目的是缺值估計(jì)和散點(diǎn)數(shù)據(jù)的格網(wǎng)化［12］，目前格網(wǎng)DEM多數(shù)由散點(diǎn)或地形特征數(shù)據(jù)的插值獲得，DEM插值過(guò)程其實(shí)是已知點(diǎn)權(quán)重分配的過(guò)程，即DEM插值的最終結(jié)果是已知點(diǎn)高程向量的一種線性組合［13］，設(shè)已知點(diǎn)高程向量z：z1，z2，…，zn（n為已知點(diǎn)個(gè)數(shù)），則待插點(diǎn)P（xP，yP）的高程zP可表示為

根據(jù)公式（1），格網(wǎng)P和Q的高程估值為

上式的k和l表示已知點(diǎn)權(quán)重向量，對(duì)式（2）寫(xiě)成矩陣形式如下

圖1 散點(diǎn)插值Fig.1 Elevation points interpolation

插值條件下，格網(wǎng)DEM的噪聲誤差是指已知點(diǎn)上的噪聲誤差在插值函數(shù)作用下傳遞到格網(wǎng)點(diǎn)上的噪聲誤差，格網(wǎng)DEM上的噪聲誤差的中誤差為格網(wǎng)DEM噪聲精度。根據(jù)方差協(xié)方差傳播率，并顧及已知點(diǎn)是獨(dú)立觀測(cè)數(shù)據(jù)，于是格網(wǎng)P、Q的方差協(xié)方差陣為

上式為2×2階矩陣，主對(duì)角線元素分別表示格網(wǎng)P、Q的指已知點(diǎn)噪聲誤差經(jīng)插值的傳遞誤差，表達(dá)了格網(wǎng)P、Q插值精度

而非對(duì)角線元素、分別格網(wǎng)P、Q之間的協(xié)方差，表達(dá)了臨近格網(wǎng)P、Q之間的誤差的相關(guān)性

由上式可知，造成臨近格網(wǎng)的噪聲誤差存在相關(guān)性的根本原因是插值過(guò)程中共用了已知點(diǎn)，若臨近格網(wǎng)在插值過(guò)程中沒(méi)有共用已知點(diǎn)，其協(xié)方差為0，表示兩者不考慮或不存在相關(guān)性。

3 格網(wǎng)DEM坡度計(jì)算模型的誤差分析

格網(wǎng)DEM坡度計(jì)算通常采用鄰域分析方法（neighborhood analysis，NA），即以固定尺寸如3×3、5×5或7×7等二維滑動(dòng)窗口逐行或逐列對(duì)格網(wǎng)DEM數(shù)據(jù)進(jìn)行掃描，每次滑動(dòng)過(guò)程中，對(duì)局部窗口里的格網(wǎng)高程值按照特定差分方式估計(jì)偏導(dǎo)數(shù)，用于估計(jì)中心格網(wǎng)的坡度值。一般的，格網(wǎng)DEM窗口分析算法可以寫(xiě)成局部窗口格網(wǎng)高程值的函數(shù)，以3×3窗口分析為例（如圖2所示），設(shè)9個(gè)格網(wǎng)對(duì)應(yīng)的高程值為z0、z1、…、z8，窗口分析模型可以寫(xiě)成如下形式

圖2 DEM 3×3鄰域窗口Fig.2 3×3DEM local moving window

根據(jù)誤差傳播定理中函數(shù)誤差的傳播公式［14］，格網(wǎng)DEM鄰域分析的中誤差為

式中，右邊括號(hào)內(nèi)的變量表示窗口分析函數(shù)對(duì)局部格網(wǎng)z＝ ［z0z1…z8］的偏導(dǎo)數(shù)，Dz為z＝［z0z1…z8］的方差協(xié)方差陣，其主對(duì)角線為格網(wǎng)的方差，而非對(duì)角線上的元素為格網(wǎng)之間的協(xié)方差。本文以坡度算法為研究對(duì)象，其他局部窗口分析算法的誤差分析與此類同。在坡度計(jì)算過(guò)程中，坡度提取的噪聲誤差是格網(wǎng)DEM上的噪聲誤差在插值函數(shù)作用下傳遞到格網(wǎng)點(diǎn)上的噪聲誤差，即已知點(diǎn)上的隨機(jī)誤差經(jīng)兩次傳遞得到的噪聲誤差，格網(wǎng)上坡度的噪聲誤差的中誤差為坡度噪聲精度。根據(jù)文獻(xiàn)［2］的推導(dǎo)，坡度中誤差為

式中，Dz為z＝［z0z1z2z3z4z5z6z7z8］T的協(xié)方差陣（公式（10）），p＝APz，q＝Aqz，而Ap和Aq是向量z的差分系數(shù)，不同的坡度差分算法有不同的系數(shù)，為方便計(jì)算，在文獻(xiàn)［2］的基礎(chǔ)上，這里對(duì)4種坡度計(jì)算模型中的差分系數(shù)作進(jìn)一步整理，得到如表1的差分系數(shù)。

表1 坡度計(jì)算模型的差分系數(shù)Tab.1 Difference coefficient of slope algorithms

需要指出的是，在坡度中誤差的計(jì)算過(guò)程中（公式（10）），對(duì)于Dz的求解，文獻(xiàn)［2］采用的是經(jīng)驗(yàn)協(xié)方差函數(shù)，本文所采用的是理論協(xié)方差模型（公式（4）），在計(jì)算鄰域窗口任意兩臨近格網(wǎng)的誤差空間相關(guān)性時(shí)，需要在插值過(guò)程中記錄插值點(diǎn)的點(diǎn)號(hào)和權(quán)重，并且在協(xié)方差陣計(jì)算過(guò)程中搜索共用點(diǎn)的點(diǎn)號(hào)和權(quán)重。

4 試驗(yàn)分析

為了對(duì)上述理論進(jìn)行分析論證，試驗(yàn)選擇了機(jī)載三維激光測(cè)量（Lidar）散點(diǎn)數(shù)據(jù)，散點(diǎn)的垂直中誤差（RMSE）為0.3～0.5m，這里取0.4m作為散點(diǎn)高程精度的試驗(yàn)值。考慮到源數(shù)據(jù)密度太大，對(duì)樣區(qū)隨機(jī)抽30%左右的散點(diǎn)，得到2512個(gè)點(diǎn)作為試驗(yàn)數(shù)據(jù)（圖3（b））。在插值方法上選擇DEM地表重建過(guò)程中3種典型的插值方法：不規(guī)則三角網(wǎng)（triangular irregular network，TIN）的線性插值（linear interpolation for TIN，LIT）、反距離權(quán)（inversed distance weighted，IDW），徑向基函數(shù)插值（racial basic function，RBF），搜索點(diǎn)數(shù)為13個(gè)，IDW選取的核函數(shù)為一次反距離函數(shù)，RBF核函數(shù)選擇二次曲面函數(shù)，并建立TIN（圖3（c）），作為L(zhǎng)IT插值的源數(shù)據(jù)。

圖3 試驗(yàn)數(shù)據(jù)Fig.3 Test data

具體的試驗(yàn)過(guò)程如下：

（1）數(shù)據(jù)預(yù)處理，散點(diǎn)構(gòu)建TIN，獲得TIN數(shù)據(jù)，獲取散點(diǎn)數(shù)據(jù)的先驗(yàn)方差D，確定插值參數(shù)（搜索點(diǎn)數(shù)和插值核函數(shù)）。

（2）插值處理，逐行逐列進(jìn)行插值，每插值一次記錄一次散點(diǎn)的權(quán)重k、l，并記錄散點(diǎn)的ID號(hào)。

（3）第1次誤差傳遞，按照公式（6），由權(quán)重向量k和先驗(yàn)方差D計(jì)算DEM精度，獲得DEM精度數(shù)據(jù)。

（4）局部窗口處理，對(duì)格網(wǎng)DEM逐行逐列進(jìn)行局部窗口處理，包括：① 公共點(diǎn)搜索，搜索窗口范圍內(nèi)任意兩格網(wǎng)對(duì)在插值過(guò)程中共用已知點(diǎn)的ID號(hào)，并獲取共用點(diǎn)權(quán)重l；② 計(jì)算協(xié)方差矩陣，根據(jù)公式7，由k、l計(jì)算局部窗口9格網(wǎng)的協(xié)方差矩陣Dz。

（5）第2次誤差傳遞，按照公式（10），計(jì)算坡度精度，獲得坡度噪聲精度數(shù)據(jù)，并進(jìn)行統(tǒng)計(jì)分析。具體計(jì)算方法的流程圖如圖4所示。

圖4 坡度精度計(jì)算過(guò)程Fig.4 Flow chart of DEM accuracy calculation

試驗(yàn)按照上述步驟分別計(jì)算3種插值模型下的4種差分算法的坡度精度，并考慮顧及和不顧及噪聲誤差的鄰域空間自相關(guān)性兩種方式。經(jīng)過(guò)計(jì)算，獲得了24幅坡度精度數(shù)據(jù)，如表2所示，表格橫向表示顧及和不顧及鄰域格網(wǎng)的誤差相關(guān)性下4種差分算法，縱向表示不同的插值模型。

對(duì)表中24幅坡度中誤差數(shù)據(jù)進(jìn)行統(tǒng)計(jì)分析，獲取每種插值方法和差分方式下坡度噪聲精度的最小值（min）、最大值（max）、均值（mean），結(jié)果如表3所示。

考察表1中坡度中誤差的均值（圖5），發(fā)現(xiàn)顧及DEM誤差的空間自相關(guān)性能明顯提高坡度計(jì)算模型的精度，而不同的插值方法和差分算法也不同程度地影響坡度噪聲精度，具體表現(xiàn)在：在格網(wǎng)DEM插值方法上，IDW方法生成的格網(wǎng)DEM在提取坡度時(shí)具有較高的噪聲精度，其次依次是RBF和LIT，通過(guò)對(duì)計(jì)算程序的跟蹤分析，造成這種現(xiàn)象的原因是IDW采用歸一化權(quán)重，而其他兩種插值方式存在權(quán)重大于1或負(fù)數(shù)的情況；在差分方式上，基于二階差分的坡度噪聲精度最高，其次依次是三階不帶權(quán)差分、三階反距離權(quán)差分和三階反距離平方權(quán)差分。這與文獻(xiàn)［2］的結(jié)論有所不同，其根本原因在于：

表2 格網(wǎng)DEM坡度中誤差Tab.2 RMSE of slope extracted from grid DEM

表3 各試驗(yàn)樣區(qū)坡度中誤差統(tǒng)計(jì)Tab.3 Statistics of root mean square error for slope （°）

圖5 坡度噪聲精度統(tǒng)計(jì)圖Fig.5 Statistics chart of slope noise accuracy

（1）協(xié)方差計(jì)算模型不同。DEM誤差自相關(guān)原因在于內(nèi)插時(shí)公共點(diǎn)的存在，與內(nèi)插鄰域大小、公共點(diǎn)數(shù)等有關(guān)。在文獻(xiàn)［2］中，協(xié)方差模型采用指數(shù)模型和線性模型，這是一種經(jīng)驗(yàn)?zāi)Ｐ停淝疤崾荄EM誤差自相關(guān)應(yīng)符合指數(shù)或者線性變化趨勢(shì)。而本文采用的是理論模型，其協(xié)方差值是通過(guò)內(nèi)插計(jì)算模型和原始數(shù)據(jù)誤差計(jì)算獲得的。通過(guò)對(duì)試驗(yàn)中的協(xié)方差數(shù)據(jù)進(jìn)一步分析，發(fā)現(xiàn)DEM噪聲誤差空間自相關(guān)分布并不符合線性或者指數(shù)模型（如圖6所示）。

圖6 DEM 3×3窗口噪聲誤差協(xié)方差分布Fig.6 Spatial distribution of noise error covariance at local windows 3×3

（2）DEM噪聲誤差自相關(guān)模型是各向異性的。文獻(xiàn)［2］中的指數(shù)模型或者線性模型，其關(guān)鍵變量是距離而與方向無(wú)關(guān)，即DEM自相關(guān)程度是各項(xiàng)同性的，本文通過(guò)計(jì)算獲取的協(xié)方差具有各項(xiàng)異性（如圖6所示），即不同方向的自相關(guān)性變異程度是不同的。

（3）在坡度精度評(píng)價(jià)上，坡度計(jì)算精度與參與計(jì)算格網(wǎng)點(diǎn)數(shù)有關(guān)。二階差分采用周圍4個(gè)格網(wǎng)點(diǎn)，而三階差分系列采用的是周圍8個(gè)格網(wǎng)點(diǎn)（如表1所示），前者相關(guān)系數(shù)為＝10，而后者為＝36，根據(jù)公式（10），參與計(jì)算的格網(wǎng)點(diǎn)的數(shù)量會(huì)影響坡度噪聲精度。

（4）通過(guò)研究，揭示了DEM噪聲誤差不是均一的，DEM噪聲誤差描述需要從全局考慮，即給出DEM誤差的空間分布，而不應(yīng)該只給出中誤差這一全局性精度指標(biāo)。

5 結(jié) 論

DEM建立和地形分析過(guò)程中鄰域操作，使得DEM誤差在空間上具有特定自相關(guān)性，顧及DEM噪聲誤差的空間自相關(guān)性這不但提高DEM數(shù)據(jù)本身的質(zhì)量描述，同時(shí)也改善了基于DEM的地形分析的精度。DEM噪聲誤差的空間自相關(guān)性源自于DEM插值過(guò)程中共用的散點(diǎn)，在DEM數(shù)字地形分析特別是鄰域分析不能忽略誤差的空間自相關(guān)性。本文通過(guò)理論推導(dǎo)和試驗(yàn)分析，以坡度為研究對(duì)象，在顧及和不顧及DEM噪聲誤差空間自相關(guān)性兩種情況下對(duì)4種常見(jiàn)坡度計(jì)算模型的噪聲精度進(jìn)行分析，本文分析方法的主要特點(diǎn)有：① 顧及了DEM噪聲誤差的空間自相關(guān)性；②采用DEM噪聲誤差理論空間自相關(guān)模型。

經(jīng)過(guò)理論推導(dǎo)和試驗(yàn)分析，有如下初步結(jié)論：①顧及DEM誤差空間自相關(guān)性能夠從整體上減少坡度中誤差，在數(shù)字地形分析精度分析中不能忽略噪聲誤差的空間自相關(guān)性，以提高精度預(yù)測(cè)準(zhǔn)確性；② 坡度計(jì)算模型的中誤差與格網(wǎng)DEM插值方法和坡度計(jì)算過(guò)程中的差分算法有關(guān)，在差分方式上，二階差分的精度最高，其次依次是三階不帶權(quán)差分、三階反距離權(quán)差分和三階反距離平方權(quán)差分，在格網(wǎng)DEM插值方法上，IDW方法生成的格網(wǎng)DEM在提取坡度時(shí)具有較高的坡度，其次依次是RBF和LIT。

以上僅探討了DEM噪聲誤差的坡度精度評(píng)價(jià)，本文進(jìn)一步的研究工作是DEM逼近誤差對(duì)地形參數(shù)的精度影響分析。

［1］ ZHOU Qiming，LIU Xuejun.Digital Terrain Analysis［M］.Beijing：Science Press，2006.（周啟鳴，劉學(xué)軍.數(shù)字地形分析［M］.北京：科學(xué)出版社，2006.）

［2］ LIU Xuejun，BIAN Lu，LU Huaxing，et al.The Accuracy Assessment on Slope Algorithms with DEM Error Spatial Autocorrelation［J］.Acta Geodaetica et Catographica Sinica，2008，37（2）：200－206.（劉學(xué)軍，卞璐，盧華興，等.顧及DEM誤差自相關(guān)的坡度計(jì)算模型精度分析［J］.測(cè)繪學(xué)報(bào)，2008，37（2）：200－206.）

［3］ STENFAN K.The Effect of DEM Raster Resolution on First Order，Second Order and Compound Terrain Derivatives［J］.Transaction in GIS，2004，8（1）：83－111.

［4］ THIERRY L，STEPHANE J，CHRISTOPHE F.Very High Resolution Digital Elevation Models：Do They Improve Models of Plant Species Distribution？ ［J］.Ecological Modeling，2006，198：139－153.

［5］ OKSANEN J，SARJAKOSHI T.Error Propagation Analysis of DEM－based Surface Derivatives ［J］.Computers and Geosciences，2005，31（2）：1015－1027.

［6］ FLORINSKY I.Accuracy of Local Topographic Variables Derived from Digital Elevation Models［J］.International Journal of Geographical Information Sciences，1998，12（1）：47－61.

［7］ JONES K.A Comparison of Algorithms Used to Compute Hill Slope as a Property of the DEM ［J］.Computer and Geosciences，1998，24（4）：315－323.

［8］ LIU Xuejun.On the Accuracy of the Algorithm for Interpretation Grid－based Digital Elevation Model［D］.Wuhan：Wuhan University，2002.（劉學(xué)軍.基于規(guī)則格網(wǎng)數(shù)字高程模型解譯算法誤差分析與評(píng)價(jià)［D］.武漢：武漢大學(xué)，2002.）

［9］ GOODECHILD M.The Spatial Data Infrastructure of Environmental Modeling［M］.London：Longman，1991.

［10］ LEE J.Digital Elevation Models：Issues of Data Accuracy and Application［EB／OL］.［2010－08－12］.http：∥www.esri.com／library／userconf／proc96.

［11］ HEUVELINK B.Error Propagation in Environmental Modeling with GIS［M］.Boca Raton：CRC Press，1998.

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

［13］ LU Huaxing.Research on DEM Error Model［D］.Nanjing：Nanjing Normal University，2008.（盧華興.DEM誤差模型研究［D］.南京：南京師范大學(xué)，2008.）

［14］ TAO Benzao，YU Zongchou.Basic Principle of Survey Adjustment［M］.Beijing：Publishing House of Surveying and Mapping，1996.（陶本藻，於宗儔.測(cè)量平差基礎(chǔ)［M］.北京：測(cè)繪出版社，1996.）

［15］ FLEMING M，HOFFER R.Machine Processing of Landsat MSS Data and DMA Topographic Data for Forest Cover Type Mapping［R］.West Lafayette：Purdue University Laboratory for Applications of Remote Sensing，1979.

［16］ SHARPNACK D，AKIN G.An Algorithm for Computing Slope and Aspect from Elevations［J］.Photogrammetric Survey，1969，35：247－248.

［17］ HORN B.Hill Shading and the Reflectance Map［C］∥Proceedings of IEEE，1981，69（1）：14－47.

［18］ UNWIN.Introductory Spatial Analysis［M］.London and New York：Methuen，1981.