基于DEM的分布式并行通視分析算法研究

張剛，湯國安＊，宋效東，楊坤

（1.南京師范大學(xué)虛擬地理環(huán)境教育部重點(diǎn)實驗室，江蘇 南京 210023；2.南京師范大學(xué)計算機(jī)科學(xué)與技術(shù)學(xué)院，江蘇 南京 210023）

0 引言

地形通視性分析（Terrain Inter－Visibility Analysis）是指基于DEM數(shù)據(jù)判斷地形上任意兩點(diǎn)之間是否可見的技術(shù)方法，通視的條件取決于視點(diǎn)與目標(biāo)點(diǎn)間是否存在妨礙視線的障礙物。通視分析實質(zhì)上屬于對地形進(jìn)行最優(yōu)化處理的范疇［1］。目前，基于規(guī)則格網(wǎng)DEM的地形通視性分析已經(jīng)廣泛應(yīng)用于通信、軍事、房地產(chǎn)、考古、景觀設(shè)計等多個領(lǐng)域。但隨著空間數(shù)據(jù)采集技術(shù)的迅速發(fā)展，海量、高精度DEM數(shù)據(jù)的應(yīng)用越來越廣泛，通視分析算法的計算量也呈現(xiàn)出指數(shù)級增長趨勢，傳統(tǒng)的計算機(jī)處理技術(shù)已經(jīng)不能有效地提高基于海量DEM數(shù)據(jù)的通視分析執(zhí)行效率［2］。因此，如何利用現(xiàn)有計算資源，提高計算效率，降低算法時間復(fù)雜度成為通視分析要解決的關(guān)鍵問題。

近年來，隨著海量空間數(shù)據(jù)分析的不斷增長以及高性能地學(xué)應(yīng)用的不斷推動，并行空間分析方法成為高性能地學(xué)計算的發(fā)展趨勢［3］。并行空間分析方法是將GIS空間分析方法與并行計算技術(shù)相融合的過程，旨在通過多種計算資源解決復(fù)雜的空間分析問題，提高空間數(shù)據(jù)處理的速度和質(zhì)量。通視分析作為地形分析的重要組成部分，也是空間分析不可或缺的內(nèi)容［4］。研究并行通視分析算法，成為提高海量數(shù)據(jù)下通視分析效率的有效途徑。

針對通視分析的高效計算問題，國內(nèi)外學(xué)者做了一定研究。Floriani等［5］在MIMD架構(gòu)基礎(chǔ)上，通過靜態(tài)和動態(tài)兩種并行數(shù)據(jù)劃分策略，提出了一種基于TIN數(shù)據(jù)的并行可視性分析算法，但是由于TIN數(shù)據(jù)存儲結(jié)構(gòu)的復(fù)雜性，算法執(zhí)行效率有待進(jìn)一步提高，使用范圍有一定的局限性。Kidner等［6］提出一種基于數(shù)據(jù)并行的反向可視性分析算法，各處理節(jié)點(diǎn)同步執(zhí)行串行通視算法，獲取由面到區(qū)域的可視域，該算法有效利用了各計算節(jié)點(diǎn)的計算能力并取得了較優(yōu)的加速比，但是計算精度相對不高。由于高分辨率海量DEM數(shù)據(jù)的點(diǎn)對點(diǎn)通視分析算法計算量巨大，Mills等［7］針對LOS視線設(shè)計了數(shù)據(jù)劃分策略，對傳統(tǒng)點(diǎn)點(diǎn)通視分析算法進(jìn)行了并行化處理，由于測試環(huán)境較為簡單，實驗結(jié)果有待進(jìn)一步驗證。在此基礎(chǔ)上，本文從負(fù)載均衡角度，通過分析DEM數(shù)據(jù)存儲的結(jié)構(gòu)特征，提出了一種面向分布并行環(huán)境的數(shù)據(jù)劃分策略，構(gòu)建了并行Bresenham地形通視分析算法，并以全國90mSRTM為數(shù)據(jù)源，分析了并行通視分析的運(yùn)行效率，為并行環(huán)境下的地形通視分析提供了新的思路。

1 并行通視分析算法設(shè)計框架

1.1 數(shù)據(jù)劃分策略

并行計算是指同時使用多種計算資源解決復(fù)雜問題的過程，其包括共享存儲模型、消息傳遞模型和數(shù)據(jù)并行模型3類［8］。由于柵格數(shù)據(jù)空間分析具有計算數(shù)據(jù)密集且數(shù)據(jù)組織規(guī)律性強(qiáng)的特點(diǎn)［9］，大部分柵格數(shù)據(jù)空間分析的并行模式均采用數(shù)據(jù)并行處理方式［10］。數(shù)據(jù)并行是將串行任務(wù)需要處理的原始大數(shù)據(jù)集劃分為從節(jié)點(diǎn)處理的子數(shù)據(jù)集，各處理節(jié)點(diǎn)獨(dú)立執(zhí)行串行程序，獲取計算結(jié)果。數(shù)據(jù)并行策略的基礎(chǔ)是數(shù)據(jù)劃分，針對DEM數(shù)據(jù)的矩陣式特征，本文從數(shù)據(jù)并行的角度對基于DEM的分布式并行通視算法進(jìn)行分析。DEM數(shù)據(jù)劃分主要包括規(guī)則劃分和不規(guī)則劃分兩種，其中規(guī)則劃分又分為一維劃分和二維劃分（圖1）。從數(shù)據(jù)結(jié)構(gòu)角度，DEM數(shù)據(jù)是一個m行×n列的矩陣，一般通過數(shù)組的形式進(jìn)行組織。因此，規(guī)則劃分相較于不規(guī)則劃分更易實現(xiàn)且有利于并行計算系統(tǒng)的負(fù)載平衡。實驗表明，在規(guī)則劃分中，一維數(shù)據(jù)劃分程序的執(zhí)行效率優(yōu)于二維劃分，而DEM數(shù)據(jù)在計算機(jī)物理內(nèi)存中是按行存儲的，故本文選擇一維行劃分作為其數(shù)據(jù)劃分方法：將數(shù)據(jù)集按行平均劃分為n等份于n個處理節(jié)點(diǎn)，如果存在剩余行，則分配給最后一個處理節(jié)點(diǎn)。

圖1 DEM數(shù)據(jù)劃分方法Fig.1 Regular partition of DEM

針對基于DEM的通視分析算法，其有效數(shù)據(jù)區(qū)為視點(diǎn)與目標(biāo)點(diǎn)所對應(yīng)的外接矩形區(qū)。如果在數(shù)據(jù)處理過程中將整個DEM數(shù)據(jù)讀取后再進(jìn)行并行通視分析操作，會大大增加程序負(fù)載量及計算冗余，降低程序的執(zhí)行效率。因此，在對DEM進(jìn)行數(shù)據(jù)劃分操作前，首先完成有效數(shù)據(jù)區(qū)的裁切處理（圖2）。在數(shù)據(jù)裁剪過程中，通過視點(diǎn)與目標(biāo)點(diǎn)位置計算外接矩形及其左上角點(diǎn)的絕對坐標(biāo)，更新DEM數(shù)據(jù)頭文件信息并獲取有效數(shù)據(jù)區(qū)，合并生成新的DEM數(shù)據(jù)，為數(shù)據(jù)劃分提供數(shù)據(jù)源。

圖2 數(shù)據(jù)裁切過程Fig.2 Process of extracting DEM

在數(shù)據(jù)并行過程中，主節(jié)點(diǎn)將數(shù)據(jù)劃分后的子數(shù)據(jù)集發(fā)送給從節(jié)點(diǎn)，從節(jié)點(diǎn)獨(dú)立完成數(shù)據(jù)處理，但針對海量DEM數(shù)據(jù)，每個從節(jié)點(diǎn)所能容納的最大計算量是有限的，為了保證每個計算節(jié)點(diǎn)處理的數(shù)據(jù)不超過其本身所能容納的最大計算量，本文還設(shè)計了內(nèi)存限制參數(shù)，以保證每個計算節(jié)點(diǎn)的最大計算效率，以獲取最優(yōu)加速比。根據(jù)數(shù)據(jù)量大小、計算節(jié)點(diǎn)數(shù)目以及內(nèi)存限制，可以得到每個處理節(jié)點(diǎn)循環(huán)計算的次數(shù)：

在理想狀態(tài)下，數(shù)據(jù)劃分后計算節(jié)點(diǎn)間的處理過程是相互獨(dú)立的，不需要進(jìn)行通信交互。然而，地形通視分析算法需要實時判斷每個子數(shù)據(jù)集中目標(biāo)點(diǎn)的可視性，以便及時結(jié)束運(yùn)算，提高計算效率，因此，子數(shù)據(jù)結(jié)果集之間需要實時通信，即結(jié)果數(shù)據(jù)通信。

1.2 并行通視處理

1.2.1 并行通視分析原理 地形通視性分析主要是研究視點(diǎn)與目標(biāo)點(diǎn)是否可視的問題，具有簡單復(fù)雜性、不可逆性以及可視不變性三大特征［4］。由于目前規(guī)則格網(wǎng)DEM的廣泛應(yīng)用，基于DEM的地形通視性算法眾多［11－15］，其基本原理為：做從視點(diǎn)V到目標(biāo)點(diǎn)T的射線，判斷視線所經(jīng)過的地面點(diǎn)高程與相應(yīng)視線上點(diǎn)的高程之間的關(guān)系，如果在視點(diǎn)與目標(biāo)點(diǎn)之間，有任意一點(diǎn)的高程高于射線高程，則兩點(diǎn)不通視，否則通視?；贒EM通視分析算法基本原理大致相同，不同之處主要集中在高程內(nèi)插方法和可視判斷原則方面，其中，最為著名的是Bresenham算法，該算法基于整數(shù)運(yùn)算，將柵格單元看做是同質(zhì)的，獲取柵格高程的內(nèi)層循環(huán)緊密而且高效，避免了其他通視算法的插值運(yùn)算過程。因此，本文選擇Bresenham算法作為并行通視分析的核心算法。

基于DEM的通視分析并行算法采用數(shù)據(jù)并行的處理方式，按照一維行數(shù)據(jù)劃分策略，根據(jù)內(nèi)存限制，將更新后的DEM的數(shù)據(jù)劃分為不同的子數(shù)據(jù)集合，各個計算機(jī)處理節(jié)點(diǎn)分別讀取子數(shù)據(jù)集數(shù)據(jù)后，同步執(zhí)行Bresenham通視分析算法，并將計算結(jié)果發(fā)送到主進(jìn)程，主進(jìn)程根據(jù)各個進(jìn)程的返回結(jié)果判斷視點(diǎn)與目標(biāo)點(diǎn)之間的可視情況。如果子進(jìn)程計算結(jié)果存在不可視情況，則視點(diǎn)與目標(biāo)點(diǎn)之間不可視；否則兩點(diǎn)可視。并行通視分析原理如圖3所示。

圖3 并行通視分析原理Fig.3 Principle of parallel visibility analysis

1.2.2 Bresenham通視分析算法 假設(shè)直線位于第一象限，從起點(diǎn)（x1，y1）到終點(diǎn)（x2，y2）。直線方程可以表示為y＝kx＋b（k∈［0，1］）。式中：

如圖4所示，當(dāng)直線y＝kx＋b沿x主軸增加1個單位，即xi＋1＝xi＋1時，在數(shù)學(xué)意義上其y副軸也會相應(yīng)的增加斜率大小k。由于計算機(jī)以柵格化的方式進(jìn)行圖像表達(dá)，實際點(diǎn)位在繪制過程中都會存在一定的誤差w（w∈［－0.5，0.5］），此時，yi＋1的真值為yi＋1＝y(tǒng)i＋k＋w，而點(diǎn)yi＋1在計算機(jī)屏幕進(jìn)行繪制時，只可能選擇yi＋1＝y(tǒng)i＋1或者yi＋1＝y(tǒng)i，選擇的原則為k＋w遵從四舍五入，分段函數(shù)表達(dá)式為：

yi＋1繪制過程中造成的誤差記為wnew，則：

在程序執(zhí)行過程中，為消除浮點(diǎn)數(shù)運(yùn)算，將上述誤差公式兩邊同乘以dx，并令φ＝dx＊w，則：

圖4 Bresenham算法原理Fig.4 Principle of Bresenham algorithm

通過Bresenham算法可以得到視線經(jīng)過的所有柵格單元的行列號，利用DEM及視線所經(jīng)過的柵格單元行列號獲取地表高程TerrainElev；如果當(dāng)前點(diǎn)高程值為NoDataValue，即空值，則默認(rèn)此點(diǎn)可視，繼續(xù)計算下一點(diǎn)，否則，通過比較視線上當(dāng)前點(diǎn)高程ViewElev與TerrainElev的大小判斷是否可視。

綜上，Bresenham通視算法思想為：1）計算dx＝x2－x1，dy＝y(tǒng)2－y1，計算誤差初始值φ＝0；2）求取直線下一點(diǎn)位置xi＋1＝xi＋1，若φ＋dy≥dx，yi＋1＝y(tǒng)i＋1，否則yi＋1＝y(tǒng)i；3）計算點(diǎn)（xi＋1，yi＋1）的地表高程TerrainElev及視線高程ViewElev，比較兩者大?。喝鬡iewElev＞TerrainElev，則視點(diǎn)與目標(biāo)點(diǎn)不可視，否者執(zhí)行步驟4；4）更新誤差值φ，若φ＋dy≥dx，φ＝φ＋dy－dx，否則φ＝φ＋dy，轉(zhuǎn)步驟2，直至點(diǎn)（x2，y2），程序結(jié)束。

2 實驗結(jié)果分析

筆者在Microsoft Visual Studio 2010環(huán)境下，利用C?、GDAL1.6開源庫及MPICH2并行計算平臺實現(xiàn)了并行通視分析算法。實驗的硬件環(huán)境為分布式并行環(huán)境，其主節(jié)點(diǎn)服務(wù)器CPU為Intel Xeon E5645，主頻2.4GHZ，內(nèi)存32G，核心數(shù)24個；從節(jié)點(diǎn)服務(wù)器CPU為Xeon E5620，主頻2.4GHZ，內(nèi)存16G，核心數(shù)8個，集群環(huán)境共包含20個從節(jié)點(diǎn)，主節(jié)點(diǎn)與從節(jié)點(diǎn)之間通過千兆網(wǎng)進(jìn)行連接。采用全國90m分辨率SRTM陸地表面地形數(shù)據(jù)作為數(shù)據(jù)源，其數(shù)據(jù)大小為13.2GB（74 869列×46 784行），從中隨機(jī)選取3條視線（長度分別為3 836km、4 473km、4 283km）分析并行通視分析代碼效率，如圖5、圖6所示。

在實驗過程中，為了對通視分析并行執(zhí)行效率進(jìn)行有效的評估，提出了性能加速比的概念模型。設(shè)N為實驗過程中總進(jìn)程數(shù)目，P為當(dāng)前進(jìn)程數(shù)目，T（N，P）為當(dāng)進(jìn)程數(shù)為P時通視分析并行執(zhí)行時間，T（N，1）為通視分析在進(jìn)程數(shù)為1時的執(zhí)行時間，也就是串行執(zhí)行時間。定義Radio為通視分析并行算法執(zhí)行加速比，則：

圖5 并行通視算法性能測試（通視情況）Fig.5 Test of parallelism of inter－visibility analysis algorithm （visible）

圖6 并行通視算法性能測試（不通視情況）Fig.6 Test of parallelism of inter－visibility analysis algorithm （invisible）

從通視分析并行化算法性能測試圖中可以看到：1）不論視點(diǎn)與目標(biāo)點(diǎn)之間是否通視，隨著處理節(jié)點(diǎn)數(shù)的增加，各類型并行化算法執(zhí)行時間的降低速率變小，直至趨向于平緩。處理節(jié)點(diǎn)數(shù)小于32時，執(zhí)行時間變化突出，變化速率較高，而當(dāng)處理節(jié)點(diǎn)數(shù)大于32時，逐漸趨向于同一水平狀態(tài)，執(zhí)行時間變化緩慢甚至存在效率降低的趨勢。因此，基于海量地形數(shù)據(jù)進(jìn)行分布式并行通視分析時，其計算效率與進(jìn)程數(shù)具有一定的關(guān)系。2）地形通視分析并行化后，不通視情況下最大加速比達(dá)到7，通視情況下最大加速比達(dá)到4，計算效率明顯高于單個處理節(jié)點(diǎn)，即高于串行算法執(zhí)行效率。在通視與不通視狀況下，并行通視算法加速比存在不一致的情況，這是因為算法的并行性能在一定程度上受到地形數(shù)據(jù)的影響。從研究區(qū)域的地形復(fù)雜性角度分析，當(dāng)研究區(qū)域地形起伏較大時，不通視狀況下的計算效率體現(xiàn)得尤為突出；從程序設(shè)計的角度分析，在程序設(shè)計過程中采用了單一進(jìn)程終止策略，通視分析并行算法在運(yùn)行過程中，每一個處理節(jié)點(diǎn)都會執(zhí)行相同的通視算法，當(dāng)任意一個節(jié)點(diǎn)獲得不通視的結(jié)果后程序就會自動結(jié)束此節(jié)點(diǎn)的運(yùn)行過程，并將不通視的結(jié)果通過MPI通信機(jī)制發(fā)送給主節(jié)點(diǎn)，由于MPI函數(shù)庫不存在同時終止所有處理節(jié)點(diǎn)運(yùn)行的機(jī)制，因此當(dāng)所有子節(jié)點(diǎn)將返回結(jié)果發(fā)送到主進(jìn)程后，主進(jìn)程對其從進(jìn)程的發(fā)送結(jié)果進(jìn)行分析，只要有一個進(jìn)程的返回結(jié)果是不通視，則視點(diǎn)與目標(biāo)點(diǎn)之間即為不通視，否則兩點(diǎn)之間通視。通過這種單一進(jìn)程終止策略可有效地提高并行算法的執(zhí)行效率。從圖中可以看出，隨著處理節(jié)點(diǎn)數(shù)目的增加，程序的執(zhí)行效率存在一定的波動，這主要是由于代碼執(zhí)行效率比較高，進(jìn)程間通信時間不穩(wěn)定直接影響著并行程序的執(zhí)行加速比。

3 結(jié)論

本文針對基于DEM的分布式并行通視分析算法進(jìn)行了分析研究，提出了一種通用的有效數(shù)據(jù)可達(dá)的DEM數(shù)據(jù)劃分策略，實現(xiàn)了Bresenham通視分析并行算法設(shè)計，經(jīng)驗證，本研究所設(shè)計實現(xiàn)的并行算法在計算時間以及加速比方面能夠體現(xiàn)出較好的優(yōu)勢，大大提高了海量數(shù)據(jù)下通視分析的執(zhí)行效率，為并行環(huán)境下的地形通視分析提供了新的思路。但是，如何進(jìn)一步提高并行I／O的讀取效率以及減少進(jìn)程間通信時間，進(jìn)一步提高并行通視分析算法加速比，還需要進(jìn)一步研究。

［1］ 李志林，朱慶.數(shù)字高程模型［M］.武漢：武漢大學(xué)出版社，2001.

［2］ 呂品，張金芳，魯敏，等.基于視線的地形可視性分析研究［J］.計算機(jī)工程與應(yīng)用，2006（8）：223－226.

［3］ 宋效東，劉學(xué)軍，湯國安，等.DEM與地形分析的并行計算［J］.地理與地理信息科學(xué)，2012（4）：1－7.

［4］ 周啟鳴，劉學(xué)軍.數(shù)字地形分析［M］.北京：科學(xué)出版社，2006.249－266.

［5］ FLORIANI L D，MONTANI C，SCOPIGNO R.Parallelizing visibility computations on triangulated terrains［J］.International Journal of Geographical Information Systems，1994，8（6）：515－531.

［6］ KIDNER D B，RALLINGS P J，WARE J A.Parallel processing for terrain analysis in GIS：Visibility as a case study［J］.GeoInformatica，1997，1（2）：183－207.

［7］ MILLS K，F(xiàn)OX1G，HEIMBACH R.Implementing an intervisibility analysis model on a parallel computing system［J］.Computers ＆ Geosciences，1992，18（8）：1047－1054.

［8］ 張林波，遲學(xué)斌，莫則堯，等.并行計算導(dǎo)論［M］.北京：清華大學(xué)出版社，2006.

［9］ 江嶺，湯國安，劉凱，等.局部型地形因子并行計算方法研究［J］.地球信息科學(xué)學(xué)報，2012，6（14）：761－767.

［10］ CLARKE K C.Geocomputation′s future at the extremes：High performance computing and nanoclients［J］.Parallel Compu－ting，2003，29（10）：1281－1295.

［11］ FLORIANI L D，MARZANO P，PUPPO E.Line of sight communication on terrain models［J］.International Journal of Geographic Information Systems，1994，8（4）：329－342.

［12］ 應(yīng)申，李霖，梅洋，等.增量法地形可視計算與分析［J］.測繪學(xué)報，2007，36（2）：192－197.

［13］ 王智杰，邱曉剛，李革.RSG地形通視性快速算法設(shè)計［J］.計算機(jī)仿真，2004，21（12）：92－95.

［14］ 魯敏，張金芳，范植華，等.基于DEM的視域分析與計算［J］.計算機(jī)仿真，2006，23（5）：171－175.

［15］ 劉旭紅，劉玉樹，張國英.利用最大仰角插值技術(shù)的通視性分析算法研究［J］.計算機(jī)輔助設(shè)計與圖形學(xué)學(xué)報，2005，17（5）：971－975.