周小華,王學(xué)志,周園春,孟 珍

(1.中國(guó)科學(xué)院計(jì)算機(jī)網(wǎng)絡(luò)信息中心,北京 100083;2.中國(guó)科學(xué)院大學(xué),北京 100049)

1 引言

因碳排放增加導(dǎo)致的全球氣候變暖問題日益突出,正對(duì)人類現(xiàn)在及未來的社會(huì)經(jīng)濟(jì)活動(dòng)產(chǎn)生深遠(yuǎn)影響。2016年全世界178個(gè)締約方共同簽署了《巴黎協(xié)定》,呼吁國(guó)際社會(huì)共同努力將全球溫度上升幅度控制在2 ℃以內(nèi)。我國(guó)也于2020年提出了2030年碳達(dá)峰與2060年碳中和的目標(biāo)[1]。通過碳衛(wèi)星檢測(cè)碳排放的方式具有連續(xù)、穩(wěn)定、大范圍、不受地形限制等優(yōu)點(diǎn),為全球各個(gè)區(qū)域在不同時(shí)期的碳排放量化評(píng)估提供了重要的技術(shù)手段。

近年來,中國(guó)、美國(guó)、加拿大及日本等國(guó)家相繼發(fā)射多顆碳衛(wèi)星,目前仍在軌運(yùn)行的有OCO-2(美國(guó),2014年)、TANSAT(中國(guó),2016年)、GOSAT-2(日本,2018年)、OCO-3(美國(guó),2018年)、DQ-1(中國(guó),2022年)等[2]。大部分碳衛(wèi)星基于自身運(yùn)行軌道以線性散點(diǎn)的形式采集并存儲(chǔ)探測(cè)數(shù)據(jù),在進(jìn)行大區(qū)域CO2濃度分布量化評(píng)估時(shí),需要通過空間插值來彌補(bǔ)缺失區(qū)域值。

在上述衛(wèi)星中,OCO-3衛(wèi)星時(shí)空覆蓋范圍更廣,數(shù)據(jù)經(jīng)過嚴(yán)格的校準(zhǔn)處理,具有更高的測(cè)碳精度,且所有數(shù)據(jù)均可公開訪問,是目前進(jìn)行全球碳排放檢測(cè)與研究的重要數(shù)據(jù)源[3-5]。因此,本文以O(shè)CO-3衛(wèi)星數(shù)據(jù)為例進(jìn)行面向碳衛(wèi)星數(shù)據(jù)的Kriging插值算法的優(yōu)化研究。

2 相關(guān)工作

空間插值是指利用周圍已知樣本點(diǎn)基于擬合函數(shù),對(duì)未知點(diǎn)的值進(jìn)行預(yù)測(cè)的方法。常用的空間插值方法從原理上可分為2大類:確定性方法和地學(xué)統(tǒng)計(jì)方法。確定性方法基于已知點(diǎn)的相似度或曲面的平滑程度來預(yù)測(cè)未知點(diǎn)的值,如反距離權(quán)重法、徑向基函數(shù)插值法等;地學(xué)統(tǒng)計(jì)方法則通過量化樣本數(shù)據(jù)點(diǎn)之間的空間自相關(guān)性,構(gòu)造空間結(jié)構(gòu)模型,進(jìn)而對(duì)未知點(diǎn)的值做出預(yù)測(cè),如Kriging插值法[6]。

Kriging插值算法可以給出最優(yōu)線性無偏估計(jì)BLUP(Best Linear Unbiased Prediction),在插值時(shí)綜合考慮了距離及空間自相關(guān)性等因素,具有較高的插值精度[7]。Kriging算法根據(jù)前提條件與應(yīng)用場(chǎng)景的不同又可分為普通Kriging算法、泛Kriging算法、協(xié)同Kriging算法及析取Kriging算法等。其中普通Kriging算法是ArcGIS推薦使用的默認(rèn)插值算法,應(yīng)用最為廣泛,為此本文將基于普通Kriging算法對(duì)OCO-3數(shù)據(jù)的大區(qū)域尺度插值進(jìn)行研究。

Kriging插值過程比較復(fù)雜:在插值前需要計(jì)算所有樣本點(diǎn)兩兩之間的距離與協(xié)方差,然后擬合得到半變異函數(shù)模型,進(jìn)而對(duì)未知點(diǎn)進(jìn)行插值預(yù)測(cè),其中包含矩陣乘積和矩陣求逆等耗時(shí)操作。假設(shè)有n個(gè)數(shù)據(jù)點(diǎn)和m個(gè)未知點(diǎn),Kriging插值的整體時(shí)間復(fù)雜度為O(m·n3),在大區(qū)域尺度上直接應(yīng)用Kriging算法進(jìn)行空間插值計(jì)算量很大[8]。

以O(shè)CO-3數(shù)據(jù)為例,覆蓋中國(guó)區(qū)域的約有150萬個(gè)數(shù)據(jù)點(diǎn),假設(shè)已知樣本點(diǎn)占一半,即75萬個(gè)點(diǎn),僅計(jì)算這75萬個(gè)點(diǎn)中任意兩點(diǎn)之間的距離與協(xié)方差需要約5 000億次運(yùn)算,基于普通配置的單臺(tái)服務(wù)器完成上述運(yùn)算需要30天左右,耗時(shí)難以接受。

目前針對(duì)Kriging插值算法的并行加速優(yōu)化方案主要有2種。一種是在高性能計(jì)算環(huán)境中改造Kriging算法,利用眾核處理器、超高速網(wǎng)絡(luò)等硬件優(yōu)勢(shì)對(duì)其進(jìn)行加速。文獻(xiàn)[9,10]在單機(jī)多核服務(wù)器上采用OpenMP實(shí)現(xiàn)了并行Kriging插值;文獻(xiàn)[11-13]基于MPI、MKL等技術(shù)結(jié)合Cholesky分解法、Gauss-Jordan消元法等矩陣優(yōu)化算法實(shí)現(xiàn)了支持多達(dá)512核的并行Kriging插值加速;文獻(xiàn)[14]解決了空間分布不均勻的待插值數(shù)據(jù)在多核CPU上的負(fù)載均衡問題。近年來,隨著GPU核數(shù)的不斷攀升以及共享顯存的持續(xù)擴(kuò)大,GPU的并行計(jì)算能力已遠(yuǎn)超主流CPU的,并在科學(xué)計(jì)算領(lǐng)域取得了廣泛應(yīng)用。文獻(xiàn)[15-17]介紹了基于OpenCL、CUDA等技術(shù)在GPU上實(shí)現(xiàn)并行Kriging插值的方法;文獻(xiàn)[18-20]進(jìn)一步研究了基于CPU-GPU協(xié)同的Kriging并行插值的實(shí)現(xiàn)途徑與方法。

以上工作均取得了較高的加速比,顯著提升了Kriging插值效率,但由于對(duì)特殊硬件設(shè)備及編程環(huán)境的依賴,存在使用成本高、開發(fā)編程難等問題,而且普遍缺乏容錯(cuò)機(jī)制,水平擴(kuò)展能力差[21]。

另一種方案是基于目前較成熟的分布式計(jì)算框架(如Spark、Hadoop等)在多臺(tái)機(jī)器上實(shí)現(xiàn)Kriging插值的并行加速。分布式計(jì)算框架采用計(jì)算向數(shù)據(jù)遷移、分布式內(nèi)存、懶計(jì)算等技術(shù)提升海量數(shù)據(jù)的處理效率,同時(shí)具有較強(qiáng)的容錯(cuò)能力,可在大規(guī)模集群上實(shí)現(xiàn)快速可靠的并行計(jì)算。文獻(xiàn)[22,23]在 Hadoop上實(shí)現(xiàn)了分布式Kriging并行插值,與基于CUDA-GPU的方式相比,雖然效率較低,但水平擴(kuò)展能力強(qiáng),且具有更好的可靠性和易用性。文獻(xiàn)[24]介紹了基于Spark實(shí)現(xiàn)分布式Kriging并行插值的技術(shù)方案,并就處理時(shí)間、CPU內(nèi)存利用率等與Hadoop進(jìn)行了對(duì)比,其整體效率達(dá)到了Hadoop的5倍左右。以上方案本質(zhì)上均是基于MapReduce架構(gòu)對(duì)數(shù)據(jù)分塊然后進(jìn)行分布式插值,數(shù)據(jù)分塊一般采用固定空間區(qū)域或滑動(dòng)窗口等方式。固定空間區(qū)域會(huì)使邊界附近的未知點(diǎn)因近鄰樣本點(diǎn)選擇不全而導(dǎo)致插值結(jié)果不準(zhǔn)?；瑒?dòng)窗口方式雖然可以在一定程度上避免這一問題,但會(huì)因不同子數(shù)據(jù)塊空間區(qū)域重疊而產(chǎn)生大量的重復(fù)計(jì)算,影響整體計(jì)算效率。

文獻(xiàn)[25]將所有待插值數(shù)據(jù)看作一個(gè)整體,優(yōu)化了Kriging插值過程中的矩陣求逆乘積運(yùn)算,在Spark上實(shí)現(xiàn)了面向大矩陣的分布式Kriging插值。OCO-3數(shù)據(jù)插值所采用的近鄰點(diǎn)一般不會(huì)超過50個(gè),對(duì)應(yīng)矩陣規(guī)模較小,求逆乘積等速度非?？?沒必要進(jìn)行矩陣運(yùn)算的分布式加速。

此外,OCO-3原始數(shù)據(jù)按照國(guó)家空間站ISS(International Space Station)航線采集,是一種線性時(shí)空數(shù)據(jù),在通過Kriging插值生成最終的CO2濃度分布數(shù)據(jù)時(shí),還需要一些特殊的空間操作,如坐標(biāo)重投影、并行時(shí)空檢索、矢量結(jié)果柵格化等,而Hadoop、Spark等針對(duì)這類空間操作的加速效果并不理想[26]。

本文綜合以上問題,提出另一種方案:將過程復(fù)雜、計(jì)算量大的Kriging插值過程重構(gòu)為可多步執(zhí)行的DAG結(jié)構(gòu)的工作流,并基于具體的數(shù)據(jù)特征對(duì)其中的復(fù)雜計(jì)算環(huán)節(jié)進(jìn)行優(yōu)化;然后將計(jì)算任務(wù)進(jìn)一步拆分,細(xì)化插值粒度,形成可在分布式環(huán)境下并行執(zhí)行的計(jì)算任務(wù),在底層通過一個(gè)雙層架構(gòu)的分布式工作流調(diào)度引擎對(duì)拆分后的插值任務(wù)進(jìn)行并行加速。與其他方案相比,本文方案對(duì)硬件環(huán)境無特殊要求,插值任務(wù)調(diào)度執(zhí)行效率高,且架構(gòu)靈活易實(shí)現(xiàn),可根據(jù)空間區(qū)域尺度和待插值數(shù)據(jù)量動(dòng)態(tài)調(diào)整集群規(guī)模,具有較強(qiáng)的可擴(kuò)展性。

3 數(shù)據(jù)統(tǒng)計(jì)特征分析

OCO-3衛(wèi)星數(shù)據(jù)中包含多項(xiàng)與CO2相關(guān)的指標(biāo),最常用的是XCO2,即二氧化碳柱濃度,單位為PPM(Parts Per Million volume),其數(shù)據(jù)覆蓋空間為[-180,-52,180,52],平均采樣分辨率約為2.5 km*2.5 km,采集軌道如圖1所示。

Figure 1 OCO-3 satellite data collection track圖1 OCO-3衛(wèi)星數(shù)據(jù)采集軌道

基于Kriging算法對(duì)OCO-3數(shù)據(jù)進(jìn)行插值包含2個(gè)步驟:(1)分析OCO-3數(shù)據(jù)的統(tǒng)計(jì)特征并對(duì)其空間相關(guān)關(guān)系進(jìn)行建模;(2)基于構(gòu)建的空間模型進(jìn)行插值。針對(duì)未知點(diǎn)p的插值公式如式(1)所示:

(1)

λ應(yīng)當(dāng)是一組使得p點(diǎn)預(yù)估值與真實(shí)值誤差最小的權(quán)重系數(shù),據(jù)此構(gòu)造目標(biāo)函數(shù)f(λ)如式(2)所示:

(2)

將式(1)代入式(2),結(jié)合普通Kriging插值的應(yīng)用條件,使用拉格朗日乘數(shù)法可得λ的求解公式[27],如式(3)所示。

(3)

其中,rij是基于距離求得的半變異函數(shù)值;樣本點(diǎn)個(gè)數(shù)n由樣本點(diǎn)密度與半變異函數(shù)模型決定,通?？刂圃?0～100內(nèi),太少會(huì)使部分相關(guān)性強(qiáng)的點(diǎn)不能參與插值,太多會(huì)讓過多無關(guān)點(diǎn)參與插值,二者都會(huì)導(dǎo)致插值結(jié)果不準(zhǔn)。

半變異函數(shù)模型基于近鄰樣本點(diǎn)的距離與協(xié)方差通過最小二乘法擬合得到。常見的半變異函數(shù)模型有球形模型、圓形模型、指數(shù)模型、高斯模型及線性模型等。本文隨機(jī)選取100個(gè)樣本點(diǎn),并分別計(jì)算與其相鄰的500個(gè)數(shù)據(jù)點(diǎn)的距離與協(xié)方差,然后以5 km為區(qū)間計(jì)算協(xié)方差的平均值,得到各個(gè)樣本點(diǎn)對(duì)線性模型的擬合度分布,如圖2所示。

Figure 2 Distribution of linear model fitting圖2 線性模型的擬合度分布

由圖2可知,擬合度超過0.5的數(shù)據(jù)點(diǎn)占比接近80%,因此本文選取線性模型作為OCO-3數(shù)據(jù)的插值函數(shù)模型。

4 整體概述

基于以上分析,本文設(shè)計(jì)并實(shí)現(xiàn)如圖3所示的分布式Kriging插值加速框架。該框架以用戶指定的時(shí)空范圍及目標(biāo)分辨率為輸入,輸出結(jié)果為柵格化的XCO2濃度分布影像,插值過程分為2個(gè)階段:Kriging插值工作流構(gòu)建和工作流任務(wù)的分布式并行加速。在工作流構(gòu)建階段,將串行的Kriging插值過程細(xì)化拆分為多個(gè)相對(duì)獨(dú)立的任務(wù)模塊,并根據(jù)其輸入輸出的依賴關(guān)系確定不同任務(wù)的執(zhí)行時(shí)序;然后根據(jù)數(shù)據(jù)特征與計(jì)算模式將每個(gè)獨(dú)立的任務(wù)模塊進(jìn)一步拆分為可并行執(zhí)行的子任務(wù);最后由底層面向DAG工作流的分布式調(diào)度引擎負(fù)責(zé)加速執(zhí)行。

Figure 3 Distributed Kriging interpolation framework for OCO-3 data圖3 面向OCO-3數(shù)據(jù)的分布式Kriging插值框架

底層的分布式調(diào)度引擎基于管理與計(jì)算分離的雙層架構(gòu)實(shí)現(xiàn)。第1層主要負(fù)責(zé)控制整個(gè)工作流任務(wù)按序向前推進(jìn)。第2層負(fù)責(zé)工作流中具體子任務(wù)在分布式環(huán)境中的調(diào)度執(zhí)行。

基于上述結(jié)構(gòu)劃分,下文將首先介紹Kriging插值工作流的構(gòu)建細(xì)節(jié)及其中關(guān)鍵步驟的優(yōu)化方法;然后介紹如何將Kriging插值工作流轉(zhuǎn)換為可在分布式環(huán)境中并行執(zhí)行的任務(wù)以及面向該任務(wù)的分布式調(diào)度引擎的設(shè)計(jì)與實(shí)現(xiàn)。

5 插值工作流的構(gòu)建與優(yōu)化

結(jié)合OCO-3衛(wèi)星數(shù)據(jù)特征構(gòu)建Kriging插值工作流,如圖4所示。

Figure 4 Workflow of Kriging interpolation for large region圖4 面向大區(qū)域尺度的Kriging插值工作流

整個(gè)工作流包含3個(gè)主要模塊:(1) 數(shù)據(jù)解析處理:通過時(shí)空檢索得到指定時(shí)空范圍內(nèi)的原始線性排列的散點(diǎn)數(shù)據(jù),將其網(wǎng)格化,然后計(jì)算并緩存所有數(shù)據(jù)網(wǎng)格之間的距離與協(xié)方差;(2) 網(wǎng)格遍歷插值:逐個(gè)遍歷所有未知網(wǎng)格,基于其近鄰網(wǎng)格之間的距離與協(xié)方差擬合半變異函數(shù)模型,以此對(duì)未知網(wǎng)格進(jìn)行插值;(3) 插值結(jié)果匯總:整合所有數(shù)據(jù)網(wǎng)格與插值網(wǎng)格,結(jié)合空間位置信息,將其柵格化為XCO2濃度分布影像。

在上述工作流中,時(shí)空檢索、近鄰點(diǎn)篩選以及全局距離協(xié)方差計(jì)算是比較耗時(shí)的3個(gè)環(huán)節(jié),為此本文基于網(wǎng)格剖分、數(shù)據(jù)分鏈等方式對(duì)其進(jìn)行了有針對(duì)性的優(yōu)化,使整體插值效率得到顯著提升,具體優(yōu)化細(xì)節(jié)如下文所述。

5.1 分布式并行時(shí)空檢索

直接基于大尺度不規(guī)則的空間區(qū)域進(jìn)行檢索十分耗時(shí),為此,本文在對(duì)原始OCO-3數(shù)據(jù)點(diǎn)構(gòu)建空間索引(RTree)的基礎(chǔ)上,進(jìn)一步將空間檢索范圍由中心點(diǎn)進(jìn)行遞歸網(wǎng)格化拆分,將單個(gè)不規(guī)則大區(qū)域的空間相交查詢轉(zhuǎn)換為多個(gè)矩形網(wǎng)格的并行相交查詢。算法偽代碼如算法1所示。

算法1基于網(wǎng)格的并行時(shí)空檢索算法

輸入：最小網(wǎng)格G(w,h),空間查詢范圍R。

輸出：與R相交的OCO-3數(shù)據(jù)點(diǎn)集合RST。

//計(jì)算R的外包矩形

x0,x1=min(R.coords[:,0]),max(R.coords[:,0]);

y0,y1=min(R.coords[:,1]),max(R.coords[:,1]);

//調(diào)整外包矩形使其為G的整數(shù)倍

a,b,c,d=align(x0,y0,x1,y1);

grids=[[a,b,c,d]];

G1=[];//位于R內(nèi)部的網(wǎng)格

G2=[];//位于R邊界的網(wǎng)格

whilenotgrids.empty():

g=grids.pop();

ifR.contains(g)://g位于R內(nèi)部

G1.append(g);

elseifR.intersects(g)://g位于R邊界

ifg.area()≤Smin:

G2.append(g);

else:

//以G為單位將g近似等分為4個(gè)子網(wǎng)格

subgrids=quad_split(g,G);

forsginsubgrids:

grids.append(sg);

endfor

endif

endif

endif

endwhile

rst1=parallel(query_by_extent,G1);

rst2=parallel(query_by_intersect,G2);

RST=merge(rst1,rst2);

returnRST;

算法1中的query_by_extent通過直接比較數(shù)值大小的方式獲得目標(biāo)數(shù)據(jù),query_by_intersect則先通過數(shù)值比較的方式快速縮小搜索空間,然后通過精確空間相交判斷獲得目標(biāo)數(shù)據(jù)。parallel是將以上2個(gè)函數(shù)分布到不同服務(wù)器節(jié)點(diǎn)并行執(zhí)行的方法。

圖5是基于以上算法對(duì)新疆省邊界進(jìn)行網(wǎng)格化拆分后的效果。

Figure 5 Illustration of grid dividing of spatial range圖5 空間范圍網(wǎng)格拆分示意圖

對(duì)空間區(qū)域網(wǎng)格化后,所有網(wǎng)格都可獨(dú)立地由不同節(jié)點(diǎn)上的計(jì)算進(jìn)程進(jìn)行查詢。大部分?jǐn)?shù)據(jù)可以基于內(nèi)部少量大網(wǎng)格通過query_by_extent快速查詢得到,剩余數(shù)據(jù)則基于邊界小網(wǎng)格通過query_by_intersect查詢得到,小網(wǎng)格雖然數(shù)量多,但網(wǎng)格面積大大縮小,查詢速度也有顯著提升。

采用上述方法分別檢索中國(guó)大陸、新疆、陜西及北京等4個(gè)不同大小區(qū)域的OCO-3數(shù)據(jù),與直接基于單個(gè)大區(qū)域的查詢方式的耗時(shí)對(duì)比如表1所示。

Table 1 Comparison of spatial retrieval time in different regions表1 不同區(qū)域尺度的空間檢索耗時(shí)對(duì)比 s

由表1可知,上述基于網(wǎng)格拆分的并行時(shí)空檢索優(yōu)化方法可以顯著提升檢索效率,且空間范圍越大,速度優(yōu)勢(shì)越明顯。

檢索到原始數(shù)據(jù)后按照目標(biāo)分辨率對(duì)其進(jìn)行網(wǎng)格化處理:取處于同一個(gè)網(wǎng)格的所有樣本點(diǎn)的均值作為該網(wǎng)格的XCO2值,網(wǎng)格位置用中心點(diǎn)表示,同時(shí)對(duì)所有網(wǎng)格按照行列進(jìn)行編號(hào)。網(wǎng)格編號(hào)規(guī)則如式(4)所示:

(4)

其中,gm為行編號(hào);gn為列編號(hào);(px,py)為網(wǎng)格的中心點(diǎn)坐標(biāo);dt為用戶指定的分辨率,即網(wǎng)格邊長(zhǎng)。

5.2 近鄰數(shù)據(jù)點(diǎn)篩選

KDTree是常用的近鄰數(shù)據(jù)點(diǎn)篩選方法,但其在查詢時(shí),每次都要回溯到根節(jié)點(diǎn),且效率會(huì)隨著樣本數(shù)據(jù)量的增加大幅下降。另一種思路是基于空間填充曲線將待索引空間分為多級(jí)網(wǎng)格,并對(duì)每個(gè)網(wǎng)格賦予一個(gè)帶有位置信息的編號(hào),進(jìn)而將二維空間的坐標(biāo)點(diǎn)映射到一維空間,基于該特性可以快速得到處于同一空間區(qū)域的數(shù)據(jù)點(diǎn)。GeoHash、Google S2等都是基于該思路的空間索引算法,空間填充曲線雖然在一維空間里盡可能地維持了空間本地性,但仍有一定的”跳躍”,即部分編號(hào)相鄰的網(wǎng)格實(shí)際相距較遠(yuǎn),直接通過編號(hào)前綴獲得的候選節(jié)點(diǎn)并不全部準(zhǔn)確。

考慮到OCO-3樣本數(shù)據(jù)網(wǎng)格化后分布比較均勻,本文設(shè)計(jì)并實(shí)現(xiàn)了更快的近鄰搜索算法,如圖6所示。

Figure 6 Illustration of neighbor search based on grids圖6 網(wǎng)格化近鄰搜索示意圖

圖6中,gx為待插值網(wǎng)格,網(wǎng)格編號(hào)為(rx,cy);sr為近鄰搜索的最大距離;dt為網(wǎng)格邊長(zhǎng);ceil為向上取整函數(shù);floor為向下取整函數(shù)。在圓形區(qū)域內(nèi)部及邊界的數(shù)據(jù)網(wǎng)格即為所求近鄰數(shù)據(jù)網(wǎng)格,具體篩選步驟如下:

(1)計(jì)算圓形區(qū)域的外包矩形范圍:Ro= [ox0,oy0,ox1,oy1],計(jì)算公式見圖6;

(2)計(jì)算圓形區(qū)域的內(nèi)接正方形范圍:Ri=[ix0,iy0,ix1,iy1],計(jì)算公式見圖6,該范圍內(nèi)的所有數(shù)據(jù)網(wǎng)格均為待插值點(diǎn)的近鄰數(shù)據(jù)網(wǎng)格;

(3)逐個(gè)遍歷介于Ro與Ri之間的數(shù)據(jù)網(wǎng)格,篩選出所有與gx距離小于sr的網(wǎng)格;

(4)匯總步驟(2)與(3)的結(jié)果,即為gx的所有近鄰數(shù)據(jù)網(wǎng)格。

步驟(1)通過外包矩形快速鎖定一個(gè)較小的候選網(wǎng)格范圍。步驟(2)則通過內(nèi)接矩形直接確定了大部分近鄰數(shù)據(jù)網(wǎng)格。以上2個(gè)步驟都是基于簡(jiǎn)單的公式計(jì)算完成,時(shí)間復(fù)雜度為O(1)。介于Ro與Ri之間的數(shù)據(jù)網(wǎng)格需要逐個(gè)遍歷,但數(shù)量較小,不會(huì)對(duì)整體耗時(shí)產(chǎn)生明顯影響。

5.3 全局距離協(xié)方差計(jì)算

Kriging插值的核心步驟是計(jì)算未知點(diǎn)近鄰數(shù)據(jù)網(wǎng)格兩兩之間的距離和協(xié)方差,進(jìn)而以此擬合得到用于插值的半變異函數(shù)模型。由于相鄰未知網(wǎng)格的近鄰數(shù)據(jù)網(wǎng)格重疊度較高,本文采用預(yù)先計(jì)算并緩存的方式來避免重疊區(qū)域的重復(fù)計(jì)算。此外,考慮到插值只用近鄰數(shù)據(jù),為此在執(zhí)行計(jì)算時(shí),只采用一定范圍內(nèi)的數(shù)據(jù)網(wǎng)格,范圍大小由目標(biāo)分辨率與數(shù)據(jù)分布特征確定。

直接計(jì)算所有數(shù)據(jù)網(wǎng)格兩兩之間的距離與協(xié)方差計(jì)算量太大(如圖7a所示),耗時(shí)無法接受,而且由于交叉太多,難以并行加速。分塊計(jì)算再合并的方式(如圖7b所示)可以實(shí)現(xiàn)并行加速,但合并過程比較復(fù)雜,且難以應(yīng)用距離限制計(jì)算量。針對(duì)以上問題,本文采用如圖8所示的分行鏈?zhǔn)浇Y(jié)構(gòu)對(duì)數(shù)據(jù)進(jìn)行重新組織,以實(shí)現(xiàn)分布式環(huán)境下的全局距離與協(xié)方差的并行計(jì)算。

Figure 7 Two modes of global distance and covariance computation圖7 2種全局距離協(xié)方差計(jì)算模式

Figure 8 Linked structure of data samples圖8 樣本數(shù)據(jù)點(diǎn)的鏈?zhǔn)浇M織結(jié)構(gòu)

圖8中,xn表示數(shù)據(jù)網(wǎng)格的橫坐標(biāo),y(n,i)表示數(shù)據(jù)網(wǎng)格的縱坐標(biāo),v(n,i)為XCO2值。距離基于位置(xn,y(n,i))計(jì)算,協(xié)方差基于v(n,i)計(jì)算,采用當(dāng)前行向右,跨行由中間向兩側(cè)的遍歷策略進(jìn)行計(jì)算(如圖9所示)。算法偽代碼如算法2所示。

Figure 9 Illustration of distance and covariance computation between neighbor data points圖9 相鄰數(shù)據(jù)點(diǎn)距離協(xié)方差計(jì)算示意圖

算法2相鄰數(shù)據(jù)網(wǎng)格距離協(xié)方差計(jì)算

輸入：行列號(hào)為(M,N)的數(shù)據(jù)網(wǎng)格G。

輸出：相鄰數(shù)據(jù)網(wǎng)格之間的距離與協(xié)方差字典DVS。

forginnext(M,N)://行內(nèi)向右遍歷

//距離、協(xié)方差計(jì)算

d,v=dist(g,G),cov(g,G);

ifd>MAXDIST:

break;

endif

DVS[(G,g)]=(d,v);

endfor

m,n=G.M+1,G.N;//開始逐行向下遍歷

whilem

forginnext(m,n)://行內(nèi)向右遍歷

d,v=dist(g,G),cov(g,G);

ifd>MAXDIST:

break;

endif

DVS[(G,g)]=(d,v);

endfor

forginpre(m,n)://行內(nèi)向左遍歷

d,v=dist(g,G),cov(g,G);

ifd>MAXDIST:

break;

endif

DVS[(G,g)]=(d,v);

endfor

m+=1;

endwhile

其中,MAXDIST為限定相鄰網(wǎng)格的距離常量,根據(jù)待插值區(qū)域的樣本點(diǎn)分布設(shè)置;DVS的主鍵為(m0,n0,m1,n1),m,n為數(shù)據(jù)網(wǎng)格編號(hào),且滿足:(m0

6 插值任務(wù)的分布式并行加速

本節(jié)介紹如何將上述Kriging插值工作流以及優(yōu)化算法轉(zhuǎn)換為可并行執(zhí)行的任務(wù)模塊,并基于不同任務(wù)模塊之間輸入輸出的依賴關(guān)系確定調(diào)度執(zhí)行的順序,最后通過一套雙層架構(gòu)的分布式調(diào)度引擎實(shí)現(xiàn)上述工作流任務(wù)的并行加速。

6.1 工作流任務(wù)模塊劃分

整個(gè)工作流分為5個(gè)主任務(wù):并行時(shí)空檢索(T1)、全局距離/協(xié)方差計(jì)算(T2)、近鄰點(diǎn)篩選(T3)、半變異函數(shù)擬合并插值(T4)和插值結(jié)果柵格化(T5),如圖10所示。

Figure 10 Kriging interpolation workflow for distributed environment圖10 面向分布式環(huán)境的Kriging插值工作流

主任務(wù)T1基于時(shí)間范圍區(qū)間劃分及空間范圍網(wǎng)格化拆分為m個(gè)可分布式并行執(zhí)行的子任務(wù),T1的輸出結(jié)果為網(wǎng)格化的數(shù)據(jù)點(diǎn)。主任務(wù)T2以T1的輸出作為輸入,基于分行鏈?zhǔn)酱鎯?chǔ)結(jié)構(gòu)拆分為n個(gè)可分布式并行執(zhí)行子任務(wù),T2的結(jié)果以鍵值對(duì)的形式緩存至Redis內(nèi)存數(shù)據(jù)庫(kù)。主任務(wù)T3并行遍歷待插值網(wǎng)格,基于5.2節(jié)描述的算法獲取每個(gè)待插值網(wǎng)格的近鄰數(shù)據(jù)網(wǎng)格集合。主任務(wù)T4以該近鄰數(shù)據(jù)網(wǎng)格集合以及Redis緩存為輸入完成最終的半邊異函數(shù)擬合與插值,T4的輸出結(jié)果最終通過主任務(wù)T5柵格化為GeoTIFF文件進(jìn)行存儲(chǔ)。主任務(wù)的執(zhí)行時(shí)序如圖11所示。

Figure 11 Executing sequence of main tasks圖11 主任務(wù)執(zhí)行時(shí)序

圖11中,主任務(wù)T2、T3為并行執(zhí)行,其他主任務(wù)均為串行執(zhí)行,執(zhí)行到具體某個(gè)主任務(wù)時(shí),主任務(wù)內(nèi)部的子任務(wù)可通過獨(dú)立調(diào)度在分布式環(huán)境下并行執(zhí)行。

6.2 分布式調(diào)度執(zhí)行與時(shí)序控制

基于6.1節(jié)所述的工作流任務(wù)特點(diǎn),本文設(shè)計(jì)并實(shí)現(xiàn)如圖12所示的雙層任務(wù)調(diào)度框架。第1層為唯一的主調(diào)度器,負(fù)責(zé)主任務(wù)的調(diào)度管理與時(shí)序控制;第2層包含多個(gè)分調(diào)度器,每個(gè)分調(diào)度器管理一個(gè)或多個(gè)計(jì)算節(jié)點(diǎn),并具體負(fù)責(zé)子任務(wù)的調(diào)度執(zhí)行,分調(diào)度器由主調(diào)度器管理,接受并執(zhí)行主調(diào)度器發(fā)送過來的調(diào)度命令并向其匯報(bào)自己的狀態(tài)。

Figure 12 Framework of double-tier scheduling engine for Kriging interpolation workflow圖12 面向Kriging插值工作流的雙層任務(wù)調(diào)度框架

所有主任務(wù)按照依賴關(guān)系與執(zhí)行順序生成一個(gè)DAG結(jié)構(gòu)的任務(wù)流,其任務(wù)節(jié)點(diǎn)以任務(wù)記錄的形式存入主任務(wù)隊(duì)列。任務(wù)記錄中包含任務(wù)標(biāo)識(shí)(id)、前序任務(wù)(pre)、后續(xù)任務(wù)(next)、狀態(tài)(status)、輸入(input)、輸出(output)、執(zhí)行時(shí)間(timestamp)等信息。主調(diào)度器遍歷該任務(wù)隊(duì)列,取出所有滿足條件(無前序依賴,且輸入條件已滿足)的主任務(wù)交給分調(diào)度器執(zhí)行,分調(diào)度器將主任務(wù)對(duì)應(yīng)的子任務(wù)通過分布式任務(wù)調(diào)度引擎將其分發(fā)到計(jì)算集群上并行執(zhí)行,并實(shí)時(shí)監(jiān)控其運(yùn)行狀態(tài)。當(dāng)前主任務(wù)執(zhí)行結(jié)束后,分調(diào)度器向主調(diào)度器匯報(bào),主調(diào)度器據(jù)此更新主任務(wù)隊(duì)列中相關(guān)的任務(wù)狀態(tài)、前序依賴、輸入輸出等信息,并繼續(xù)調(diào)度執(zhí)行符合條件的主任務(wù)。

為了實(shí)現(xiàn)任務(wù)狀態(tài)快速更新,降低任務(wù)調(diào)度開銷,主調(diào)度器與分調(diào)度器均采用C++開發(fā)實(shí)現(xiàn),并對(duì)圖10所示的工作流進(jìn)行拓?fù)渑判蛐纬扇蝿?wù)隊(duì)列,然后將其以鍵值對(duì)的形式存儲(chǔ)至有序字典中。其中key為任務(wù)標(biāo)識(shí),value為包含執(zhí)行該任務(wù)的函數(shù)名稱、參數(shù)以及前后任務(wù)等信息。

位于底層計(jì)算節(jié)點(diǎn)的任務(wù)執(zhí)行模塊基于Dask實(shí)現(xiàn)。Dask打破了Python中GIL(Global Interpreter Lock)的限制,支持多核并行計(jì)算與核外計(jì)算,且具有不亞于OpenBLAS的矩陣運(yùn)算效率[28]。

考慮到OCO-3樣本點(diǎn)總數(shù)量在10億以內(nèi),內(nèi)存占用不超過50 GB,為了提升數(shù)據(jù)訪問效率,不同計(jì)算節(jié)點(diǎn)之間的數(shù)據(jù)共享通過內(nèi)存數(shù)據(jù)庫(kù)Redis實(shí)現(xiàn)。

分調(diào)度器在執(zhí)行任務(wù)分配時(shí),優(yōu)先將同一主任務(wù)的子任務(wù)分配到相同計(jì)算節(jié)點(diǎn)上,以盡可能降低任務(wù)調(diào)度管理以及跨節(jié)點(diǎn)數(shù)據(jù)傳輸?shù)葞淼念~外開銷。計(jì)算節(jié)點(diǎn)根據(jù)自身可用的計(jì)算資源設(shè)置任務(wù)并發(fā)數(shù)上限,當(dāng)某個(gè)節(jié)點(diǎn)負(fù)載接近上限時(shí),再由分調(diào)度器將待執(zhí)行任務(wù)分配到其他節(jié)點(diǎn)。

在上述框架中,主調(diào)度器、分調(diào)度器與任務(wù)執(zhí)行等各模塊相對(duì)獨(dú)立,各個(gè)模塊之間的任務(wù)調(diào)度、運(yùn)行狀態(tài)、啟停命令等消息傳遞則基于ProtoBuf序列化協(xié)議實(shí)現(xiàn),該協(xié)議與XML、JSON相比,序列化速度更快、數(shù)據(jù)體積更小[29]。

7 實(shí)驗(yàn)與結(jié)果分析

根據(jù)處理級(jí)別與數(shù)據(jù)內(nèi)容的不同,OCO-3衛(wèi)星數(shù)據(jù)產(chǎn)品分為18種類型。其中,OCO3_L2_Standard V10r為官方推薦數(shù)據(jù)集,自2019/08/05持續(xù)更新至今。截至2021/12/31,已發(fā)布了11 380個(gè)數(shù)據(jù)文件,共包含109 315 194數(shù)據(jù)點(diǎn),發(fā)布的數(shù)據(jù)產(chǎn)品采用NetCDF-4/HDF5格式存儲(chǔ)。為方便后續(xù)分析,本文將原始數(shù)據(jù)中的時(shí)間、位置、XCO2等信息以四元組(timestamp,longitude,latitude,xco2)形式進(jìn)行重新組織,重組后的數(shù)據(jù)存入關(guān)系型數(shù)據(jù)庫(kù)PostgreSQL中。

為測(cè)試和評(píng)估本文所提出的優(yōu)化算法與計(jì)算框架的性能效率,選取了2021年的OCO-3數(shù)據(jù)進(jìn)行實(shí)驗(yàn)分析,測(cè)試環(huán)境為6臺(tái)相同配置(如表2所示)的服務(wù)器搭建的集群,服務(wù)器之間的網(wǎng)絡(luò)帶寬為1 Gbps。其中,5臺(tái)服務(wù)器用于搭建計(jì)算集群,每臺(tái)服務(wù)器的計(jì)算任務(wù)并發(fā)上限數(shù)為16,另外1臺(tái)用作Redis數(shù)據(jù)庫(kù)服務(wù)器,所有服務(wù)器掛載均基于MooseFS搭建的分布式文件系統(tǒng)。

為測(cè)試不同區(qū)域尺度下的插值效率,本文選取了中國(guó)大陸、新疆、陜西、北京等4個(gè)邊界不規(guī)則且大小不同的測(cè)試區(qū)域,如圖13所示。

Figure 13 Interpolation regions of different sizes圖13 不同大小的插值區(qū)域

基于5.1節(jié)介紹的并行時(shí)空檢索方法得到以上4個(gè)區(qū)域的原始數(shù)據(jù)后,按0.1度分辨率對(duì)原始樣本點(diǎn)數(shù)據(jù)進(jìn)行網(wǎng)格化,并取平均值作為該網(wǎng)格的XCO2數(shù)值。以上4個(gè)區(qū)域包含的原始樣本數(shù)據(jù)點(diǎn)個(gè)數(shù)、數(shù)據(jù)網(wǎng)格數(shù)、未知網(wǎng)格數(shù)如表3所示。

Table 3 Size of test data表3 測(cè)試數(shù)據(jù)量

由表3可知,OCO-3數(shù)據(jù)在不同區(qū)域的分布并不均勻,新疆的數(shù)據(jù)覆蓋比例比陜西的高20%左右,以上4個(gè)區(qū)域的數(shù)據(jù)分布情況如圖14所示。

Figure 14 Data coverage of 4 different regions圖14 4個(gè)不同區(qū)域的數(shù)據(jù)覆蓋情況

Figure 15 Images rendered from interpolation results of four regions圖15 4個(gè)區(qū)域的插值結(jié)果渲染圖

基于上述環(huán)境,采用本文所述方法框架得到的插值結(jié)果如圖15所示。完成以上4個(gè)區(qū)域插值的整體耗時(shí)如圖16所示。

Figure 16 Comparison of time consumed by interpolation in four different regions圖16 4個(gè)不同區(qū)域的插值耗時(shí)對(duì)比

由圖16可知,雖然陜西與北京的面積、待插值數(shù)據(jù)量相差20倍左右,但耗時(shí)僅相差3倍左右,數(shù)據(jù)量的增加并未導(dǎo)致耗時(shí)等比例增加。主要是因?yàn)楸本﹨^(qū)域的數(shù)據(jù)規(guī)模較小,插值時(shí)不能完全利用已有的計(jì)算資源,當(dāng)有少量的數(shù)據(jù)增加時(shí),更多未被利用的資源開始參與計(jì)算,因此陜西相對(duì)北京的插值耗時(shí)增加不明顯。但是,當(dāng)網(wǎng)格數(shù)據(jù)量增加到15 638且待插值網(wǎng)格達(dá)到2 368時(shí),單臺(tái)服務(wù)器的計(jì)算資源已難以支撐,更多的計(jì)算節(jié)點(diǎn)開始參與插值,耗時(shí)顯著增加,但增加幅度仍遠(yuǎn)小于數(shù)據(jù)量的增長(zhǎng)幅度。當(dāng)區(qū)域尺度由新疆?dāng)U大到全國(guó)時(shí),網(wǎng)格數(shù)據(jù)量增加了4倍左右,待插值網(wǎng)格增加了10倍左右,但耗時(shí)僅增加了5倍左右,說明本文所設(shè)計(jì)的方法框架針對(duì)大區(qū)域尺度的插值具有良好的并行加速效果。

以上4個(gè)地區(qū)在不同任務(wù)階段的耗時(shí)對(duì)比如圖17所示。4個(gè)區(qū)域最終的柵格化(T5)都沒有超過0.1 s,與其他階段耗時(shí)差異太大,故未在圖17上展示該階段的耗時(shí)。由圖17可知,4個(gè)階段耗時(shí)最多的是T2階段,該階段需要計(jì)算相鄰數(shù)據(jù)點(diǎn)兩兩之間的距離和協(xié)方差并緩存至Redis,時(shí)間復(fù)雜度約為O(n2),計(jì)算量大且有較多的I/O開銷;T3與T4階段主要操作對(duì)象為待插值網(wǎng)格,其數(shù)量少于數(shù)據(jù)網(wǎng)格,且其中近鄰篩選與插值可以通過簡(jiǎn)單計(jì)算或緩存查詢直接完成,時(shí)間復(fù)雜度約為O(n),因此耗時(shí)均小于相同區(qū)域的T2階段的。

Figure 17 Comparison of time consumed at different stages in four test regions圖17 4個(gè)測(cè)試區(qū)域在不同階段的耗時(shí)對(duì)比

本文選取Spark作為對(duì)比框架。Spark是一種基于分布式內(nèi)存的大數(shù)據(jù)分析處理框架,支持交互式的數(shù)據(jù)查詢與大批量任務(wù)的并行執(zhí)行,GeoSpark對(duì)Spark進(jìn)行了擴(kuò)展,使其具有了大規(guī)模時(shí)空數(shù)據(jù)處理能力[30]。圖18展示了GeoSpark與本文所用網(wǎng)格化方法在時(shí)空數(shù)據(jù)檢索階段(T1)的效率對(duì)比。由圖18可知,本文所采用的方法在4個(gè)區(qū)域上的檢索效率均高于GeoSpark的,而且區(qū)域越大,本文所述方法的效率優(yōu)勢(shì)越明顯。

Figure 18 Comparison of time consumed between GeoSpark and the proposed method圖18 GeoSpark與本文并行時(shí)空檢索方法耗時(shí)對(duì)比

考慮到基于矩陣的批量插值模式需要大內(nèi)存的支持且不適合不規(guī)則的大區(qū)域,所以基于Spark的插值依然采用先計(jì)算距離/協(xié)方差并緩存然后逐點(diǎn)插值的流程。全局距離/協(xié)方差計(jì)算(T2)、近鄰點(diǎn)篩選(T3)和半變異函數(shù)擬合并插值(T4)3個(gè)階段在Spark上的實(shí)現(xiàn)方式如圖19所示。

Figure 19 Flow of Kriging interpolation based on Spark圖19 基于Spark的Kriging插值流程

Spark主要通過數(shù)據(jù)分區(qū)后獨(dú)立并行計(jì)算來提升計(jì)算效率,而距離/協(xié)方差計(jì)算需要大量交叉運(yùn)算,難以直接分成獨(dú)立的數(shù)據(jù)塊。為此在計(jì)算之前,本文將所有數(shù)據(jù)網(wǎng)格通過兩兩組合的形式形成圖19b所示的四元組;然后再由Spark進(jìn)行并行計(jì)算,后續(xù)的近鄰查詢、半變異函數(shù)擬合與插值仍采用本文所述的方法,具體的子任務(wù)由Spark負(fù)責(zé)調(diào)度執(zhí)行?；赟park完成以上操作的耗時(shí)如圖20所示。

Figure 20 Comparison of time consumed between Spark and the proposed method圖20 本文方法與Spark插值耗時(shí)對(duì)比

由圖20可知,在大區(qū)域尺度上,本文方法相對(duì)Spark具有明顯的效率優(yōu)勢(shì),原因主要有2個(gè):(1)在計(jì)算距離協(xié)方差時(shí),本文采用的分行鏈?zhǔn)酱鎯?chǔ)結(jié)構(gòu)可以有效地結(jié)合距離范圍限制參與計(jì)算的數(shù)據(jù)量,而Spark的分塊計(jì)算模式難以利用該條件;(2)在后續(xù)Kriging插值階段,Spark自身的內(nèi)存計(jì)算、彈性數(shù)據(jù)集等技術(shù)難以對(duì)逐點(diǎn)插值的計(jì)算流程進(jìn)行加速,而其內(nèi)部復(fù)雜的任務(wù)調(diào)度、監(jiān)控等操作增加了任務(wù)之外的額外開銷。

8 結(jié)束語

碳衛(wèi)星數(shù)據(jù)是進(jìn)行碳排放和碳追蹤研究的重要數(shù)據(jù)源,但在大區(qū)域尺度上基于Kriging算法進(jìn)行插值時(shí)計(jì)算量很大,而且難以直接進(jìn)行并行加速。為此,本文以常用OCO-3衛(wèi)星數(shù)據(jù)為例,設(shè)計(jì)了一種DAG結(jié)構(gòu)的Kriging插值工作流,并實(shí)現(xiàn)了可對(duì)其并行加速的雙層分布式調(diào)度引擎,通過調(diào)整插值過程、細(xì)化插值粒度以及優(yōu)化其中的近鄰查詢、距離/協(xié)方差計(jì)算等關(guān)鍵環(huán)節(jié)使該工作流可以在分布式環(huán)境下快速執(zhí)行。實(shí)驗(yàn)結(jié)果表明,本文方法框架可以高效完成不同區(qū)域尺度上的OCO-3快速插值,在大區(qū)域尺度上,相對(duì)Spark也具有明顯的速度優(yōu)勢(shì)。后續(xù)將進(jìn)一步探索更多碳衛(wèi)星觀測(cè)數(shù)據(jù),研究不同碳衛(wèi)星數(shù)據(jù)的融合方法,形成一套更加準(zhǔn)確高效的碳衛(wèi)星觀測(cè)數(shù)據(jù)插值融合框架。