徐志敏,王利雙
(1.長江勘測規劃設計研究有限責任公司,武漢 430010;2.長江空間信息技術工程有限公司(武漢),武漢 430010;3.湖北省水利信息感知與大數據工程技術研究中心,武漢 430010;4.國家知識產權局專利局專利審查協作湖北中心,武漢 430205)
數字高程模型(digital elevation model,DEM)是國家基礎地理數據的重要組成部分,也是進行三維空間分析和應用的基礎。對DEM進行全面的精度分析和質量評價是保證其正確應用的基本前提[1]。目前主要的DEM精度評價模型很多,其中高程中誤差模型[2]已經成為國家測繪行業的采納標準(CH/T 1015.2—2007、CH/T 9009.2—2010)。然而研究表明,高程中誤差只能評價DEM采樣點的高程誤差,但無法評價DEM模擬表面對真實地表的逼近誤差,即使高程中誤差為0,DEM模擬表面也不一定能真實表達地表起伏。也正因此,學界提出了多種DEM誤差的評價模型,旨在對DEM誤差進行更加科學的評價,包括重構等高線模型[3]、地形信息熵模型[4]、Strahler積分模型[5]、填挖方模型[6]與地形描述精度模型[7]等。
其中,DEM地形描述精度模型由于其科學性和實用性,引起了學界的關注和研究。湯國安等[8]最早提出了地形描述誤差的概念和提取方法,并對地形描述誤差的空間結構和定量計算模型做了研究。張勇等[9]以黃土丘陵溝壑區為例,采用統計和比較分析的方法,研究了不同空間尺度下的DEM地形描述誤差的數學轉換模型。劉春等[10]基于地形復雜因子概念,利用線性回歸分析方法探討了柵格DEM地形描述精度與地形復雜因子和平均高程之間的線性關系。賈敦新等[11]以模擬實驗獲得的高精度DEM為數據源,分析坡度計算精度與DEM數據誤差和地形描述誤差的耦合關系。齊曉飛等[12]提出了一種基于散點圖矩陣建立多元回歸模型的方法,并依此構建了DEM地形描述誤差與分辨率、相對高差、凹凸系數的關系公式。王春等[13-14]借鑒誤差橢圓和ε誤差帶建模的思想,提出了地形描述誤差的中誤差場和極值場,并研究了這些場的基本特性。可見,DEM地形描述精度的定量模型一直是學界關注的熱點。
本文對地形描述誤差(以下簡稱Et)與DEM空間分辨率和平均剖面曲率之間的定量關系進行了重新研究。從地形描述精度的概念出發,分析了現有的地形描述精度定量模型的缺陷,提出了優化的冪函數定量模型,并通過實驗對比了本文所提出定量模型的擬合精度和泛化精度的優勢。
湯國安等首次提出了地形描述誤差的概念和計算方法,并以地形描述誤差概念和計算方法為基礎,選取了六種不同地貌類型(平原、低丘、丘陵、中山、高山以及混合類型)的DEM作為實驗樣區,分別計算了各樣區在不同DEM空間分辨率下的地形描述誤差的均方差。首先,對各樣區的DEM地形描述誤差均方差和DEM空間分辨率之間的關系采用y=ax+b線性模型進行了線性回歸;然后,對方程系數b與對應樣區的平均剖面曲率V分別再次進行線性回歸;最后,得到DEM地形描述誤差均方差與空間分辨率R和平均剖面曲率V的定量計算公式,如式(1)所示。
RMSEt=(0.006 3V+0.006 6)R-0.022V+0.241 5
(1)
對此,王光霞等[15]指出:上述研究中系數b與平均剖面曲率的線性相關性并不好,利用直線擬合勢必產生較大的誤差。因此,對系數b與平均剖面曲率采用了二次函數模型進行了擬合,將參數a與平均剖面曲率的相關系數從0.989 9提高到0.998 9,將系數b與平均剖面曲率的相關系數從0.842 5提高到0.988 2,最終得到擬合精度更高的定量模型,計算方法如式(2)所示。該方法得到的定量模型與式(1)相比較具有更好的擬合效果,將在地勢較平緩的地貌類型(如平原、低丘和丘陵地區)上的擬合中誤差提高了一個數量級。
RMSEt=(0.006 1V+0.002 7R)+0.001V2-0.064 9V+0.569 5
(2)
事實上,根據湯國安等對地形描述誤差的定義,DEM地形描述誤差Et是假定DEM高程采樣誤差為0的條件下,模擬地面與實際地面之間的差異。
如圖1所示,弧面AB代表實際地形表面,面AB代表離散采樣點所表達的模擬表面,AB點在水平面上的投影點之間的距離為DEM的分辨率。在點C處的DEM模擬地面高程與實際地面高程差為C點的DEM地形描述誤差Etc。實際計算中,實際地面高程值一般采用4鄰域的DEM柵格點高程的平均值進行替代。

圖1 地形描述誤差示意圖
那么,根據DEM地形描述誤差的定義和計算方法,在理想情況下,當DEM平均剖面曲率趨向于0時,DEM模擬地面為近似平面,任意柵格點四鄰域的高程平均值與其本身高程差為0,地形描述誤差為0;而當DEM分辨率趨向于0時,DEM模擬地面能夠完全表達實際地面,地形描述誤差也為0。因此,對于地形描述誤差Et與DEM空間分辨率R和平均剖面曲率V之間的定量模型RMSEt=f(R,V)應該滿足式(3)。
(3)
然而,式(1)和式(2)均為帶常數項的多項式模型,在DEM分辨率R和平均剖面曲率V分別趨向于0時的極限均為帶常數的不定式,因而均不滿足式(3)。
圖2顯示了當平均剖面曲率為0,即地形絕對平坦的情況下,式(1)、式(2)對應模型的地形描述精度與DEM分辨率之間的函數關系。根據式(3),當地形絕對平坦時,DEM地形描述誤差與DEM分辨率不相關,且恒為0。但式(1)、式(2)在平均剖面曲率為0時,地形描述誤差與DEM分辨率為線性關系,隨著地形分辨率的提升,地形描述誤差逐漸下降,且恒為某一大于0的值。

圖2 V為0時式(1)、式(2)中RMS Et與R的關系
這顯然于地形描述誤差的固有性質不符。因此,式(1)、式(2)對應的模型均欠缺一定的科學性。事實上,文獻[9]和文獻[11]研究所得的定量模型也存在類似問題。這些模型在對DEM誤差進行定量建模的時候都沒有顧及到DEM誤差與DEM地形特征之間的內在關系。因此,地形描述誤差與DEM空間分辨率和平均剖面曲率之間的定量模型需要基于式(3)這一前提進行重新設計和研究。
考慮到上述DEM地形描述誤差的重要性質,本文基于式(3)的前提,設計了基于冪函數的地形描述誤差Et與DEM空間分辨率R和平均剖面曲率V的定量關系模型,如式(4)所示。
RMSEt=aRbVc
(4)
式中:a、b、c為模型系數。
將(R,V,RMSEt)看作是三維空間中的離散點,將式(4)看作空間曲面方程。模型求解過程即可看作三維離散點的曲面擬合。由于式(4)與自變量的關系是非線性的,不能用線性最小二乘法求多元函數極值的方法獲得參數估計值,因此需要采取復雜的優化算法進行求解。常用的優化求解算法包括最速下降算法、高斯-牛頓法和列文伯格-馬夸爾特算法。本文采用基于高斯-牛頓迭代的非線性最小二乘法,對離散三維點進行曲面擬合。將三維離散點和三維曲面的擬合效果在MATLAB中顯示,結果如圖3所示。

圖3 曲面擬合結果
擬合得到的模型系數a、b、c的值分別為:0.005 713、1.056和0.943 1。因此,基于冪函數的地形描述誤差與DEM空間分辨率和平均剖面曲率之間的定量模型為式(5)所示。
RMSEt=0.005 713R1.056V0.943 1
(5)
進一步地,計算并比較式(1)、式(2)和式(5)對樣本數據的擬合優度相關指標,包括誤差平方和SSE、確定系數R-square以及標準差RMSE,結果如表1所示。可見,本文所提出定量模型的擬合效果很好,擬合精度與式(2)基本相當,且各項指標明顯優于式(1)。

表1 擬合優度
對式(5)進行改寫,可得DEM分辨率是地形描述誤差和平均剖面曲率的函數,即式(6),以及平均剖面曲率是地形描述誤差和DEM分辨率的函數,即式(7)。DEM地形描述誤差定量模型在生產實踐中的實用價值在于:基于給定的DEM的限差,式(5)可以直接根據分辨率和剖面曲率檢驗DEM是否符合精度要求;式(6)可以推算一定地貌類型范圍內符合精度要求的最低DEM分辨率;式(7)可以推算一定分辨率的DEM產品中符合精度要求的地貌類型區域。
(6)
(7)
擬合精度評價的是定量模型公式對已知的實驗樣本數據的逼近程度,但是模型的擬合精度并不能保證模型在生產實踐中的實用性,還需要評價模型對未知樣本數據的估算效果,即泛化精度。
本文分別選取了平原、丘陵和中山地帶的0.5 m高分辨率DEM進行了實驗。數據來源于美國圣地亞哥超級計算機中心的開放地形數據OpenTopography網站,這個網站開放了高分辨率地形數據及處理工具(http://opentopo.sdsc.edu/lidar?format=sd)。數據區域為美國加尼福尼亞州Napa Valley,地理范圍為1×104m2,高程分布在海拔-1~634 m之間,從西南到東北方向依次分布著平原丘陵和山地,包含了比較豐富的地形種類。本文選取了該DEM范圍內的平原、丘陵和中山的三塊地形作為實驗樣區,源數據和實驗樣區的分布如圖4所示。

圖4 源數據與實驗樣區
采用文獻[7]關于地形描述誤差的計算方法,基于0.5 m分辨率的原始DEM,將地形描述精度計算窗口依次按照3×3、5×5、7×7擴大,分別計算得到3種地形在分辨率為1 m、2 m、4 m三種情況下的地形描述誤差,總計九組數據。同時,采用式(1)、式(2)、式(5)分別計算九組不同DEM分辨率和平均剖面曲率情況下的地形描述誤差估算值,結果如表2所示。將地形描述精度的計算值與式(1)、式(2)和式(5)的估算值進行統計比較。可以看出,式(5)的定量模型對于高分辨率DEM的地形描述誤差估算精度都在83.3%以上,同時隨著DEM分辨率的提高,模型的估算精度也呈現出增加的趨勢,可見本文的定量模型對于高分辨率DEM具有很好的泛化性能。而對于式(1)和式(2),其定量模型的估算精度不太穩定,式(1)的估算精度在45.5%和84.7%之間波動,式(2)的估算精度在49.4%和98.9%之間波動,甚至隨著DEM分辨率的提高,一些樣區地形描述誤差的均方差估算值產生了不符合邏輯的負值。

表2 泛化精度對比結果
因此,實驗結果印證了本文指出的式(1)和式(2)的問題:該定量模型在設計的過程中只追求對現有樣本數據的擬合精度,沒有考慮到地形描述精度的固有性質并根據對應自變量和因變量的固有數學關系來選擇合適的擬合模型,導致模型對原有樣本數據的擬合精度較好,但是對于更高分辨率DEM樣本的泛化性能不理想。
本文對地形描述精度與DEM分辨率和地形平均剖面曲率之間的定量關系模型進行了研究。首先,從DEM地形描述精度的概念出發,分析了現有的帶常數項的多項式模型的缺陷;然后,提出了基于冪函數的定量模型,計算并對比了模型對樣本數據擬合精度;最后,基于 4 m、2 m和1 m的高分辨率DEM,比較了模型在新的樣區的泛化精度。分析結果表明,本文的冪函數模型在保證良好的擬合精度的同時,對高分辨率的DEM具有較好的泛化性能。
隨著測繪遙感技術的發展和計算機軟硬件性能的提升,高分辨率DEM生產和應用由于能夠滿足用戶對三維地形精度和體驗的要求,得到越來越廣泛的應用。在這種趨勢下,本文的定量模型對于在生產實踐中,根據限差要求推算適宜的DEM分辨率具有更好的指導意義和實用價值。
然而,還需注意到,雖然本文基于冪函數的定量模型具有較好的擬合精度與泛化性能,但是滿足“當DEM空間分辨率和平均剖面曲率各自趨向于0,地形描述誤差亦趨向于0這一客觀事實”的冪函數模型形式卻不限于本文這一種。事實上,將本文的冪函數定量模型進行線性組合,所得到新模型都同樣滿足上述條件。因此,要找到更加精確、實用性更好的模型仍然需要進一步的探索。