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

盧華興，劉學軍，王永君，田梓文

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

1 引 言

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

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

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

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

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

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

圖1 散點插值Fig.1 Elevation points interpolation

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

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

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

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

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

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

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

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

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

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

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

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

4 試驗分析

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

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

具體的試驗過程如下：

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

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

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

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

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

圖4 坡度精度計算過程Fig.4 Flow chart of DEM accuracy calculation

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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

［1］ ZHOU Qiming，LIU Xuejun.Digital Terrain Analysis［M］.Beijing：Science Press，2006.（周啟鳴，劉學軍.數(shù)字地形分析［M］.北京：科學出版社，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.（劉學軍，卞璐，盧華興，等.顧及DEM誤差自相關的坡度計算模型精度分析［J］.測繪學報，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.（劉學軍.基于規(guī)則格網(wǎng)數(shù)字高程模型解譯算法誤差分析與評價［D］.武漢：武漢大學，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.（李新，程國棟，盧玲.空間插值方法比較［J］.地球科學進展，2000，15（3）：260－265.）

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

［14］ TAO Benzao，YU Zongchou.Basic Principle of Survey Adjustment［M］.Beijing：Publishing House of Surveying and Mapping，1996.（陶本藻，於宗儔.測量平差基礎［M］.北京：測繪出版社，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.