高 昂，何青松

(1. 黃岡中學(xué)，湖北 黃州 438000；2. 武漢大學(xué) 資源與環(huán)境科學(xué)學(xué)院，湖北 武漢 430079)

湖北省地類(lèi)表面積柵格與矢量疊加統(tǒng)計(jì)算法探究

高 昂1，何青松2

(1. 黃岡中學(xué)，湖北 黃州 438000；2. 武漢大學(xué) 資源與環(huán)境科學(xué)學(xué)院，湖北 武漢 430079)

本文通過(guò)DEM生成表面積柵格，然后在表面積柵格的基礎(chǔ)上提出了一種顧及多邊形邊界精確裁剪的柵格與矢量疊加統(tǒng)計(jì)算法CVOA。與傳統(tǒng)的做法相比，在處理時(shí)間上極大縮短，能夠滿(mǎn)足大規(guī)模的工程應(yīng)用。本文以湖北省具有典型平原地貌特征的仙桃市，丘陵地貌的羅田縣以及高山地貌的神農(nóng)架林區(qū)為研究區(qū)，計(jì)算各個(gè)研究區(qū)內(nèi)的各個(gè)地類(lèi)的表面積取得了很好的效果。

投影平面面積；地表表面面積；DEM；邊界；工程應(yīng)用；地貌特征

地表表面面積是指結(jié)合地形條件的地球表面的面積。在地理信息技術(shù)領(lǐng)域，對(duì)地球進(jìn)行建模時(shí)通常將地球簡(jiǎn)化為一個(gè)標(biāo)準(zhǔn)的橢球體，而不考慮地球表面實(shí)際的起伏情況。在資源調(diào)查中，基于橢球面積而不是表面積測(cè)算得到的面積常常被用來(lái)評(píng)價(jià)生物資源、土地資源規(guī)模、生態(tài)服務(wù)價(jià)值等[1-4]。然而不論是人類(lèi)還是自然界中的其他各種生物都生活在有地形起伏的地形上[5-6]。由于地球表面高地起伏較大，導(dǎo)致地表表面面積與橢球面積會(huì)在山區(qū)和丘陵地區(qū)會(huì)存在較大差異。因此，僅僅以地物的橢球面積來(lái)表示地物的面積是不能滿(mǎn)足資源調(diào)查(國(guó)土調(diào)查、林業(yè)調(diào)查)等實(shí)際應(yīng)用的需要[7]。在生物地理學(xué)領(lǐng)域，地形通常被認(rèn)為通過(guò)異質(zhì)性和分離性顯著的增加了生物多樣性[8-9]；文獻(xiàn)[10]在估算水稻種植面積時(shí)注意到了水稻種植易受地表起伏影響的特點(diǎn),在分析時(shí)加入了坡度信息[10]；在風(fēng)化作用中，地表表面積被認(rèn)為是至關(guān)重要的決定因素[11]。

隨著GIS技術(shù)以及Digital Elevation Models(DEM)模型發(fā)展，衍生出越來(lái)越多的基于DEM求算地表表面積的算法[12-16]。然而這些算法多只存在于理論或小尺度的應(yīng)用上?；贒EM的統(tǒng)計(jì)單元表面面積統(tǒng)計(jì)主要包括以下4種：(1)利用DEM構(gòu)建區(qū)域的不規(guī)則三角網(wǎng)(TIN)[17-18]，在此基礎(chǔ)上疊加需要計(jì)算表面面積的矢量多邊形，由不規(guī)則三角網(wǎng)與多邊形的疊加計(jì)算得到各多邊形的表面面積[19-20]； (2)直接基于數(shù)字高程模型數(shù)據(jù)構(gòu)建規(guī)則三角網(wǎng)，在此基礎(chǔ)上疊加需要計(jì)算表面面積的矢量多邊形，由規(guī)則三角網(wǎng)與多邊形的空間位置關(guān)系，匯總落入多邊形內(nèi)的空間三角形的面積得到各多邊形的表面面積[16]；(3)基于相鄰的4個(gè)DEM柵格像元，利用積分公式推算出由4個(gè)相鄰像元構(gòu)成的區(qū)域的曲面面積，進(jìn)而依據(jù)該曲面與多邊形的空間位置關(guān)系，匯總落入多邊形內(nèi)的曲面面積，從而得到多邊形的表面面積[21]。(4)利用DEM上3*3的像元窗口構(gòu)建8個(gè)空間三角形，計(jì)算3*3區(qū)域內(nèi)的中心像元區(qū)域的表面面積，進(jìn)而利用柵格像元與多邊形進(jìn)行疊加，匯總落入多邊形內(nèi)柵格像元值作為多邊形的表面面積[7]。然而上述相關(guān)算法模型都沒(méi)有經(jīng)過(guò)大規(guī)模工程應(yīng)用的檢驗(yàn)，且回避了怎樣計(jì)算目標(biāo)區(qū)域邊界附近面積的問(wèn)題。

中國(guó)幅員廣闊，橢球投影面積就約有960萬(wàn)平方千米，且中國(guó)是個(gè)多山地地貌的國(guó)家，山地面積占到全國(guó)總面積的2/3以上[22]。國(guó)務(wù)院第一次全國(guó)地理國(guó)情普查領(lǐng)導(dǎo)小組辦公室規(guī)定縣級(jí)行政單元各個(gè)用地類(lèi)型的表面積統(tǒng)計(jì)是國(guó)家的一項(xiàng)關(guān)鍵的基本國(guó)情統(tǒng)計(jì)指標(biāo)。 因此如何快速且準(zhǔn)確的計(jì)算出各個(gè)縣級(jí)行政單元各個(gè)地類(lèi)的表面面積對(duì)科研人員提出了巨大的挑戰(zhàn)。本文正是作者及其合作者在思考解決這個(gè)工程級(jí)應(yīng)用問(wèn)題時(shí)的一些思考和實(shí)驗(yàn)嘗試。本文主要解決兩個(gè)問(wèn)題：(1)在DEM基礎(chǔ)上如何高效率、高精度的獲取表面積柵格文件；(2)柵格與矢量疊加統(tǒng)計(jì)時(shí)如何在顧及矢量多邊形邊界基礎(chǔ)高效率的計(jì)算多邊形的表面積。

1 數(shù)據(jù)及方法

1.1 研究區(qū)及使用數(shù)據(jù)

仙桃市:位于湖北中南部，江漢平原腹地。地勢(shì)平坦，起伏甚微，是湖北省江漢平原中心城市, 面積2 538平方千米,見(jiàn)圖1(a)。

羅田縣：位于湖北省東北部、大別山南麓，地處東經(jīng)115°06′至115°46′，北緯30°35′至31°16′之間。境內(nèi)多200米以下小山丘，丘陵廣布。面積2 145平方千米，見(jiàn)圖1(b)。

神農(nóng)架林區(qū)：位于湖北省西部。境內(nèi)多高山，以神農(nóng)頂最高，高程為3 105.4米，是華中地區(qū)最高點(diǎn)，呈現(xiàn)典型的山地地貌。面積3 253平方千米，見(jiàn)圖1(c)。

本文使用的數(shù)據(jù)來(lái)自于全國(guó)第二次土地利用調(diào)查數(shù)據(jù)中用地類(lèi)型數(shù)據(jù)以及研究區(qū)的30m分辨率的DEM。

圖1 研究區(qū)

1.2 方法

1.2.1 基于DEM生成表面積柵格

(1)坐標(biāo)之間的轉(zhuǎn)換。為了在三維空間中精確的描述地球表面有個(gè)點(diǎn)的位置，通常需要借助某種特定的坐標(biāo)系統(tǒng)得以完成。常用的坐標(biāo)系統(tǒng)主要可以分為3種：1)大地坐標(biāo)系，2)平面投影坐標(biāo)系， 3)空間直角坐標(biāo)系。在大地坐標(biāo)系下，各個(gè)點(diǎn)的位置是直接基于經(jīng)緯度坐標(biāo)進(jìn)行表示，點(diǎn)的位置是精確、沒(méi)有偏移和變形的；在投影坐標(biāo)系下，各點(diǎn)的位置是一種平面位置，依據(jù)投影方法的不同，投影前后的幾何形狀可能在角度、長(zhǎng)度和面積等方面都產(chǎn)生一定的變形和偏移；空間直角坐標(biāo)系是一種三維坐標(biāo)系，各點(diǎn)的位置由三維坐標(biāo)系下的真實(shí)的X,Y,Z坐標(biāo)描述，沒(méi)有因?yàn)橥队皫?lái)的變化。三種坐標(biāo)系中，大地坐標(biāo)系和空間直角坐標(biāo)系表示的點(diǎn)位和幾何形狀是沒(méi)有變形的，而基于平面直角坐標(biāo)系的點(diǎn)位和幾何形狀是有變形或偏移的。

在計(jì)算三角形面積時(shí)，為了避免平面投影產(chǎn)生的變形，需要以空間直角坐標(biāo)系作為三角形頂點(diǎn)坐標(biāo)存儲(chǔ)單位。在大地坐標(biāo)系下任意一點(diǎn)(B,L,H)可以采用以下計(jì)算公式轉(zhuǎn)換為空間直角坐標(biāo)系下的(X,Y,Z)，計(jì)算公式如下：

(1)

(2)

(3)

式中：B,L,H分別表示某點(diǎn)的緯度、經(jīng)度和高程值；X,Y,Z分別表示該點(diǎn)在三維空間直角坐標(biāo)系中的坐標(biāo)值；Cos(B)為緯度的余弦值，Sin(B)是緯度的正弦值，Sin(L)是經(jīng)度的正弦值，Cos(L)是經(jīng)度的余弦值。e2為地球橢球體的第一偏心率的平方，作為地球橢球體的已知常數(shù)輸入。N的計(jì)算公式如公式(4)：

(4)

式中：a為地球橢球體的長(zhǎng)半軸。a和e2均為常數(shù)，可在地球參考橢球體的定義中查找得到。

(2)地表表面面積計(jì)算原理。計(jì)算地表表面面積的前提是對(duì)地表表面進(jìn)行數(shù)值模擬和逼近。數(shù)字高程模型DEM，通常利用柵格圖像來(lái)表示地球表面高程(也稱(chēng)為“海拔”)的高低起伏情況。DEM以離散的、規(guī)則的網(wǎng)格對(duì)地表進(jìn)行了模擬[7]。然而，地表是一個(gè)連續(xù)的表面，基于DEM無(wú)法直接得到連續(xù)的表面。因此，在計(jì)算地表表面時(shí)，通常基于DEM上離散的柵格點(diǎn)構(gòu)建相互連接的三角網(wǎng)(在空間直角坐標(biāo)系中，任意3個(gè)點(diǎn)都能保證在同一個(gè)平面上，而3個(gè)以上的點(diǎn)無(wú)法保證在同一個(gè)平面上，也無(wú)法計(jì)算面積)。以相互連接的多個(gè)三角形構(gòu)成的三角網(wǎng)來(lái)逼近地表真實(shí)的連續(xù)表面[19]。給定空間三角形頂點(diǎn)的A,B，C及其坐標(biāo)(XA,YA,ZA)、(XB,YB,ZB)和(XC,YC,ZC)，其面積S可利用海倫公式(5)得出：

(5)

式中：p=(a+b+c)/2。a,b,c分別表示三角形的邊長(zhǎng)?？臻g三角形的邊長(zhǎng)可以采用下式公式(6)計(jì)算得到：

(6)

式中：lij表示任意兩點(diǎn)i,j的三維空間距離，(xi,yi,zi)和(xj,yj,zj)分別表示兩點(diǎn)的三維坐標(biāo)。

設(shè)DEM采用大地坐標(biāo)系進(jìn)行存儲(chǔ)，柵格像元的行數(shù)為M,列數(shù)為N，柵格分辨率為C。則生成表面積柵格的具體步驟如下： 首先創(chuàng)建一個(gè)(M-1)*2行、(N-1)*2列的空的柵格文件，柵格文件中各柵格像元的值默認(rèn)為0。設(shè)原始DEM的起始點(diǎn)經(jīng)緯度坐標(biāo)為(B1,L1)，則新柵格文件的起始點(diǎn)經(jīng)緯度坐標(biāo)為(B1+C,L1+C)。

接著按行遍歷DEM柵格，按照以下方法構(gòu)建規(guī)則三角形并計(jì)算柵格像元區(qū)域的表面面積：

(1)設(shè)i表示當(dāng)前遍歷的行號(hào)、j表示當(dāng)前遍歷的列號(hào)，i初始值為2，j初始值為2,從DEM的第二行第二列開(kāi)始遍歷。則當(dāng)前遍歷的柵格像元行列坐標(biāo)為(i,j)。按照?qǐng)D2中的方法，分別獲得行列號(hào)為(i-1,j-1),(i,j-1)和(i-1,j)柵格像元點(diǎn)，構(gòu)成相鄰的4個(gè)點(diǎn)并分別進(jìn)行編號(hào)。當(dāng)前點(diǎn)為3號(hào)點(diǎn)。

(2)由步驟(1)相鄰的4個(gè)點(diǎn)可在平面上得到一個(gè)四邊形。按照?qǐng)D3的示意圖，內(nèi)插出四邊形四條邊的中點(diǎn)高程值和對(duì)角線(xiàn)交點(diǎn)的高程值。圖3中，任意一條邊的中點(diǎn)，如6、7、8、9號(hào)點(diǎn)的高程值和經(jīng)緯度坐標(biāo)由其邊上兩點(diǎn)的高程值的平均值計(jì)算得到。

(3)由步驟(1)和(2)可以得到空間上相鄰的8個(gè)空間三角形，并按照?qǐng)D3的方式進(jìn)行編號(hào)：I、II、III、IV、V、VI、VII、VIII。

(4)利用海倫公式分別計(jì)算8個(gè)空間三角形的面積。

(5)由三角形I、II可得到一個(gè)新的平面四邊形編號(hào)為(一)，三角形III、IV可得到一個(gè)新的平面四邊形編號(hào)為(二)；三角形V、VI可得到一個(gè)新的平面四邊形編號(hào)為(三)；三角形VII、VIII可得到一個(gè)新的平面四邊形編號(hào)為(四)，如圖4所示。

(6)根據(jù)8個(gè)空間三角形與4個(gè)新的四邊形的組成關(guān)系，分別匯總?cè)切蔚拿娣e得到4個(gè)新的四邊形的面積，分別記為S1、S2、S3、S4。

(7)將4個(gè)新的四邊形的面積寫(xiě)入到新的柵格文件的對(duì)應(yīng)行列號(hào)中。其中，S1寫(xiě)入新柵格的行列坐標(biāo)(2i-1,2j-1)，S2寫(xiě)入新柵格的行列坐標(biāo)(2i-1,2j)，S3寫(xiě)入新柵格的行列坐標(biāo)(2i,2j)，S4寫(xiě)入新柵格的行列坐標(biāo)(2i,2j-1)。重復(fù)執(zhí)行步驟(1)～(7)直至遍歷完所有的行列。

由圖4可知，新的柵格文件的空間分辨率是原來(lái)的一半，空間精度更高。新柵格文件中，柵格像元的值表示了該區(qū)像元占據(jù)空間區(qū)域的地表表面面積，該柵格文件定義為數(shù)字地表面積模型，標(biāo)記為SAM。

圖2 相鄰4個(gè)DEM柵格點(diǎn)示意圖

圖3 DEM內(nèi)插原理示意圖

圖4 四邊形構(gòu)建原理示意圖

1.2.2 矢量與柵格疊加運(yùn)算

柵格與矢量空間疊加時(shí)(見(jiàn)圖5)，只有柵格像元的中心點(diǎn)落在多邊形內(nèi)時(shí)，該像元被認(rèn)為落在多邊形內(nèi)，否則不在多邊形內(nèi)。如圖5所示，一個(gè)矩形表示一個(gè)柵格像元，矩形內(nèi)的圓點(diǎn)表示柵格像元的中心點(diǎn)。在基于柵格做空間運(yùn)算時(shí)，通常以柵格像元的中心點(diǎn)代替整個(gè)柵格以簡(jiǎn)化計(jì)算提高計(jì)算效率。圖5(a)中的多邊形的面積在ArcGIS中統(tǒng)計(jì)的值為圖5(b)中灰色填充柵格值的求和。不難發(fā)現(xiàn)這個(gè)值與該多邊形的真實(shí)面積存在誤差。當(dāng)一種地類(lèi)由很多這樣的多邊形構(gòu)成時(shí)累積的誤差更大。ArcGIS這種矢量柵格疊加統(tǒng)計(jì)方式對(duì)于摸清地類(lèi)表面積的大概十分有用，然后卻無(wú)法精確的掌握類(lèi)型的表面積。很明顯該方法的精度與DEM的空間分辨率有關(guān)，DEM分辨率越低,面積統(tǒng)計(jì)精度也越低。這種計(jì)算方法在處理比較小的圖斑的時(shí)候會(huì)對(duì)結(jié)果產(chǎn)生較大的影響。

圖5 矢量柵格疊加統(tǒng)計(jì)示意圖

對(duì)于只有部分覆蓋的柵格，想直接精確求出多邊形與該柵格相交部分的表面面積是不可能的。但是通過(guò)求該相交部分的平面面積與該柵格的平面面積得到一個(gè)比值，以該比值作為該相交部分應(yīng)該分配到的當(dāng)前柵格值的比例是個(gè)可行的辦法，他能最大程度的保證表面積的精確程度。因此此時(shí)的問(wèn)題轉(zhuǎn)到了如何快速求多邊形與所有相交柵格相交部分的平面面積。矢量數(shù)據(jù)中的多邊形由各個(gè)弧段組成，弧段則由兩點(diǎn)連線(xiàn)構(gòu)成。因此矢量多邊形可以看成是由多個(gè)點(diǎn)連線(xiàn)構(gòu)成。任意多邊形的平面面積可由任意一點(diǎn)與多邊形上依次兩點(diǎn)連線(xiàn)構(gòu)成的三角形矢量面積求和得出。而矢量面積等于三角形兩邊矢量的叉乘[23]。假設(shè)構(gòu)成多邊形P的點(diǎn)坐標(biāo)依次為(x1,y1),(x2,y2),…,(xm,ym)(圖6)。則該多邊形P的面積可由公式(7)計(jì)算：

(7)

式中：Sp為多邊形P的面積，為以P0、Pk、Pk+1為頂點(diǎn)構(gòu)成的三角形的面積。

使用該公式計(jì)算的面積頂點(diǎn)坐標(biāo)順序?yàn)槟鏁r(shí)針時(shí)候?yàn)檎槙r(shí)針時(shí)為負(fù)。但是無(wú)論正負(fù)，最終的結(jié)果取絕對(duì)值即為該多邊形的面積。

圖6 矢量面積法示意圖

計(jì)算速度是柵格與矢量疊加統(tǒng)計(jì)時(shí)需要解決的一個(gè)重要問(wèn)題。在ArcGIS實(shí)現(xiàn)中通過(guò)“掃面線(xiàn)算法(Sweep Line Algorithm)，逐個(gè)像元進(jìn)行掃描判斷柵格與多邊形的交點(diǎn)。當(dāng)像元數(shù)量以及多邊形的數(shù)量較多的情況下算法時(shí)間過(guò)長(zhǎng)，不適合大規(guī)模的工程應(yīng)用。不難發(fā)現(xiàn)盡可能的減少空間操作次數(shù)是提高統(tǒng)計(jì)速度的一個(gè)有效捷徑。因此本文作者提出了顧及邊界精確裁剪的柵格與矢量疊加統(tǒng)計(jì)算法 (CVOA)，它通過(guò)對(duì)表面積柵格進(jìn)行整行的邊框提取構(gòu)建矩形，利用該矩形裁切多邊形，得到多邊形與矩形的相交區(qū)域。我們已經(jīng)知道矢量多邊形其實(shí)是由點(diǎn)構(gòu)成，通過(guò)點(diǎn)的坐標(biāo)就可以知道點(diǎn)所在柵格中的具體位置。將一個(gè)柵格內(nèi)的所有點(diǎn)圍成的多邊形按照公式7求出平面面積，再將該面積與該柵格的平面面積得到一個(gè)比值，以該比值作為該相交部分應(yīng)該分配到的當(dāng)前柵格表面積值的比例就可以知道當(dāng)前相交部分的表面面積。

2 結(jié)果

統(tǒng)計(jì)的表面積結(jié)果如表1所示。其中表頭中Area 表示平面投影面積；ArcGIS_SA 表示用ArcGIS統(tǒng)計(jì)得到的表面積；CVOA_SA表示CVOA算法計(jì)算得到的表面積； D1=ArcGIS_SA-Area；D2=CVOA_SA-Area；D3=CVOA_SA-ArcGIS_SA。

3 討論

經(jīng)過(guò)上述三個(gè)縣為例進(jìn)行的算法實(shí)驗(yàn)，可以明顯看出本文提出的柵格與矢量疊加統(tǒng)計(jì)算法在統(tǒng)計(jì)結(jié)果和時(shí)間效率上與傳統(tǒng)方法有明顯的不同。

3.1 表面積柵格生成

本文提出的表面積柵格生成方法主要優(yōu)勢(shì)在于：(1)直接基于大地坐標(biāo)系下的DEM進(jìn)行計(jì)算，避免了構(gòu)建不規(guī)則三角網(wǎng)時(shí)需要進(jìn)行平面投影和DEM柵格點(diǎn)抽稀的精度損失；(2)避免了矢量多邊形和規(guī)則三角網(wǎng)(也是矢量)疊加時(shí)產(chǎn)生的巨大計(jì)算量，效率優(yōu)勢(shì)明顯。(3)計(jì)算簡(jiǎn)單，計(jì)算量小，無(wú)需復(fù)雜的積分運(yùn)算，適合大規(guī)模計(jì)算任務(wù)的求解。(4)生成得到的新的柵格文件的空間分辨率是原有DEM的一半，使得矢量與柵格疊加運(yùn)算時(shí)的空間精度得到進(jìn)一步提升。以(10*10)平方米分辨率為例，按照文獻(xiàn)[7]進(jìn)行柵格判斷時(shí)，每個(gè)面積柵格的數(shù)值至少為100平方米；按照本方法進(jìn)行判斷時(shí)，原有柵格被一分為四，每個(gè)柵格的數(shù)值至少為25平方米。柵格劃分越細(xì)則越能精確的描述幾何對(duì)象的邊界。

表1 表面積統(tǒng)計(jì)結(jié)果

3.2 CVOA在柵格與矢量疊加統(tǒng)計(jì)時(shí)的優(yōu)勢(shì)

CVOA顧及了柵格與矢量多邊形疊加統(tǒng)計(jì)時(shí)多邊形邊界對(duì)柵格的精確裁剪，它解決了傳統(tǒng)GIS軟件中對(duì)于柵格與多邊形部分相交時(shí)柵格要么完全屬于要么完全不屬于多邊形的“二值化”弊端，如果柵格與多邊形數(shù)量越多則這種“二值化”處理的誤差累積也有可能越來(lái)越大，使得最終計(jì)算的表面面積結(jié)果與實(shí)際數(shù)據(jù)之間存在相當(dāng)大的差異。從結(jié)果中也看出在仙桃市沙地地類(lèi)統(tǒng)計(jì)時(shí)甚至出現(xiàn)了表面積比平面投影面積少的情況，而在DEM為正的情況下這種情況理論上是不可能存在的。此外CVOA盡量減少了空間相交的次數(shù)，它用柵格的一整行的外接矩形框(也是矢量)，而非逐個(gè)柵格邊框與多邊形相交，再按照相交多邊形各個(gè)節(jié)點(diǎn)的坐標(biāo)推算點(diǎn)在柵格中的位置信息。與逐像元相交動(dòng)輒幾個(gè)小時(shí)的處理時(shí)間相比，CVOA完成統(tǒng)計(jì)的時(shí)間均在10分鐘以?xún)?nèi)。由此可見(jiàn)CVOA適合于大規(guī)模的工程應(yīng)用中。

3.3 表面積與投影面積

本文選取的三個(gè)地區(qū)仙桃市、羅田縣以及神農(nóng)架林區(qū)代表了三類(lèi)典型的地形地貌區(qū)：平原、丘陵和山地。我們定義一個(gè)平整系數(shù)來(lái)評(píng)估表面積與平面投影面積之間的相對(duì)差異,如公式(8):

(8)

式中：fc為平整系數(shù)，SA為表面面積，PA為投影平面面積。

計(jì)算表明仙桃市、羅田縣以及神農(nóng)架林區(qū)的平整系數(shù)分別是1.004 5,1.033 4,1.115 0。由此可知地形起伏越大的地方地表真實(shí)面積與投影面積的相對(duì)差異越大。從第三部分的三個(gè)縣的統(tǒng)計(jì)表格中D2結(jié)果來(lái)看，絕對(duì)數(shù)量上，三個(gè)縣的面積差異也分別達(dá)到了11.43 km2, 71.27 km2, 371.58 km2。這樣的差異對(duì)于評(píng)估研究區(qū)的固碳能力、生物資源、土地資源規(guī)模帶來(lái)了較大的誤差。換言之,在遙感數(shù)據(jù)的面積概念選擇中,不應(yīng)一味地使用垂直投影面積,而應(yīng)根據(jù)自己所需,在條件允許的情況下適當(dāng)選擇地表真實(shí)面積的概念,否則將會(huì)引起較大的誤差,甚至得出錯(cuò)誤的結(jié)論[24]。

4 結(jié)論

本文中，作者提出了一種顧及邊界精確裁剪的柵格與矢量疊加統(tǒng)計(jì)算法(CVOA)并將它應(yīng)用在了湖北省不同地貌類(lèi)型區(qū)縣的地類(lèi)表面積統(tǒng)計(jì)中。作者從兩方面來(lái)描述CVOA:基于DEM構(gòu)建表面積柵格以及用生成的表面積柵格與用地類(lèi)型矢量進(jìn)行疊加統(tǒng)計(jì)。該文主要得到以下三個(gè)結(jié)論：

(1)在不同的地貌區(qū)，投影面積并不能替代真實(shí)的表面積，兩者之間存在著較大的差異，且隨著地形起伏加大，這種差異也變得越來(lái)越大。

(2)柵格與矢量疊加統(tǒng)計(jì)時(shí)，CVOA保證了矢量多邊形邊界對(duì)柵格像元的精確“裁剪”，統(tǒng)計(jì)結(jié)果更真實(shí)可靠。同時(shí)CVOA極大地減少了空間相交的次數(shù)，計(jì)算高效，適合于大規(guī)模的工程應(yīng)用中。

(3)真實(shí)地表面積具有廣泛的應(yīng)用前景。在區(qū)域資源評(píng)價(jià)(如林業(yè)資源評(píng)估，糧食產(chǎn)量預(yù)測(cè))以及生態(tài)地理學(xué)研究(如固碳能力評(píng)價(jià)、生態(tài)服務(wù)價(jià)值計(jì)算)中，這些都需要表面積這一基礎(chǔ)數(shù)據(jù)的支持。

作者貢獻(xiàn)說(shuō)明：本文為第一作者與通訊作者共同完成。第一作者于2015年12月～2016年08月在通訊作者何青松的指導(dǎo)下進(jìn)行了相關(guān)領(lǐng)域的學(xué)習(xí)并起草了該文的寫(xiě)作。通訊作者在原文的基礎(chǔ)上進(jìn)行了整體的修改。實(shí)驗(yàn)部分(數(shù)據(jù)預(yù)處理、算法的編程以及柵格與矢量統(tǒng)計(jì))由通訊作者完成。

[1] Hoechstetter S, Thinh N X, Walz U. 3D-Indices for the Analysis of Spatial Patterns of Landscape Structure[M]. Intercarto-Intergis. 2006.

[2] Gao X, Chun L U, Yun L,etal. Ecological assets valuation of the Tibetan Plateau[J]. Journal of Natural Resources, 2003，18(2):189-196.

[3] Li J, Wang W, Hu G,etal. Changes in ecosystem service values in Zoige Plateau, China[J]. Agriculture Ecosystems & Environment, 2010, 139(4):766-770.

[4] Liu Y, Li J, Zhang H. An ecosystem service valuation of land use change in Taiyuan City, China[J]. Ecological Modelling, 2012, 225(3):127-132.

[5] Zhou T, Chen B M, Liu G,etal. Biodiversity of Jinggangshan Mountain: the importance of topography and geographical location in supporting higher biodiversity[J]. Plos One, 2015, 10(3)：e0120208.

[6] Eisenlohr P V, Alves L F, Bernacci L C,etal. Disturbances, elevation, topography and spatial proximity drive vegetation patterns along an altitudinal gradient of a top biodiversity hotspot[J]. Biodiversity & Conservation, 2013, 22(12):2767-2783.

[7] Jenness J S. Calculating landscape surface area from digital elevation models[J]. Wildlife Society Bulletin, 2004, 32(3), 829-839.

[8] Wakelyn L A. Changing habitat conditions on bighorn sheep ranges in Colorado[J]. Journal of Wildlife Management, 1987, 51(51):904-912.

[9] Rosenzweig M L. Species diversity in space and time. Species diversity in space and time[M]. Cambridge University Press, 1995.

[10] 程乾, 王人潮. 數(shù)字高程模型和多時(shí)相MODIS數(shù)據(jù)復(fù)合的水稻種植面積遙感估算方法研究[J]. 農(nóng)業(yè)工程學(xué)報(bào), 2005, 21(5):89-92.

[11] Gregory A. Pope. Interal weathering in quartz grains[J]. Physical Geography, 2013, 16(4):315-338.

[12] Beasom S L, Wiggers E P, Giardino J R. A technique for assessing land surface ruggedness[J]. Journal of Wildlife Management, 1983, 47(4):1163-1166.

[13] Hodgson M E. What cell size does the computed slope/aspect angle represent?[J]. Photogrammetric engineering and remote sensing, 1995, 61:513-517.

[14] Kundu S N, Pradhan B. Surface area processing in GIS. GIS Development. 2004. (https://www.researchgate.net/publication/233390687_Surface_area_processing_in_GIS).

[15] Zhang Y, Zhang L N, Yang C D,etal. Surface area processing in GIS for different mountain regions[J]. Forestry Studies in China, 2011, 13(4):311-314.

[16] Berry J K. Use surface area for realistic calculations[J]. Geoworld, 2002, 9(2):11-14.

[17] Song G, Chen Y, Zhou Y. Algorithm for Estimating Mountainous Surface Area Based on GPS Data. Proceedings of the 2nd International Conference on Green Communications and Networks 2012 (GCN 2012): Volume 2[R]. Springer Berlin Heidelberg, 2013:679-687.

[18] 張偉, 李?lèi)?ài)農(nóng). 基于SRTM DEM數(shù)據(jù)的全國(guó)地表面積計(jì)算研究[J]. 地理與地理信息科學(xué), 2014, 30(3):51-55.

[19] Wang K, Lo C P. An assessment of the accuracy of triangulated irregular networks (tins) and lattices in arc/info[J]. Transactions in GIS,1999, 3(2), 161-174.

[20] 謝成磊, 趙榮, 梁勇. 基于地理坐標(biāo)的地理國(guó)情統(tǒng)計(jì)單元表面面積精確計(jì)算[J]. 遙感信息, 2014(4):47-51.

[21] 江帆, 呂曉華, 王仲蘭. 基于復(fù)化公式的DEM表面積算法分析[J]. 測(cè)繪學(xué)院學(xué)報(bào), 2005, 22(4):263-265.

[22] 伍光和. 自然地理學(xué)[M].北京：高等教育出版社，2000.

[23] O'Rourke, J. Computational Geometry in C. Computational geometry in C[M]. Cambridge University Press, 1998.

[24] 李文航, 龔建華. 遙感地表真實(shí)面積和垂直投影面積差異的定量化模擬[J]. 計(jì)算機(jī)應(yīng)用研究, 2008, 25(4):983-985.

責(zé)任編輯 喻曉敏

P3;TP3

A

1003-8078(2016)06-0048-07

2016-09-19 doi 10.3969/j.issn.1003-8078.2016.06.14

高昂，男，湖北黃岡人，黃岡中學(xué)學(xué)生。

何青松，男，安徽合肥人，博士，主要研究方向?yàn)橹茍D與地理信息。