由等高線重構曲面的HASMOC方法應用研究

宋敦江，畢 誠，岳天祥

1.中國科學院 科技政策與管理科學研究所，北京 100190

2.中國兵器工業計算機應用技術研究所，北京 100089

3.中國科學院 地理科學與資源研究所，北京 100101

由等高線重構曲面的HASMOC方法應用研究

宋敦江1，畢 誠2，岳天祥3

1.中國科學院 科技政策與管理科學研究所，北京 100190

2.中國兵器工業計算機應用技術研究所，北京 100089

3.中國科學院 地理科學與資源研究所，北京 100101

1 引言

雖然新一代地球空間信息技術如InSAR、LiDAR和數字攝影測量等技術可快速獲取DEM數據，但是傳統的由人工解譯地形得到的等高線仍然是一種十分重要的用于建立DEM的數據源[1]。廣泛應用于科研與人們日常生活中的中國1∶25萬、1∶5萬及1∶1萬DEM數據就是由等高線和離散點通過約束TIN（Triangulated Irregular Networks）方法構建得到[2-3]。等高線是從攝影測量立體相對模型中通過人工判斷提取獲得，是人對地形地貌（Landforms）抽象理解和概括而獲得的一種數據，它蘊含著大量的地形特征信息[3]。等高線表示地形時有很多類似TIN模型之處，如在地形變化平緩的地方等高線稀疏，在地形變化急促的地方等高線密集。為了充分利用等高線蘊涵的大量地形特征信息，在建立DEM時應選擇針對等高線的方法，而不是選擇通用的空間插值方法，如樣條法（Spline）、反距離權重（IDW）及克里金（Kriging）方法。早在1728年人們就開始應用等高線[4]，近20多年來，國內外對于等高線構建DEM方法的研究卻一直沒間斷過，由等高線建立地形曲面的方法主要有TIN方法[5]、薄板樣條法（Thin Plate Spline，TPS）[6]、等高線膨脹法（Contour Dilation）[7]、最大中間等高線法（Max Intermediate Contours，MIC）[8]、最陡坡度法[9]和地形骨架線法（Skeleton）[10]及綜合方法（綜合使用中間等高線、最陡方向及地形骨架線等方法）[11]。這些方法各有優劣，如TIN方法在谷地和山峰生成“平三角”，兩等高線間利用線性插值，與實際地形的連續過渡的現象不符；最陡坡度法考慮到兩條等高線之間的線性關系，而實際上兩條等高線之間的最陡方向線常常難以求得，這使得最陡坡法成為一個思想方法，而不是一個具體的算法，限于篇幅，其他方法不一一論述。

HASM（High Accuracy Surface Modeling）是我國學者21世紀初提出的基于微分幾何曲面論的曲面建模新方法[12]。曲面論的基本思想是，一個曲面除了它的空間位置外，它的形狀由它的第一類基本式和第二類基本式決定[13]。HASM通過迭代計算獲得曲面的兩類基本式系數，從而模擬得到一個曲面的形狀。HASM方法目前主要用于由離散點模擬曲面的研究。HASMOC方法也被稱為基于曲面論和優化控制理論的曲面建模方法，它是在HASM方法基礎上的擴展，增加了更多的約束條件，目前主要用于由等高線（也可以包含離散點）重構地形曲面[14-15]。HASMOC方法的應用最早發表于2010年的科學計算與優化國際會議[16]，當時用于數學曲面的等高線重構曲面的研究，得到較理想的結果，詳細結果可參見《Surface Modeling：High Accuracy and High Speed Methods》第7章，2011年HASMOC方法申請了國家專利，比較了HASMOC方法與薄板樣條的模擬結果的區別[17]，但是以上相關文獻中仍然缺少本文的很多細節的描述，特別是對于不同方法的模擬結果DEM高程頻率分布特性的統計比較、回放等高線的細節分析。

多條等高線形成的“等高線群”能較好地概括和表現地形地貌特征，等高線既蘊含著定量的信息，如等高線所經過位置點的高程信息；還蘊含著定性信息，如兩條等高線間的網格點的高程上下界，等高線還蘊含著山峰、山谷和鞍部以及山脊線、溝谷線等定性信息[1]。等高線樹是用于表示等高線空間拓撲關系的一種數據結構，等高線拓撲關系主要包括父子包含關系，以及兄弟并列關系。等高線樹的研究方興未艾，自1963年“等高線樹”這一概念出現以后[18]，涌現了很多構建等高線樹的方法[19]。通過等高線樹，可以快速地計算出研究區域內任意一點與等高線的空間關系，如該點被哪條等高線包含，不被哪條等高線包含，從而獲得區域內任意一個格網點高程的上下界（為包圍該格網點的所有等高線高程的最小值與最大值，開區間）。

本文分別應用HASMOC方法和TIN方法根據等高線重構地形曲面，比較分析了這兩種方法得到的重構曲面的一些特性，包括曲面的光滑性、回放的等高線、高程值的分布頻率等。由于TIN方法的介紹說明非常常見，這里不再贅述。

2 HASMOC方法

設{(xi，yj)|xi=i×h，yj=j×h，0≤i≤I+1，0≤j≤J+1}是計算區域Ω進行均勻正交剖分的網格，其中h為網格分辨率，則HASM主方程Dirichlet邊界問題的有限差分迭代形式可表達為：

表示迭代過程中網格點(x，y)上第n次迭代的模擬ij值，等變量的具體含義，可參見文獻[12]。

設

以上兩個方程可以用矩陣形式表達為：

其中，A和bn分別為方程組（1）中第1個方程的系數矩陣和右端項向量；B和cn分別為方程組（1）中第2個方程的系數矩陣和右端項向量。設：

則HASM方法可以表示為如下等式約束最小二乘問題：

其中S∈RT×(I×J)和t∈RT×1分別為采樣矩陣（布爾矩陣）和采樣向量，T為采樣點個數。建立等式約束最小二乘方法的目的是為了在保證采樣點處模擬值接近采樣值的條件下，曲面的整體模擬誤差最小。

等高線（結合圖廓邊界），將研究區域分割成多個子區域。對于每個子區域內的格網點，根據包圍它的等高線，可以確定格網點高程值的范圍。如圖1所示，方格代表一個柵格單元，中間的點表示該柵格的中心點，標注的數字21～88表示格網點編號，等高線間距為2.5 m，等高線的高程值用紅色數字標注。

對于每個網格點，可以根據其所在的子區域，獲得其高程范圍。l、u分別表示網格點高程值上下界的列向量，f表示待求的模擬變量。網格點51被等高值為117.5和120的等高線所包圍，對網格點51的高程，應該滿足不等式l(51)＜f(51)＜u(51)，其中l(51)=117.5和u(51)=120。網格點46，被一條高程為132.5的山頂等高線包圍，同時其附近還存在一個采樣點（高程為134），這種情況，可以近似認為，這個采樣點是山峰點，故可得不等式l(46)＜f(46)＜u(46)，其中l(46)=132.5，u(46)=134。若網格點46所在的這個山峰沒有采樣點，則其高程所滿足的不等式中，上界值u(46)為135，即為圍成這片區域等高線的等高值加上1倍的等高距。對于洼地區域，可做類似處理。

圖1 根據等高線確定離散格網點上下界約束示意圖

對于每個網格點，其高程都可以建立不等式方程，從而得到不等式組，故曲面應該一直滿足不等式約束優化模型：

方程（5）即為HASMOC方法，它可以使用Matlab 7.8的lsqlin函數進行求解。lsqlin主要用于解決不等式約束的最小二乘的曲線擬合問題：

x=lsqlin(C，d，A，b，Aeq，beq，lb，ub，x0，οptiοns)，這里的C是一個m×n的稀疏矩陣，它是方程組（1）的矩陣系數，lb和ub表示的從等高線中求得的變量的上下界，x0是待求量的初始值，οptiοns是關于優化求解參數的選項，本文的約束優化求解的參數皆為默認參數，初始值為元素全部為0的列向量。

3 應用案例

中國科學院千煙洲生態試驗站位于江西省中部，隸屬于吉安市泰和縣，處于115°04′13″E，26°44′48″N，是一個由約80個小山丘，9條溝谷組成的小山村，總面積約為2.04 km2，海拔高程60～150 m左右，屬典型的紅壤丘陵區。地貌類型主要是低丘，丘頂渾圓，海拔多在100 m以下，最高達147 m，相對高度20～50 m，丘坡坡度以10°～30°居多。本實驗數據是在江西省吉安市泰和縣紙質地形圖的基礎上，針對千煙洲生態試驗站附近的地形圖掃描然后人工矢量化得到的數據，1954年北京坐標系，1956年黃海高程系，等高距2.5 m，199?年（原圖看不清），該區域內多植被，且大部分為水稻田。為驗證HASMOC方法在稀疏等高線情況下的有效性，從2.5 m等高距的等高線中抽取了20 m等高距的等高線（見圖2）（共計76條等高線）進行重構地形曲面的實驗。

圖2 千煙洲掃描矢量化原始等高線數據（20 m等高距）

圖3 DEM結果比較圖（高程單位：m）

利用圖2中的等高線數據，分別使用TIN和HASMOC方法建立5 m分辨率DEM（圖3），TIN方法模擬得到的結果為等高線的凸包，HASMOC模擬方法返回的是外接矩形區域（521行×546列），為了統一起見，只顯示了與TIN相同范圍的DEM。從圖4中可以看出，TIN在地形中的山頂和谷地出現了很多“磨平”現象，遺棄了很多地形細節，而HASMOC方法則基本保留了地形細節，地形曲面連續光滑、自然。在多處“彎月”型山谷中，TIN方法構建的DEM嚴重失真，未保證建立的DEM忠實于原始等高線；在多個馬鞍處，TIN方法也出現磨平現象。通過回放等高線（圖4）可看出，TIN方法回放的等高線（共計35條）與原始等高線的相差非常大，而HASMOC方法回放的等高線與原始等高線的形狀基本一致（都是76條等高線）。類似地，還進行了從2.5 m等高距的等高線中抽取了10 m等高距實驗，得到類似的結果，TIN方法回放的等高線（110條）與原始等高線（152條）相差非常大，而HASMOC方法則與原始等高線基本吻合，等高線的條數相等。

圖4 回放等高線（藍色+橫條）與原始等高線（黑色）疊加分析圖

理論上，通過對等高線間各子區域的約束優化控制，HASMOC方法能保證回放等高線的最大偏離距離不超過一個柵格單元，這點可以通過圖5得到佐證。圖5中的數字編號是等高線的標識號，等高線69的回放等高線與原始等高線的差別較大，但是兩條等高線的最大距離沒有超過一個柵格單元寬度。等高線53也存在3處相差較大的情況，但是都沒有超過理論上的最大距離。

圖5 HASMOC方法回放等高線局部放大圖

根據文獻[20]的研究，由等高線建立的DEM容易受到等高線數據本身的影響，如DEM中等高線高程值附近的柵格單元數相對較多，這是DEM構建方法無法避免的，這與實際地形嚴重不符，于是DEM的高程頻率分布柱狀圖也是檢驗構建DEM方法優劣的一個標準。從圖6中可以看出，TIN方法構建的DEM受原始等高線數據的影響很大，分別在高程值為80、100和120處的高程點分布較多，HASMOC方法模擬結果也受到了一定的影響，但是相對來說影響很小。

圖6 高程頻率統計柱狀圖

4 結束語

HASMOC方法較TIN方法模擬得到了更加符合實際的DEM：（1）從DEM結果來看，TIN產生大量的磨平現象，丟失了很多地形細節，HASMOC保留了大量的地形細節；（2）HASMOC回放等高線與原始等高線基本吻合，TIN的回放等高線與原始等高線差別很大；（3）TIN方法構建的DEM高程值分布頻率受等高線數據影響大，HASMOC受到的影響小。

由等高線（或等值線）重構曲面方法在機械鍛造等工業加工過程中應用廣泛，HASMOC方法采用曲面離散化方法，最后的求解過程計算量大，本文運用Matlab中的lsqlin函數進行求解，可能速度上無法滿足實際企業生產的實時需求，運算速度有待提高。運用GPU并行計算技術實現大規模優化控制問題的快速求解，是下一步的研究工作。

[1]Gallant C，Hutchinson M F.Digital elevation models and representation of terrain shape[M]//Terrain Analysis：Principles and Applications.[S.l.]：John Wiley and Sons，2000.

[2]Li Zhilin，Zhu Qing，Gold C.Digital terrain modeling：principles and methodology[M].[S.l.]：Taylor&Francis，2005.

[3]El-Sheimy N，Valeo C.Digital terrain modelling：acquisition，manipulation and applications[M].Norwood：Artech House，Inc，2005.

[4]Cayley A.On contour and slope lines[J].Philosoph Mag，1859，18：264-268.

[5]Garcia，Bello A.A contour line based triangulating algorithm[C]// Bresnahan P，Corwin E，Cowen D.Proceedings of the 5th International Symposium on Spatial Data Handling，1992：411-423.

[6]Hutchinson，Michael.Calculation of hydrologically sound digital elevation models[C]//Proceedings of the Third International Symposium on Spatial Data Handling.Columbus，Ohio：International Geographical Union，1988.

[7]Taud H，Jean-Francois P，Alvarez R.DEM generation by contour line dilation[J].Computers&Geosciences，1999，25（7）：775-783.

[8]Gousie M B，Franklin W R.Constructing a DEM from gridbased data by computing intermediate contours[C]//GIS 2003：Proceedings of the Eleventh ACM International Symposium on Advances in Geographic Information Systems，2003：71-77.

[9]Ardiansyah P O D，Yokoyama R.DEM generation method from contour lines based on the steepest slope segment chain and a monotone interpolation function[J].ISPRS Journal of Photogrammetry and Remote Sensing，2002，57（1/2）：86-101.

[10]Thibault D，Gold C M.Terrain reconstruction from contours by skeleton construction[J].GeoInformatica，2000，4：349-373.

[11]胡鵬.新數字高程模型[M].北京：測繪出版社，2007.

[12]Yue Tianxiang，Du Zhengping，Song Dunjiang.A new method of surface modeling and its application to DEM construction[J].Geomorphology，2007，91（12）：161-172.

[13]Ciarlet P G，Larsonneur F.On the recovery of a surface with prescribed first and second fundamental forms[J].J Math Pures Appl，2002，81：167-185.

[14]Yue Tianxiang，Song Dunjiang，Du Zhengping.Chapter 7“An optimal control method of HASM for DEM construction and its validation”[M]//Surface Modeling：High Accuracy and High Speed Methods.[S.l.]：CRC Press，2011.

[15]岳天祥，杜正平，宋敦江.基于曲面論和優化控制理論的曲面建模方法：中國，201110021504.8[P].2011-05-11.

[16]Song Dunjiang，Yue Tianxiang，Du Zhengping.DEM construction from contour lines based on regional optimum control[C]// Proceedings of the Third International Joint Conference on Computational Sciences and Optimization，2010：162-165.

[17]宋敦江，岳天祥，杜正平.一種由等高線構建DEM的新方法[J].武漢大學學報：信息科學版，2012，37（4）：472-476.

[18]Boyell R，Reston H.Hybrid techniques for real-time radar simulation[C]//Proceedings of the Fall Joint Computer Conference，Las Vegas，1963：445-458.

[19]Chen Jun，Qiao Chaofe，Zhao Renliang.A voronoi interior adjacency-based approach for generating a contour tree[J]. Computers&Geosciences，2004，30：355-367.

[20]Carla R，Carrara A，Bitelli G.Comparison of techniques for generating digital terrain models from contour lines[J].International Journal of Geographic Information Science，1997，11（5）：451-473.

SONG Dunjiang1,BI Cheng2,YUE Tianxiang3

1.Institute of Policy and Management,Chinese Academy of Sciences,Beijing 100190,China
2.Computer Application Technology Institute of China North Industries Group Corporation,Beijing 100089,China
3.Institute of Geographical Sciences and Natural Resources Research,Chinese Academy of Sciences,Beijing 100101,China

The HASMOC（High Accuracy Surface Modeling-Optimal Control）method is based on HASM,on which more constraint equations are added.The latter is mainly used to reconstruct surface from discrete points,while the former can be also used to reconstruct terrain surface from contours.By minimizing the norm of the combining equation of the two basic equations of the HASM method,subject to the equality constraints from sample points and bound constraints for each grid points which are surrounded by contours,HASMOC method can preserve the smoothness of the terrain surface,and the high fidelity to the original contours.A real contour lines example is given,and results from HASMOC are compared with that from TIN（Triangulated Irregular Network）method,as for the respects of retrieved contour,details of surface and histogram of surface height distribution,the former is superior to the latter in terrain surface reconstruction from real contour lines.

High Accuracy Surface Modeling（HASM）;optimal control;Triangulated Irregular Network（TIN）;retrieved contours;histogram of height distribution

HASM優化控制方法（High Accuracy Surface Modeling-Optimal Control，HASMOC）是在高精度曲面建模（HASM）方法的基礎上，增加更多約束條件方程后形成的一種方法。通過對等高線間格網點高程范圍的約束優化控制，最小化HASM基本方程的模，HASMOC方法既能保證地形曲面的整體光滑性，又保證地形曲面對于原始等高線數據的忠實性。實際案例表明，HASMOC方法得到的地形曲面結果優于TIN方法的地形曲面模擬結果；比較分析地形曲面的回放等高線、地形光滑程度和地形曲面的高程分布頻率等，可以看出，HASMOC方法能較好地克服TIN的缺點。

高精度曲面建模（HASM）；優化控制；不規則三角網（TIN）；回放等高線；高程分布頻率

A

TP391

10.3778/j.issn.1002-8331.1304-0349

SONG Dunjiang,BI Cheng,YUE Tianxiang.Application of HASMOC method in terrain surface reconstruction from contours.Computer Engineering and Applications,2013,49（18）：171-175.

中科院135創新項目（No.Y201131Z06）；國家自然科學基金青年科學基金（No.40801187）；國家自然科學基金杰出青年科學基金（No.40825003）。

宋敦江（1979—），男，博士，副研究員，研究領域為計算管理學；畢誠（1975—），男，高級工程師，研究領域為機電一體化與可視化；岳天祥（1963—），男，博士，研究員，研究領域為資源環境模型與系統模擬。E-mail：songdj@casipm.ac.cn

2013-04-24

2013-07-23

1002-8331（2013）18-0171-05