謝成磊，趙榮，梁勇

(1.山東農(nóng)業(yè)大學(xué) 信息科學(xué)與工程學(xué)院，山東 泰安 271018；2.中國測(cè)繪科學(xué)研究院，北京 100830)

1 引 言

基本統(tǒng)計(jì)是地理國情普查工作的重要內(nèi)容，其根據(jù)地理國情普查采集的幾何特征類型和地理實(shí)體對(duì)象，對(duì)統(tǒng)計(jì)單元內(nèi)地形地貌、植被覆蓋、荒漠與裸露地表、水域、交通網(wǎng)絡(luò)、居民地與設(shè)施等內(nèi)容進(jìn)行統(tǒng)計(jì)[1]。地表覆蓋作為面狀要素是地理國情監(jiān)測(cè)的基礎(chǔ)，其表面面積是幾何特征類型中面狀要素在某時(shí)刻的一項(xiàng)基本統(tǒng)計(jì)指標(biāo)，不同時(shí)間狀態(tài)下的表面面積能夠反映地表覆蓋面積的動(dòng)態(tài)變化。目前有很多方法來統(tǒng)計(jì)表面積，李軼平等(2007)運(yùn)用投影面積和坡度的正割值計(jì)算表面積[2]，王長海(2010)用逐步分割多變形的方法獲得多邊形的表面積[3]，江帆(2005)將1個(gè)像素格分成2個(gè)三角形來統(tǒng)計(jì)DEM的表面積[4]，Jeff S.Jenness(2004)用8個(gè)周圍毗鄰點(diǎn)將1個(gè)像素格分成8個(gè)空間三角形計(jì)算DEM的表面積[5]。相比之下，Jeff S.Jenness的工作在處理柵格數(shù)據(jù)上更具有借鑒意義。這些算法只是提供了基本的計(jì)算原理，并沒有考慮到地球曲率的影響，且回避了怎樣計(jì)算目標(biāo)區(qū)域邊界附近面積的問題，因此，其并不太適合計(jì)算地表面積。本文的算法在前人的基礎(chǔ)上做出了適合計(jì)算地理表面積的改進(jìn)。

2 方 法

DEM影像數(shù)據(jù)是將一系列高程點(diǎn)按規(guī)則組織起來的矩陣格網(wǎng)，可用這些高程點(diǎn)來模擬地表的形態(tài)特征，計(jì)算表面積基本的原理是將DEM組織成三角網(wǎng)，網(wǎng)中的三角形可以近似的代表該三角形區(qū)域的面積，所有相鄰的三角形組成的三角網(wǎng)的面積可以表示這個(gè)DEM的表面積。本文DEM的坐標(biāo)系統(tǒng)是地理坐標(biāo)系統(tǒng)。

2.1 原 理

以4個(gè)相毗鄰的像素格為計(jì)算單元，將4個(gè)像素的中心點(diǎn)與這4個(gè)像素連接起來，組成4個(gè)三角形，通過三角形(Ⅰ、Ⅱ、Ⅲ、Ⅳ)的面積計(jì)算得到單元的面積。

圖1 4個(gè)相毗鄰的像素格組成一個(gè)計(jì)算單元示意圖

三角形(Ⅰ、Ⅱ、Ⅲ和Ⅳ)的面積可用海倫-秦九韶公式計(jì)算。

(1)

其中，a、b、c代表三角形的邊長，S是三角形的面積。三角形的邊長就是三角形的頂點(diǎn)之間的歐式距離，DEM高程點(diǎn)在地理坐標(biāo)下，需要轉(zhuǎn)為空間直角坐標(biāo)才能計(jì)算各邊長的長度。

圖2 過地球旋轉(zhuǎn)軸一個(gè)截面的各項(xiàng)參數(shù)示意圖

圖2中O為橢球中心點(diǎn)，a是橢球的長半軸，b為橢球的短半軸，μ為橢圓離心角，θ為橢圓極角，β為地面點(diǎn)Q的緯度，P為過Q點(diǎn)的法線與橢球的交點(diǎn)。

橢圓的參數(shù)方程：

(2)

橢圓中各角度參數(shù)之間的關(guān)系[6]：

(3)

計(jì)算地心O到橢球面一點(diǎn)P的距離：

(4)

計(jì)算地心O到地表點(diǎn)Q的長，令R=OQ：

(5)

其中，PQ=H=h+ζ，H為大地高，h為正常高，ζ為高程異常。此時(shí)可以認(rèn)為H為DEM的值。

?為地表點(diǎn)Q和地心連線與赤道面的夾角，δ=∠POQ，

(6)

計(jì)算地表點(diǎn)Q的空間直角坐標(biāo)(x，y，z)：

(7)

其中，L為經(jīng)度。將式(2)到式(6)帶入式(7)便可將DEM數(shù)據(jù)某點(diǎn)的地理坐標(biāo)轉(zhuǎn)化為空間直角坐標(biāo)。

空間直角坐標(biāo)系中任意兩點(diǎn)A(xA，yA，zA)和點(diǎn)B(xB，yB，zB)的歐式距離ΔDAB：

(8)

圖1中4個(gè)柵格值是已知的，而中心點(diǎn)的高程值是不存在的，此時(shí)將其轉(zhuǎn)換成空間直角坐標(biāo)系會(huì)發(fā)現(xiàn)式(5)、式(6)中的H是未知的，其值需要通過內(nèi)插得到，可以用反距離加權(quán)法[7-8]對(duì)中心點(diǎn)的高程值進(jìn)行估計(jì)：

(9)

2.2 流 程

本文的目的是計(jì)算多邊形圈出范圍的表面積，因此在計(jì)算時(shí)要考慮多邊形與柵格點(diǎn)的關(guān)系。計(jì)算流程如圖3所示：

圖3 表面積計(jì)算的基本流程圖

(1)根據(jù)多邊形的四至范圍獲取恰好覆蓋多邊形的DEM數(shù)據(jù)，DEM的四至范圍作為第一個(gè)矩形。

(2)對(duì)現(xiàn)有的矩形進(jìn)行分割，矩形的邊長長度單位是DEM柵格的長度，在分割時(shí)不能分成同樣大小的矩形，只能分成4個(gè)大小盡可能相同的矩形。

(3)判斷所有的矩形與多邊形的位置關(guān)系。若矩形完全在多邊形內(nèi)，以上文中提到的以單元為單位計(jì)算矩形內(nèi)DEM數(shù)據(jù)，累加各個(gè)計(jì)算單元的面積，計(jì)算完成后刪除矩形；若矩形完全在多邊形外，直接刪除矩形；若矩形與多邊形邊界有接觸且矩形邊長不為1個(gè)柵格長度，返回第2步，直至所有的矩形成為邊長都為1個(gè)像素長度的矩形。

(4)當(dāng)剩下的矩形的4條邊長為一個(gè)單位長度時(shí)，提取矩形與多邊形重合的部分，這部分重新組成小于像素方格的小多邊形。對(duì)各個(gè)小多邊形沒有高程的邊界點(diǎn)用式(9)賦予高程值，用多邊形上的點(diǎn)組成帶有約束條件的Delaunay三角網(wǎng)[10]，計(jì)算并累加每個(gè)三角形的面積。得到各個(gè)小多邊形的面積，累加到多邊形內(nèi)部面積。

計(jì)算步驟中的矩形與多邊形的位置關(guān)系可簡單分為三大類：多邊形內(nèi)，多邊形外，多邊形邊界上，根據(jù)矩形頂點(diǎn)在多邊形內(nèi)的個(gè)數(shù)這三類可延伸為6種類型，如圖4所示(實(shí)心方塊代表在多邊形內(nèi)的網(wǎng)格點(diǎn)，灰色部分為多邊形)：

圖4 矩形與多邊形關(guān)系示意圖

(1)1個(gè)網(wǎng)格點(diǎn)在多邊形內(nèi)：矩形只有1個(gè)點(diǎn)在多變形內(nèi)(圖4(a))；

(2)2個(gè)網(wǎng)格點(diǎn)在多邊形內(nèi)：矩形只有2個(gè)相鄰點(diǎn)在多邊形內(nèi)(圖4(b))；

(3)2個(gè)網(wǎng)格點(diǎn)在多邊形內(nèi)：矩形只有2個(gè)對(duì)角點(diǎn)在多變形內(nèi)(圖4(c))；

(4)3個(gè)網(wǎng)格點(diǎn)在多邊形內(nèi)：矩形只有3個(gè)點(diǎn)在多邊形內(nèi)(圖4(d))；

(5)4個(gè)網(wǎng)格點(diǎn)在多邊形內(nèi)：矩形4個(gè)點(diǎn)在多邊形內(nèi)，矩形的某條邊有可能被多邊形邊界穿過偶數(shù)次，矩形內(nèi)部可能有多邊形邊界點(diǎn)(圖4(e))；

(6)0個(gè)網(wǎng)格點(diǎn)在多邊形內(nèi)：矩形4個(gè)點(diǎn)在多邊形外，矩形的某條邊有可能被多邊形邊界穿過偶數(shù)次，矩形內(nèi)部可能有多邊形的點(diǎn)(圖4(f))。

(1)、(2)、(3)和(4)4種情況下矩形與多邊形有交集，表示其在多邊形的邊界上，(5)和(6)情況下矩形有在多邊形邊界上的可能性，但(5)在一般情況下表示矩形在多變形內(nèi)，(6)在一般情況下表示矩形在多變形外，但在某些情況下(5)、(6)矩形與多邊形是有交集的，此時(shí)其在多邊形邊界上。

矩形與多邊形位置關(guān)系的判斷最基本的是判斷矩形頂點(diǎn)與多變形的關(guān)系，其基本原理是以影像點(diǎn)為端點(diǎn)向右做出射線，若射線與多邊形的交點(diǎn)個(gè)數(shù)為奇數(shù)，則該點(diǎn)在多邊形內(nèi)，否則該點(diǎn)在多邊形外。射線法也存在缺點(diǎn)，在使用時(shí)需要考慮以下幾種特殊情況：

(1)射線恰好穿過多邊形最高點(diǎn)或最低點(diǎn)；

(2)射線穿過多邊形某條邊后，又恰好穿過某個(gè)極值點(diǎn)；

(3)射線端點(diǎn)恰好在多邊形上；

(4)射線恰好與某一條邊重合。

其中，灰色表示點(diǎn)在多邊形內(nèi)，黑色表示點(diǎn)在多邊形外。圖5右圓形窗口表示某邊界處的局部放大圖，白線為多邊形的邊線，黑點(diǎn)為在多變形內(nèi)部的高程點(diǎn)。

圖5 某縣縣界為對(duì)象驗(yàn)證射線法結(jié)果圖

結(jié)合各矩形頂點(diǎn)與多邊形的關(guān)系以及矩形內(nèi)是否有多邊形的點(diǎn)，便可將矩形分成圖4中的6種類型。

3 實(shí)驗(yàn)結(jié)果

選取中國四川省及附近的20縣市為實(shí)驗(yàn)區(qū)域，分別用本文算法(算法A)和ArcGIS(算法B)計(jì)算這些地方的表面面積，結(jié)果如表1所示：

表1 算法A與算法B表面積計(jì)算結(jié)果比較

表1表明算法A所計(jì)算的試驗(yàn)區(qū)域表面面積的數(shù)值平均值比算法B所得數(shù)值平均值大約0.61%。

表2表明同一區(qū)域DEM分辨率的增加會(huì)使該區(qū)域的表面積也隨之增加。

表3 當(dāng)圖斑尺寸接近DEM分辨率時(shí)兩種算法計(jì)算結(jié)果比較

表3和圖6表示ArcGIS在處理圖斑邊界時(shí)不精確，當(dāng)圖斑的尺寸達(dá)到DEM分辨率的級(jí)別時(shí)，ArcGIS并不能統(tǒng)計(jì)出圖斑面積。

ArcGIS計(jì)算各縣面積需要平面投影后的DEM數(shù)據(jù)，然后用shp文件中的多邊形進(jìn)行Mask處理，才能統(tǒng)計(jì)得到縣行政區(qū)劃的面積。GIS軟件通常只有像素方格中心位于多邊形內(nèi)才會(huì)將此像素格進(jìn)行表面積統(tǒng)計(jì)，所以Mask處理后的影像邊界是鋸齒狀的，在處理邊界問題顯得粗糙。用這種方法統(tǒng)計(jì)出的地表面積精度和DEM的空間分辨率是息息相關(guān)的，DEM分辨率越低，面積統(tǒng)計(jì)精度也越低。這種計(jì)算方法在處理比較小的圖斑的時(shí)候會(huì)對(duì)結(jié)果產(chǎn)生較大的影響。圖6和表3表示ArcGIS在處理某些圖斑邊界時(shí)不夠精確，當(dāng)圖斑的尺寸達(dá)到DEM分辨率的級(jí)別時(shí)，ArcGIS并不能統(tǒng)計(jì)出圖斑面積。

圖6 ArcGIS處理圖斑效果圖

圖6中方塊組成的圖形是GIS軟件根據(jù)多邊形截取的DEM，其覆蓋的范圍與多邊形有差別。當(dāng)多邊形非常小時(shí)，會(huì)被GIS軟件忽略掉(圖6(b)，虛線為1個(gè)DEM影像像素的范圍，其內(nèi)部有1個(gè)多邊形)。

算法A是直接基于地理坐標(biāo)數(shù)據(jù)，省去了平面投影過程，避免了投影變形，且算法在邊界處完全按照實(shí)際邊界進(jìn)行處理，即考慮多邊形邊界切割DEM影像像素的實(shí)際情況進(jìn)行計(jì)算，可以避免投影誤差、邊界處面積計(jì)算誤差，并且本算法在圖斑面積非常小的情況下也能夠計(jì)算該圖斑的面積。

4 結(jié)束語

本文提出的基于地理坐標(biāo)的表面積算法直接采用大地坐標(biāo)進(jìn)行計(jì)算，避免了傳統(tǒng)軟件計(jì)算時(shí)因平面投影產(chǎn)生的誤差傳遞到表面積統(tǒng)計(jì)結(jié)果。在處理計(jì)算區(qū)域邊界時(shí)，按照邊界的實(shí)際位置進(jìn)行處理，使計(jì)算結(jié)果更加貼近真實(shí)值。然而算法也存在不足之處，如目前沒有找到更有效率的方法處理多邊形內(nèi)有蟲洞的情況，這是以后需要解決的問題。

參考文獻(xiàn)：

[1] 國務(wù)院第一次全國地理國情普查領(lǐng)導(dǎo)小組辦公室.地理國情普查基本統(tǒng)計(jì)[M].北京：測(cè)繪出版社，2013.

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

[3] 王長海.DEM模型的三維地表面積算法研究[J].紅水河，2010，29(3)：153-154.

[4] 苗春葆.點(diǎn)與多邊形關(guān)系的射線法[J].電腦編程技巧與維護(hù)，2008(3)：56-58.

[5] 張錦明，郭麗萍，張小丹.反距離加權(quán)插值算法中插值參數(shù)對(duì)DEM插值誤差的影響[J].測(cè)繪科學(xué)技術(shù)學(xué)報(bào)，2012，29(1)：51-56.

[6] JENNESS J S.Calculating landscape surface area from digital elevation models[J].Wildlife Society Bulletin，2004，32(3)：829-839.

[7] AGUILAR F J，AGüERA F，AGUILAR M A，et al.Effects of terrain morphology，sampling density，and interpolation methods on grid DEM accuracy[J].Photogrmmetric Engineering and Remote Sensing，2005：805-816.

[8] WEINTRIT A.So，What is actually the distance from the equator to the pole?—Overview of the meridian distance approximations[J].The International Journal on Marine Navigation and Safety of Sea Transportation，2013，2(7)：259-272.

[9] TSENG W K，LEE H S.Navigation on a great ellipse[J].Journal of Marine Science and Technology，2010，18(3)：369-375.

[10] SHEWCHUK J R.Delaunay refinement algorithms for triangular mesh generation[J].Computational Geometry，2002，22：21-74.