基于地理坐標的地理國情統計單元表面面積精確計算

謝成磊，趙榮，梁勇

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

1 引 言

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

2 方 法

DEM影像數據是將一系列高程點按規則組織起來的矩陣格網，可用這些高程點來模擬地表的形態特征，計算表面積基本的原理是將DEM組織成三角網，網中的三角形可以近似的代表該三角形區域的面積，所有相鄰的三角形組成的三角網的面積可以表示這個DEM的表面積。本文DEM的坐標系統是地理坐標系統。

2.1 原 理

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

圖1 4個相毗鄰的像素格組成一個計算單元示意圖

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

(1)

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

圖2 過地球旋轉軸一個截面的各項參數示意圖

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

橢圓的參數方程：

(2)

橢圓中各角度參數之間的關系[6]：

(3)

計算地心O到橢球面一點P的距離：

(4)

計算地心O到地表點Q的長，令R=OQ：

(5)

其中，PQ=H=h+ζ，H為大地高，h為正常高，ζ為高程異常。此時可以認為H為DEM的值。

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

(6)

計算地表點Q的空間直角坐標(x，y，z)：

(7)

其中，L為經度。將式(2)到式(6)帶入式(7)便可將DEM數據某點的地理坐標轉化為空間直角坐標。

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

(8)

圖1中4個柵格值是已知的，而中心點的高程值是不存在的，此時將其轉換成空間直角坐標系會發現式(5)、式(6)中的H是未知的，其值需要通過內插得到，可以用反距離加權法[7-8]對中心點的高程值進行估計：

(9)

2.2 流 程

本文的目的是計算多邊形圈出范圍的表面積，因此在計算時要考慮多邊形與柵格點的關系。計算流程如圖3所示：

圖3 表面積計算的基本流程圖

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

(2)對現有的矩形進行分割，矩形的邊長長度單位是DEM柵格的長度，在分割時不能分成同樣大小的矩形，只能分成4個大小盡可能相同的矩形。

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

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

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

圖4 矩形與多邊形關系示意圖

(1)1個網格點在多邊形內：矩形只有1個點在多變形內(圖4(a))；

(2)2個網格點在多邊形內：矩形只有2個相鄰點在多邊形內(圖4(b))；

(3)2個網格點在多邊形內：矩形只有2個對角點在多變形內(圖4(c))；

(4)3個網格點在多邊形內：矩形只有3個點在多邊形內(圖4(d))；

(5)4個網格點在多邊形內：矩形4個點在多邊形內，矩形的某條邊有可能被多邊形邊界穿過偶數次，矩形內部可能有多邊形邊界點(圖4(e))；

(6)0個網格點在多邊形內：矩形4個點在多邊形外，矩形的某條邊有可能被多邊形邊界穿過偶數次，矩形內部可能有多邊形的點(圖4(f))。

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

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

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

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

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

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

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

圖5 某縣縣界為對象驗證射線法結果圖

結合各矩形頂點與多邊形的關系以及矩形內是否有多邊形的點，便可將矩形分成圖4中的6種類型。

3 實驗結果

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

表1 算法A與算法B表面積計算結果比較

表1表明算法A所計算的試驗區域表面面積的數值平均值比算法B所得數值平均值大約0.61%。

表2表明同一區域DEM分辨率的增加會使該區域的表面積也隨之增加。

表3 當圖斑尺寸接近DEM分辨率時兩種算法計算結果比較

表3和圖6表示ArcGIS在處理圖斑邊界時不精確，當圖斑的尺寸達到DEM分辨率的級別時，ArcGIS并不能統計出圖斑面積。

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

圖6 ArcGIS處理圖斑效果圖

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

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

4 結束語

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

參考文獻：

[1] 國務院第一次全國地理國情普查領導小組辦公室.地理國情普查基本統計[M].北京：測繪出版社，2013.

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

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

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

[5] 張錦明，郭麗萍，張小丹.反距離加權插值算法中插值參數對DEM插值誤差的影響[J].測繪科學技術學報，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.