鄭茂騰，周順平，熊小東，朱俊峰

1. 中國(guó)地質(zhì)大學(xué)(武漢)國(guó)家地理信息系統(tǒng)工程技術(shù)研究中心，湖北 武漢 430074； 2. 北京中測(cè)智繪有限責(zé)任公司，北京 100830

鄭茂騰，周順平，熊小東，等.GPU并行區(qū)域網(wǎng)平差[J].測(cè)繪學(xué)報(bào)，2017,46(9)：1193-1201.

10.11947/j.AGCS.2017.20160636.

ZHENG Maoteng,ZHOU Shunping,XIONG Xiaodong,et al.GPU Parallel Bundle Block Adjustment[J]. Acta Geodaetica et Cartographica Sinica,2017,46(9):1193-1201. DOI:10.11947/j.AGCS.2017.20160636.

GPU并行區(qū)域網(wǎng)平差

鄭茂騰1，周順平1，熊小東2，朱俊峰2

1. 中國(guó)地質(zhì)大學(xué)(武漢)國(guó)家地理信息系統(tǒng)工程技術(shù)研究中心，湖北 武漢 430074； 2. 北京中測(cè)智繪有限責(zé)任公司，北京 100830

為了進(jìn)一步解決大數(shù)據(jù)量帶來的平差效率低下的問題，引入GPU并行計(jì)算技術(shù)，同時(shí)使用預(yù)條件共軛梯度法以及不精確牛頓解法求解區(qū)域網(wǎng)平差過程中的法方程，構(gòu)建了適用于GPU并行計(jì)算的全新的區(qū)域網(wǎng)平差技術(shù)流程。本文方法避免了存儲(chǔ)法方程系數(shù)矩陣，而是在需要的時(shí)候?qū)崟r(shí)的計(jì)算該矩陣，使得本文算法相較于傳統(tǒng)的算法所需的計(jì)算機(jī)內(nèi)存空間大幅減少(僅需要存儲(chǔ)平差原始數(shù)據(jù)即可)，平差計(jì)算速度明顯提升，同時(shí)計(jì)算精度與傳統(tǒng)方法相當(dāng)。初步試驗(yàn)證明，本文的方法在普通電腦上僅需要約1.5 min即可完成對(duì)4500張影像、近900萬像點(diǎn)數(shù)據(jù)的平差計(jì)算，且計(jì)算精度達(dá)到子像素級(jí)。

GPU并行計(jì)算；區(qū)域網(wǎng)平差；預(yù)條件共軛梯度；不精確牛頓解；大數(shù)據(jù)

隨著傾斜攝影、無人機(jī)攝影、車載相機(jī)攝影、普通手持相機(jī)攝影等新的數(shù)據(jù)獲取方法的引入，航空攝影測(cè)量的數(shù)據(jù)源日趨多源化，同時(shí)，地理空間信息也逐步進(jìn)入了大數(shù)據(jù)時(shí)代[1]，傳統(tǒng)的航空攝影測(cè)量技術(shù)手段已經(jīng)難以處理多源化、大規(guī)模的影像數(shù)據(jù)，學(xué)者們針對(duì)這些新的數(shù)據(jù)也做了大量的研究[2-5]，但區(qū)域網(wǎng)平差的效率問題始終沒有得到很好的解決。有的學(xué)者對(duì)法方程矩陣的帶寬進(jìn)行了優(yōu)化[6-9]，還有學(xué)者則直接提出了一種塊狀稀疏矩陣壓縮格式，并采用預(yù)條件共軛梯度法(preconditioned conjugate gradient，PCG)以及不精確牛頓解法來求解法方程，大幅減少了平差所需的內(nèi)存，提高了平差效率，同時(shí)保證平差精度與傳統(tǒng)方法相當(dāng)[10-11]，但也同時(shí)指出，如果影像數(shù)持續(xù)增加(超過2萬張影像)，該矩陣壓縮方法仍然需要耗費(fèi)大量?jī)?nèi)存，將可能導(dǎo)致內(nèi)存無法開辟而使得平差無法進(jìn)行。本文是文獻(xiàn)[10]的后續(xù)研究，為了解決影像數(shù)據(jù)持續(xù)增加帶來的問題，本文將不存儲(chǔ)法方程，而是引入圖形處理器(graphics processing unit,GPU)并行計(jì)算，實(shí)時(shí)的計(jì)算法方程，法方程的求解方法仍然采用PCG及不精確牛頓解法，該方法在提高平差效率的同時(shí)，還使得法方程求解過程易于并行化實(shí)現(xiàn)，為GPU并行計(jì)算的引入提供了必要條件。

近年來，隨著GPU通用計(jì)算的發(fā)展，GPU的應(yīng)用已經(jīng)不僅僅局限于圖形顯示領(lǐng)域，很多CPU的計(jì)算任務(wù)也可以由GPU并行計(jì)算來完成，而且GPU提供了數(shù)十倍乃至于上百倍于CPU的計(jì)算性能，使得GPU在非圖形顯示的高性能計(jì)算領(lǐng)域得到了廣泛的應(yīng)用。GPU的高性能計(jì)算能力主要得益于其遠(yuǎn)多于CPU的計(jì)算核心數(shù)，但是其數(shù)據(jù)邏輯運(yùn)算單元較少，因此GPU更加擅長(zhǎng)于大量簡(jiǎn)單的運(yùn)算工作，而CPU則擅長(zhǎng)于處理復(fù)雜的計(jì)算任務(wù)(如大量的邏輯判斷，復(fù)雜的數(shù)據(jù)依賴等情況)。

本文研究的區(qū)域網(wǎng)平差問題中，法方程的構(gòu)建以及求解過程就需要高性能計(jì)算能力的支撐，特別是隨著大數(shù)據(jù)時(shí)代的到來，區(qū)域網(wǎng)平差需要求解的法方程數(shù)據(jù)量更大，計(jì)算能力需求也大幅增加。盡管有很多有效的矩陣分解、求逆等大型矩陣的計(jì)算工具，但是矩陣的存儲(chǔ)問題、區(qū)域網(wǎng)平差流程的并行化問題仍然沒有得到解決。使用GPU并行計(jì)算是解決大數(shù)據(jù)量問題的重要手段。在此之前，也有很多學(xué)者在GPU并行區(qū)域網(wǎng)平差方面做了很多研究工作[12-28]。在算法方面，有些方法應(yīng)用了共軛梯度法(CG)，但是沒有使用預(yù)條件矩陣，也沒有使用不精確牛頓解來降低求解法方程時(shí)的迭代次數(shù)[12-14]。有些方法仍然存儲(chǔ)了法方程，處理大數(shù)據(jù)時(shí)，仍需要大量的內(nèi)存[13-17]。在技術(shù)流程方面，傳統(tǒng)區(qū)域網(wǎng)平差流程需要進(jìn)行全面的修改，以便于引進(jìn)GPU并行計(jì)算，而且在GPU內(nèi)存使用、任務(wù)分配及GPU代碼設(shè)計(jì)等方面都需要進(jìn)行具體的優(yōu)化，才能達(dá)到最優(yōu)的并行效率，但是上述方法中，文獻(xiàn)[17]和文獻(xiàn)[29]最大化地利用了GPU高性能計(jì)算的特點(diǎn)，沒有存儲(chǔ)法方程系數(shù)矩陣，而是實(shí)時(shí)地計(jì)算法方程，但是二者在具體的GPU內(nèi)存使用和任務(wù)分配等方面的描述較少，且其處理精度也不高。本文所做的工作就是結(jié)合預(yù)條件共軛梯度法及不精確牛頓解法，引入GPU并行計(jì)算，通過對(duì)傳統(tǒng)的區(qū)域網(wǎng)平差流程進(jìn)行全面改進(jìn)，使其能夠進(jìn)行并行化設(shè)計(jì)，為GPU并行計(jì)算的引入提供必要條件。并行任務(wù)分配必須滿足負(fù)載均衡的原則，才能使得并行效率得到最大化；同時(shí)GPU內(nèi)存的協(xié)調(diào)調(diào)度對(duì)GPU并行計(jì)算也是至關(guān)重要的，有大量試驗(yàn)表明，如果內(nèi)存處理不得當(dāng)，其運(yùn)算效率會(huì)受到嚴(yán)重影響，甚至?xí)陀谄胀℅PU并行算法。因此，為了使得平差效率最大化，本文重點(diǎn)闡述了GPU并行計(jì)算的任務(wù)分配及內(nèi)存調(diào)度方法。

綜上所述，本文引入GPU并行計(jì)算，使用預(yù)條件共軛梯度法及其不精確牛頓解法求解區(qū)域網(wǎng)平差流程中，建立適用于GPU并行計(jì)算的全新的區(qū)域網(wǎng)平差技術(shù)流程，同時(shí)對(duì)GPU并行計(jì)算中的內(nèi)存使用，任務(wù)分配等具體環(huán)節(jié)進(jìn)行優(yōu)化，以最大化GPU并行計(jì)算的效率。

1 理論與方法

1.1 基于預(yù)條件共軛梯度法的區(qū)域網(wǎng)平差

區(qū)域網(wǎng)平差的數(shù)學(xué)模型在文獻(xiàn)[10]中已有詳細(xì)的描述，本文僅給出簡(jiǎn)單的介紹。傳統(tǒng)的法方程表達(dá)形式如式(1)所示

(1)

式中，U對(duì)應(yīng)于法方程中外方位元素未知數(shù)部分；V對(duì)應(yīng)于法方程中地面點(diǎn)未知數(shù)部分；W表示二者的協(xié)方差矩陣；Δxc表示外方位元素未知數(shù)向量；lc表示法方程常數(shù)項(xiàng)中外方位元素部分；lp表示法方程常數(shù)項(xiàng)中地面點(diǎn)未知數(shù)部分。

此時(shí)，通過解算式(1)可以得到外方位元素未知數(shù)改正數(shù)，地面點(diǎn)未知數(shù)則通過利用新的外方位元素進(jìn)行前方交會(huì)得到。為了避免對(duì)法方程直接求逆，本文引入了預(yù)條件共軛梯度法對(duì)法方程式(1)進(jìn)行求解。同時(shí)為了進(jìn)一步降低迭代次數(shù)，還引入了不精確牛頓解法、預(yù)條件共軛梯度法和不精確牛頓解法在文獻(xiàn)[10]中有詳細(xì)的介紹，本文不再贅述。

1.2 GPU并行區(qū)域網(wǎng)平差技術(shù)流程

引入GPU并行計(jì)算、預(yù)條件共軛梯度法及不精確牛頓解法后，全新的區(qū)域網(wǎng)平差技術(shù)流程如圖1所示。與文獻(xiàn)[10]中的流程類似，不同點(diǎn)是本文的流程中有兩個(gè)步驟是在GPU中運(yùn)行的(圖1中，被填充的步驟)，其計(jì)算結(jié)果被拷貝回CPU內(nèi)存，然后繼續(xù)下一步流程。

圖1 GPU并行區(qū)域網(wǎng)平差技術(shù)流程Fig.1 The whole workflow of the GPU parallel bundle adjustment

如圖1所示，新的區(qū)域網(wǎng)平差流程可以分為兩個(gè)迭代循環(huán)，一個(gè)是區(qū)域網(wǎng)平差迭代，一個(gè)是預(yù)條件共軛梯度法迭代，后者包含于前者，每一次區(qū)域網(wǎng)平差迭代，都包含一個(gè)預(yù)條件共軛梯度法迭代循環(huán)，該迭代循環(huán)是用來求解法方程的。只有當(dāng)區(qū)域網(wǎng)平差迭代收斂后，整個(gè)區(qū)域網(wǎng)平差流程結(jié)束。

1.3 GPU并行區(qū)域網(wǎng)平差優(yōu)化

GPU并行計(jì)算的效率容易受到GPU并行任務(wù)分配以及GPU內(nèi)存的使用方案的影響。各線程的計(jì)算任務(wù)應(yīng)當(dāng)滿足均衡性才能最大地發(fā)揮并行計(jì)算的高效性。另外，GPU全局內(nèi)存容量大，一般可達(dá)到2～4 GB，且所有的線程均可以訪問，但是其存取效率與存取的順序有著密切的聯(lián)系，如果存取的順序正好滿足全局內(nèi)存讀取的一致性，則存取效率高，否則其存取效率變得十分低下。共享內(nèi)存的存取效率高，不受存取順序的影響，但是同一個(gè)塊(block)內(nèi)的線程才能訪問該塊的共享內(nèi)存，不同的塊之間的線程無法訪問對(duì)方的共享內(nèi)存。而且共享內(nèi)存的存儲(chǔ)容量非常小(本文所使用GPU的共享內(nèi)存僅有48 KB)。因此合理地使用各種內(nèi)存，會(huì)直接影響到GPU并行計(jì)算的效率。綜上所述，GPU并行計(jì)算任務(wù)的優(yōu)化，主要包括GPU并行任務(wù)分配優(yōu)化以及內(nèi)存的使用優(yōu)化。本文中的區(qū)域網(wǎng)平差技術(shù)流程有兩個(gè)步驟使用了GPU并行計(jì)算，下文將對(duì)這兩個(gè)步驟的并行任務(wù)分配和內(nèi)存使用策略進(jìn)行詳細(xì)討論。

1.3.1 GPU并行區(qū)域網(wǎng)平差任務(wù)分配

針對(duì)圖1中使用GPU并行計(jì)算的兩個(gè)步驟:①計(jì)算常數(shù)項(xiàng)以及預(yù)條件矩陣；②計(jì)算矩陣向量積。其任務(wù)分配情況各有不同。

1.3.1.1 常數(shù)項(xiàng)向量以及預(yù)條件矩陣的并行計(jì)算任務(wù)分配

常數(shù)項(xiàng)以及預(yù)條件矩陣的計(jì)算方法在文獻(xiàn)[10]中有明確的描述，利用每一個(gè)像點(diǎn)及其對(duì)應(yīng)的物方點(diǎn)坐標(biāo)，外方位元素可以計(jì)算得到常數(shù)項(xiàng)向量的一個(gè)分量，將所有像點(diǎn)計(jì)算得到的常數(shù)項(xiàng)向量的分量累加起來，即可得到完整的常數(shù)項(xiàng)向量。預(yù)條件矩陣是一個(gè)塊狀對(duì)角矩陣，由對(duì)角線上的塊狀子矩陣組成，每一個(gè)像點(diǎn)在計(jì)算常數(shù)項(xiàng)向量時(shí)，同時(shí)也可以計(jì)算得到其對(duì)應(yīng)的預(yù)條件矩陣分量，將所有像點(diǎn)計(jì)算得到的預(yù)條件矩陣的分量累加起來，即可得到完整的預(yù)條件矩陣。因此常數(shù)項(xiàng)向量以及預(yù)條件矩陣的計(jì)算可以放在同一個(gè)并行任務(wù)中進(jìn)行。其任務(wù)分配方法有兩種，以物方點(diǎn)為單位進(jìn)行任務(wù)分配和以像方點(diǎn)為單位進(jìn)行任務(wù)分配。在平差原始數(shù)據(jù)中，每一個(gè)連接點(diǎn)是由一個(gè)物方點(diǎn)，以及該物方點(diǎn)在若干影像上成像得到的像點(diǎn)坐標(biāo)組成，一個(gè)物方點(diǎn)對(duì)應(yīng)至少2個(gè)像方點(diǎn)，不同的物方點(diǎn)對(duì)應(yīng)的像點(diǎn)數(shù)存在較大差別，最少2個(gè)，最大可達(dá)到數(shù)百個(gè)(與測(cè)區(qū)影像的重疊度有關(guān))。如圖2所示，從向下的箭頭代表以物方點(diǎn)為單位分配的任務(wù)，向上的箭頭代表以像方點(diǎn)為單位分配的任務(wù)，可見以物方點(diǎn)為單位分配的任務(wù)計(jì)算量各不相同，任務(wù)量根據(jù)該物方點(diǎn)所包含的像點(diǎn)數(shù)有關(guān)，顯然該任務(wù)分配方法不滿足并行任務(wù)負(fù)載均衡的原則。而以像點(diǎn)為單位分配的任務(wù)計(jì)算量均為一個(gè)像點(diǎn)的計(jì)算量，各任務(wù)之間負(fù)載十分均衡，有利于并行計(jì)算效率達(dá)到最大化。因此本文對(duì)該步驟采用以像點(diǎn)為單位分配計(jì)算任務(wù)，即每個(gè)像點(diǎn)的計(jì)算任務(wù)交給一個(gè)GPU線程去計(jì)算，然后累加每個(gè)線程的結(jié)果，即可得到完整的常數(shù)項(xiàng)向量和預(yù)條件矩陣。

圖2 以物方點(diǎn)為單位以及以像方點(diǎn)為單位的任務(wù)分配方式Fig.2 The ground-point based task assignment

1.3.1.2 法方程系數(shù)矩陣與殘差向量(預(yù)條件共軛梯度法)的矩陣向量積的并行計(jì)算任務(wù)分配

法方程系數(shù)矩陣一般是指改化后的法方程系數(shù)矩陣，即法方程里面僅含有外方位元素未知數(shù)相關(guān)的部分，地面點(diǎn)未知數(shù)部分已經(jīng)被消去了。改化的法方程系數(shù)矩陣中每一組外方位元素在對(duì)角線上存在一個(gè)矩陣，每?jī)山M外方位元素(對(duì)應(yīng)于具有重疊度的兩張影像)在非對(duì)角線上存在一個(gè)協(xié)方差矩陣，這個(gè)協(xié)方差矩陣的構(gòu)成必須是依賴于兩個(gè)像點(diǎn)，如果仍然以像點(diǎn)為單位進(jìn)行任務(wù)分配，在計(jì)算非對(duì)角線上的協(xié)方差矩陣時(shí)，不同像點(diǎn)對(duì)應(yīng)的線程需要進(jìn)行通信，此時(shí)便不滿足各線程之間獨(dú)立性要求。但如果以物方點(diǎn)為單位進(jìn)行任務(wù)分配，則各線程之間負(fù)載不均衡。針對(duì)上述問題，本文采用一種基于協(xié)方差的任務(wù)分配方法，如圖3所示，同一個(gè)物方點(diǎn)所對(duì)應(yīng)的每?jī)蓚€(gè)像點(diǎn)構(gòu)成一個(gè)計(jì)算任務(wù)，假設(shè)一個(gè)物方點(diǎn)對(duì)應(yīng)3個(gè)像點(diǎn)，則該物方點(diǎn)需要分配9個(gè)線程(如圖3中，3個(gè)編號(hào)分別為5、6、7的像點(diǎn),組成9個(gè)綠色方格，每個(gè)綠色方格代表一個(gè)線程)，每個(gè)線程處理其中每?jī)蓚€(gè)像點(diǎn)的計(jì)算任務(wù)。這樣，各線程之間的任務(wù)既滿足負(fù)載均衡的要求也滿足獨(dú)立運(yùn)算的要求了。

注：每一個(gè)小方格代表一個(gè)線程(為了簡(jiǎn)化圖片，僅用箭頭標(biāo)出了前面6個(gè)線程)，每個(gè)線程處理每?jī)蓚€(gè)像點(diǎn)的計(jì)算任務(wù)，方格中的編號(hào)代表該線程處理的兩個(gè)像點(diǎn)的編號(hào)圖3 基于協(xié)方差的任務(wù)分配方式Fig.3 Parallel task assignment based on covariance task; a rectangle represents a thread

1.3.2 GPU并行區(qū)域網(wǎng)平差內(nèi)存使用優(yōu)化

區(qū)域網(wǎng)平差的輸入數(shù)據(jù)主要包括，連接點(diǎn)數(shù)據(jù)(包括控制點(diǎn)數(shù)據(jù))和外方位元素?cái)?shù)據(jù)。相對(duì)于連接點(diǎn)數(shù)據(jù)，外方位元素?cái)?shù)據(jù)較小，而且在一次區(qū)域網(wǎng)平差迭代中，外方位元素的值不會(huì)發(fā)生改變，因此可以將其存儲(chǔ)在只讀的紋理內(nèi)存中。連接點(diǎn)數(shù)據(jù)則存儲(chǔ)在全局內(nèi)存中。連接點(diǎn)數(shù)據(jù)包括物方點(diǎn)坐標(biāo)數(shù)據(jù)及其對(duì)應(yīng)的若干個(gè)像點(diǎn)坐標(biāo)數(shù)據(jù)。如前所述，GPU中的全局內(nèi)存容量最大，但是其存取速度受到存取順序的影響，為了使每一個(gè)線程按順序讀取每一個(gè)像點(diǎn)，可以將像點(diǎn)數(shù)據(jù)和物方點(diǎn)數(shù)據(jù)分開存儲(chǔ)，所有的像點(diǎn)數(shù)據(jù)按照線程讀取順序來存儲(chǔ)，物方點(diǎn)數(shù)據(jù)則按照像點(diǎn)順序存儲(chǔ)。

按照前一節(jié)的任務(wù)分配，每一個(gè)線程都會(huì)生成一個(gè)臨時(shí)向量，如果將這些臨時(shí)向量直接累加至全局內(nèi)存中，則會(huì)導(dǎo)致大量線程同時(shí)訪問全局內(nèi)存中同一個(gè)地址，因此會(huì)嚴(yán)重影響累加的速度。但如果將每一個(gè)臨時(shí)向量都存儲(chǔ)在全局內(nèi)存中，則需要大量的全局內(nèi)存空間。為了解決這一矛盾，可考慮充分使用GPU的塊內(nèi)共享內(nèi)存。假設(shè)GPU核函數(shù)被分配了64個(gè)塊，每個(gè)塊256個(gè)線程，則總線程數(shù)為64×256=16 384。本文的內(nèi)存使用策略是將每一個(gè)塊中的所有線程(256個(gè)線程)的臨時(shí)向量存儲(chǔ)在該塊的共享內(nèi)存中。在全局內(nèi)存中開辟64個(gè)結(jié)果向量的內(nèi)存空間，每一個(gè)塊對(duì)應(yīng)一個(gè)全局內(nèi)存位置，用于累加該塊中所有線程的運(yùn)算結(jié)果。當(dāng)每個(gè)塊的所有線程都運(yùn)行結(jié)束時(shí)，使用該塊中的一個(gè)線程(一般為0號(hào)線程)將共享內(nèi)存中的臨時(shí)向量，累加至全局內(nèi)存中對(duì)應(yīng)于該塊的內(nèi)存位置，所有塊都運(yùn)行結(jié)束并累加結(jié)束后，將每個(gè)塊計(jì)算得到的結(jié)果(存儲(chǔ)至全局內(nèi)存中的64個(gè)向量)累加起來，即可得到最終完整的結(jié)果向量。這樣需要的全局內(nèi)存僅為64個(gè)結(jié)果向量的大小，需要的共享內(nèi)存僅為256個(gè)臨時(shí)向量的大小。以計(jì)算常數(shù)項(xiàng)向量為例，假設(shè)測(cè)區(qū)影像數(shù)為1024，且未知數(shù)僅包含外方位元素，則常數(shù)項(xiàng)向量所需內(nèi)存空間為1024×6×4 bytes=24 kB，而每個(gè)線程計(jì)算得到的臨時(shí)向量的大小為6×4 Bytes，使用上述共享內(nèi)存方法所需的全局內(nèi)存空間僅為64×24 KB=1.5 MB，共享內(nèi)存空間僅為256×6×4 Bytes=6 KB。而且各線程直接訪問共享內(nèi)存的效率比訪問全局內(nèi)存的效率要高得多，因此采用上述共享內(nèi)存加全局內(nèi)存的數(shù)據(jù)存儲(chǔ)方法既有利于節(jié)省內(nèi)存空間，同時(shí)也提高了內(nèi)存訪問效率。

2 試驗(yàn)結(jié)果與分析

2.1 數(shù)據(jù)及試驗(yàn)環(huán)境

本文直接使用文獻(xiàn)[10]中試驗(yàn)數(shù)據(jù)來進(jìn)行對(duì)比試驗(yàn)，共10組數(shù)據(jù)，分別是真實(shí)的無人機(jī)數(shù)據(jù)，手機(jī)拍攝的影像數(shù)據(jù)以及網(wǎng)絡(luò)上下載的圖片數(shù)據(jù)。測(cè)區(qū)的影像數(shù)最小為40，最大的測(cè)區(qū)影像數(shù)達(dá)到13 682，具體的數(shù)據(jù)信息見文獻(xiàn)[10]。本文中所有試驗(yàn)使用的主機(jī)配置為Inter(R) Core(TM) i7-6700HQ CPU 2.60 GHz,8.00 GB RAM，顯卡配置為NVIDIA GeForce GTX970M，顯存3 GB，軟件環(huán)境為Windows 10操作系統(tǒng)，GPU程序開發(fā)使用Microsoft Visual Studio 2013以及CUDA 8.0平臺(tái)。

2.2 內(nèi)存使用對(duì)比

分別采用傳統(tǒng)的方法，文獻(xiàn)[10]的方法以及本文方法對(duì)10組試驗(yàn)數(shù)據(jù)進(jìn)行處理，統(tǒng)計(jì)程序運(yùn)行過程中的實(shí)際內(nèi)存使用情況，其結(jié)果如表1所示，不同方法的內(nèi)存消耗如圖4所示。

表1不同方法的內(nèi)存消耗情況對(duì)比

Tab.1Comparisonofmemoryrequirementofconventionalmethodandproposedmethod

測(cè)區(qū)影像數(shù)傳統(tǒng)方法/MB文獻(xiàn)[10]的方法/MB本文方法/MB1 40 2 2 59252161365384161471488171574531458318863946420787961314135120819361214510249945855960136342910136825293932031477

圖4 不同方法的內(nèi)存消耗對(duì)比Fig.4 Comparison of memory requirement of different methods

從表1和圖4可以看出，當(dāng)影像數(shù)較少時(shí)，本文方法所需內(nèi)存較大，這主要是由于使用GPU計(jì)算時(shí)，CUDA配置占用了一部分CPU內(nèi)存，當(dāng)影像數(shù)增加至1000左右時(shí)，本文方法對(duì)內(nèi)存的需求明顯低于傳統(tǒng)辦法的內(nèi)存需求，文獻(xiàn)[10]的方法雖然在一定程度上減少了內(nèi)存需求，但是由于該方法還是存儲(chǔ)了壓縮后的法方程系數(shù)矩陣，而本文方法則完全不存儲(chǔ)法方程系數(shù)矩陣，因此本文方法則需要更少的內(nèi)存。傳統(tǒng)的方法對(duì)內(nèi)存的需求隨著影像數(shù)的增加呈現(xiàn)近似指數(shù)式增長(zhǎng)，當(dāng)影像數(shù)逐漸增大至10 000張以上時(shí)，傳統(tǒng)的辦法需要至少53 GB內(nèi)存，文獻(xiàn)[10]的方法經(jīng)過對(duì)法方程系數(shù)矩陣壓縮后僅需要約3 GB內(nèi)存，而本文方法不存儲(chǔ)法方程系數(shù)矩陣，僅需要約1.5 GB內(nèi)存用于存儲(chǔ)平差必需的原始數(shù)據(jù)(連接點(diǎn)，內(nèi)外方位元素?cái)?shù)據(jù)等)?？梢酝茢啵粲跋駭?shù)繼續(xù)增加，文獻(xiàn)[10]的方法所需內(nèi)存還會(huì)繼續(xù)增加，并最終導(dǎo)致普通計(jì)算機(jī)無法開辟如此大內(nèi)存而使得平差失敗，本文的GPU并行算法所需的存儲(chǔ)空間僅與平差原始數(shù)據(jù)的大小相關(guān)，若可以將這些原始數(shù)據(jù)存放在磁盤或者其他外部存儲(chǔ)設(shè)備中，理論上本文方法可以處理的影像數(shù)可以達(dá)到系統(tǒng)的磁盤或者外部存儲(chǔ)設(shè)備的容量上限，但實(shí)際情況中仍需要考慮這些磁盤或者外部存儲(chǔ)設(shè)備的讀取效率。

2.3 平差效率對(duì)比

3種方法對(duì)10組數(shù)據(jù)處理所消耗的時(shí)間對(duì)比如表2所示，不同算法的計(jì)算耗時(shí)如圖5所示。

表2 采用不同算法的所耗費(fèi)的時(shí)間統(tǒng)計(jì)

圖5 不同算法的計(jì)算耗時(shí)對(duì)比Fig.5 The computation time for different methods

從表2和圖5可以看出，傳統(tǒng)方法平差的計(jì)算耗時(shí)也隨著影像數(shù)的增加呈現(xiàn)指數(shù)式增加，文獻(xiàn)[10]的方法效率比傳統(tǒng)方法有所提升，而本文算法的計(jì)算效率最高，且隨著影像數(shù)的增加，僅僅呈現(xiàn)出線性的增長(zhǎng)趨勢(shì)。文獻(xiàn)[10]的方法的計(jì)算速度比傳統(tǒng)方法提升了3～10倍，而本文方法則大幅提升了10～25倍。且隨著數(shù)據(jù)量的增大，本文方法效率的提升幅度越大。值得指出的是，當(dāng)影像數(shù)達(dá)到4585及以上時(shí)，傳統(tǒng)的算法由于受到內(nèi)存大小的限制，無法進(jìn)行平差，文獻(xiàn)[10]的方法和本文方法仍然可以正常進(jìn)行平差處理。同時(shí)也可以推斷，當(dāng)影像數(shù)繼續(xù)增加時(shí)，即便是在配備大內(nèi)存的圖形工作站上，文獻(xiàn)[10]的方法也必將受到內(nèi)存的限制而無法進(jìn)行，若可以將原始數(shù)據(jù)存放在磁盤或者其他外部存儲(chǔ)設(shè)備中，理論上本文方法可以處理的影像數(shù)沒有上限。

2.4 平差精度對(duì)比

本文算法在內(nèi)存使用和計(jì)算效率方面均優(yōu)于傳統(tǒng)的算法和文獻(xiàn)[10]中的方法，為了進(jìn)一步驗(yàn)證本文方法在平差精度方面的表現(xiàn)，本節(jié)將比較3種方法的平差精度。還是對(duì)上述10組數(shù)據(jù)分別采用3種方法進(jìn)行處理，并將其精度指標(biāo)進(jìn)行對(duì)比分析(本文由于沒有使用控制點(diǎn)和檢查點(diǎn)，因此精度指標(biāo)用相對(duì)精度來衡量，即用像點(diǎn)殘差來表示)。具體精度統(tǒng)計(jì)信息如表3所示，3種方法平差后在x、y方向的殘差如圖6和圖7所示。

表3 不同算法結(jié)果的精度對(duì)比

圖6 3種方法平差后在x方向的殘差對(duì)比Fig.6 The comparison of reprojection error of the conventional method and proposed method in x direction

圖7 3種方法平差后在y方向的殘差對(duì)比Fig.7 The comparison of reprojection error of the conventional method and proposed method in y direction

如表3和圖6、圖7所示，文獻(xiàn)[10]的方法在平差精度上與傳統(tǒng)算法相當(dāng)，本文方法的精度略低于另外兩種方法，可能是由于GPU并行計(jì)算僅支持單精度浮點(diǎn)計(jì)算，但是這種精度差異在實(shí)際應(yīng)用中可忽略不計(jì)。

2.5 與同類方法對(duì)比

文獻(xiàn)[17]也是用GPU并行計(jì)算以及預(yù)條件共軛梯度法進(jìn)行區(qū)域網(wǎng)平差，其使用的數(shù)據(jù)大部分也來自于本文從網(wǎng)上下載的數(shù)據(jù)(http:∥grail.cs.washington.edu/projects/bal/),從該網(wǎng)址中找到了與本文中相同的3組試驗(yàn)數(shù)據(jù)，并將其試驗(yàn)精度結(jié)果與本文的精度結(jié)果進(jìn)行對(duì)比，如表4所示。

表4 與國(guó)際上同類方法的對(duì)比

該方法與本文的方法均使用了GPU并行計(jì)算，且均不存儲(chǔ)法方程系數(shù)矩陣，因此在內(nèi)存消耗方面是一致的，都只需要將連接點(diǎn)以及外方位元素?cái)?shù)據(jù)載入內(nèi)存即可；速度方面，文獻(xiàn)[17]的方法的效率略高，但是差距不大，可以說計(jì)算效率是屬于同級(jí)別的，引起這種實(shí)際效率差別的原因是多方面的，可能是使用的軟硬件環(huán)境有所區(qū)別(文獻(xiàn)[17]使用的是兩個(gè)Telsa C1060系列的GPU，而本文使用的是單個(gè)NVIDIA GeForce GTX970M GPU；文獻(xiàn)[17]中試驗(yàn)的軟件運(yùn)行環(huán)境是Linux操作系統(tǒng)，而本文是在Windows 10操作系統(tǒng)下展開)，也有可能是由于文獻(xiàn)[17]的區(qū)域網(wǎng)平差流程與本文有所差異。另外，值得注意的是，本文的方法的精度要比文獻(xiàn)[17]的方法好，本文的精度可達(dá)到子像素級(jí)，而文獻(xiàn)[17]的方法精度均大于一個(gè)像素，這其中的原因可能是由于文獻(xiàn)[17]在平差迭代過程中設(shè)置了較低的退出門檻(閾值)，迭代次數(shù)較少，這同時(shí)也可能是文獻(xiàn)[17]的速度比本文要快的原因之一。

綜上所示，本文算法不僅節(jié)省了內(nèi)存使用量，提高了計(jì)算速度，同時(shí)還保證了平差的最終精度與傳統(tǒng)算法相當(dāng)。

3 結(jié) 論

本文是文獻(xiàn)[10]的后續(xù)研究，針對(duì)文獻(xiàn)[10]中遇到的內(nèi)存瓶頸問題，本文在文獻(xiàn)[10]的研究基礎(chǔ)上，引入了計(jì)算速度更加高效的GPU并行計(jì)算，在平差的過程中不存儲(chǔ)法方程系數(shù)矩陣，而且在需要的時(shí)候，利用GPU實(shí)時(shí)的計(jì)算法方程系數(shù)矩陣中的相關(guān)位置的值，并在此基礎(chǔ)之上，建立了全新的GPU并行區(qū)域網(wǎng)平差方法流程，使之在保證計(jì)算精度的同時(shí)，可以節(jié)省內(nèi)存使用量，并提高計(jì)算速度。本文共使用了10組真實(shí)數(shù)據(jù)，包括無人機(jī)影像、傾斜影像、手機(jī)影像，以及網(wǎng)絡(luò)下載的影像分別進(jìn)行對(duì)比試驗(yàn)。通過對(duì)試驗(yàn)數(shù)據(jù)的對(duì)比分析，可以得出以下結(jié)論：①當(dāng)影像數(shù)較少時(shí)，本文算法由于使用GPU計(jì)算時(shí)CUDA配置占用了部分CPU內(nèi)存，因此導(dǎo)致總占用內(nèi)存較大，而當(dāng)影像數(shù)增大至1000以上時(shí)，本文算法比文獻(xiàn)[10]的方法進(jìn)一步節(jié)省了2～3倍的內(nèi)存空間，且隨著影像數(shù)的繼續(xù)增加，本文算法將更具優(yōu)勢(shì)；②本文算法的計(jì)算速度比文獻(xiàn)[10]的方法也提高了10～20倍；③本文算法比傳統(tǒng)算法以及文獻(xiàn)[10]的方法的精度略有降低，但該精度差異在實(shí)際應(yīng)用中處于可接受范圍。

本算法雖然節(jié)省了內(nèi)存空間，提高了計(jì)算速度，但是由于本算法使用預(yù)條件共軛梯度法以及不精確牛頓解法迭代求解法方程，同時(shí)又使用了僅有單浮點(diǎn)精度的GPU并行計(jì)算，因此本算法的穩(wěn)定性仍需要大量真實(shí)數(shù)據(jù)進(jìn)行測(cè)試驗(yàn)證。另外本文方法雖然不存儲(chǔ)法方程系數(shù)矩陣，但是內(nèi)存中仍然存儲(chǔ)了平差原始數(shù)據(jù)(連接點(diǎn)數(shù)據(jù)，內(nèi)外方位元素?cái)?shù)據(jù)等)，而隨著影像數(shù)增加，這些平差原始數(shù)據(jù)所需的內(nèi)存也會(huì)隨之增加，最后也必將導(dǎo)致內(nèi)存不足，若可以將這些平差原始數(shù)據(jù)存儲(chǔ)在磁盤或者外部存儲(chǔ)設(shè)備中，理論上本算法能夠處理的影像數(shù)沒有上限。

[1] 李德仁. 展望大數(shù)據(jù)時(shí)代的地球空間信息學(xué)[J]. 測(cè)繪學(xué)報(bào), 2016, 45(4): 379-384. DOI: 10.11947/j.AGCS.2016.20160057.

LI Deren. Towards Geo-spatial Information Science in Big Data Era[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(4): 379-384. DOI: 10.11947/j.AGCS.2016.20160057.

[2] MARQUARDT D W. An Algorithm for Least-squares Estimation of Nonlinear Parameters[J]. Journal of the Society for Industrial and Applied Mathematics, 1963, 11(2): 431-441.

[3] 陳馳, 楊必勝, 彭向陽. 低空UAV激光點(diǎn)云和序列影像的自動(dòng)配準(zhǔn)方法[J]. 測(cè)繪學(xué)報(bào), 2015, 44(5): 518-525. DOI: 10.11947/j.AGCS.2015.20130558.

CHEN Chi, YANG Bisheng, PENG Xiangyang. Automatic Registration of Low Altitude UAV Sequent Images and Laser Point Clouds[J]. Acta Geodaetica et Cartographica Sinica, 2015, 44(5): 518-525. DOI: 10.11947/j.AGCS.2015.20130558.

[4] 閆利, 費(fèi)亮, 葉志云, 等. 大范圍傾斜多視影像連接點(diǎn)自動(dòng)提取的區(qū)域網(wǎng)平差法[J]. 測(cè)繪學(xué)報(bào), 2016, 45(3): 310-317, 338. DOI: 10.11947/j.AGCS.2016.20140673.

YAN Li, FEI Liang, YE Zhiyun, et al. Automatic Tie-points Extraction for Triangulation of Large-scale Oblique Multi-view Images[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(3): 310-317, 338. DOI: 10.11947/j.AGCS.2016.20140673.

[5] 季順平, 史云. 車載全景相機(jī)的影像匹配和光束法平差[J]. 測(cè)繪學(xué)報(bào), 2013, 42(1): 94-100, 107.

JI Shunping, SHI Yun. Image Matching and Bundle Adjustment Using Vehicle-based Panoramic Camera[J]. Acta Geodaetica et Cartographica Sinica, 2013, 42(1): 94-100, 107.

[6] 王祥, 張永軍, 黃山, 等. 旋轉(zhuǎn)多基線攝影光束法平差法方程矩陣帶寬優(yōu)化[J]. 測(cè)繪學(xué)報(bào), 2016, 45(2): 170-177. DOI: 10.11947/j.AGCS.2016.20150282.

WANG Xiang, ZHANG Yongjun, HUANG Shan, et al. Bandwidth Optimization of Normal Equation Matrix in Bundle Block Adjustment in Multi-baseline Rotational Photography[J]. Acta Geodaetica et Cartographica Sinica, 2016, 45(2): 170-177. DOI: 10.11947/j.AGCS.2016.20150282.

[7] 林詒勛. 稀疏矩陣計(jì)算中的帶寬最小化問題[J]. 運(yùn)籌學(xué)學(xué)報(bào), 1983, 2(1): 20-27.

LIN Yixun. Band Width Minimization Problem in Sparse Matrix Computations[J]. Chinese Journal of Operations Research, 1983, 2(1): 20-27.

[8] GIBBS N E, POOLE JR W G, STOCKMEYER P K. An Algorithm for Reducing the Bandwidth and Profile of a Sparse Matrix[J]. SIAM Journal on Numerical Analysis, 1976, 13(2): 236-250.

[9] 鄭志鎮(zhèn), 李尚健, 李志剛. 稀疏矩陣帶寬減小的一種算法[J]. 華中理工大學(xué)學(xué)報(bào), 1998, 26(12): 43-45.

ZHENG Zhizhen, LI Shangjian, LI Zhigang. A New Algorithm for Reducing Bandwidth of Sparse Matrix[J]. Journal of Huazhong University of Science & Technology, 1998, 26(12): 43-45.

[10] 鄭茂騰, 張永軍, 朱俊峰, 等. 一種快速有效的大數(shù)據(jù)區(qū)域網(wǎng)平差方法[J]. 測(cè)繪學(xué)報(bào), 2017, 46(2): 188-197. DOI: 10.11947/j.AGCS.2017.20160293.

ZHENG Maoteng, ZHANG Yongjun, ZHU Junfeng, et al. A Fast and Effective Block Adjustment Method with Big Data[J]. Acta Geodaetica et Cartographica Sinica, 2017, 46(2): 188-197. DOI: 10.11947/j.AGCS.2017.20160293.

[11] ZHENG Maoteng, ZHANG Yongjun, ZHOU Shunping, et al. Bundle Block Adjustment of Large-scale Remote Sensing Data with Block-based Sparse Matrix Compression Combined with Preconditioned Conjugate Gradient[J]. Computers & Geosciences, 2016(92): 70-78.

[12] LIU Xin, GAO Wei, HU Zhanyi. Hybrid Parallel Bundle Adjustment for 3D Scene Reconstruction with Massive Points[J]. Journal of Computer Science and Technology, 2012, 27(6): 1269-1280.

[13] 佟國(guó)峰, 蔣昭炎, 葉檸, 等. 大場(chǎng)景三維重建中多核并行捆集調(diào)整算法[J]. 控制與決策, 2013, 28(29): 1403-1408.

TONG Guofeng, JIANG Zhaoyan, YE Ning. Multi-core Bundle Adjustment Algorithm Using Parallel Processing in Large-scale 3D Scene Reconstruction[J]. Control and Decision, 2013, 28(29): 1403-1408.

[14] AGARWAL S, SNAVELY N, SEITZ S M, et. al. Bundle Adjustment in the Large[C]∥DANIILIDIS K, MARAGOS P, PARAGIOS N. Computer Vision-ECCV 2010. Berlin Heidelberg: Springer, 2010: 29-42.

[15] AGARWAL S, FURUKAWA Y, SNAVELY N, et al. Building Rome in a Day[J]. Communications of the ACM, 2011, 54(10): 105-112.

[16] LI Ruipeng, SAAD Y. GPU-accelerated Preconditioned Iterative Linear Solvers[J]. The Journal of Supercomputing, 2013, 63(2): 443-466.

[17] WU Changchang. Multicore Bundle Adjustment[C]∥Proceedings of 2011 IEEE Conference on Computer Vision and Pattern Recognition. Colorado Springs, CO: IEEE, 2011: 3057-3064.

[18] CHOUDHARY S, GUPTA S, NARAYANAN P J. Practical Time Bundle Adjustment for 3D Reconstruction on the GPU[C]∥Proceedings of the 11th European Conference on Trends and Topics in Computer Vision-Volume Part II. Heraklion, Crete, Greece: Springer, 2010: 423-435.

[20] BENNER P, EZZATTI P, KRESSNER D, et al. A Mixed-Precision Algorithm for the Solution of Lyapunov Equations on Hybrid CPU-GPU Platforms[J]. Parallel Computing Archive, 2011, 37(8): 439-450.

[21] BHASKARAN-NAIR K, MA Wenjing, KRISHNAMOORTHY S, et al. Noniterative Multireference Coupled Cluster Methods on Heterogeneous CPU-GPU Systems[J]. Journal of Chemical Theory and Computation, 2013, 9(4): 1949-1957.

[22] CHAI Jun, SU Huayou, WEN Mei, et al. Resource-efficient Utilization of CPU/GPU-based Heterogeneous Supercomputers for Bayesian Phylogenetic Inference[J]. The Journal of Supercomputing, 2013, 66(1): 364-380.

[23] TAN Y S, LEE B S, HE Bingsheng, et al. A Map-reduce Based Framework for Heterogeneous Processing Element Cluster Environments[C]∥Proceedings of the 12th IEEE/ACM International Symposium on Cluster, Cloud and Grid Computing (CCGrid). Ottawa, ON, Canada: IEEE, 2012: 57-64.

[24] WEN Mei, SU Huayou, WEI Wenjie, et al. Using 1000+ GPUs and 10000+ CPUs for Sedimentary Basin Simulations[C]∥Proceedings of 2012 IEEE International Conference on Cluster Computing (CLUSTER). Beijing, China: IEEE, 2012: 27-35.

[25] NEWCOMBE R A, LOVEGROVE S J, DAVISON A J. DTAM: Dense Tracking and Mapping in Real-time[C]∥Proceedings of 2011 IEEE International Conference on Computer Vision. Barcelona, Spain: IEEE, 2011: 2320-2327.

[26] ACOSTA A, CORUJO R, BLANCO V, et al. Dynamic Load Balancing on Heterogeneous Multicore/MultiGPU Systems[C]∥Proceedings of 2010 International Conference on High Performance Computing and Simulation (HPCS). Caen, France: IEEE, 2010: 467-476.

[27] AGULLEIRO J I, VZQUEZ F, GARZN E M, et al. Hybrid Computing: CPU+ GPU Co-processing and Its Application to Tomographic Reconstruction[J]. Ultramicroscopy, 2012(115): 109-114.

[28] AGULLO E, AUGONNET C, DONGARRA J, et al. QR Factorization on a Multicore Node Enhanced with Multiple GPU Accelerators[C]∥Proceedings of 2011 IEEE International Parallel & Distributed Processing Symposium. Anchorage, AK: IEEE, 2011: 932-943.

[29] HNSCH R, DRUDE I, HELLWICH O. Modern Methods of Bundle Adjustment on the GPU[C]∥ISPRS Annals of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume III-3, 2016 XXIII ISPRS Congress. Prague, Czech Republic: Copernicus Publications,

GPU Parallel Bundle Block Adjustment

ZHENG Maoteng1,ZHOU Shunping1,XIONG Xiaodong2,ZHU Junfeng2

1. National Engineering Research Center for Geographic Information System, China University of Geosciences, Wuhan 430074, China； 2. Smart Mapping Technology Lnc., Beijing 100830, China

To deal with massive data in photogrammetry, we introduce the GPU parallel computing technology. The preconditioned conjugate gradient and inexact Newton method are also applied to decrease the iteration times while solving the normal equation. A brand new workflow of bundle adjustment is developed to utilize GPU parallel computing technology. Our method can avoid the storage and inversion of the big normal matrix, and compute the normal matrix in real time. The proposed method can not only largely decrease the memory requirement of normal matrix, but also largely improve the efficiency of bundle adjustment. It also achieves the same accuracy as the conventional method. Preliminary experiment results show that the bundle adjustment of a dataset with about 4500 images and 9 million image points can be done in only 1.5 minutes while achieving sub-pixel accuracy.

GPU parallel computing; bundle adjustment; preconditioned conjugate gradient; inexact Newton method; big data

The National Natural Science Foundation of China (No. 41601502); The China Postdoctoral Science Foundation (No. 2015M572224); The Fundamental Research Funds for the Central Universities(Nos. CUG160838；CUG170664)

ZHENG Maoteng(1987—)，male, PhD, associate researcher, majors in photogrammetry and computer vision.

ZHOU Shunping

P237

A

1001-1595(2017)09-1193-09

國(guó)家自然科學(xué)基金(41601502)；中國(guó)博士后科學(xué)基金面上項(xiàng)目(2015M572224)；中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)資金(CUG160838;CUG170664)

(責(zé)任編輯：張艷玲)

2016-12-22

修回日期： 2017-08-14

鄭茂騰(1987—)，男，博士，副研究員，研究方向?yàn)楹娇蘸教鞌z影測(cè)量與計(jì)算機(jī)視覺。

E-mail： tengve@163.com

周順平

E-mail： zhoushunping@mapgis.com