王朔,修佳鵬,楊正球

(北京郵電大學(xué) 計(jì)算機(jī)學(xué)院(國(guó)家示范性軟件學(xué)院),北京 100089)

0 引言

數(shù)字高程模型(digital elevation model,DEM)是存儲(chǔ)為浮點(diǎn)或整數(shù)值的矩形陣列,用來(lái)表示地面高程的一種實(shí)體地面模型。越來(lái)越多學(xué)者將DEM用于提取某一地區(qū)的水文和地貌特征。洼地指近似封閉的比周圍地面低的地形,在DEM中體現(xiàn)為一塊區(qū)域的高程完全被其相鄰的較高高程包圍,如圖1(a)所示,黃色的單元即為洼地。有些洼地是真實(shí)的自然地形,但大多數(shù)的洼地是偽洼地,因DEM收集和處理中的技術(shù)問題產(chǎn)生[1]。

注:流向同一邊緣點(diǎn)的單元label一致,label 1為環(huán)繞DEM的低地勢(shì)點(diǎn),連接label 1的都是DEM邊緣點(diǎn)label。圖1 平行PF算法示意圖

因DEM中存在洼地或偽洼地,在進(jìn)行水文分析時(shí),二者都會(huì)使DEM中的水流方向不能判斷或水流逆流,導(dǎo)致大部分水文提取算法無(wú)法直接依據(jù)DEM進(jìn)行計(jì)算[2]。所以首先要對(duì)DEM進(jìn)行填洼,這是水文分析的必要預(yù)處理[3]。

填洼指填充洼地,在DEM中通過(guò)提高單元的高程值使得洼地得以填充,從而獲得平整的土地,使得DEM中任何單元都有一條通向DEM邊界的非上升路徑[4],圖1(g)即為填洼后的DEM。近年來(lái),許多填洼算法被提出并應(yīng)用,其中有3類最具代表性[5]。第一類是J&D算法[6],先掃描 DEM識(shí)別出其中的洼地,然后將洼地提升相鄰單元中高程值最低者,但該類算法的時(shí)間復(fù)雜度過(guò)高,為O(N2)。第二類是P&D算法[7]。該算法以一極大高程值水面淹沒整個(gè)DEM,在水流可流出的前提下,不斷迭代去除多余的水,直到每個(gè)單元的高程都具有到相鄰單元的下坡。理論上,P&D算法最壞情況的時(shí)間復(fù)雜度為O(N1.5),相比J&D算法,降低了時(shí)間成本。第三類算法[8-11],即Priority-Flood(PF)算法。該算法引入優(yōu)先隊(duì)列這一數(shù)據(jù)結(jié)構(gòu),隊(duì)列內(nèi)高程值遞增排列。首先將DEM的所有邊緣單元入隊(duì),再將它們依次出隊(duì),對(duì)出隊(duì)單元的相鄰單元進(jìn)行填洼,然后將相鄰單元入隊(duì),直到隊(duì)列為空,算法結(jié)束。對(duì)于浮點(diǎn)數(shù)據(jù),該算法的時(shí)間復(fù)雜度為O(NlogN),對(duì)于整形數(shù)據(jù),該算法的時(shí)間復(fù)雜度為O(N),相較于J&D和P&D算法,大大縮短了時(shí)間成本。

如果單張DEM可以完全適應(yīng)計(jì)算機(jī)的內(nèi)存(random access memory,RAM),那么上述3種算法都可以有效地執(zhí)行填洼操作。但隨著機(jī)載激光雷達(dá)技術(shù)的發(fā)展,DEM的分辨率不斷提高,分辨率高達(dá)1 m或萬(wàn)億單元的超大規(guī)模DEM已可以獲取。雖然計(jì)算機(jī)的處理性能和內(nèi)存顯著增加,但增速遠(yuǎn)遠(yuǎn)低于DEM數(shù)據(jù)規(guī)模增速,DEM通常不能適應(yīng)一臺(tái)普通計(jì)算機(jī)的RAM[12],則現(xiàn)有的3類單機(jī)串行填洼算法無(wú)法實(shí)施填洼。

1 并行PF算法

針對(duì)單機(jī)串行填洼的缺陷,許多學(xué)者提出了J&D算法、P&D算法的并行化。但對(duì)于效率最高的PF算法,受限于洼地可能橫跨切分區(qū)域,直到2016年才由Richard Barnes提出了PF算法的并行實(shí)現(xiàn)。該并行算法的執(zhí)行效率高于原PF算法,同時(shí)實(shí)現(xiàn)了大尺度DEM的填洼,該算法成功對(duì)2萬(wàn)億單元DEM進(jìn)行填洼,這個(gè)測(cè)試集比以前文獻(xiàn)中的任何測(cè)試都要大3個(gè)數(shù)量級(jí),但計(jì)算時(shí)間只增加了1～2個(gè)數(shù)量級(jí),是目前效率最高的平行填洼算法。圖1描述了該算法的全過(guò)程。首先,將單張DEM按指定長(zhǎng)高切分為多張小尺度的、可以適應(yīng)內(nèi)存的虛擬圖塊(tile)。如圖1(b)所示,切分時(shí)遇到了并行填洼最為棘手的情況:一塊洼地被切分到了多張tile里。該算法分為3個(gè)階段進(jìn)行計(jì)算。

第一階段,處理單張瓦片(tile)。生產(chǎn)者(producer)將tile分配給各個(gè)消費(fèi)者(consumers)節(jié)點(diǎn),各節(jié)點(diǎn)依據(jù)文獻(xiàn)[12]的algorithm 1對(duì)所有tile并行執(zhí)行基礎(chǔ)PF算法,得到每個(gè)tile填洼后的高程值(非全局的填洼結(jié)果)。將所有流向該tile同一邊緣點(diǎn)的單元都標(biāo)記相同的標(biāo)簽(label),同時(shí)得到label之間的最小溢出高程值映射(assMap)。最小溢出高程值的定義見式(1)。

assMapLabelA_LabelB=min{max(demc,demn)|c∈LabelA,n∈LabelB且n∈c的鄰點(diǎn)}

(1)

式中:c和n是屬于dem的兩個(gè)柵格點(diǎn);demc和demn分別指c點(diǎn)和n點(diǎn)在dem的高程值;assMapLabelA_LabelB指水流如果在LabelA與LbaelB間流動(dòng)應(yīng)滿足的最小高程值。

第二階段,處理邊緣。依據(jù)文獻(xiàn)[12]的algorithm 2,通過(guò)tile之間的相鄰邊,關(guān)聯(lián)起tile與tile之間label的最小溢出高程值,由此得到了全局label間的最小溢出高程(globalAssMap)。

第三階段,對(duì)globalAssMap執(zhí)行PF算法,得到每個(gè)label應(yīng)填洼到的最小高程值labelWithMin,并將tile中所有l(wèi)abel所屬單元提升到該label對(duì)應(yīng)的labelWithMin值。此時(shí)的tile即為填洼后的最終結(jié)果?？梢愿鶕?jù)用戶需求將tile合并成一張總的DEM,或者直接將tile分塊存放。

2 Barnes并行PF算法的優(yōu)化

本文基于Barnes的并行PF算法(簡(jiǎn)稱原算法),針對(duì)輸入為一張大尺度DEM的場(chǎng)景,在切分策略及l(fā)abel的標(biāo)記策略提出優(yōu)化(簡(jiǎn)稱優(yōu)化算法)。將原算法第一、第二階段的所有操作都限制在一張tile,而無(wú)需等待其他tile構(gòu)建完畢,節(jié)省了原算法需要跨8張tile處理邊緣邊和構(gòu)建全局映射所消耗的時(shí)間,從而提高算法的整體效率。

2.1 帶halo的切分策略

原算法在DEM切分tile的時(shí)候,嚴(yán)格按tile指定的長(zhǎng)和高進(jìn)行切分,見圖1(b)。本文提出一種帶光環(huán)(halo)的切分方法。halo是以tile為中心,來(lái)自其相鄰的8個(gè)tile的相鄰單元格,值為高程值,這些額外的單元像光環(huán)一樣圍繞在原始tile周圍。halo和tile的邊緣邊(edge)各為4條邊。以圖2為例,以4×4為尺寸進(jìn)行DEM切分,切分為了9張tile,以tile4為例,黃色區(qū)域即tile4的halo,紫色區(qū)域?yàn)閠ile4的edge,紫色區(qū)域加藍(lán)色區(qū)域即為tile,tile加黃色區(qū)域即為tileWithHalo。tileWithHalo里存在12個(gè)四角點(diǎn)(corner),corner是指tile中的左上、右上、左下、右下4個(gè)頂點(diǎn),依次定義為corner0、corner1、corner2、corner3。如圖2所示,包括tile的4個(gè)corner以及相鄰tile的8個(gè)corner,分別位于edge與halo。若指定tile的尺寸為a×a,經(jīng)此切分策略切分后,tileWithHalo的尺寸為(a+2)×(a+2)。

圖2 帶halo的切分策略

2.2 label的標(biāo)記及globalAssMap構(gòu)建

原算法要等所有tile都執(zhí)行完畢PF算法,再去其他節(jié)點(diǎn)尋找該tile的相鄰8個(gè)tile來(lái)構(gòu)建全局的globalAssMap。是同步操作,因等待其他tile執(zhí)行完畢,會(huì)消耗過(guò)多時(shí)間。

優(yōu)化算法經(jīng)過(guò)2.1節(jié)的切分策略得到tileWithHalo,tileWithHalo已經(jīng)帶有相鄰8個(gè)tile的edge,所以在構(gòu)建globalAssMap的過(guò)程中,只要按照一定的label的標(biāo)記策略,根據(jù)halo計(jì)算出相鄰tile的edge的label值,就可以將原算法的第一、第二階段整合,將所有操作限制在一張tileWithHalo中進(jìn)行,不涉及其他tile或其他節(jié)點(diǎn),也不需要等待所有tile構(gòu)建完畢再去構(gòu)建globalAssMap。具體的標(biāo)記策略如下。

步驟1:計(jì)算tile的label起始值(startNum)。對(duì)每個(gè)tile分配唯一值num,每個(gè)tile最大label數(shù)量為tileSize×4個(gè)。當(dāng)num等于0時(shí),若label從0開始標(biāo)記,會(huì)與1沖突,因?yàn)?已經(jīng)置為了DEM邊緣低地勢(shì)點(diǎn)label,為保證label唯一性,統(tǒng)一定義為偏移2個(gè)數(shù)值,每個(gè)tile的startNum為num×tileSize×4+2。

步驟2:計(jì)算corner的label值。corner的label值是固定的,要保證edge的corner label值和相鄰tileWithHalo中通過(guò)halo算出的label值一一對(duì)應(yīng)。本文賦予corner一個(gè)固定值(式(2)),labelcornern為該label對(duì)應(yīng)4個(gè)corner的值,可以推出相鄰tile的8個(gè)corner的label值。以tile4為例,計(jì)算出tile4及相鄰tile的startNum后,其corner值見式(3),startNumtilen對(duì)應(yīng)tile下的startNum。

(2)

(3)

步驟3:計(jì)算edge非corner的label值。要保證tileWithHalo通過(guò)edge計(jì)算的label值和相鄰tileWithHalo中通過(guò)halo算出的label值一一對(duì)應(yīng)。以tile4為例,tile4的labeledge_top每個(gè)單元和tile2的labelhalo_bottom都相等。若以計(jì)算corner的方式給edge中每個(gè)非corner一個(gè)固定的label值,這樣雖然簡(jiǎn)單易行,但會(huì)產(chǎn)生過(guò)多的label值,影響算法的效率。本文通過(guò)優(yōu)先隊(duì)列,設(shè)定startNum+4、startNum+4+(edgeSize-2)、startNum+4+(edgeSize-2)×2、startNum+4+(edgeSize-2)×3分別為上、下、左、右4條edge的非corner的起始label(startLabel),如式(4)所示,startLabeln表示該tile對(duì)應(yīng)邊緣邊的起始label,為減少label值,且保證流向tile同一邊緣點(diǎn)的單元label值相同,本文設(shè)計(jì)PFBORD算法,將edges所有單元分別進(jìn)入一個(gè)最小優(yōu)先隊(duì)列,見algorithum 3。

(4)

algorithm 3 PFBORDInput: Q, tileWithHalo, labels, startLabelOutput: labels1:whileQ 不為空 do2:c ← Q.pop3: iflabels(c) == 0then4:foreach n ∈ c 的鄰點(diǎn)且同為 edge 或同為 halo do

步驟4:計(jì)算halo非corner的label值。要保證tileWithHalo通過(guò)halo計(jì)算的label值和相鄰tileWithHalo中通過(guò)edge算出的label值一一對(duì)應(yīng)。流程同步驟3,但startLabel與edge不同,以tile4為例,如式(5)所示,需要與相鄰tileWithHalo的edge對(duì)應(yīng)。將4條halo分別進(jìn)入最小優(yōu)先隊(duì)列,注意此時(shí)的halo要剔除左上、右上、左下、右下tile的corner點(diǎn),然后執(zhí)行PFBORD算法。

(5)

步驟5:構(gòu)建globalAssMap。首先對(duì)tile執(zhí)行基礎(chǔ)PF算法,得到:tile填洼后的高程值(非全局填洼結(jié)果),存放在tileWithHalo中;tile所有單元的label;tile的assMap。由于halo的label相同于相鄰tile edge的label,所以通過(guò)tileWithHalo構(gòu)建edge與halo label的assMap,即可把相鄰tile相關(guān)聯(lián),與原算法的第二階段達(dá)到的結(jié)果一致。最后,將所有tileWithHalo的assMap進(jìn)行合并即為globalAssMap。

經(jīng)過(guò)以上5步,得到關(guān)鍵的globalAssMap,接下來(lái)的步驟與原算法一致,本過(guò)程的偽代碼見algorithm 4。圖3描述了優(yōu)化算法的整體流程。

algorithm 4 label 的標(biāo)記及 globalAssMap 構(gòu)建輸入: tileWithHalo,tileSize,tileOff,offStartNum輸出: assMap,labels,tileWithHalo1: 定義二維數(shù)組 labels 且與 tileWithHalo 維度相同,初始值置為 02: 定義 assMap 為空 Map3: 定義 Q1,Q2,Q3,Q4 為最小優(yōu)先隊(duì)列4: 根據(jù) tileOff 找到該tile相鄰的8個(gè)tile的tileOff,再根據(jù)offStartNum找到相鄰的8個(gè) tile。

圖3 優(yōu)化算法流程

2.3 優(yōu)化算法的Spark實(shí)現(xiàn)

原算法與優(yōu)化算法,實(shí)際上是兩個(gè)MapReduce過(guò)程,可以通過(guò)OpenMP、MPI、Spark等技術(shù)實(shí)現(xiàn),原文選擇了MPI進(jìn)行實(shí)現(xiàn)。MPI的語(yǔ)言更原生,速度更快,但是代碼量巨大,非常繁瑣,在接口的易用性上,也低于spark。用spark來(lái)實(shí)現(xiàn)原算法及優(yōu)化算法,可為并行填洼算法提供一個(gè)新的解決方案。

spark是大數(shù)據(jù)領(lǐng)域基于內(nèi)存的快速而通用的計(jì)算引擎。彈性分布式數(shù)據(jù)集(resilient distributed datasets,RDD)是spark的基本數(shù)據(jù)結(jié)構(gòu),是位于內(nèi)存中的對(duì)象集合,同時(shí)可以將每個(gè)RDD分成多個(gè)分區(qū),不同的分區(qū)保存在不同的集群節(jié)點(diǎn)上。同時(shí),spark為用戶提供了豐富的算子對(duì)RDD進(jìn)行操作。優(yōu)化算法的spark實(shí)現(xiàn),就是將數(shù)據(jù)抽象成RDD并進(jìn)行操作,在由一個(gè)主節(jié)點(diǎn)(master)和多個(gè)從節(jié)點(diǎn)(slave)的集群中,以經(jīng)緯1°作為DEM劃分tile的范圍,tileSize為tile的尺寸,整個(gè)優(yōu)化算法的偽代碼見algorithm5。

algorithm 5 優(yōu)化算法主流程輸入: DEM,tileSize輸出: 填洼后的 tile1: 定義 numOff為 Map //tile 與經(jīng)緯度的映射2: 定義 offStartNum 為 Map //經(jīng)緯度與 Label 起始值的映射3: 獲取 DEM 的維度范圍 (xStart, xend)4: 獲取 DEM 的經(jīng)度范圍 (yStart, yend)5: num ← -16: forx = xStart to xenddo7:fory = yStart to yenddo8:num ← num + 19:numOff.put(num,(x, y))10:offStartNum.put((x, y), num * 4 * tile-Size + 2)//每個(gè) tile 的 label 起始值11: end for12: end for

3 實(shí)驗(yàn)與分析

本文在spark3.1.1集群上,驗(yàn)證優(yōu)化算法的正確性及高效性,集群配置如表1所示,采用GDAL2.4開源庫(kù),選取了3個(gè)數(shù)據(jù)集,分別在正確性、執(zhí)行時(shí)間上將原算法與本算法進(jìn)行對(duì)比。

表1 spark集群配置

3.1 測(cè)試數(shù)據(jù)集

測(cè)試數(shù)據(jù)集覆蓋了中國(guó)經(jīng)度73°E～135°E,緯度18°N～53°N的區(qū)域,在此區(qū)域分別選取10 m、30 m、90 m分辨率的DEM數(shù)據(jù),3個(gè)數(shù)據(jù)集的信息見表2。

表2 測(cè)試數(shù)據(jù)集

3.2 正確性實(shí)驗(yàn)

將優(yōu)化算法分別與原算法、單機(jī)PF算法對(duì)比,來(lái)驗(yàn)證優(yōu)化算法填洼結(jié)果的正確性。首先,分別用3種算法進(jìn)行填洼,得到每個(gè)tile填洼后的結(jié)果,再將tile合并成整張DEM,分別記為DEMoptimize_n、DEMoriginal_n、DEMsingle_n,其中n為DEM里所有單元數(shù)。將DEMoptimize_n分別與DEMoriginal_n、DEMsingle_n對(duì)應(yīng)相等的單元高程值一一相減,若差值和為0,則說(shuō)明兩種算法得到的結(jié)果一致,優(yōu)化算法能得到正確的結(jié)果。因?yàn)閱螜C(jī)PF算法不支持大分辨率的數(shù)據(jù)集,所以只能在small數(shù)據(jù)集上對(duì)比,n為1.3×107,見式(6),result表示結(jié)果差值,tileoptimize_n、tileoriginal_n、tilesingle_n分別表示n點(diǎn)在優(yōu)化算法、原算法、單機(jī)串行算法的高程值。對(duì)比結(jié)果見表3、表4,可以得到result=0,從而驗(yàn)證優(yōu)化算法的填洼結(jié)果正確。

表3 優(yōu)化算法與原算法的正確性對(duì)比

表4 優(yōu)化算法與單機(jī)PF算法的正確性對(duì)比

(tileoptimize_n-tilesingle_n)] (6)

(6)

3.3 執(zhí)行時(shí)間

選取10 m、30 m、90 m數(shù)據(jù)集,在集群分別執(zhí)行優(yōu)化算法與原算法進(jìn)行填洼,各執(zhí)行3次,取3次的平均值作為執(zhí)行時(shí)間。將執(zhí)行時(shí)間進(jìn)行對(duì)比,來(lái)驗(yàn)證優(yōu)化算法效率優(yōu)于原算法,結(jié)果如圖4所示。由分析可知,優(yōu)化算法在3個(gè)數(shù)據(jù)集都至少減少14%的執(zhí)行時(shí)間,在10 m數(shù)據(jù)集優(yōu)化效果最明顯,減少了約37%。效率優(yōu)化體現(xiàn)在label的標(biāo)記及globalAssMap構(gòu)建上,原算法是同步的,而優(yōu)化算法是異步的,優(yōu)化算法將所有的操作限制在一張tileWithHalo上,無(wú)需等待相鄰tile構(gòu)建完畢再去遍歷。

圖4 執(zhí)行時(shí)間對(duì)比結(jié)果

優(yōu)化效果在小數(shù)據(jù)集體現(xiàn)并不明顯,這是因?yàn)閮?yōu)化算法halo的切分策略及l(fā)abel標(biāo)記也會(huì)增加少量耗時(shí)。在數(shù)據(jù)集小的情況下,單張tile僅有144萬(wàn)的單元點(diǎn),同步異步構(gòu)建并不會(huì)出現(xiàn)較大時(shí)間差。優(yōu)化算法更適用于需要大量label標(biāo)記或切分tile較多、構(gòu)建globalAssMap需要等待相鄰tile的執(zhí)行結(jié)果的場(chǎng)景,優(yōu)化算法異步的優(yōu)勢(shì)才會(huì)體現(xiàn)得更為明顯。

4 結(jié)束語(yǔ)

隨著機(jī)載激光雷達(dá)技術(shù)的發(fā)展,DEM規(guī)模的增速遠(yuǎn)遠(yuǎn)大于普通計(jì)算機(jī)RAM的增速,若單張DEM不能完全適應(yīng)RAM,則傳統(tǒng)的單機(jī)PF算法無(wú)法進(jìn)行填洼。Barnes提出的并行PF算法解決了上述問題。本文基于此并行PF算法,針對(duì)輸入一張大規(guī)模DEM且未進(jìn)行劃分的情況,提出了帶光環(huán)的劃分策略,將相鄰tile的邊緣邊劃入到一張tileWithHalo中,由此可以優(yōu)化原算法的第一、第二階段,將label的標(biāo)記及globalAssMap構(gòu)建完全限制在一張tileWithHalo進(jìn)行操作,節(jié)約了等待相鄰tile構(gòu)建完畢及遍歷相鄰tile的時(shí)間。在填洼實(shí)驗(yàn)中,3個(gè)測(cè)試數(shù)據(jù)集都至少減少14%的執(zhí)行時(shí)間。通過(guò)本文的研究,提高了并行填洼的計(jì)算效率,為并行填洼提供了新的視角和方法。