潘雅輝,相方莉
(浙江長征職業(yè)技術(shù)學(xué)院,浙江杭州310023)
建立高質(zhì)量的數(shù)字高程模型(Digital Elevation Model,DEM)在土壤侵蝕研究、土地利用規(guī)劃與評價(jià)、景觀設(shè)計(jì)等領(lǐng)域具有重要意義[1-2].目前,利用等高線生成DEM的方法較多,但在精度、效率等方面仍存在一些問題.如當(dāng)?shù)雀呔€數(shù)據(jù)不足時(shí),等高線直接內(nèi)插法在確定內(nèi)插所需要的參考點(diǎn)時(shí)會(huì)出現(xiàn)問題;變換函數(shù)法精度較低;趨勢面分析法需要的樣點(diǎn)數(shù)較多等.此外,上述DEM生成算法在計(jì)算機(jī)實(shí)現(xiàn)及集成方面也存在一定難度[3-6].因此,本文設(shè)計(jì)開發(fā)了一種數(shù)字高程模型生成算法:動(dòng)態(tài)網(wǎng)格模板法(DGT);并采用OpenGL圖像庫和組件技術(shù),開發(fā)了一個(gè)三維可視化插件,實(shí)現(xiàn)了原始等高線數(shù)據(jù)、遙感數(shù)據(jù)及DEM數(shù)據(jù)的多源數(shù)據(jù)集成和一體化顯示.
DEM內(nèi)插算法的實(shí)現(xiàn)是以待定點(diǎn)(未知高程值點(diǎn))的周邊相鄰已知高程值的點(diǎn)為參考,求解待定點(diǎn)高程值的過程.任何內(nèi)插算法都有一定的假設(shè)前提:原始地形的起伏變化是連續(xù)和平滑的;相鄰點(diǎn)之間,無論是未知值的點(diǎn)還是已知值的點(diǎn),它們之間都有很大的關(guān)聯(lián)性或相關(guān)性,只有這樣才有可能通過相鄰的已知值的點(diǎn)來求解出或內(nèi)插出待定點(diǎn)的高程值.本文所采用的DGT算法屬于逐點(diǎn)內(nèi)插.
在移動(dòng)擬合法中,一般可以選取與待定點(diǎn)鄰近的、高程值已知的數(shù)據(jù)點(diǎn)作為參考點(diǎn),來擬合一個(gè)多項(xiàng)式曲面,曲面擬合公式如下:
式中 X、Y、Z 是各參考點(diǎn)的坐標(biāo)值,A、B、C、D、E、F 為待求解的參數(shù)[7]134.移動(dòng)擬合法主要涉及以下兩方面的問題:(1)如果參考點(diǎn)太多,會(huì)影響計(jì)算效率;而太少,則會(huì)影響計(jì)算精度甚至根本無法求解.為了保證有充足的參考點(diǎn),就必須確定合理的最小鄰域.(2)參考點(diǎn)距離待定點(diǎn)有近有遠(yuǎn),與待定點(diǎn)的相關(guān)性會(huì)有所不同,因此確定各參考點(diǎn)的權(quán)重也很重要.
最小鄰域的確定有搜索范圍和參考點(diǎn)數(shù)量兩個(gè)影響因素.圖1為基于范圍選點(diǎn)的例子,其選中的點(diǎn)(參考點(diǎn))都位于以待定點(diǎn)為圓心、R為半徑的圓形范圍之內(nèi).R值的大小取決于原始數(shù)據(jù)點(diǎn)密度情況和可能影響的范圍.根據(jù)二次曲面方程求解的需要,同時(shí)保證插值的效率,可采用動(dòng)態(tài)搜索圓解決這個(gè)問題,即根據(jù)數(shù)據(jù)點(diǎn)平均密度確定圓內(nèi)數(shù)據(jù)點(diǎn)(平均要有10個(gè)),從而求得搜索圓的半徑,其公式如下:
式中N為總點(diǎn)數(shù),A為總面積.
圖1 基于范圍選參考點(diǎn)
移動(dòng)擬合法需要求解復(fù)雜的誤差方程,為了運(yùn)算簡便快捷可將加權(quán)平均法看作移動(dòng)擬合法的特例,其參考點(diǎn)范圍的確定與移動(dòng)擬合法相同[7]136.求解待定點(diǎn)x的高程值的公式如下:
式中Zi是第i個(gè)參考點(diǎn)的高程值,n為參考點(diǎn)數(shù)量,xi是第i個(gè)參考點(diǎn)的權(quán)重,Zx是待定點(diǎn)x的高程值.對于權(quán)重的確定,一般采取與距離相關(guān)的權(quán)函數(shù),常用的權(quán)函數(shù)有:
式中R是搜索圓的半徑,r是待定點(diǎn)到參考點(diǎn)的距離,p是參考點(diǎn)的權(quán)重.
在本研究中,等高線生成DEM模型的基本步驟為:(1)將地形圖中的等高線要素層數(shù)字化生成含高程屬性值的ESRI Shape file文件格式.(2)確定X、Y方向的網(wǎng)格數(shù),并生成網(wǎng)格矩陣.在本算法中,網(wǎng)格不一定是規(guī)則正方形,X值和Y值可以不同.確定X、Y方向網(wǎng)格數(shù)量的公式為:X方向網(wǎng)格數(shù)量 =等高線要素層X方向空間范圍/網(wǎng)格X方向?qū)挾?Y方向網(wǎng)格數(shù)量 =等高線要素層Y方向空間范圍/網(wǎng)格Y方向?qū)挾?(3)初始化參考網(wǎng)格值,即給參考點(diǎn)賦初始值.(4)根據(jù)已賦值網(wǎng)格點(diǎn)內(nèi)插出未知高程的網(wǎng)格點(diǎn)的高程值.
等高線柵格化是指構(gòu)建一個(gè)N*M的網(wǎng)格矩陣,并初始化每個(gè)網(wǎng)格的高程值.具體步驟如下:(1)為每個(gè)網(wǎng)格賦以-999999的初始高程值.若某網(wǎng)格高程值為-999999,則表示該網(wǎng)格高程值還未求解.(2)遍歷每個(gè)網(wǎng)格,判斷是否有等高線經(jīng)過網(wǎng)格,如有1條等高線經(jīng)過,則采用直線內(nèi)插算法把等高線高程值賦給該網(wǎng)格;若大于1條等高線經(jīng)過網(wǎng)格,則取距離加權(quán)平均值,距離值以等高線到網(wǎng)格中心的垂直距離計(jì)算.(3)對無等高線經(jīng)過的網(wǎng)格,如圖2(a)中的g1和g2網(wǎng)格,柵格化后的高程值仍為-999999,圖2(b)中如黑色網(wǎng)格所示.對于這些未獲取高程值的網(wǎng)格,將采用DGT算法通過內(nèi)插計(jì)算得到高程值.
圖2 等高線柵格化
DGT算法的核心是確定動(dòng)態(tài)網(wǎng)格模板,選擇合適的權(quán)函數(shù)與插值函數(shù)來計(jì)算待定點(diǎn)高程值.
3.3.1 參考網(wǎng)格數(shù)量的確定
相對于政府一元獨(dú)大的管理模式,多元主體合作共治的最大特點(diǎn),在于它可以把多元主體置于社會(huì)治理系統(tǒng)之中,實(shí)現(xiàn)多元主體的優(yōu)勢互補(bǔ),克服多元主體各自的局限,實(shí)現(xiàn)社會(huì)治理系統(tǒng)的整體優(yōu)化。而從新時(shí)代構(gòu)建共建共治共享格局的要求講,它可以解決社會(huì)治理領(lǐng)域社會(huì)組織和公眾參與發(fā)展不平衡、不充分問題,最大限度地激發(fā)社會(huì)組織活力和公眾參與熱情,促進(jìn)社會(huì)進(jìn)步,減少社會(huì)治理成本。
參考網(wǎng)格的數(shù)量通常是取4或8個(gè),經(jīng)過對內(nèi)插效果的反復(fù)測試對比分析,DGT算法默認(rèn)取8個(gè)參考網(wǎng)格進(jìn)行計(jì)算.
3.3.2 插值函數(shù)與權(quán)函數(shù)的確定
DGT算法中插值函數(shù)采用式(3),權(quán)函數(shù)采用與距離相關(guān)的權(quán)函數(shù):
式中pi是參考點(diǎn)i的權(quán),(X0,Y0)是待定點(diǎn)的X坐標(biāo)和Y坐標(biāo),(Xi,Yi)是參考點(diǎn)i的X坐標(biāo)和Y坐標(biāo).
3.3.3 確定動(dòng)態(tài)網(wǎng)格模板
DGT算法中,動(dòng)態(tài)網(wǎng)格模板需在一開始確定.圖3為一個(gè)5×5的模板,黑色網(wǎng)格為待插值網(wǎng)格,其值通過搜索周圍的24個(gè)網(wǎng)格內(nèi)插計(jì)算而得(其中可能含有待插值網(wǎng)格).為了避免所選取的數(shù)據(jù)點(diǎn)集中在某個(gè)方向而造成參考點(diǎn)分布不均,DGT算法以待插網(wǎng)格為中心,將5×5平面分為1、2、3、4四個(gè)區(qū),每個(gè)區(qū)有6個(gè)網(wǎng)格.內(nèi)插計(jì)算時(shí),DGT算法將從每個(gè)區(qū)取2個(gè)已知高程的網(wǎng)格作加權(quán)平均,這就克服了傳統(tǒng)移動(dòng)擬和數(shù)據(jù)點(diǎn)偏向的缺點(diǎn).模板矩陣的大小是通過區(qū)的大小確定的.圖4中區(qū)的大小是3×2,即Matrix(區(qū)的長軸)=3,所以模板的大小是5×5.如果區(qū)的大小是4×3,則模板大小是7×7,依此類推.
圖3 動(dòng)態(tài)網(wǎng)格模板
3.3.4 參考網(wǎng)格的搜索
在求解網(wǎng)格的高程值前,需要搜索到8個(gè)參考點(diǎn).DGT算法首先從待插網(wǎng)格周圍的(Matrix-1)圈的網(wǎng)格開始搜索,搜索的順序是:區(qū)1,區(qū)2,區(qū)3,區(qū)4.每個(gè)區(qū)的搜索方式相同,其流程如圖4所示.DGT算法在每個(gè)區(qū)取兩個(gè)已知高程值的網(wǎng)格作加權(quán)平均,若在某個(gè)區(qū)域搜索不到2個(gè)已知高程值的網(wǎng)格做參考點(diǎn),則模板在該區(qū)向外擴(kuò)張1個(gè)網(wǎng)格,即在第Matrix++圈的網(wǎng)格中搜索參考點(diǎn),直到取到2個(gè)已知高程值的參考點(diǎn)為止.當(dāng)搜索超界即超出數(shù)據(jù)邊界時(shí),若仍未找到2個(gè)參考點(diǎn),結(jié)束搜索.
圖4 參考網(wǎng)格的搜索
dtmDX:X方向的網(wǎng)格數(shù)
dtmDY:Y方向的網(wǎng)格數(shù)
minDGXx:高程模型原點(diǎn)即左下角X坐標(biāo)
minDGXy:高程模型原點(diǎn)即左下角Y坐標(biāo)
CellSize:格網(wǎng)間距值
NODATA_value:格網(wǎng)沒有數(shù)值時(shí)的標(biāo)記,用-999999表示
……
為了能更直觀地顯示插值效果,本研究基于OpenGL圖形庫開發(fā)了一個(gè)三維可視化組件,將DGT算法采用C++語言進(jìn)行程序?qū)崿F(xiàn)后,集成到該三維可視化組件中.此處選擇浙江省安溪流域1:10000比例尺地形圖12幅,經(jīng)掃描跟蹤數(shù)字化生成Shapefile文件格式的等高線數(shù)據(jù).圖5中等高線數(shù)據(jù)經(jīng)DGT算法內(nèi)插,生成ARC/INFO ASCII GRID格式的DEM數(shù)據(jù),經(jīng)簡單著色處理,以柵格圖像方式顯示如圖6所示.圖7為DEM數(shù)據(jù)在三維可視化組件中以網(wǎng)格方式顯示的三維效果.
DGT算法采用動(dòng)態(tài)網(wǎng)格模板的方式進(jìn)行參考點(diǎn)搜索,不用進(jìn)行離散點(diǎn)值的求解,也不用采用搜索圓來搜索內(nèi)插所需要的參考點(diǎn),提高了插值的運(yùn)算效率.該算法已被成功運(yùn)用于安溪流域,測試結(jié)果顯示:對于較為連續(xù)和平滑的地形,DGT算法計(jì)算效率較為理想;但對于尖銳變化或極不規(guī)則的地形區(qū)域,該算法會(huì)出現(xiàn)數(shù)據(jù)平滑性“低頻濾波”現(xiàn)象,有待以后不斷改進(jìn)和完善.
[1]葉愛中,夏軍,王綱勝,等.基于數(shù)字高程模型的河網(wǎng)提取及子流域生成[J].水利學(xué)報(bào),2005,36(5):531-537.
[2]張彩霞,楊勤科,段建軍.高分辨率數(shù)字高程模型的構(gòu)建方法[J].水利學(xué)報(bào),2006,37(8):1009-1014.
[3]王家耀.空間信息系統(tǒng)原理[M].北京:科學(xué)出版社,2001:158-160
[4]張超.地理信息系統(tǒng)實(shí)習(xí)教程[M].北京:高等教育出版社,2000:190-193
[5]邱衛(wèi)寧.根據(jù)等高線建立數(shù)字高程模型[J].武漢測繪科技大學(xué)學(xué)報(bào),1994,19(3):199-203.
[6]余鵬,劉麗芳.利用地形圖生產(chǎn)DEM數(shù)據(jù)的研究[J].測繪通報(bào),1998(10):16-18.
[7]李志林,朱慶.數(shù)字高程模型[M].武漢:武漢大學(xué)出版社,2001.