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

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

鄭茂騰，周順平，熊小東，等.GPU并行區(qū)域網(wǎng)平差[J].測繪學(xué)報，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. 中國地質(zhì)大學(xué)(武漢)國家地理信息系統(tǒng)工程技術(shù)研究中心，湖北 武漢 430074； 2. 北京中測智繪有限責(zé)任公司，北京 100830

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

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

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

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

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

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

1 理論與方法

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

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

(1)

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

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

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

引入GPU并行計算、預(yù)條件共軛梯度法及不精確牛頓解法后，全新的區(qū)域網(wǎng)平差技術(shù)流程如圖1所示。與文獻[10]中的流程類似，不同點是本文的流程中有兩個步驟是在GPU中運行的(圖1中，被填充的步驟)，其計算結(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)平差流程可以分為兩個迭代循環(huán)，一個是區(qū)域網(wǎng)平差迭代，一個是預(yù)條件共軛梯度法迭代，后者包含于前者，每一次區(qū)域網(wǎng)平差迭代，都包含一個預(yù)條件共軛梯度法迭代循環(huán)，該迭代循環(huán)是用來求解法方程的。只有當(dāng)區(qū)域網(wǎng)平差迭代收斂后，整個區(qū)域網(wǎng)平差流程結(jié)束。

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

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

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

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

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

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

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

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

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

注：每一個小方格代表一個線程(為了簡化圖片，僅用箭頭標(biāo)出了前面6個線程)，每個線程處理每兩個像點的計算任務(wù)，方格中的編號代表該線程處理的兩個像點的編號圖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ù)主要包括，連接點數(shù)據(jù)(包括控制點數(shù)據(jù))和外方位元素數(shù)據(jù)。相對于連接點數(shù)據(jù)，外方位元素數(shù)據(jù)較小，而且在一次區(qū)域網(wǎng)平差迭代中，外方位元素的值不會發(fā)生改變，因此可以將其存儲在只讀的紋理內(nèi)存中。連接點數(shù)據(jù)則存儲在全局內(nèi)存中。連接點數(shù)據(jù)包括物方點坐標(biāo)數(shù)據(jù)及其對應(yīng)的若干個像點坐標(biāo)數(shù)據(jù)。如前所述，GPU中的全局內(nèi)存容量最大，但是其存取速度受到存取順序的影響，為了使每一個線程按順序讀取每一個像點，可以將像點數(shù)據(jù)和物方點數(shù)據(jù)分開存儲，所有的像點數(shù)據(jù)按照線程讀取順序來存儲，物方點數(shù)據(jù)則按照像點順序存儲。

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

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

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

本文直接使用文獻[10]中試驗數(shù)據(jù)來進行對比試驗，共10組數(shù)據(jù)，分別是真實的無人機數(shù)據(jù)，手機拍攝的影像數(shù)據(jù)以及網(wǎng)絡(luò)上下載的圖片數(shù)據(jù)。測區(qū)的影像數(shù)最小為40，最大的測區(qū)影像數(shù)達到13 682，具體的數(shù)據(jù)信息見文獻[10]。本文中所有試驗使用的主機配置為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平臺。

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

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

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

Tab.1Comparisonofmemoryrequirementofconventionalmethodandproposedmethod

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

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

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

2.3 平差效率對比

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

表2 采用不同算法的所耗費的時間統(tǒng)計

圖5 不同算法的計算耗時對比Fig.5 The computation time for different methods

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

2.4 平差精度對比

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

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

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

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

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

2.5 與同類方法對比

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

表4 與國際上同類方法的對比

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

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

3 結(jié) 論

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

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

[1] 李德仁. 展望大數(shù)據(jù)時代的地球空間信息學(xué)[J]. 測繪學(xué)報, 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激光點云和序列影像的自動配準(zhǔn)方法[J]. 測繪學(xué)報, 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] 閆利, 費亮, 葉志云, 等. 大范圍傾斜多視影像連接點自動提取的區(qū)域網(wǎng)平差法[J]. 測繪學(xué)報, 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]. 測繪學(xué)報, 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]. 測繪學(xué)報, 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]. 運籌學(xué)學(xué)報, 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é)報, 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]. 測繪學(xué)報, 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] 佟國峰, 蔣昭炎, 葉檸, 等. 大場景三維重建中多核并行捆集調(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

國家自然科學(xué)基金(41601502)；中國博士后科學(xué)基金面上項目(2015M572224)；中央高校基本科研業(yè)務(wù)費專項資金(CUG160838;CUG170664)

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

2016-12-22

修回日期： 2017-08-14

鄭茂騰(1987—)，男，博士，副研究員，研究方向為航空航天攝影測量與計算機視覺。

E-mail： tengve@163.com

周順平

E-mail： zhoushunping@mapgis.com