地形特征約束下的失真DEM修復方法

周 波，徐俊波

(1. 安徽省智慧城市與地理國情監測重點實驗室，安徽 合肥 230031； 2. 合肥工業大學，安徽 宣城 242000)

地形特征約束下的失真DEM修復方法

周 波1,2，徐俊波1,2

(1. 安徽省智慧城市與地理國情監測重點實驗室，安徽 合肥 230031； 2. 合肥工業大學，安徽 宣城 242000)

數字高程模型(DEM)在描述地形時會存在一定程度的失真現象，對失真DEM進行修復可以有效地提高其精度及保真度。本文針對DEM在描述規則地形時存在的失真現象，提出了用地形特征約束失真DEM形態結構的失真修復方法，以規則地形中的平地結構地形為例，進行了失真修復試驗。首先根據地形特征確定平地地形的DEM應為平面結構，再采用最小二乘法(LS)及隨機抽樣一致性算法(RANSAC)對其平面方程進行參數估計，進而得到修復的DEM數據。試驗結果表明，本文所提的修復方法有效可行，修復后的DEM數據不僅提高了精度，還能充分地反映原地形的地形特征。兩種修復算法的修復能力相接近，當存在大量格網點時，都能達到很好的修復效果，但RANSAC更能適應高程異常的情況，魯棒性更好。

數字高程模型；地形特征；最小二乘法；隨機抽樣一致性算法；失真修復

數字高程模型(digital elevation model，DEM)從1958年提出至今，已經成為地形的數字描述和空間模擬的數據基礎，為數字地球的發展提供了數據支持[1]。目前獲取DEM方法多樣，可以從航空影像、遙感圖像及現有地形圖中得到。由于DEM的數據來源、生產方式、內插公式等存在差異，以及操作流程中可能存在的不規范操作，會導致DEM與地理對象的地形特征不符現象，如面狀水域、人工平地等這類屬于平面結構地形的地理對象出現“平地不平”的失真現象[2-5]。

為了提高DEM精度，眾學者作了很多此類研究，提出過不同的減少誤差的方法[6-13]，但大多是從DEM生產過程中提高DEM精度，對已經失真的DEM卻沒有一個很好的修復策略。本文針對諸如“平地不平”這類規則地形的失真現象，結合規則地形的地形特征與失真DEM的數據分布特點，提出用地形特征約束DEM形態結構，并用參數估計算法求得符合地形起伏狀態的結構表達式，進而修復失真數據的方法，以求得到更高精度的DEM數據。

1 問題解析

失真的DEM通常不能充分且確切地反映原地形的整體特征，即地理對象的地形特征，包括地形的高低起伏狀態、陡緩程度等。如圖1所示，所選樣區為一大型水庫，其地形特征為平面結構，水域高程應相等或線性變化，但從圖1失真高程數據中可以看出，樣區高程并不相等且無規律變化，該處數據明顯失真。

圖1 樣區地形圖與失真DEM

對此類地形失真現象進行總結，具有以下特點：①失真區域為規則地形，地形特征可用數學公式表示，如圖1中所選的水庫樣區所屬的平面結構地形，在數學中可以用三維平面方程表示；②失真數據雖然存在誤差，但數據整體依然在真值附近的可信范圍內波動，可能存在極少數異常值，但所占比例可以忽略不計。圖1中DEM高程數據雖然存在誤差，但與等高線標識的高程值進行比較，其依然在合理范圍內波動。

結合以上兩個特點，對規格地形失真DEM進行修復的思路如圖2所示，規則地形特征所對應的數學公式描述的是DEM真值應具有的形態結構，失真格網點分布在真值附近的可信范圍內，對其進行修復就是要使其整體符合地形特征且圖中標識的誤差盡可能小。運用參數估計算法(修復算法)求取失真DEM的高程真值近似表達式，再根據此表達式求得修復格網點，進而可得更高精度的DEM數據。

圖2 失真DEM修復示意圖

地形特征約束了失真DEM整體形態結構，修復算法則可求出此結構的具體表達式，但不同修復算法求取的表達式不盡相同，修復后的數據精度也不同。為了得到精度更高，更能反映地形特征的DEM數據，本文選取最小二乘法(least square method，LS)與隨機抽樣一致性算法(random sample consensus，RANSAC)進行修復試驗，并對修復結果進行比較。在處理失真問題時，LS是運用最廣泛的參數估計算法之一[14-16]，其原理為利用觀測向量的殘差平方和最小來求取模型參數，進而求得修復值。但LS雖然考慮了誤差因素，卻不能規避異常數據對參數的影響(DEM中常表現在部分格網點的高程突變)，為此，本文選取了能夠提高數據精度、減小異常數據對模型參數干擾的RANSAC與LS進行對比試驗[17-19]，并對修復結果進行了質量檢測，比較出最優的修復算法。

2 修復方法

2.1 修復流程

由于樣區高程真值不可得，且抽樣具有隨機性，小樣本不具有普適性，為了方便對修復后的DEM數據進行質量檢測，以及避免采樣過少可能導致的結論偏差，本文首先用模擬數據代替采樣數據進行修復試驗，找到更好的修復算法，再用采樣數據對結論進行驗證。

對失真DEM數據進行數學模擬時，針對平地地形，在三維平面上模擬了一組平面規則格網點，將此三維平面視為DEM真值平面，再為這組格網點高程添加誤差向量，將未添加誤差的格網數據視為DEM真值數據，含有誤差的視為失真DEM數據，采用LS與RANSAC對此模擬數據進行修復。運用兩種算法進行修復試驗，可得兩個不同修復平面方程，根據這兩個平面方程可以得到不同的修復DEM數據。

為了對這兩種算法的修復效果進行評估，分別求出所得的兩組平面與DEM真值平面之間的夾角及修復后所得DEM的精度(本文通過修復前后數據高程中誤差大小進行判斷)，根據夾角的大小及每種算法修復后所得數據精度，比較兩種算法的優劣性，找出最優修復算法。

2.2 修復算法

2.2.1 采用最小二乘法進行修復

最小二乘法成立的前提是自變量為真值矢量[20]，對于DEM數據而言，當其平面向量的平面位置呈規則格網排列時，常將其平面坐標省略，只保留含有一維向量序列的高程數據(格網DEM)[1]。因此，用LS修復失真DEM數據時可以將其平面位置向量視為真值矢量，只需考慮高程Z方向上的誤差情況，解算相應參數。

假設平面樣區對應的三維平面方程為[20]

Z=aX+bY+c

(1)

式中，a、b、c為未知參數。

對失真DEM進行處理，得到n個格網點，可表示為

Vi={(xi,yi,zi)i=1,2,…,n}

(2)

根據最小二乘原理，修復失真DEM應使S值最小，即

(3)

需滿足式(4)，即使式(5)成立。

(4)

(5)

整理可得

(6)

解方程組得到參數a、b、c。

2.2.2 采用隨機抽樣一致性算法進行修復

隨機抽樣一致性算法由Fishchler和Bolles于1981年提出，是用于從觀測數據中估計出數學模型參數的迭代算法。RANSAC的思想是假定觀測數據中含有內點(inliers，對估計模型參數有用的點)，也包含異常數據外點(outliers，偏離正常范圍很遠，不符合數學模型的數據)，并且這組數據受噪聲影響，RANSAC認為只要存在inliers數據點，就能估計出數學模型。其不同于LS等參數估計方法利用所有觀測數據估計模型參數后再舍去誤差較大的點，而是利用滿足條件的盡可能少的數據并用一致性數據集將其逐漸擴大，摒棄異常數據。

用RANSAC對平面結構地形的失真DEM數據進行修復的基本步驟如下[6]：

(1) 在失真格網點中先任取3個格網點構成一個平面，并遍歷其他格網點，計算每個格網點到此平面的距離。

(2) 設定閾值，當小于閾值時，認為此格網點為內點，所有內點構成一個一致集(consensus set)。

(3) 重復以上過程，重新隨機抽取新的一致集。在完成一定的抽樣次數后，若未找到一致集則算法失敗，否則在所有一致集中選取最大一致集。

(4) 根據判斷內外點并采用LS對最大一致集 中所含的內點進行參數估計得到模型 ，算法結束。

3 試驗與分析

3.1 用模擬數據進行試驗

在一般情況下，若隨機采樣點超過30個，可以認為誤差服從正態分布[1]，樣區格網點數量遠大于30個，故可認為樣區DEM誤差服從正態分布，對樣區DEM進行數學模擬時可進行如下操作。

(1) 在平面Z=X+Y上模擬出100個規則格網點({(xi,yi,zi) i=1,2,…,500}，如圖3(a)所示)，將此平面DEM視為真值平面。

(2) 為這組格網點的Z值添加一組服從N(0,1)標準正態分布的誤差向量(如圖3(b)所示)，將未添加誤差的格網點視為DEM真值數據，添加誤差的格網點視為失真DEM數據。

(3) 重復步驟(1)、(2)，再添加3組不同誤差向量，且格網點數量分別為500、1000、1500，形成4組數據，用LS和RANSAC對其修復。

為了比較LS與RANSAC在處理格網點出現明顯異常情況下的處理能力，重復上述3個步驟重新模擬4組數據并在每組數據中添加3%的異常點(如圖3(c)所示，圓圈標注的即為異常點)，用此含有異常點的模擬數據再次進行試驗。

給高程Z(Z=[z1,z2,…,zn]T)添加服從N(0,1)正態分布的誤差向量后，數據整體呈現出有規律狀態，但是數據不處于同一平面，部分數值見表1，可以模擬出DEM格網數據的分布情況。

圖3 模擬格網點

XY原Z值添加誤差后Z值1121.332317080487221234.753613980245871346.467044741743381454.267128688964841564.085975998275272133.495818268900712243.295755411228772354.922712760197142465.900135725502732577.322427621851123142.691606423194513255.237142922925953365.582581650806693475.510174415820733587.858084950512644153.954015534123144268.864112237669824377.05020539890223

3.2 試驗結果

3.2.1 不含異常的模擬數據試驗結果分析

采用LS進行修復，試驗結果如圖4、圖5所示，可以看出4組試驗中采用LS修復后所得平面與真值平面之間的夾角(如圖4所示)及高程中誤差(如圖5所示)都比RANSAC小?？傮w上看，兩種算法都提高了數據精度，與真值平面的夾角和高程中誤差都會隨著格網點的增加而減小。

圖4 不含異常點數據修復前后的平面間夾角

圖5 不含異常點數據修復前后的高程中誤差

3.2.2 含有異常的模擬數據試驗結果分析

用含有異常的模擬數據進行試驗，試驗結果如圖6、圖7所示，可以看出每組試驗中采用RANSAC進行修復，所得平面與真值平面的夾角(如圖6所示)和高程中誤差(如圖7所示)都小于LS。總體上格網點數目越多修復效果越好，且修復能力越來越接近。

圖6 含有異常點數據修復前后的平面夾角

圖7 含有異常點數據修復前后的高程中誤差

3.3 結果分析

試驗采用的模擬數據一方面能夠避免抽樣的隨機性，可以模擬出不同誤差的失真數據，另一方面方便對修復后所得數據進行質量檢測，試驗結論總結如下：

(1) 當不存在異常點時，用LS修復所得數據相比RANSAC更加貼近真值數據，但隨著數據量增大，兩者的修復效果越來越接近。

(2) 當存在明顯異常點時，RANSAC的數據修復能力要比LS強，所得平面更接近原平面，所得數據精度也更高。這也說明RANSAC能夠更好地適應高程異常情況，比LS更加穩定，魯棒性更好。

(3) 兩種算法的修復效果都隨著格網點的增加越來越好，說明當存在大量格網點時兩種算法都能達到很好的修復效果。

3.4 用樣區規則格網DEM數據進行檢驗

為了驗證上述結論是否適用于實際采樣數據，本文從1∶10 000等高線地形圖中隨機選取了4個水庫樣區，并獲取它們的DEM數據，由于數據的誤差情況未知，在試驗前對LS與RANSAC的修復效果保有未知性。樣區DEM格網點數目是遞進的，分別為691、1210、2381、3288，由于格網點所屬真實平面不可求，樣區為水庫，其面積遠小于湖泊、河流，可視為靜水狀態，故本文將其視為水平面，其高程為等高線標注的值。求出修復后的數據所處平面與水平面的夾角(如圖8所示)，以及修正前后DEM高程中誤差(如圖9所示)結果。

圖8 采樣數據修復后所得平面與水平面夾角

圖9 采樣數據修復前后高程中誤差

用實際采樣中的失真數據進行結論驗證試驗后，試驗結果與采用模擬試驗時基本保持一致，具體可總結為如下幾個方面：

(1) 當存在大量格網點時，采用LS與RANSAC進行失真修復所得平面與水平面之間的夾角精度可達0.001 arc(如區域1—3)甚至0.000 1 arc(如區域4)，兩種算法都能降低高程中誤差的大小，修復所得DEM的高程中誤差都遠遠小于國家標準[22]。

(2) 由于樣區數目較少，且數據中不存在異常數據，從試驗結果上看LS的修復效果要稍優于RANSAC。但當采樣較多時，存在異常數據的樣區不可避免，當存在異常數據時，根據采用模擬數據試驗時RANSAC具有更好的魯棒性，本文推薦采用RANSAC進行失真修復。

綜合來看，兩種算法的修復效果會隨著格網點增多而變好，且修復能力會越來越接近，試驗結果基本與用模擬數據時保持一致。兩種算法都能夠達到很好的修復效果，但與LS相比，RANSAC更能夠適應異常點的干擾，具有很好的魯棒性。

4 結 語

本文用地形特征約束失真DEM的形態結構，并采用了兩種算法估算其真值表達式，進而展開修復工作。結果表明，采用LS與RANSAC都能達到很好的修復效果，且隨著格網點數目的增多，修復效果都越來越好，但RANSAC比LS更加穩定，魯棒性更好。

本文只是針對規則地形中的平面結構地形進行了失真修復試驗，但文中的失真修復方法具有普適性，所提出的兩種算法也提供了選擇思路。其他規則地形如半球面狀、錐面狀、拋物面狀地形等，都可以進行失真修復，只要確定了地形特征模型再選擇相應算法，即可運用本文所提的相應技術方法達到修復效果。

[1] 李志林，朱慶.數字高程模型[M].武漢：武漢大學出版社，2000.

[2] 周波.語義約束的地形建模方法研究[D].南京：南京師范大學，2012.

[3] 盧華興.DEM誤差模型研究[D].南京：南京師范大學，2008.

[4] 吳艷蘭，胡海，胡鵬，等.數字高程模型誤差及其評價的問題綜述[J].武漢大學學報(信息科學版)，2011，36(5)：568-574.

[5] AGUILAR F J，MILLS J P，DELGADO J，et al. Modeling Vertical Error in LiDAR-derived Digital Elevation Models [J]. ISPRS Journal or Photogrammetry and Remote Sensing，2010，65(1)：103-110.

[6] 胡鵬，吳艷蘭，胡海.數字高程模型精度評定的基本理論[J].地球信息科學，2003，5(3)：64-70.

[7] 胡鵬，吳艷蘭，胡海.再論DEM精度評定的基本理論問題[J].地球信息科學，2005，7(3)：28-33.

[8] HU Peng，LIU Xiaohang，HU Hai. Accuracy Assessment of Digital Elevation Models Based on Approximation Theory [J]. Photogrammetric Engineering and Remote Sensing, 2009, 75(1)：49-56.

[9] GUAN Yunlan，CHENG Xiaojun,SHI Guigang.A Robust Method for Fitting a Plane to Point Cloud [J]. Journal of Tongji University(Natural Science)，2008，36(7)：981-984.

[10] 趙爭. 地形復雜區域InSAR高精度DEM提取方法[J]. 測繪學報, 2016, 45(11)：1385.

[11] 張錦明, 游雄, 萬剛. DEM插值參數優選的試驗研究[J]. 測繪學報, 2014(2):178-185.

[12] 王琴, 陳蜜, 劉書軍,等. 利用升降軌道SAR數據獲取DEM的試驗研究[J]. 測繪通報, 2015(6):39-43.

[13] 范強, 朱添翼, 李永化. DEM曲面分區內插方法研究[J]. 測繪通報, 2016(11)：64-66.

[14] 王峰，丘廣新，程效軍.改進的魯棒迭代最小二乘平面擬合算法[J].同濟大學學報(自然科學版)，2011，39(9)：1350-1354.

[15] 魯鐵定.總體最小二乘平差理論及其在測繪數據處理中的應用[D].武漢：武漢大學，2010.

[16] 章嘉人.語義規則約束下的DEM與二維矢量數據融合方法研究[D].南京：南京師范大學，2011.

[17] SCHNABEL R，WAHL R，KLEIN R. Efficient RANSAC for Point-cloud Shape Detection [J]. Computer Graphics Forum，2007，26(2)：214-226.

[18] 宋衛燕.RANSAC算法及其在遙感圖像處理中的應用[D].華北電力大學，2011.

[19] ZHOU Chunlin，ZHU Hehua，LI Xiaojun. Research and Application of Robust Plane Fitting Algorithm with RANSAC [J]. Computer Engineering and Applications，2011，47(7)：177-179.

[20] 陶本藻，邱衛寧.誤差理論與測量平差[M].武漢：武漢大學出版社，2012.

[21] 管興胤，李真富，宋朝暉.貝葉斯統計在測量數據擬合處理中的應用[J].核電子學與探測技術，2010，30(4)：503-505，523.

[22] 中華人民共和國測繪行業標準.基礎地理信息數字產品1∶10000 1∶50000數字線劃圖：CH/T 1011—2005[S].北京：測繪出版社，2005.

Recovery of Distorted DEM under Constraint of Terrain Feature

ZHOU Bo1,2,XU Junbo1,2

(1. Anhui Key Laboratory of Smart City and Geographical Condition Monitoring，Hefei 230031， China；2. HeFei University of Technology, Xuancheng 242000, China)

There is a certain degree of distortion in Digital Elevation Model (DEM) depicting terrain. It is necessary of recovering distorted DEM to improve the accuracy and fidelity. In this paper, a method of recovering distorted DEM of regular terrain under constraint of terrain feature is proposed by taking a flat terrain object as an example. Firstly, the DEM of flat terrain should be a flat structure according to terrain feature, and it can be described by a plane equation. Then the expression of plane can be estimated by parameter estimation algorithms, this paper selects least square method (LS) and random sample consensus algorithm (RANSAC) to estimate it. And the recovered DEM will be obtained based on the expression. The recovery results reveal that the method is effective and feasible, and the recovered DEM not only has a higher accuracy, but also can reflect the terrain feature fully and exactly. With massive grid points, both algorithms can achieve a great recovery effect. However, RANSAC algorithm has a better robustness than LS and RANSAC can better adapt to the mutations of elevation.

digital elevation model; terrain feature; least square method; random sample consensus algorithm; recovery

周波，徐俊波.地形特征約束下的失真DEM修復方法[J].測繪通報，2017(8)：56-61.

10.13474/j.cnki.11-2246.2017.0254.

2016-12-26；

2017-02-20

安徽省智慧城市與地理國情監測重點實驗室開放性課題基金(2016-K-05Z)

周 波(1981—)，男，博士，講師，研究方向為GIS建模與應用。E-mail：zhoubo810707@163.com

P237

A

0494-0911(2017)08-0056-06