曹振宇，譚明建，沈艷敏

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

一、引 言

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

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

二、相關(guān)定義

1.多邊形

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

2.橢球面積

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

圖1 多邊形要素

3.表面積

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

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

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

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

圖2 橢球面上任意多邊形面積計(jì)算

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

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

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

圖3表示的復(fù)雜多邊形要素包含P1和P2兩個(gè)部分，由弧段 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(diǎn)構(gòu)成、圖3中弧段 e包含頂點(diǎn)(V1、V2、V3、V4、V5)。

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

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

四、多邊形表面積計(jì)算方法

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

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

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

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

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

圖6 多邊形表面積計(jì)算流程圖

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

五、試驗(yàn)和結(jié)論

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

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

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

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

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

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

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

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

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

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