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

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

0 引言

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

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

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

1 并行通視分析算法設計框架

1.1 數(shù)據劃分策略

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

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

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

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

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

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

1.2 并行通視處理

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

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

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

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

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

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

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

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

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

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

2 實驗結果分析

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

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

3 結論

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

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

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

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

［4］ 周啟鳴，劉學軍.數(shù)字地形分析［M］.北京：科學出版社，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］ 張林波，遲學斌，莫則堯，等.并行計算導論［M］.北京：清華大學出版社，2006.

［9］ 江嶺，湯國安，劉凱，等.局部型地形因子并行計算方法研究［J］.地球信息科學學報，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］ 應申，李霖，梅洋，等.增量法地形可視計算與分析［J］.測繪學報，2007，36（2）：192－197.

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

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

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