謝成磊,趙榮,梁勇
(1.山東農(nóng)業(yè)大學(xué) 信息科學(xué)與工程學(xué)院,山東 泰安 271018;2.中國測(cè)繪科學(xué)研究院,北京 100830)
基本統(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)。
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)。
以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)
本文的目的是計(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種類型。
選取中國四川省及附近的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ì)算該圖斑的面積。
本文提出的基于地理坐標(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.