廖洋洋, 陳建華, 張 倩, 劉 帥, 粟 解

(1.成都理工大學(xué) 地球物理學(xué)院，成都 610059；2.阿壩州自然資源與科技信息研究所，汶川 623000)

0 引言

傳統(tǒng)的地形測量由于其自身的局限且存在高危險(xiǎn)性的隱患，在一些人不能到達(dá)的區(qū)域無法進(jìn)行測量。隨著地理信息技術(shù)與互聯(lián)網(wǎng)技術(shù)的快速發(fā)展，三維實(shí)景地形模型已經(jīng)能夠在瀏覽器中呈現(xiàn)。在瀏覽器中最受歡迎的三維地理信息庫是Cesium，它通過JavaScript封裝的三維地理信息可視化的庫，內(nèi)部有大量的接口方便使用。三維WebGIS受到了越來越多的關(guān)注，其中對(duì)地形測量是其中的一個(gè)熱點(diǎn)。Michael等[1]使用新的Web技術(shù)(如HTML5、WebGL)，實(shí)現(xiàn)高效的客戶端分析，旨在推動(dòng)基于Web的3D地理空間分析發(fā)展；Leila等[2]考慮到計(jì)算地形可見性問題，概述了在規(guī)則網(wǎng)格和不規(guī)則網(wǎng)格解決地形表面的算法；王恒[3]基于Cesium對(duì)地方大地測量平臺(tái)數(shù)據(jù)庫進(jìn)行設(shè)計(jì)，提供了三維大地測量數(shù)據(jù)的展示、定位和查詢；范俊甫[4]對(duì)Cesium框架多源電子地圖瓦片方案進(jìn)行設(shè)計(jì)，為地理信息系統(tǒng)與遙感技術(shù)項(xiàng)目應(yīng)用提供了一種高靈活、低成本的解決方案；朱栩逸等[5]在GIS框架和Web服務(wù)的基礎(chǔ)上，提出了一種基于Cesium的三維WebGIS搭建方案，并簡單地實(shí)現(xiàn)了繪制功能；李偉等[6]使用ArcGIS軟件的三維分析功能，通過建立數(shù)字高程模型，使用不規(guī)則三角網(wǎng)(TIN)，對(duì)鄱陽湖洪水淹沒面積進(jìn)行測量。

上述研究取得了較好的效果，但是目前的方法還無法滿足Web環(huán)境下三維地形面積測量的現(xiàn)實(shí)需求。這里以三維開源地理信息庫(Cesium)和地形切片技術(shù)為基礎(chǔ)，重點(diǎn)針對(duì)真實(shí)三維地形模型的表面積測量進(jìn)行研究，為地形測量提供一種有效、可靠的測量手段。

1 系統(tǒng)架構(gòu)

系統(tǒng)架構(gòu)圖見圖1。采用瀏覽器到服務(wù)器再到數(shù)據(jù)層的方式獲取基本數(shù)據(jù)。根據(jù)Cesium官方文檔說明，加載的地形格式數(shù)據(jù)可以為.terrain格式，使用CesiumLab開源工具，可以將數(shù)字高程數(shù)據(jù)格式轉(zhuǎn)換為Cesium能識(shí)別的.terrain數(shù)據(jù)格式。使用IIS發(fā)布地形數(shù)據(jù)，Cesium內(nèi)置有讀取加載.terrain的方法，使用CesiumTerrainProvider方法可以加載已經(jīng)發(fā)布的地形數(shù)據(jù)。將遙感影像數(shù)據(jù)使用Geoserver數(shù)據(jù)進(jìn)行發(fā)布，使用Cesium內(nèi)置imageryProvider方法加載發(fā)布好的影像數(shù)據(jù)，以便在有起伏的影像數(shù)據(jù)上繪制需要測量的區(qū)域。Cesium內(nèi)置有許多方便的接口供我們使用，其中提供的接口方法sampleTerrainMostDetailed，可以快速地獲取任何分辨率下地形切片任意一點(diǎn)最大精細(xì)程度的高度，而不必重新寫一種前后端交互的新方法，也不需要后端操作DEM方法，這一種專門讀取.terrain的方法，考慮到實(shí)現(xiàn)的復(fù)雜程度，所以采用前端計(jì)算的方式來實(shí)現(xiàn)整個(gè)算法，該方法主要針對(duì)Cesium提供一種計(jì)算思路。

圖1 系統(tǒng)架構(gòu)

2 具體算法

2.1 算法流程

算法的整體流程如圖2所示，通過繪制需要測量的區(qū)域，可以得到一個(gè)任意形狀多邊形，將多邊形分解成為多個(gè)三角形，使用雙線性插值的方法分別在每個(gè)三角形的最小邊界矩形中插值，每個(gè)三角形各自判斷插值過后的哪些點(diǎn)在自己內(nèi)部，排除掉不在自己內(nèi)部點(diǎn)，然后獲取各自內(nèi)部點(diǎn)的高度，將這些在內(nèi)部的點(diǎn)重新連接起來構(gòu)成新的三角網(wǎng)，從而計(jì)算每個(gè)小三角形的面積，以此類推，直到分割過后所有的三角形都計(jì)算完面積，最后累加求和，得到整個(gè)區(qū)域的面積。

圖2 算法流程

2.2 區(qū)域切割

利用Cesium的ScreenSpaceEventHandler方法，用來捕捉我們繪制區(qū)域的事件?？紤]到繪制的區(qū)域可能是任意形狀的多邊形，這里使用耳切法(Ear Clipping)對(duì)區(qū)域進(jìn)行切割。

耳切法：連續(xù)頂點(diǎn)V1、V2、V3組成的內(nèi)部不包含其他任意頂點(diǎn)的三角形，將其稱為簡單多邊形的耳朵(內(nèi)角小于180°)[7]。V1與V3之間連接的線段稱之為多邊形對(duì)角線，點(diǎn)V2稱為耳朵頂點(diǎn)。雖然可以認(rèn)為三角形任何一個(gè)頂點(diǎn)都能是耳朵，但是原理上認(rèn)為一個(gè)三角形只有一個(gè)耳朵。一個(gè)由4個(gè)以上頂點(diǎn)組成的多邊形至少有2個(gè)不重疊的耳朵。根據(jù)這個(gè)理論特性，可以利用遞歸思想來解決三角形切割的方法。針對(duì)包含由n個(gè)頂點(diǎn)組成的簡單多邊形，找到其第一個(gè)耳朵，移除該耳朵頂點(diǎn)，此時(shí)剩余頂點(diǎn)將組成一個(gè)n-1個(gè)點(diǎn)的多邊形。重復(fù)這個(gè)操作，直到剩余三個(gè)頂點(diǎn)，遞歸結(jié)束。

2.3 區(qū)域插值

若繪制區(qū)域不是三角形則區(qū)域切割后會(huì)形成至少兩個(gè)以上的三角形，對(duì)每個(gè)三角形求取最小邊界矩形。最小邊界矩形是以三角形三個(gè)頂點(diǎn)(xk，yk)為基準(zhǔn)，尋找出xmin、ymin、xmax、ymax4個(gè)值，然后以(xmin，ymin)、(xmin，ymax)、(xmax,ymax)、(xmax,ymin)4個(gè)點(diǎn)構(gòu)成最小邊界矩形。

當(dāng)獲取到每個(gè)三角形最小邊界矩形的時(shí)候，使用線性插值方法，對(duì)每個(gè)最小邊界矩形進(jìn)行等間距插值。線性方程即線性插值的核心，使用線性插值的原因是因?yàn)槠洳僮骱唵危⑶夷芫鶆虻貙?duì)地區(qū)進(jìn)行插值覆蓋。其已知兩個(gè)點(diǎn)的坐標(biāo)(x0，y0)與(x1，y1)，假設(shè)存在未知函數(shù)f，則函數(shù)上的一點(diǎn)可以表示為式(1)。

f(x,y)=f(x0)=w(f(x1)-f(x0))

(1)

其中w的比率為式(2)。

(2)

按照自定義精度在矩形范圍內(nèi)進(jìn)行線性插值，精度越高，計(jì)算時(shí)間越長，然后將每個(gè)點(diǎn)與三角形進(jìn)行判斷，判斷點(diǎn)是否在三角形內(nèi)部。

獲取到所有在三角形內(nèi)部的點(diǎn)后，使用Cesium提供Sample Terrain Most Detailed方法就能獲取每個(gè)內(nèi)部的點(diǎn)的最大細(xì)化程度的高程值，。然后使用哈希表來存儲(chǔ)每個(gè)點(diǎn)高程信息。

2.4 點(diǎn)集的三角剖分

將任意多邊形區(qū)域切割成互相獨(dú)立的三角形，然后再對(duì)每個(gè)三角形的內(nèi)部進(jìn)行插值，由此每一個(gè)三角形對(duì)應(yīng)了插值后一個(gè)點(diǎn)的集合。每一個(gè)集合要做一個(gè)三角格網(wǎng)化，筆者使用Delaunay三角剖分算法對(duì)點(diǎn)集合進(jìn)行三角格網(wǎng)創(chuàng)建。

實(shí)現(xiàn)Delaunay三角格網(wǎng)有Lawson算法和Bowyer-Watson算法等。Lawson算法是由Lawson在1977年提出，該算法具有思路簡單明了、編程簡單并且容易實(shí)現(xiàn)的優(yōu)點(diǎn)，同時(shí)被稱為逐點(diǎn)插入法[8]。Lawson的基本原理為首先創(chuàng)建一個(gè)包含散點(diǎn)集合中所有散點(diǎn)的簡單多邊形，將散點(diǎn)集合中未插入的散點(diǎn)逐一插入到三角格網(wǎng)中，并對(duì)新插入的點(diǎn)做空外接圓檢測和局部調(diào)整，要求重新生成的三角格網(wǎng)均符合Delaunay準(zhǔn)則，最后移除包含初始化的多邊形頂點(diǎn)的三角形。

筆者使用改進(jìn)的Lawson算法Bowyer-Watson算法實(shí)現(xiàn)地形三角格網(wǎng)剖分。Bowyer-Watson算法描述如下：

1)創(chuàng)建一個(gè)包含所有散點(diǎn)多邊形或者三角形，放入鏈表中。

2)按順序?qū)⒓现械纳Ⅻc(diǎn)插入到第一步創(chuàng)建的多邊形或者三角形中，在鏈表中找出外接圓包含插入散點(diǎn)的三角形，刪除外接圓的公共邊，形成Delaunay空腔。

3)插入的散點(diǎn)連接Delaunay空腔上的頂點(diǎn)，形成新的三角格網(wǎng)。

4)重復(fù)步驟2)、步驟3)，直到所有散點(diǎn)插入完畢。

2.5 區(qū)域面積求和

根據(jù)前面步驟能獲取到每個(gè)插值點(diǎn)的高度的信息，可以利用海倫公式(3)對(duì)三角格網(wǎng)中的每個(gè)三角形進(jìn)行面積的計(jì)算。

(3)

其中：a、b、c分別代表三角形三個(gè)邊的長度；p代表的是三角形的周長。每一個(gè)邊的邊長可以使用歐式距離公式來計(jì)算，n維歐式距離表示為式(4)。

(4)

其中：n表示維度；x、y為在同一維度上的值。最后將三角格網(wǎng)中三角形的面積逐一相加，得到三維地形表面積測量的結(jié)果。

3 實(shí)例分析

筆者以Cesium.js作為3D WebGIS前端框架，IIS服務(wù)器為地形服務(wù)提供數(shù)據(jù)響應(yīng)，用CesiumLab處理30 m分辨率的DEM數(shù)據(jù)產(chǎn)生.terrain格式的地形切片數(shù)據(jù)，以GeoServer為影像地圖服務(wù)器，本文所有的計(jì)算均是由瀏覽器內(nèi)核計(jì)算，無后臺(tái)服務(wù)器提供計(jì)算。實(shí)驗(yàn)環(huán)境是采用的win10操作系統(tǒng)，Cpu為i5-3470，瀏覽器是谷歌瀏覽器。

Cesium是由JavaScript腳本語言編寫的前端輕量級(jí)開源WebGIS框架，也是開源三維地圖引擎。Cesium不需要其他任何瀏覽器插件支持，只需要瀏覽器支持WebGL功能即可，并且具有跨平臺(tái)、跨瀏覽器的優(yōu)點(diǎn)[9]。GeoServer是開源地圖服務(wù)器，能方便用戶快速的發(fā)布地圖影像數(shù)據(jù)。IIS是由微軟公司開發(fā)的Web服務(wù)程序。

在筆者研發(fā)的三維平臺(tái)上，分別選取地勢較為平坦、地勢較為陡峭、大面積區(qū)域、小面積區(qū)域、凸多邊形陡峭區(qū)域、凹多邊形陡峭區(qū)域、凸多邊形平坦區(qū)域、凹多邊形平坦區(qū)域共計(jì)8組相同控制點(diǎn)區(qū)域，利用傳統(tǒng)的投影方式與本文提出的三維地形面積測量方法來進(jìn)行實(shí)驗(yàn)，以15 m為插值點(diǎn)之間的間距來計(jì)算貼地面積。最終實(shí)驗(yàn)對(duì)比結(jié)果見表1。

表1 面積測量結(jié)果對(duì)比分析

實(shí)驗(yàn)在地勢陡峭區(qū)域的投影面積效果圖如圖3所示，貼地表面積如圖4所示，投影面積為4 805 292.3 m2，貼地表面積為5 696 831.239 m2。

圖3 投影面積

圖4 貼地面積

由表1可知，在選擇同樣的控制點(diǎn)的情況下，投影面積始終小于貼地表面積，符合實(shí)驗(yàn)的預(yù)期結(jié)果。在地勢較為平坦的情況下，投影面積和貼地表面積相差不是很大，在地勢較為陡峭的區(qū)域，投影面積和貼地表面積相差較大，這兩者從側(cè)面驗(yàn)證了本文提出方法的可靠性。

4 結(jié)論

筆者提出了一種基于Web的三維地形表面積測量方法，其具有無需后臺(tái)服務(wù)器計(jì)算，依靠瀏覽器內(nèi)核計(jì)算，同時(shí)能夠快速嵌入前端的特點(diǎn)。該方法在Cesium庫的基礎(chǔ)上結(jié)合耳切法、Delaunay三角網(wǎng)算法與線性插值構(gòu)建區(qū)域三角網(wǎng)，然后利用海倫公式來求取三維地形表面積。其中耳切法、Delaunay三角網(wǎng)算法具有比較成熟且應(yīng)用較廣的特點(diǎn)，能有效解決三維三角形格網(wǎng)面積的計(jì)算。該方法可以輔助人們?cè)诘貏荻盖碗y以測量的區(qū)域進(jìn)行測量，例如滑坡面積測量以及露頭巖層面積的測量，一定程度上減少了人們到高危險(xiǎn)區(qū)域的危險(xiǎn)性以及費(fèi)時(shí)、費(fèi)力的局限性，同時(shí)也具有高度靈活測量的優(yōu)點(diǎn)。