劉思成，蔣志強

（四川大學(xué) 計算機學(xué)院，四川 成都 610065）

近年來，隨著計算機軟件與硬件水平的提高，特別是GPU 處理性能的顯著提高，在三維地球仿真領(lǐng)域，人們對于數(shù)字地形地貌渲染的細節(jié)要求越來越高。在地形仿真中，通常采用高精細度的數(shù)字高程模型（Digital Elevation Model）結(jié)合衛(wèi)星影像貼圖的方法對地球進行仿真，精度越高的衛(wèi)星貼圖與數(shù)字高程越能提高最終的渲染效果。然而，高清衛(wèi)星影像貼圖需求的存儲空間極大，加載高清的衛(wèi)星影像對計算機硬件有更高的要求，高精度的影像衛(wèi)片以TB 為存儲單位，過于精細的影像衛(wèi)片在加載速度、實時性上難以滿足三維仿真場景的需求。此外，通過數(shù)字高程模型和衛(wèi)星影像組合的方式所展示的地形地貌無法清晰、準確地呈現(xiàn)，如山脊、山谷等復(fù)雜的地貌，損失相關(guān)地形細節(jié)。在進行地形地貌相關(guān)方面的研究時，多數(shù)研究重點放在數(shù)字高程的加載以及地形的曲面細分，例如王宇等[1]提出加載地形的細化算法。此外通常采用地形暈渲技術(shù)對地理特征進行提取，然后進行著色處理，例如胡丹丹等[2]提出的數(shù)字地貌暈渲的特征增強。但上述方法通常對于三維地形的實際渲染效果提升有限，特別是在數(shù)字高程與衛(wèi)星影像貼圖結(jié)合時難以清晰直觀地展示地勢的起伏。

為此本文提出了一種基于提取地形特征的地形增強方法。通過特征提取算法實時提取出山脊、山谷等地貌特征，再通過著色器結(jié)合實時光照與法線在原有衛(wèi)星影像基礎(chǔ)上進行融合，從而達到以較小的性能代價來增強地形地貌特征的效果，使得地形更加具有層次感，保留更多的地形細節(jié)。

1 相關(guān)理論與算法

1.1 數(shù)字高程模型

在進行三維地形渲染需要用到相應(yīng)的地形數(shù)據(jù)，地形數(shù)據(jù)通常是以二維結(jié)構(gòu)儲存，是帶有空間位置特征和地形屬性特征的數(shù)字描述。目前較為常用的做法是采用數(shù)字高程模型[3]。數(shù)字高程模型是一種儲存了地形表面高度數(shù)據(jù)的幾何模型，同時還存儲了對應(yīng)的地理位置信息。通過數(shù)字高程模型可以計算出地形的坡向、坡度等數(shù)據(jù)。

數(shù)字高程模型按照數(shù)據(jù)組織表達形式分類，其中常見的主要有兩類：規(guī)則矩形格網(wǎng)和不規(guī)則三角網(wǎng)，如圖1 所示。

圖1 數(shù)字高程模型

規(guī)則網(wǎng)格通常是使用相同的三角形、正方形、矩形等圖元將空間區(qū)域切割為規(guī)則的格柵單元，其中每一個切割出來的單元對應(yīng)一個點位的高度數(shù)據(jù)，在格柵間隔較大時容易丟失地形細節(jié)，比如山川的起伏、河流分布，但規(guī)則網(wǎng)格制作簡單使用方便，是較為通用的做法。不規(guī)則網(wǎng)格是使用有限個采樣點獲取離散數(shù)據(jù)，將這些離散點通過一定的規(guī)則組合成連續(xù)的三角面，能夠很好地反映地形細節(jié)，構(gòu)建過程復(fù)雜，且使用時需要花費大量時間去處理各個三角面之間的位置關(guān)系，因此本文為了高效地進行三維地形可視化采用了以規(guī)則網(wǎng)格為組織形式的數(shù)字高程模型。

綜合本文需求，選取90 M 分辨率的數(shù)字高程模型進行研究。

1.2 特征提取區(qū)域

地形特征點指地表上各類剛性物體起伏所組成獨特性象征和標志，主要包括山頂點（peak）、凹陷點（pit）、脊點（ridge）、谷點（channel）、鞍點（pass），平地點（plane）等6 種，見表1。

表1 地形特征點

山脊是由兩個坡向相反坡度不一的斜坡相遇組合而成條形脊狀延伸的凸形地貌形態(tài)。山脊線是由地形的最高點連線而成，山脊線能夠很好地反映出山脈的整體走勢。山脊線的數(shù)學(xué)描述[4]為：設(shè)D1 為正向上升類型的反向斜坡，D2 為負向下降類型的反向斜坡，設(shè)地面任意一點為P，點P 的高程為Hp，設(shè)它們的交線R=D1∩D2，若存在點P∈R，使得對于任意一點I∈R，都有Hp-Hi>0，則稱R 為山脊線，D=D1∪D2 為山脊。

山谷是指地貌中等高線的彎曲部分向高處凸出處。山谷線與山脊線相反，地形中最低點的連線稱為山谷線。山谷線的數(shù)學(xué)描述[4]為：設(shè)D1 為正向上升類型的反向斜坡，D2 為負向下降類型的反向斜坡，設(shè)地面任意一點為P，點P 的高程為Hp，設(shè)它們的交線R=D1∩D2，若存在點P∈R，使得對于任意一點I∈R，都有Hp-Hi＜0，則稱R 為山谷線，D=D1∪D2 為山谷。

在地形特征提取的相關(guān)研究中，通常著重處理山脊線和山谷線，這兩類特征辨識度高效果明顯，對于突出地勢起伏有著重要作用，而其他的地形特征處理起來復(fù)雜，且不易直觀顯示，因此根據(jù)本文的研究需求也同樣選擇處理山脊線和山谷線這兩個地形特征。

1.3 數(shù)字地形分析

數(shù)字高程模型缺乏普通紙質(zhì)地圖具有的直觀性，所含有的地形信息具有一定的隱蔽性。從數(shù)據(jù)結(jié)構(gòu)上來看，數(shù)字高程模型是一組離散的呈規(guī)則或不規(guī)則三角網(wǎng)的高程數(shù)據(jù)，其數(shù)據(jù)本身無法反映地形地貌的特征。因此利用數(shù)字高程模型進行地形地貌分析需要一種基于離散數(shù)據(jù)的信息分析方法，即數(shù)字地形分析（Digital Ter-rain Analysis，DTA），定義為在數(shù)字高程模型上進行地形屬性計算和特征提取的數(shù)字信息處理技術(shù)，而用于進行地形屬性計算的數(shù)學(xué)模型和地形特征提取算法稱為DEM 解譯算法（DEM Interpretation Algorithms）[5]。DTA 技術(shù)是各種與地形因素相關(guān)的空間模擬技術(shù)的基礎(chǔ)，其結(jié)果通常被理解為更為廣義上的DTM，即數(shù)字地面模型。DEM、DTA、DTM 的關(guān)系如圖2所示。地形的山脊線、山谷線、河流等地形特征線是整個數(shù)字地貌的核心，因此合理地運用數(shù)字地形分析從數(shù)字高程中提取出地形特征在三維地形渲染和水文分析等方面都具有重要意義。

圖2 DEM、DTA、DTM 關(guān)系圖

1.4 法線計算

為了模擬出地形在光照下的陰影以及凹凸不平的效果，傳統(tǒng)方法通常根據(jù)模型單獨計算出一張法線貼圖，即保存了法線信息的一組二維數(shù)據(jù)。而在三維地形渲染中，數(shù)字高程模型數(shù)據(jù)量巨大，單獨渲染法線貼圖以及處理法線貼圖的時間成本過大，因此本文通過著色器實時計算法線（圖3）。

圖3 法線計算示意圖

2 算法流程

基于特征提取的地形增強算法描述如下：首先通過平面曲率法對數(shù)字高程模型進行處理，采用正反地形平面曲率的歸一化值作為權(quán)重來處理山脊和山谷區(qū)域，使其突出顯示地形特征區(qū)域。然后通過濾波對提取過后的特征進行處理，減少噪點，最后通過著色器計算地形法線調(diào)整光照，將提取出的地形特征設(shè)置適當?shù)念伾?、透明度與衛(wèi)星影像貼圖融合。流程如圖4 所示。

圖4 算法流程

2.1 特征提取

目前常見的地形特征提取方法有基于圖像處理技術(shù)的方法、基于地形表面幾何形態(tài)分析的方法、流水物理分析相結(jié)合的方法以及平面曲率法等。優(yōu)劣分析見表2。

表2 特征提取方法分析

綜上比較，根據(jù)本文需要，為了保證著色器實時處理的效率選擇采用效率最高、處理最簡便的平面曲率法。采用平面曲率法所獲取的地形地貌特征線較真實地形的準確度不如基于幾何形態(tài)與地表流水模擬法，但本文重在反映地形的大致走勢，因此可以滿足本文的研究需求。

平面曲率是指地形表面上任意一點，用過該點的水平面沿水平方向切地形表面所得到的曲線在該點的曲率值。它描述的是地表曲面沿水平方向的彎曲、變化情況。如果平面曲率為正值，其在水平方向上導(dǎo)數(shù)為0且二次偏導(dǎo)數(shù)大于0，則柵格屬于山脊點；若平面曲率為負值，其在水平方向上導(dǎo)數(shù)為0 且二次偏導(dǎo)數(shù)小于0，則柵格點屬于山谷點。即山脊是平面曲率為正值且在局部區(qū)域方向軸上取最大值的區(qū)域；山谷是平面曲率為負值且在局部區(qū)域方向軸上為最小值的區(qū)域。

本文采用平面曲率法提取山脊線、山谷線。地形曲率是指地形曲面在各個界面方向上的形狀，凹凸變化的反映，是平面點位的函數(shù)。反映了地形結(jié)構(gòu)和形態(tài)，在地表過程模擬、水文、土壤等領(lǐng)域有著重要的應(yīng)用價值和意義。

平面曲率又稱水平凸度，是等高線的方向變化率，可以反映地形的高低起伏[6]。當平面曲率點大于0時為山脊線點，小于0 時為山谷線點，通過計算得到山脊線點、山谷線點連接起來即可得到相應(yīng)的山脊線、山谷線。

假設(shè)地形表面S 的曲面方程z=f（x，y），點P 為曲面S 上任意一點，將曲面S：z=f（x，y）的各階偏導(dǎo)數(shù)用下列符號代替：

根據(jù)平面曲率的推導(dǎo)過程及其數(shù)學(xué)意義，曲面S與XOY 水平面相交線在P 點處的二階導(dǎo)數(shù)即是平面曲率，計算公式如下：

提取結(jié)果如圖5 所示。

圖5 處理效果對比圖

2.2 特征處理

根據(jù)平面曲率法獲取的地形特征數(shù)據(jù)有大量噪點，因此通過過濾來進行處理。

地形特征提取后獲得的是一組二維數(shù)據(jù)，因此對這組數(shù)據(jù)的處理，本文采用圖像濾波的方法。常見的圖像濾波方法如下：

（1）均值濾波[7]，也稱為線性濾波，屬于低通濾波。具體流程是選取一個大小適當?shù)膎*n 的模板，采樣這個模板中的每一個點，取其平均值作為該點的像素值。這種方法簡單易懂，但不能徹底去除噪點，同時還會降低圖像的清晰度。

（2）中值濾波[8]是一種非線性濾波，具體流程是選取一個大小適當?shù)膎*n 的模板，采樣這個模板中的每一個點，取其中間值作為該點的像素值。這種方法對于校驗噪聲有著顯著的效果，但是對圖像的破壞較大，易出現(xiàn)大量斷裂，導(dǎo)致圖像不連續(xù)，且無法處理圖像中存在的跳躍變化。

（3）高斯濾波是一種線性平滑濾波，屬于低通濾波。具體流程是選取一個大小適當?shù)膎*n 的模板，采樣這個模板中的每一個點，并計算模板內(nèi)每一個加權(quán)平均值，用得到的值去替代該點的像素值，該方法能夠最大程度地保留圖像的總體灰度分布。

對比幾種濾波方法，通過平面曲率提取出的地形特征數(shù)據(jù)存在大量的噪聲，而提取出的地形地貌特征通常屬于高頻分量，因此在對初步提取出的地形特征數(shù)據(jù)處理中，應(yīng)采用低通濾波進行處理。考慮到濾波時計算效率以及本文對于地形特征提取的精度要求不高。本文采用3*3 模板的均值濾波進行處理。結(jié)果如圖6 所示。

圖6 過濾效果對比圖

2.3 衛(wèi)星影像融合

通過前面章節(jié)提取出的地形特征經(jīng)過處理后與衛(wèi)星影像進行融合，同時實時計算法線使地形與太陽光所生成的陰影更加立體、真實，可以以較小的性能開銷獲得較大的地形地貌特征細節(jié)的提升。首先著色器（shader）讀取地形高程，通過平面曲率法提取出地形特征，然后通過濾波再處理得到過濾后的地形特征，再為地形特征中的山脊、山谷設(shè)定適當?shù)念伾屯该鞫龋ū疚膶⑸郊诡伾O(shè)定為白色，透明度0.6，山谷設(shè)定為黑色，透明度0.2），得到最終的地形特征圖，設(shè)置光線參數(shù)，再根據(jù)計算出的法線將地形特征圖與衛(wèi)星影像混合后輸出到數(shù)字高程模型上。通過透明度適當?shù)那驶旌显谛l(wèi)星影像地圖上增加三維地形的整體立體感，同時能夠較好地顯示山脊山谷的特征，展示地勢起伏。如圖7 所示。

圖7 融合流程圖

核心代碼如下：

3 結(jié)果與分析

仿真實驗相關(guān)參數(shù)見表3。

表3 算法參數(shù)

渲染效果對比如圖8 所示。

圖8 渲染效果對比

關(guān)閉垂直同步后在三維地球中進行漫游實時幀率，如圖9 所示。

圖9 實時幀率對比

運行1 000 s 平均幀率，結(jié)果如圖10 所示。

圖10 幀率對比

從對比圖中可以看到通過提取地形特征進行加強過后的三維地形效果有了顯著提升，相比原始衛(wèi)星影像地形的立體感提高了很多，其次地勢起伏、高低落差更加明顯，有利于準確展示地形地貌的細節(jié)。在渲染效率方面影響較小。

4 結(jié)束語

開啟地形增強后對于整個系統(tǒng)的效率影響較小，但對整體的地形仿真效果有了很大提升，通過此方法提取的地形特征適用于提高地形渲染效果，對展示地形地勢起伏有較大作用，提高了地形的立體感，對于作戰(zhàn)指揮或是地理信息展示具有重要意義。目前本方法著重對山脊、山谷進行了強化，對于其他的地形特征表現(xiàn)的準確度和細節(jié)保留程度還有提高的空間，其次選取濾波方法取決于如何平衡效率與渲染效果，本文選擇了較為簡單高效的濾波方法，但是過濾效果比較一般，此外針對不同地區(qū)所需的特征線透明度、顏色等參數(shù)設(shè)置還有待進一步研究。