大地坐標(biāo)系中多邊形面積計算方法研究

曹振宇，譚明建，沈艷敏

(1.武漢大學(xué)測繪遙感信息工程國家重點實驗室，湖北武漢430079;2.四川省基礎(chǔ)地理信息中心，四川成都610041)

一、引 言

地理國情監(jiān)測是測繪工作者面臨的重要任務(wù)和責(zé)任工程，是地理信息更深更廣層次的應(yīng)用。重要地理信息是重要的人文或自然地理實體的位置、高程、面積、長度等信息。從地理信息數(shù)據(jù)挖掘、提取重要地理信息已經(jīng)成為研究與應(yīng)用的熱點領(lǐng)域［1］。其中，面狀分布的地理要素在GIS中常用多邊形來表示。本文系統(tǒng)地研究了地理坐標(biāo)系中多邊形要素橢球面積和表面積計算的方法。

目前商業(yè)GIS(如ArcGIS)軟件提供的多邊形面積統(tǒng)計功能一般是針對平面投影坐標(biāo)系數(shù)據(jù)，不適用于大地坐標(biāo)系數(shù)據(jù)多邊形面積計算?，F(xiàn)有的多邊形橢球面積計算方法一般基于辛普森積分，依次按照多邊形頂點順序?qū)⒍噙呅芜?、通過多邊形邊端點的緯線以及給定經(jīng)線組成梯形的面積，然后求其代數(shù)和［2-3］。尚沒有文獻(xiàn)討論復(fù)雜多邊形如有島多邊形的橢球面積計算。多邊形表面積一般基于格網(wǎng)DEM表面面積計算，采用多邊形包含柵格單元面積進(jìn)行累加得到，但大家一般只討論了投影平面直角坐標(biāo)系中表面積的具體計算方法。本文首先討論了大地坐標(biāo)系中復(fù)雜多邊形橢球面積計算方法，然后討論了大地坐標(biāo)系多邊形表面積計算方法。

二、相關(guān)定義

1.多邊形

多邊形由首尾相連的弧段構(gòu)成。如圖1所示，多邊形包含內(nèi)外邊界弧段，按左手法則，其中邊界弧段的正向?qū)?yīng)多邊形內(nèi)部，負(fù)向?qū)?yīng)的區(qū)域成為島。圖中P1為多邊形內(nèi)部，P2對應(yīng)的區(qū)域即為島。

2.橢球面積

橢球面積是指多邊形在參考橢球面上的面積。

圖1 多邊形要素

3.表面積

本文表面積是指多邊形在地形表面上的曲面面積。

三、復(fù)雜多邊形橢球面積計算方法

1.基于辛普森積分多邊形橢球面積計算方法

如圖2所示，多邊形ABCD的面積就等于4個梯形圖塊(ABB1A1、BCC1B1、CDD1C1、DAA1D1)面積的代數(shù)和。L0為任意經(jīng)線。任意多邊形面積計算按照多邊形頂點順序逐個計算相鄰兩點連線與任意經(jīng)線構(gòu)成的梯形面積Si然后累加。頂點按逆時針順序得到的結(jié)果為正，按順時針順序得到的結(jié)果為負(fù)。

圖2 橢球面上任意多邊形面積計算

橢球面上任一梯形圖塊面積計算公式為

式中，A、B、C、D、E 為常數(shù)。

2.復(fù)雜多邊形要素橢球面積計算方法

圖3表示的復(fù)雜多邊形要素包含P1和P2兩個部分，由弧段 a、b、c、d、e、f組成。a、e 為外邊界弧段，稱為外環(huán)。b、c、d為a包含的內(nèi)邊界弧段，稱為內(nèi)環(huán)。圖3中f為e包含的內(nèi)環(huán)?；《斡梢幌盗许旤c構(gòu)成、圖3中弧段 e包含頂點(V1、V2、V3、V4、V5)。

圖3 復(fù)雜多邊形區(qū)域

對圖3所示的復(fù)雜多邊形要素，采用圖4所示的流程，首先遍歷多邊形要素的外環(huán)，按頂點順時針順序利用式(2)計算外環(huán)面積Si。然后遍歷當(dāng)前外環(huán)包含的內(nèi)環(huán)，按頂點逆時針順序利用式(2)計算內(nèi)環(huán)包含的面積并求和則當(dāng)前外環(huán)所在區(qū)域的面積為依次計算每個外環(huán)所在區(qū)域的面積累加即為多邊形要素的橢球面積。

四、多邊形表面積計算方法

多邊形表面積是指多邊形在地形表面上的空間曲面面積。如圖5所示，規(guī)則格網(wǎng)表示多邊形所在區(qū)域的DEM，多邊形表面積S可以表示為以數(shù)字高程模型DEM、多邊形Polygon為參數(shù)的函數(shù)：S=f(polygon，DEM)。

圖4 復(fù)雜多邊形橢球面積計算流程圖

圖5 基于格網(wǎng)DEM計算多邊形表面積

多邊形P的表面積計算流程如圖6所示。先求P的最小外切矩形，然后將P的最小外切矩形按d x，d y大小分為n×m個格網(wǎng)，從DEM數(shù)據(jù)查詢每個格網(wǎng)節(jié)點高程。然后遍歷格網(wǎng)中所有矩形，判斷矩形和P的位置關(guān)系，格網(wǎng)矩形和多邊形的關(guān)系有分離、包含和相交3種。對包含在多邊形P的矩形格網(wǎng)分割成兩個三角形，按式(4)計算矩形格網(wǎng)的表面積。

和多邊形相交的矩形格網(wǎng)先計算矩形格網(wǎng)的粗糙度 CZ=S表面積/S橢球面積，然后根據(jù)文中復(fù)雜多邊形橢球面積的計算方法計算相交部分的橢球面積S相交橢球面積，則相交部分的表面積為：S=CZ×S相交橢球面積。

圖6 多邊形表面積計算流程圖

多邊形表面積是所有包含在多邊形中格網(wǎng)表面積和所有與多邊形相交的格網(wǎng)相交部分表面積的累加。

五、試驗和結(jié)論

按照本文提出的橢球面積和表面積的計算方法，利用1∶50 000數(shù)字高程模型和行政境界數(shù)據(jù)，計算了縣級行政區(qū)域的橢球面積和表面積。試驗結(jié)果見表1。

試驗結(jié)果表明，本文提出的多邊形橢球面積和表面積計算方法，滿足任意復(fù)雜多邊形要素橢球面積和表面積計算的需要，并可以由此計算出如地形粗糙度等地形因子，因此可以用于重要地理信息統(tǒng)計工作。

表1 縣級行政區(qū)域面積計算結(jié)果

［1］蔣捷，黃蔚，王茜，等.基于數(shù)字地形圖數(shù)據(jù)的重要地理信息統(tǒng)計分析方法初探［J］.測繪科學(xué)，2008，33(S1)：93-95.

［2］衛(wèi)東.山西省重要地理信息統(tǒng)計分析系統(tǒng)設(shè)計與實現(xiàn)［J］.華北國土資源，2010(4)：57-60.

［3］陳婭婷，嚴(yán)泰來，朱德海.基于辛普森面積的多邊形凹凸性識別算法［J］.地理與地理信息科學(xué)，2010，26(6)：28-30.

［4］王秀云，陳曄，周厚華，等.DEM在林地資源表面積調(diào)查中的應(yīng)用［J］.南京師范大學(xué)學(xué)報：工程技術(shù)版，2006，6(1)：10-16.

［5］趙強，胡娟.MAPGIS軟件中土地利用現(xiàn)狀面積量算方法探討［J］.國土資源導(dǎo)刊：湖南，2007，4(5)：58-60.

［6］李志林，朱慶.數(shù)字高程模型［M］.武漢：武漢大學(xué)出版社，2007.

［7］吳立新，史文中.地理信息系統(tǒng)原理與算法［M］.北京：科學(xué)出版社，2003.