本發(fā)明屬于地理信息系統(tǒng)(GIS)和計算機圖學(xué)領(lǐng)域，涉及地形的分類方法和地形參數(shù)提取及簡化模型的建立方法，具體地說就是數(shù)字地圖中典型地形幾何參數(shù)獲取方法，可用于野外大型地形場景的電波傳播近似計算和快速實時顯示。

背景技術(shù)：

隨著近年來計算機、衛(wèi)星遙感等科技的迅猛發(fā)展，數(shù)字化浪潮正在席卷全球。地球作為人類活動和生活的首要載體，人們在認識自然過程中以及對自然進行改造過程中，對周圍地形地貌的環(huán)境信息一直不斷試著利用不同的方法進行描述及表達。在眾多的方法中，數(shù)字地形圖這種表達方法能夠較好對地形表面形態(tài)狀況進行描述和表達。通過地理測繪等相關(guān)手段得到數(shù)字化坐標數(shù)據(jù)，地貌仿真后產(chǎn)生三維可視化效果，地形地貌結(jié)構(gòu)可以很好被表達，這種方法在二十世紀以來應(yīng)用較為廣泛。上個世紀四十年代計算機技術(shù)誕生及其不斷蓬勃發(fā)展，以及伴隨著計算機圖學(xué)(Computer Graphics)、現(xiàn)代數(shù)學(xué)理論、數(shù)據(jù)挖掘、模式識別等相關(guān)技術(shù)理論完善和應(yīng)用，都極大地促進了數(shù)字地圖的發(fā)展及應(yīng)用。除了在互聯(lián)網(wǎng)和手機等行業(yè)數(shù)字地圖有著廣泛的應(yīng)用之外，數(shù)字地圖還有許多其他的重要應(yīng)用前景，例如我們?nèi)绻馨迅鞣N物體、個人、建筑、交通、景觀等動態(tài)的和靜態(tài)的數(shù)字信息和三維數(shù)字地圖有機結(jié)合成一個整體，就會形成一個數(shù)字化的城市應(yīng)用服務(wù)場景，將對我們的生活產(chǎn)生無可比擬的便捷和享受。對于數(shù)字地圖，常常是通過讀取各種不同精度的高程原始坐標數(shù)據(jù)構(gòu)成多邊形(常用三角形)面片，然后運用多邊形面片無限逼近真實地貌。數(shù)字高程是利用有序的數(shù)據(jù)陣列來存儲地面高程，其模型(DigitalElevation Model，簡稱DEM)是歸屬于數(shù)字地形模型(Digital TerrainModel，簡稱DTM)的一個分支。由于DEM數(shù)據(jù)是一個DTM分析研究所需的基本數(shù)據(jù)源，所以我們希望從中提取各種各樣有用的信息，可為國民經(jīng)濟發(fā)展、生態(tài)農(nóng)業(yè)、資源利用和作戰(zhàn)指揮等相關(guān)工作提供參考。DEM就是用數(shù)字建模的形式來表達現(xiàn)實地球表面環(huán)境的方法，自從二十世紀五十年代被采用，在各種應(yīng)用領(lǐng)域得到了極為廣泛的應(yīng)用。隨著計算機軟硬件水平以及計算機圖學(xué)等技術(shù)爆發(fā)式的發(fā)展，DEM從獲取方式、存儲結(jié)構(gòu)以及處理速度等方面都獲得了質(zhì)的飛躍。如今，伴隨大量可靠高精度DEM數(shù)據(jù)的產(chǎn)生，亟需根據(jù)各種現(xiàn)實需求場景開展DEM數(shù)據(jù)應(yīng)用研究。

如今的數(shù)字地圖的發(fā)展促使我們通過各種途徑獲得越來越精確的大量地圖信息，伴隨著虛擬數(shù)字地球概念的提出，數(shù)字地圖的用途不斷從軍用逐漸延伸到民用，日常的生活場景更是逐漸離不開其相關(guān)的應(yīng)用。三維地圖數(shù)據(jù)挖掘就是要以實際需求為出發(fā)點，從三維地圖數(shù)據(jù)中挖掘出可行的、規(guī)律的、有價值的信息，使之服務(wù)于實際研究和現(xiàn)實生產(chǎn)。三維地圖數(shù)據(jù)的挖掘可以通過采用諸多挖掘方法，例如通過統(tǒng)計、專家系統(tǒng)(經(jīng)驗體系)、機器學(xué)習、歸納法和模式識別等。地形地貌分類的傳統(tǒng)獲取方式主要利用人工手段辨識地形信息，其工作量大、周期長、難以完成大范圍地形信息的快速低成本獲取，針對遙感地圖自動識別的問題，目前也尚無完全自動化的針對性的實用自動分類方法?，F(xiàn)代社會生活中通信電子設(shè)備越來越多，不斷面臨著電子設(shè)備之間的耦合問題，也就是電磁干擾與兼容的問題。電磁兼容(ElectromagneticCompatibility，EMC)作為現(xiàn)代熱門的重要前言研究，主要包括了兩部分：EMI(電磁干擾)和EMS(電磁耐受性)。隨著研究不斷深入，各種電磁分析軟件不斷專業(yè)化、自動化，其中關(guān)于電波傳播計算的場景建模主要有兩種方式：一種是依靠人機交互的方式手工輸入復(fù)雜地形環(huán)境數(shù)據(jù)實現(xiàn)建模；另一種是基于數(shù)字地圖的自動建模，但其場景限于城市及其周圍范圍，場景中待計算的物及地形較為單一、尺寸不是很大。

數(shù)字地圖的地形信息提取研究內(nèi)容包括對山頂點的提取，對山底面信息的提取，對地形的分類，以及對所提取的地形進行簡化。文獻“地理信息系統(tǒng)原理及應(yīng)用”(電子工業(yè)出版社，2011，作者：胡祥培等)從理論的角度，詳細闡述了空間地理數(shù)據(jù)的獲取以及處理，并論述了空間數(shù)據(jù)的管理和地理空間分析；碩士論文“基于DEM的地形特征提取算法研究及應(yīng)用”(西安建筑科技大學(xué)，2012，作者：易煒)研究了目前現(xiàn)有的各類山頂點、鞍部點和山脊線的提取算法，對各種地形進行了定性和定量的分析，并通過對DEM格網(wǎng)數(shù)據(jù)進行等高距分層以模擬其在等高線地圖中呈現(xiàn)的特性和拓撲結(jié)構(gòu)，提出了相應(yīng)的地形特征提取算法；碩士論文“利用圖像處理技術(shù)提取地形特征線的方法研究”(西安建筑科技大學(xué)，2011，作者：方莉)在認真分析地形特征線的物理形態(tài)和結(jié)構(gòu)特征的基礎(chǔ)上，提出了一種利用灰度形態(tài)學(xué)算子提取山脊線和山谷線的新方法，另外針對山腳線的提取利用其處于山體和平地交界處這一典型物理特征，考慮將灰度圖像的統(tǒng)計原理與地形坡度相結(jié)合設(shè)計一種將平地區(qū)域分割出來的方法，最后通過提取平地的邊緣得到山腳線；論文“地貌認知及空間剖分的山頂點提取”(測繪科學(xué)，2010年9月，126-127，作者：羅明良等)通過對山頂局部形態(tài)特征及空間定量特征分析，結(jié)合傳統(tǒng)地貌學(xué)認知思路，給出了基于空間剖分的山頂點快速提取方法，空間剖分標準對應(yīng)于地貌學(xué)山地分類，通過高差限制剖分實現(xiàn)DEM分塊，并在此基礎(chǔ)上擬合函數(shù)，給出符合地貌認知的山頂點；論文“山頂點類型及其形態(tài)特征數(shù)字表達”(南京師范大學(xué)學(xué)報，2010年3月，136-130，作者：蒼學(xué)智等)在分析了山頂點所具有的空間特征、成因特征、尺度特征及點群特征的基礎(chǔ)上，依照科學(xué)性、系統(tǒng)性、實用性、可實現(xiàn)性的原則，對山頂點進行了系統(tǒng)的分類，并闡述了其定量描述方法，并且還深入探討了山頂點群的基本類型、定量描述指標與地學(xué)意義?？梢钥闯觯延醒芯砍晒嬖谝韵聠栴}：(1)研究僅涉及到地形特殊點(如山頂點、鞍部點)、特殊線(山脊線、山谷線、山腳線)等局部幾何特征的識別，缺乏對地形整體幾何形狀識別、簡化和幾何參數(shù)提取的研究；(2)研究方法較為單一，主要基于幾何圖形學(xué)，未充分利用圖像學(xué)方法進行研究；(3)未充分利用地形學(xué)中的坡度、剖面曲率、正切曲率、地形起伏度等地形因子開展研究，存在一定的局限性。

技術(shù)實現(xiàn)要素：

本發(fā)明的目的在于克服上述現(xiàn)有研究中存在的問題，提供一種數(shù)字地圖中典型地形幾何參數(shù)獲取方法，以滿足野外大型地形場景的電波傳播近似計算和快速實時建模及顯示的需要。

本發(fā)明的目的是這樣實現(xiàn)的：數(shù)字地圖中典型地形幾何參數(shù)獲取方法，其特征是：至少包括如下步驟：

步驟101：使用ArcGIS工具，打開一個ASTER GDEM類型的GeoTIFF地圖文件；

步驟102：將所有柵格點的高程值作歸一化處理，結(jié)果保存在DEM_1.gif文件中；歸一化處理的方法為：設(shè)x為任一柵格點處的量值，該量值是高程、坡度、坡度變率、曲率、地形起伏度中的任一種，x_max是所有柵格點量值的最大值，x_min是所有柵格點量值的最小值，x′是x的歸一化值，則x′＝(x-x_min)×255/(x_max-x_min)；

步驟103：計算所有柵格點的坡度值，再作歸一化處理，歸一化的坡度值保存在Slope_1.gif文件中；

步驟104：計算所有柵格點的坡度變率，再作歸一化處理，歸一化的坡度變率值保存在PoDuBianLv_1.gif文件中；

步驟105：計算所有柵格點的曲率，并對曲率求絕對值，再作歸一化處理，歸一化的曲率值保存在Curvature_1.gif文件中；

步驟106：計算所有柵格點的地形起伏度，結(jié)果保存在QiFuDu.gif文件中，再對其作歸一化處理，歸一化的地形起伏度值保存在QiFuDu_1.gif文件中；

步驟107：將歸一化的高程、坡度、坡度變率、曲率和地形起伏度做為特性指標，采用ISO(Iterative SelfOrganization)非監(jiān)督聚類得到地形的初步分類結(jié)果，每類地形有一個類別編號，初步分類結(jié)果保存在JuLei.gif文件中；

步驟108：使用GDAL庫函數(shù)打開步驟101的地圖文件，打開

QiFuDu.gif文件和JuLei.gif文件；

步驟109：依據(jù)地形平均起伏度對步驟107初步分類結(jié)果進行地形判定分類；

步驟110：用四位正整數(shù)標志地形類別及區(qū)域，前一位表示地形類別，平原、丘陵、低山、中山、高山分別用1、2、3、4、5表示；后三位表示此種類別地形的區(qū)域個數(shù)編號，從001開始累積編號，地形判定分類的結(jié)果用此種方式表示成數(shù)值之后，保存在地形分類結(jié)果文件FenLei.gif之中；

步驟111：讀取分類結(jié)果文件FenLei.gif，根據(jù)地形類別及區(qū)域標志數(shù)字選擇待研究地形區(qū)域；

步驟112：判斷地形類別標志是否為1，若是，轉(zhuǎn)至步驟113；若不是，轉(zhuǎn)至步驟115；

步驟113：將平原地形簡化為平面，求取地形類別標志為平原的所有柵格點高程的平均值，記其為平均高程Ea；

步驟114：求取地形類別標志為平原的所有柵格點地形起伏度的平均值，記其為平均起伏度Ha，轉(zhuǎn)至步驟127；

步驟115：對待研究的地形區(qū)域進行八鄰域邊界跟蹤，提取出邊界點序列存入鏈表Point_List中；

步驟116：求地形區(qū)域的面積S、地形區(qū)域的中心C(x0,y0)和地形的體積V：由鏈表Point_List存儲的邊界點序列連線形成一個邊界多邊形，用圖形學(xué)中的“交點計數(shù)法”判斷柵格點是否在邊界多邊形內(nèi)，統(tǒng)計處在邊界多邊形內(nèi)和位于邊界多邊形邊上的柵格點數(shù)目num，設(shè)單個柵格點所占面積為a，則地形區(qū)域的面積S＝num×a；地形區(qū)域的中心點C的坐標x_i、y_i是邊界多邊形內(nèi)及邊界多邊形邊上的任一柵格點的坐標；地形的體積h_i是邊界多邊形內(nèi)及邊界多邊形邊上的任一柵格點的高程；

步驟117：地形類別標志是否為2，若是，轉(zhuǎn)至步驟118；若不是，轉(zhuǎn)至步驟120；

步驟118：將丘陵地形簡化為球缺，求取球缺的底面圓半徑

步驟119：求取球缺的高度h1：解方程(h1)³+3(r1)²(h1)-6V/π＝0即得球缺的高度轉(zhuǎn)至步驟127；

步驟120：求地形區(qū)域的周長L：由鏈表Point_List存儲的邊界點序列連線形成一個邊界多邊形，邊界點即為邊界多邊形的頂點，設(shè)多邊形頂點個數(shù)為n，令第n+1頂點等于第1個頂點，邊界多邊形任一邊的邊長|P_iP_i+1|由邊的兩個頂點P_i、P_i+1通過兩點距離公式求得，則邊界多邊形的周長即為地形區(qū)域的周長L，

步驟121：求地形區(qū)域的形狀因子f：形狀因子f＝4πS/L²；

步驟122：判斷形狀因子f是否大于或等于0.7，若是，轉(zhuǎn)至步驟123；若不是，轉(zhuǎn)至步驟125；

步驟123：將地形簡化為錐形山，求錐形山底面圓的半徑

步驟124：求取錐形山的高度轉(zhuǎn)至步驟127；

步驟125：將地形簡化為楔形山，求取楔形山底面區(qū)域的最小外接矩形；

步驟126：求取楔形山的高度其中a、b為步驟125求取的外接矩形的兩個邊長，轉(zhuǎn)至步驟127；

步驟127：輸出地形幾何參數(shù)：若地形為平原，輸出平均高程Ea、平均起伏度Ha；若地形為丘陵，輸出球缺底面圓中心C(x0,y0)、半徑r1和球缺高度h1；若地形為錐形山，輸出圓錐底面圓中心C(x0,y0)、半徑r1和圓錐高度h2；若地形為楔形山，輸出楔形體底面矩形中心C(x0,y0)，邊長a、b，矩形方向角an和楔形體高度h3。

所述的步驟109中依據(jù)地形平均起伏度對步驟107初步分類結(jié)果進行地形判定分類，包括以下步驟：

步驟201：讀取JuLei.gif文件，遍歷地形類別編號，記該編號為n；

步驟202：遍歷類別n的地形范圍，從QiFuDu.gif文件中讀取相對應(yīng)位置的地形起伏度數(shù)據(jù)；

步驟203：求取該類別地形區(qū)域的地形起伏度平均值U；

步驟204：判斷U是否小于或等于20米，若是，轉(zhuǎn)至步驟205；若不是，轉(zhuǎn)至步驟206；

步驟205：該類別地形判定為平原；

步驟206：判斷U是否小于或等于200米，若是，轉(zhuǎn)至步驟207；若不是，轉(zhuǎn)至步驟208；

步驟207：該類別地形判定為丘陵；

步驟208：判斷U是否小于或等于500米，若是，轉(zhuǎn)至步驟209；若不是，轉(zhuǎn)至步驟210；

步驟209：該類別地形判定為低山；

步驟210：判斷U是否小于或等于1500米，若是，轉(zhuǎn)至步驟211；若不是，轉(zhuǎn)至步驟212；

步驟211：該類別地形判定為中山；

步驟212：該類別地形判定為高山。

所述的步驟115中對待研究的地形區(qū)域進行八鄰域邊界跟蹤，提取出邊界點序列存入鏈表Point_List中，包括以下步驟：

步驟301：基于GDAL讀取地形分類結(jié)果文件FeiLei.gif；

步驟302：根據(jù)待研究的地形區(qū)域的四位正整數(shù)標志，將其地形類別標志數(shù)字記作type；

步驟303：按照從上至下、由左到右的順序，從文件FeiLei.gif中圖像左上角依次讀取柵格點的地形類別標志數(shù)字，當?shù)谝淮巫x到的柵格點的地形類別標志數(shù)字等于type時，則將該點存入邊界點序列鏈表Point_List中，并記錄該點為邊界起點；

步驟304：令當前邊界點等于邊界起點，令下一個邊界點的搜索方向碼Directcode＝0；邊界點搜索采用八鄰域方法，每個柵格點有八個其他的柵格點與其相鄰，該柵格點正上方柵格點的搜索方向碼記為0，該柵格點右上方、正右方、右下方、正下方、左下方、正左方、左上方柵格點的搜索方向碼依次記為1、2、3、4、5、6、7；

步驟305：讀取當前邊界點沿Directcode方向的搜索點的地形類別標志數(shù)字，將其記作value；

步驟306：判斷value是否等于type，若是，轉(zhuǎn)至步驟308；若不是，轉(zhuǎn)至步驟307；

步驟307：令Directcode＝Directcode+1，轉(zhuǎn)至步驟305；

步驟308：令當前邊界點等于此搜索點；

步驟309：判斷當前邊界點是否等于邊界起點，若是，則邊界跟蹤結(jié)束；若不是，轉(zhuǎn)至步驟310；

步驟310：將當前邊界點存入鏈表Point_List中；

步驟311：令Directcode＝Directcode-2；

步驟312：判斷Directcode是否大于或等于0，若是，轉(zhuǎn)至步驟305；若不是，轉(zhuǎn)至步驟313；

步驟313：令Directcode＝Directcode+8，轉(zhuǎn)至步驟305。

所述的步驟125中將地形簡化為楔形山，求取楔形山底面區(qū)域的最小外接矩形，包括以下步驟：

步驟401：讀取地形區(qū)域邊界點序列鏈表Point_List，令鏈表Rot_List等于鏈表Point_List，ang＝3°,定義二維數(shù)組Data[30][4]；

步驟402：對鏈表Rot_List中的多邊形各頂點做以地形區(qū)域中心C(x0,y0)為旋轉(zhuǎn)中心、3°為轉(zhuǎn)角的旋轉(zhuǎn)變換，旋轉(zhuǎn)后的多邊形各頂點保存在鏈表Rot_List中；

步驟403：求鏈表Rot_List中的多邊形的軸向包圍盒的坐標x_max、x_min、y_max、y_min，計算此包圍盒邊長L1＝|x_max-x_min|、L2＝|y_max-y_min|，計算此包圍盒的面積s0＝L1×L2；

步驟404：將ang、L1、L2、s0四個值作為一行依次存入數(shù)組Data的第1列、第2列、第3列、第4列，令ang＝ang+3°；

步驟405：將步驟402～404循環(huán)29次；

步驟406：求數(shù)組Data中第4列包圍盒面積的最小值，記錄該面積最小值所在的數(shù)組行號為m；

步驟407：讀取數(shù)組Data中第m行的前3列值，分別用an、a、b記錄；令an＝90°-an，定義an為矩形方向角，它是長度為a的矩形邊與正東方向的夾角；

步驟408：輸出地形底面區(qū)域的最小外接矩形邊長a、b及矩形方向角an。

本發(fā)明有如下優(yōu)點：

(1)提出了圖像處理與圖形處理相結(jié)合的典型地形整體幾何形狀識別的方法；

(2)通過地形簡化，解決了基于數(shù)字地圖的典型地形快速幾何建模問題；

(3)所提方法可行、實用。

附圖說明

圖1是本發(fā)明的總流程圖；

圖2是地形判定分類的流程圖；

圖3是邊界跟蹤提取的流程圖；

圖4是楔形山底面區(qū)域最小外接矩形求取的流程圖；

圖5是一座山的實際形狀；

圖6是圖5所示的山進行地形分類和邊界跟蹤提取后圖示化的結(jié)果；

圖7是圖5所示山簡化并提取幾何參數(shù)后再顯示的結(jié)果。

具體實施方式

本發(fā)明所讀取的數(shù)據(jù)源，是一種數(shù)字高程模型(DEM)，數(shù)據(jù)格式是ASTER GDEM(Advanced Spaceborn Thermal Emission and Reflection Radiometer Global Digital Elevation Model)類型的GeoTIFF遙感影像數(shù)據(jù)(.gif文件)，該圖像柵格點對應(yīng)的空間分辨率約為30m×30m。借助ArcGIS工具，計算歸一化的高程、坡度、坡度變率、曲率、地形起伏度，以其為特性指標采用ISO(Iterative SelfOrganization)非監(jiān)督聚類得到地形的初步分類結(jié)果；通過GDAL(Geospatial Data Abstraction Library)這種API讀取初步分類結(jié)果文件并處理，按照各類地形平均起伏度的大小，將地形分為平原、丘陵、低山、中山、高山等幾種地形；之后，對每一類地形，通過八鄰域邊界跟蹤和地形區(qū)域底面的面積、中心位置、周長、形狀因子、最小外接矩形的計算以及區(qū)域的體積計算，將地形簡化為丘陵、錐形山、楔形山等典型地形并得到底面的中心位置和幾何尺寸，再根據(jù)體積和底面尺寸求出地形的高度尺寸，最終得到這幾種典型地形的位置、底面及高度尺寸等幾何參數(shù)。

參照圖1，本發(fā)明的數(shù)字地圖中典型地形幾何參數(shù)獲取方法包括如下步驟：

步驟101：使用ArcGIS工具，打開一個ASTER GDEM類型的GeoTIFF地圖文件，如打開ASTGTM2_N32E115_dem.gif文件；

步驟102：將所有柵格點的高程值作歸一化處理，結(jié)果保存在DEM_1.gif文件中；歸一化處理的方法為：設(shè)x為任一柵格點處的量值，該量值是高程、坡度、坡度變率、曲率、地形起伏度中的任一種，x_max是所有柵格點量值的最大值，x_min是所有柵格點量值的最小值，x′是x的歸一化值，則x′＝(x-x_min)×255/(x_max-x_min)；歸一化處理的目的是將有量綱的數(shù)據(jù)轉(zhuǎn)化為無量綱的數(shù)據(jù)并使轉(zhuǎn)化后的數(shù)據(jù)值范圍為0～255,為下面地形的模糊聚類分類做準備；

步驟103：計算所有柵格點的坡度值，再作歸一化處理，歸一化的坡度值保存在Slope_1.gif文件中；

步驟104：計算所有柵格點的坡度變率，再作歸一化處理，歸一化的坡度變率值保存在PoDuBianLv_1.gif文件中；

步驟105：計算所有柵格點的曲率，并對曲率求絕對值，再作歸一化處理，歸一化的曲率值保存在Curvature_1.gif文件中；

步驟106：計算所有柵格點的地形起伏度，結(jié)果保存在QiFuDu.gif文件中，再對其作歸一化處理，歸一化的地形起伏度值保存在QiFuDu_1.gif文件中；地形起伏度統(tǒng)計范圍取60×60柵格點，即統(tǒng)計范圍為1.8km×1.8km；

步驟107：將歸一化的高程DEM_1.gif、坡度Slope_1.gif、坡度變率PoDuBianLv_1.gif、曲率Curvature_1.gif和地形起伏度QiFuDu_1.gif做為特性指標，設(shè)定類數(shù)目如25，采用ISO(Iterative SelfOrganization)非監(jiān)督聚類得到地形的初步分類結(jié)果，共有17類，每類地形有一個類別編號，初步分類結(jié)果保存在JuLei.gif文件中；

步驟108：使用GDAL庫函數(shù)打開步驟101的地圖文件ASTGTM2_N32E115_dem.gif，打開QiFuDu.gif文件和JuLei.gif文件；

步驟109：依據(jù)地形平均起伏度對步驟107初步分類結(jié)果進行地形判定分類；

步驟109中依據(jù)地形平均起伏度對步驟107初步分類結(jié)果進行地形判定分類，參照圖2，包括以下步驟：

步驟201：讀取JuLei.gif文件，遍歷地形類別編號，記該編號為n；

步驟202：遍歷類別n的地形范圍，從QiFuDu.gif文件中讀取相對應(yīng)位置的地形起伏度數(shù)據(jù)；

步驟203：求取該類別地形區(qū)域的地形起伏度平均值U；

步驟204：判斷U是否小于或等于20米，若是，轉(zhuǎn)至步驟205；若不是，轉(zhuǎn)至步驟206；

步驟205：該類別地形判定為平原；

步驟206：判斷U是否小于或等于200米，若是，轉(zhuǎn)至步驟207；若不是，轉(zhuǎn)至步驟208；

步驟207：該類別地形判定為丘陵；

步驟208：判斷U是否小于或等于500米，若是，轉(zhuǎn)至步驟209；若不是，轉(zhuǎn)至步驟210；

步驟209：該類別地形判定為低山；

步驟210：判斷U是否小于或等于1500米，若是，轉(zhuǎn)至步驟211；若不是，轉(zhuǎn)至步驟212；

步驟211：該類別地形判定為中山；

步驟212：該類別地形判定為高山。

步驟110：用四位正整數(shù)標志地形類別及區(qū)域，前一位表示地形類別，平原、丘陵、低山、中山、高山分別用1、2、3、4、5表示；后三位表示此種類別地形的區(qū)域個數(shù)編號，從001開始累積編號，地形判定分類的結(jié)果用此種方式表示成數(shù)值之后，保存在地形分類結(jié)果文件FenLei.gif之中；

步驟111：讀取分類結(jié)果文件FenLei.gif，根據(jù)地形類別及區(qū)域標志數(shù)字選擇待研究地形區(qū)域；

步驟112：判斷地形類別標志是否為1，若是，轉(zhuǎn)至步驟113；若不是，轉(zhuǎn)至步驟115；

步驟113：將平原地形簡化為平面，求取地形類別標志為平原的所有柵格點高程的平均值，記其為平均高程Ea；

步驟114：求取地形類別標志為平原的所有柵格點地形起伏度的平均值，記其為平均起伏度Ha，轉(zhuǎn)至步驟127；

步驟115：對待研究的地形區(qū)域進行八鄰域邊界跟蹤，提取出邊界點序列存入鏈表Point_List中；

步驟115中對待研究的地形區(qū)域進行八鄰域邊界跟蹤，提取出邊界點序列存入鏈表Point_List中，參照圖3，包括以下步驟：

步驟301：基于GDAL讀取地形分類結(jié)果文件FeiLei.gif；

步驟302：根據(jù)待研究的地形區(qū)域的四位正整數(shù)標志，將其地形類別標志數(shù)字記作type；

步驟303：按照從上至下、由左到右的順序，從文件FeiLei.gif中圖像左上角依次讀取柵格點的地形類別標志數(shù)字，當?shù)谝淮巫x到的柵格點的地形類別標志數(shù)字等于type時，則將該點存入邊界點序列鏈表Point_List中，并記錄該點為邊界起點；

步驟304：令當前邊界點等于邊界起點，令下一個邊界點的搜索方向碼Directcode＝0；邊界點搜索采用八鄰域方法，每個柵格點有八個其他的柵格點與其相鄰，該柵格點正上方柵格點的搜索方向碼記為0，該柵格點右上方、正右方、右下方、正下方、左下方、正左方、左上方柵格點的搜索方向碼依次記為1、2、3、4、5、6、7；

步驟305：讀取當前邊界點沿Directcode方向的搜索點的地形類別標志數(shù)字，將其記作value；

步驟306：判斷value是否等于type，若是，轉(zhuǎn)至步驟308；若不是，轉(zhuǎn)至步驟307；

步驟307：令Directcode＝Directcode+1，轉(zhuǎn)至步驟305；

步驟308：令當前邊界點等于此搜索點；

步驟309：判斷當前邊界點是否等于邊界起點，若是，則邊界跟蹤結(jié)束；若不是，轉(zhuǎn)至步驟310；

步驟310：將當前邊界點存入鏈表Point_List中；

步驟311：令Directcode＝Directcode-2；

步驟312：判斷Directcode是否大于或等于0，若是，轉(zhuǎn)至步驟305；若不是，轉(zhuǎn)至步驟313；

步驟313：令Directcode＝Directcode+8，轉(zhuǎn)至步驟305。

步驟116：求地形區(qū)域的面積S、地形區(qū)域的中心C(x0,y0)和地形的體積V：由鏈表Point_List存儲的邊界點序列連線形成一個邊界多邊形，用圖形學(xué)中的“交點計數(shù)法”判斷柵格點是否在邊界多邊形內(nèi)，統(tǒng)計處在邊界多邊形內(nèi)和位于邊界多邊形邊上的柵格點數(shù)目num，設(shè)單個柵格點所占面積為a，則地形區(qū)域的面積S＝num×a；地形區(qū)域的中心點C的坐標x_i、y_i是邊界多邊形內(nèi)及邊界多邊形邊上的任一柵格點的坐標；地形的體積h_i是邊界多邊形內(nèi)及邊界多邊形邊上的任一柵格點的高程；

步驟117：地形類別標志是否為2，若是，轉(zhuǎn)至步驟118；若不是，轉(zhuǎn)至步驟120；

步驟118：將丘陵地形簡化為球缺，求取球缺的底面圓半徑

步驟119：求取球缺的高度h1：解方程(h1)³+3(r1)²(h1)-6V/π＝0即得球缺的高度轉(zhuǎn)至步驟127；

步驟120：求地形區(qū)域的周長L：由鏈表Point_List存儲的邊界點序列連線形成一個邊界多邊形，邊界點即為邊界多邊形的頂點，設(shè)多邊形頂點個數(shù)為n，令第n+1頂點等于第1個頂點，邊界多邊形任一邊的邊長|P_iP_i+1|由邊的兩個頂點P_i、P_i+1通過兩點距離公式求得，則邊界多邊形的周長即為地形區(qū)域的周長L，

步驟121：求地形區(qū)域的形狀因子f：形狀因子f＝4πS/L²；

步驟122：判斷形狀因子f是否大于或等于0.7，若是，轉(zhuǎn)至步驟123；若不是，轉(zhuǎn)至步驟125；

步驟123：將地形簡化為錐形山，求錐形山底面圓的半徑

步驟124：求取錐形山的高度轉(zhuǎn)至步驟127；

步驟125：將地形簡化為楔形山，求取楔形山底面區(qū)域的最小外接矩形；

步驟125中將地形簡化為楔形山，求取楔形山底面區(qū)域的最小外接矩形，參照圖4，包括以下步驟：

步驟401：讀取地形區(qū)域邊界點序列鏈表Point_List，令鏈表Rot_List等于鏈表Point_List，ang＝3°,定義二維數(shù)組Data[30][4]；

步驟402：對鏈表Rot_List中的多邊形各頂點做以地形區(qū)域中心C(x0,y0)為旋轉(zhuǎn)中心、3°為轉(zhuǎn)角的旋轉(zhuǎn)變換，旋轉(zhuǎn)后的多邊形各頂點保存在鏈表Rot_List中；

步驟403：求鏈表Rot_List中的多邊形的軸向包圍盒的坐標x_max、x_min、y_max、y_min，計算此包圍盒邊長L1＝|x_max-x_min|、L2＝|y_max-y_min|，計算此包圍盒的面積s0＝L1×L2；

步驟404：將ang、L1、L2、s0四個值作為一行依次存入數(shù)組Data的第1列、第2列、第3列、第4列，令ang＝ang+3°；

步驟405：將步驟402～404循環(huán)29次；

步驟406：求數(shù)組Data中第4列包圍盒面積的最小值，記錄該面積最小值所在的數(shù)組行號為m；

步驟407：讀取數(shù)組Data中第m行的前3列值，分別用an、a、b記錄；令an＝90°-an，定義an為矩形方向角，它是長度為a的矩形邊與正東方向的夾角；

步驟408：輸出地形底面區(qū)域的最小外接矩形邊長a、b及矩形方向角an。

步驟126：求取楔形山的高度其中a、b為步驟125求取的外接矩形的兩個邊長，轉(zhuǎn)至步驟127；

步驟127：輸出地形幾何參數(shù)：若地形為平原，輸出平均高程Ea、平均起伏度Ha；若地形為丘陵，輸出球缺底面圓中心C(x0,y0)、半徑r1和球缺高度h1；若地形為錐形山，輸出圓錐底面圓中心C(x0,y0)、半徑r1和圓錐高度h2；若地形為楔形山，輸出楔形體底面矩形中心C(x0,y0)，邊長a、b，矩形方向角an和楔形體高度h3。

仿真實例

利用本發(fā)明對國內(nèi)某座名為孤峰山的地形進行幾何參數(shù)提取。圖5是讀取處理孤峰山的DEM數(shù)據(jù)并以網(wǎng)格形式顯示的孤峰山的實際形狀。通過對孤峰山所在的ASTER GDEM類型的GeoTIFF格式的數(shù)字地圖文件用ArcGIS工具計算歸一化的高程、坡度、坡度變率、曲率、地形起伏度并以其為特性指標采用ISO非監(jiān)督聚類得到地形初步分類；再對初步分類結(jié)果按照各類地形平均起伏度的大小進行地形判定分類，選取孤峰山作為待研究地形區(qū)域，圖6是該地形區(qū)域分類和邊界跟蹤提取后圖示化的結(jié)果，有平原和中山兩種地形類別，品紅色表示中山，綠色表示平原；通過對提取的邊界求面積、中心坐標、地形體積、邊界周長和計算形狀因子，將此中山簡化為錐形山，再求出底面圓半徑和錐形山高度，最終得到錐形山底面圓中心、半徑和圓錐高度幾何參數(shù)，圖7是此孤峰山簡化并按提取的幾何參數(shù)以網(wǎng)格形式繪制的圖形。