宋敦江，畢 誠，岳天祥

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

2.中國兵器工業(yè)計算機應(yīng)用技術(shù)研究所，北京 100089

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

由等高線重構(gòu)曲面的HASMOC方法應(yīng)用研究

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

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

2.中國兵器工業(yè)計算機應(yīng)用技術(shù)研究所，北京 100089

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

1 引言

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

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

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

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

2 HASMOC方法

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

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

設(shè)

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

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

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

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

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

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

圖1 根據(jù)等高線確定離散格網(wǎng)點上下界約束示意圖

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

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

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

3 應(yīng)用案例

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

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

圖3 DEM結(jié)果比較圖（高程單位：m）

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

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

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

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

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

圖6 高程頻率統(tǒng)計柱狀圖

4 結(jié)束語

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

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

[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，F(xiàn)ranklin 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]胡鵬.新數(shù)字高程模型[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]岳天祥，杜正平，宋敦江.基于曲面論和優(yōu)化控制理論的曲面建模方法：中國，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]宋敦江，岳天祥，杜正平.一種由等高線構(gòu)建DEM的新方法[J].武漢大學(xué)學(xué)報：信息科學(xué)版，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優(yōu)化控制方法（High Accuracy Surface Modeling-Optimal Control，HASMOC）是在高精度曲面建模（HASM）方法的基礎(chǔ)上，增加更多約束條件方程后形成的一種方法。通過對等高線間格網(wǎng)點高程范圍的約束優(yōu)化控制，最小化HASM基本方程的模，HASMOC方法既能保證地形曲面的整體光滑性，又保證地形曲面對于原始等高線數(shù)據(jù)的忠實性。實際案例表明，HASMOC方法得到的地形曲面結(jié)果優(yōu)于TIN方法的地形曲面模擬結(jié)果；比較分析地形曲面的回放等高線、地形光滑程度和地形曲面的高程分布頻率等，可以看出，HASMOC方法能較好地克服TIN的缺點。

高精度曲面建模（HASM）；優(yōu)化控制；不規(guī)則三角網(wǎng)（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創(chuàng)新項目（No.Y201131Z06）；國家自然科學(xué)基金青年科學(xué)基金（No.40801187）；國家自然科學(xué)基金杰出青年科學(xué)基金（No.40825003）。

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

2013-04-24

2013-07-23

1002-8331（2013）18-0171-05