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

高 昂，何青松

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

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

高 昂1，何青松2

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

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

投影平面面積；地表表面面積；DEM；邊界；工程應用；地貌特征

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

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

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

1 數據及方法

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

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

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

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

本文使用的數據來自于全國第二次土地利用調查數據中用地類型數據以及研究區(qū)的30m分辨率的DEM。

圖1 研究區(qū)

1.2 方法

1.2.1 基于DEM生成表面積柵格

(1)坐標之間的轉換。為了在三維空間中精確的描述地球表面有個點的位置，通常需要借助某種特定的坐標系統(tǒng)得以完成。常用的坐標系統(tǒng)主要可以分為3種：1)大地坐標系，2)平面投影坐標系， 3)空間直角坐標系。在大地坐標系下，各個點的位置是直接基于經緯度坐標進行表示，點的位置是精確、沒有偏移和變形的；在投影坐標系下，各點的位置是一種平面位置，依據投影方法的不同，投影前后的幾何形狀可能在角度、長度和面積等方面都產生一定的變形和偏移；空間直角坐標系是一種三維坐標系，各點的位置由三維坐標系下的真實的X,Y,Z坐標描述，沒有因為投影帶來的變化。三種坐標系中，大地坐標系和空間直角坐標系表示的點位和幾何形狀是沒有變形的，而基于平面直角坐標系的點位和幾何形狀是有變形或偏移的。

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

(1)

(2)

(3)

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

(4)

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

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

(5)

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

(6)

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

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

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

(1)設i表示當前遍歷的行號、j表示當前遍歷的列號，i初始值為2，j初始值為2,從DEM的第二行第二列開始遍歷。則當前遍歷的柵格像元行列坐標為(i,j)。按照圖2中的方法，分別獲得行列號為(i-1,j-1),(i,j-1)和(i-1,j)柵格像元點，構成相鄰的4個點并分別進行編號。當前點為3號點。

(2)由步驟(1)相鄰的4個點可在平面上得到一個四邊形。按照圖3的示意圖，內插出四邊形四條邊的中點高程值和對角線交點的高程值。圖3中，任意一條邊的中點，如6、7、8、9號點的高程值和經緯度坐標由其邊上兩點的高程值的平均值計算得到。

(3)由步驟(1)和(2)可以得到空間上相鄰的8個空間三角形，并按照圖3的方式進行編號：I、II、III、IV、V、VI、VII、VIII。

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

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

(6)根據8個空間三角形與4個新的四邊形的組成關系，分別匯總三角形的面積得到4個新的四邊形的面積，分別記為S1、S2、S3、S4。

(7)將4個新的四邊形的面積寫入到新的柵格文件的對應行列號中。其中，S1寫入新柵格的行列坐標(2i-1,2j-1)，S2寫入新柵格的行列坐標(2i-1,2j)，S3寫入新柵格的行列坐標(2i,2j)，S4寫入新柵格的行列坐標(2i,2j-1)。重復執(zhí)行步驟(1)～(7)直至遍歷完所有的行列。

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

圖2 相鄰4個DEM柵格點示意圖

圖3 DEM內插原理示意圖

圖4 四邊形構建原理示意圖

1.2.2 矢量與柵格疊加運算

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

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

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

(7)

式中：Sp為多邊形P的面積，為以P0、Pk、Pk+1為頂點構成的三角形的面積。

使用該公式計算的面積頂點坐標順序為逆時針時候為正，順時針時為負。但是無論正負，最終的結果取絕對值即為該多邊形的面積。

圖6 矢量面積法示意圖

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

2 結果

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

3 討論

經過上述三個縣為例進行的算法實驗，可以明顯看出本文提出的柵格與矢量疊加統(tǒng)計算法在統(tǒng)計結果和時間效率上與傳統(tǒng)方法有明顯的不同。

3.1 表面積柵格生成

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

表1 表面積統(tǒng)計結果

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

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

3.3 表面積與投影面積

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

(8)

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

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

4 結論

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

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

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

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

作者貢獻說明：本文為第一作者與通訊作者共同完成。第一作者于2015年12月～2016年08月在通訊作者何青松的指導下進行了相關領域的學習并起草了該文的寫作。通訊作者在原文的基礎上進行了整體的修改。實驗部分(數據預處理、算法的編程以及柵格與矢量統(tǒng)計)由通訊作者完成。

[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] 程乾, 王人潮. 數字高程模型和多時相MODIS數據復合的水稻種植面積遙感估算方法研究[J]. 農業(yè)工程學報, 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] 張偉, 李愛農. 基于SRTM DEM數據的全國地表面積計算研究[J]. 地理與地理信息科學, 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] 謝成磊, 趙榮, 梁勇. 基于地理坐標的地理國情統(tǒng)計單元表面面積精確計算[J]. 遙感信息, 2014(4):47-51.

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

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

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

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

責任編輯 喻曉敏

P3;TP3

A

1003-8078(2016)06-0048-07

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

高昂，男，湖北黃岡人，黃岡中學學生。

何青松，男，安徽合肥人，博士，主要研究方向為制圖與地理信息。