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

背景技術(shù)：

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

步驟113：將平原地形簡(jiǎn)化為平面，求取地形類別標(biāo)志為平原的所有柵格點(diǎn)高程的平均值，記其為平均高程Ea；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

步驟202：遍歷類別n的地形范圍，從QiFuDu.gif文件中讀取相對(duì)應(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中對(duì)待研究的地形區(qū)域進(jìn)行八鄰域邊界跟蹤，提取出邊界點(diǎn)序列存入鏈表Point_List中，包括以下步驟：

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

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

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

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

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

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

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

步驟308：令當(dāng)前邊界點(diǎn)等于此搜索點(diǎn)；

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

步驟310：將當(dāng)前邊界點(diǎn)存入鏈表Point_List中；

步驟311：令Directcode＝Directcode-2；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

附圖說(shuō)明

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

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

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

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

圖5是一座山的實(shí)際形狀；

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

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

具體實(shí)施方式

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

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

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

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

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

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

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

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

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

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

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

步驟109中依據(jù)地形平均起伏度對(duì)步驟107初步分類結(jié)果進(jìn)行地形判定分類，參照?qǐng)D2，包括以下步驟：

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

步驟202：遍歷類別n的地形范圍，從QiFuDu.gif文件中讀取相對(duì)應(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ù)標(biāo)志地形類別及區(qū)域，前一位表示地形類別，平原、丘陵、低山、中山、高山分別用1、2、3、4、5表示；后三位表示此種類別地形的區(qū)域個(gè)數(shù)編號(hào)，從001開始累積編號(hào)，地形判定分類的結(jié)果用此種方式表示成數(shù)值之后，保存在地形分類結(jié)果文件FenLei.gif之中；

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

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

步驟113：將平原地形簡(jiǎn)化為平面，求取地形類別標(biāo)志為平原的所有柵格點(diǎn)高程的平均值，記其為平均高程Ea；

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

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

步驟115中對(duì)待研究的地形區(qū)域進(jìn)行八鄰域邊界跟蹤，提取出邊界點(diǎn)序列存入鏈表Point_List中，參照?qǐng)D3，包括以下步驟：

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

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

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

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

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

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

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

步驟308：令當(dāng)前邊界點(diǎn)等于此搜索點(diǎn)；

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

步驟310：將當(dāng)前邊界點(diǎn)存入鏈表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存儲(chǔ)的邊界點(diǎn)序列連線形成一個(gè)邊界多邊形，用圖形學(xué)中的“交點(diǎn)計(jì)數(shù)法”判斷柵格點(diǎn)是否在邊界多邊形內(nèi)，統(tǒng)計(jì)處在邊界多邊形內(nèi)和位于邊界多邊形邊上的柵格點(diǎn)數(shù)目num，設(shè)單個(gè)柵格點(diǎn)所占面積為a，則地形區(qū)域的面積S＝num×a；地形區(qū)域的中心點(diǎn)C的坐標(biāo)x_i、y_i是邊界多邊形內(nèi)及邊界多邊形邊上的任一柵格點(diǎn)的坐標(biāo)；地形的體積h_i是邊界多邊形內(nèi)及邊界多邊形邊上的任一柵格點(diǎn)的高程；

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

仿真實(shí)例

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