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

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

一種快速有效的大數(shù)據(jù)區(qū)域網(wǎng)平差方法

鄭茂騰1，張永軍2，朱俊峰3，熊小東3，周順平1

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

針對攝影測量影像來源多樣化、復(fù)雜化、大數(shù)據(jù)化等趨勢，傳統(tǒng)區(qū)域網(wǎng)平差算法在應(yīng)對當(dāng)前復(fù)雜多變的數(shù)據(jù)來源，矩陣排列毫無規(guī)律的法方程結(jié)構(gòu)以及大數(shù)據(jù)量帶來的高內(nèi)存需求和低計算效率等問題上，遇到了前所未有的挑戰(zhàn)，為了解決上述難題，本文引入了預(yù)條件共軛梯度法以及不精確牛頓解法求解區(qū)域網(wǎng)平差過程中的法方程，同時使用一種塊狀法方程系數(shù)矩陣壓縮存儲格式，構(gòu)建了全新的區(qū)域網(wǎng)平差技術(shù)流程。本文方法避免了直接對法方程系數(shù)矩陣的求逆，壓縮了法方程系數(shù)矩陣所需的內(nèi)存空間，使得本文算法比傳統(tǒng)算法所需計算機(jī)內(nèi)存空間大幅減少，平差計算速度明顯提升，同時保證了計算精度與傳統(tǒng)方法相當(dāng)。初步試驗證明，本文方法對4500張影像、近900萬像點數(shù)據(jù)的平差計算在普通電腦上僅需要約15 min，且計算精度達(dá)到子像素級。

區(qū)域網(wǎng)平差；預(yù)條件共軛梯度；不精確牛頓解；稀疏矩陣壓縮；大數(shù)據(jù)

區(qū)域網(wǎng)平差是遙感數(shù)據(jù)幾何處理，特別是解析空中三角測量中的關(guān)鍵技術(shù)之一。經(jīng)過幾十年的發(fā)展，其理論方法已較為成熟，應(yīng)用也十分廣泛，幾乎涉及與遙感數(shù)據(jù)幾何處理相關(guān)的各個應(yīng)用領(lǐng)域，如航空航天攝影測量，基于影像的城市三維建模等；在計算機(jī)視覺領(lǐng)域，區(qū)域網(wǎng)平差同樣是三維重建，虛擬現(xiàn)實中的十分重要的技術(shù)。傳統(tǒng)的區(qū)域網(wǎng)平差方法主要是采用列文伯格-馬夸爾特法(levenberg marquardt,LM)[1]模型，逐點構(gòu)建誤差方程，構(gòu)建法方程以及改化法方程，通過對法方程系數(shù)矩陣求逆，計算法方程的解(未知數(shù)改正數(shù))，更新未知數(shù)，通過反復(fù)迭代直到收斂，即可最終得到未知數(shù)(內(nèi)外方位元素)的解。上述模型對于處理傳統(tǒng)的規(guī)則航空攝影影像十分有效，然而隨著科技的發(fā)展，大量新興傳感器不斷涌現(xiàn)，數(shù)據(jù)獲取方式和來源發(fā)生了巨大的變化，如傾斜攝影、無人機(jī)攝影、車載相機(jī)攝影、普通手持相機(jī)攝影等，相應(yīng)的，也有很多學(xué)者對這些新數(shù)據(jù)進(jìn)行了大量的研究[1-4]，提出了一些新的方法，但平差效率問題仍然存在。與傳統(tǒng)的規(guī)則航空攝影相比，新的數(shù)據(jù)呈現(xiàn)出不同的特點，如排列不規(guī)則、重疊度的方向和大小不一等。而隨著地理空間信息逐步進(jìn)入大數(shù)據(jù)時代[5]，平差的大數(shù)據(jù)量問題也日益凸顯，新的數(shù)據(jù)源給區(qū)域網(wǎng)平差帶來了新的難題和挑戰(zhàn)。

針對傳統(tǒng)的LM區(qū)域網(wǎng)平差模型，很多學(xué)者研究了法方程的存儲和運算問題，例如，如何優(yōu)化法方程系數(shù)矩陣的最小帶寬，以減少法方程系數(shù)矩陣的內(nèi)存需求[6-9]，這些研究均是基于影像排列的規(guī)律性，如果影像的排列呈隨機(jī)性，法方程系數(shù)矩陣的結(jié)構(gòu)不滿足帶狀條件，則必須存儲整個法方程。此時法方程系數(shù)矩陣將會耗費大量內(nèi)存，即便現(xiàn)代計算機(jī)性能已經(jīng)得到大幅提升，但是當(dāng)影像數(shù)超過5000張時，普通計算機(jī)仍然會面臨內(nèi)存不足的問題，因此必須尋找其他替代方法來解決上述難題。

如前所述，法方程(改化法方程，對法方程的改化是常用辦法，為了便于表達(dá)，若無特殊說明，下文中法方程均是僅含有外方位元素未知數(shù)的改化法方程)的存儲和求逆是兩大難題，為了解決這兩大難題，可以引入一種迭代求解大型線性方程組的算法，用于替代傳統(tǒng)的直接求逆方法。最為典型的迭代求解方法包括最速下降法(steepest descent method)[10]、共軛梯度算法等。共軛梯度算法是由文獻(xiàn)[1]在1952年提出[11]，該算法采用迭代逼近的方法求解法方程，迭代次數(shù)與法方程系數(shù)矩陣的條件數(shù)密切相關(guān)。法方程系數(shù)矩陣的條件數(shù)越大，所需的迭代次數(shù)越大。為了減少迭代次數(shù)，可引入預(yù)條件矩陣，通過將該矩陣左乘于法方程等式兩邊，減少法方程系數(shù)矩陣的條件數(shù)，從而加快迭代收斂速度。預(yù)條件矩陣的選擇必須是構(gòu)造簡單、易于求逆，且左乘法方程系數(shù)矩陣之后，可以降低法方程系數(shù)矩陣的條件數(shù)。常用的預(yù)條件矩陣有雅可比預(yù)條件矩陣(Jacobi)[12-15]、對稱逐次超松弛(symmetric successive over-relaxation，SSOR)預(yù)條件矩陣[14-15]、基于QR分解的預(yù)條件矩陣[14]、平衡不完全分解(balanced incomplete factorization，BIF)預(yù)條件矩陣[16]、多尺度預(yù)條件矩陣[17]等。其中Jacobi預(yù)條件矩陣是最常用且最穩(wěn)定的預(yù)條件矩陣。引入預(yù)條件矩陣在一定程度上減少了共軛梯度法迭代次數(shù)，但當(dāng)法方程系數(shù)矩陣較大時，其迭代次數(shù)仍可達(dá)到數(shù)百次[18]，為了進(jìn)一步減少迭代次數(shù)，本文引入一種基于共軛梯度的不精確牛頓解法(inexact/truncated newton method)，其基本思想是在共軛梯度法迭代求解法方程過程中，改變迭代收斂條件，提前終止迭代。該方法用一組強(qiáng)制序列系數(shù)(forcing sequence)來替代共軛梯度法迭代中的殘差向量，每次迭代分別計算一次該系數(shù)，當(dāng)強(qiáng)制序列系數(shù)小于某個閾值時即終止迭代，此時得到的法方程解稱為它的一組不精確牛頓解[12,14]。由于整個區(qū)域網(wǎng)平差本身也是一個迭代近似逼近的過程，因此用不精確牛頓解替代法方程的精確解，在一定程度上是可行的，經(jīng)過測試表明，在保證最終精度的同時，可以大幅減少共預(yù)條件軛梯度法迭代次數(shù)。

預(yù)條件共軛梯度法避免了對法方程系數(shù)矩陣的直接求逆，在迭代求解的過程中，每次只需要計算一次法方程系數(shù)矩陣與向量的乘積即可。但該方法仍面臨法方程的存儲問題，如果不存儲法方程，則必須在每一次共軛梯度法迭代中對法方程進(jìn)行實時計算一次，這將是一個非常耗時的過程，此時可考慮引入GPU高性能并行計算方法。文獻(xiàn)[19]介紹了一種GPU并行區(qū)域網(wǎng)平差方法，使得平差的效率得到了提升，但受限于GPU的單精度運算能力，其平差精度還是不如CPU方法。本文選擇CPU方法進(jìn)行研究，提出一種對法方程系數(shù)矩陣壓縮存儲格式。常用的矩陣壓縮格式是稀疏行壓縮法(compressed sparse row，CSR)[20]，它是一種以行為單位來進(jìn)行壓縮存儲的格式，對普通稀疏矩陣較為適用。其他的壓縮方法還有如對角壓縮法(diagonal，DIA)[21]、ELLPACK/ITPACK(ELL)[22]等，但這些壓縮方法均破壞了法方程系數(shù)矩陣的塊狀結(jié)構(gòu)(每組外方位元素未知數(shù)對應(yīng)著一個6×6大小的塊)，因此實際的壓縮和恢復(fù)操作較為復(fù)雜，本文使用一種分塊壓縮方法對法方程系數(shù)矩陣進(jìn)行壓縮存儲，在法方程系數(shù)矩陣中，以每一組外方位元素未知數(shù)所構(gòu)成的6×6小區(qū)域為一個分塊，僅存儲非零塊，并建立索引，且由于法方程的對稱性，僅需要存儲上三角矩陣中的非零塊即可。在進(jìn)行矩陣-向量積運算時，當(dāng)需要某個分塊的數(shù)據(jù)時，則通過索引從內(nèi)存中取出對應(yīng)的分塊，然后乘以向量中對應(yīng)的位置，將結(jié)果累加至結(jié)果向量中，在所有分塊處理完畢之后，即可得到矩陣-向量積的結(jié)果向量。

綜上所述，本文引入預(yù)條件共軛梯度法及其不精確牛頓解法至區(qū)域網(wǎng)平差流程中，建立全新的區(qū)域網(wǎng)平差技術(shù)流程，同時使用一種分塊壓縮方法對法方程系數(shù)矩陣進(jìn)行壓縮存儲和運算[18]，以應(yīng)對區(qū)域網(wǎng)平差中大數(shù)據(jù)帶來的內(nèi)存和計算效率問題。通過分別采用本方法、傳統(tǒng)辦法以及國際上其他同類方法對總共10組真實數(shù)據(jù)進(jìn)行試驗，驗證本文方法相較于傳統(tǒng)發(fā)方法在內(nèi)存需求以及計算效率方面的優(yōu)越性，并檢查本文提出方法的平差精度。

1 理論與方法

1.1 區(qū)域網(wǎng)平差

區(qū)域網(wǎng)平差的主要目標(biāo)是恢復(fù)每張影像的內(nèi)外方位元素，從而得到像點坐標(biāo)與物方點坐標(biāo)之間的相互轉(zhuǎn)換關(guān)系，其基本理論方法已經(jīng)十分成熟，在共線條件方程的基礎(chǔ)上，對所有像點觀測值列出一個誤差方程組，然后法化得到法方程，通過求解法方程得到未知數(shù)改正數(shù)，更新未知數(shù)并重新進(jìn)入下一次迭代，直到滿足退出條件為止。以傳統(tǒng)的區(qū)域網(wǎng)平差為例，若不考慮內(nèi)方位元素作為未知數(shù)的情況，即內(nèi)方位元素為給定值，在平差之前已經(jīng)做好了檢校工作，誤差方程表達(dá)形式如式(1)所示

(1)

式中，JC、JP分別對應(yīng)誤差方程系數(shù)中的外方位元素部分和地面點未知數(shù)部分；lc、lp分別代表誤差方程常數(shù)項向量中外方位元素部分和地面點未知數(shù)部分；Δxc、Δxp分別對應(yīng)未知數(shù)改正數(shù)向量中外方位元素部分和地面點未知數(shù)部分。

將上述誤差方程法化，從而得到法方程，如式(2)所示

(2)

式中，λDC、λDP分別對應(yīng)外方位元素未知數(shù)和地面點未知數(shù)的阻尼系數(shù)(damping value)，對于沒有控制點的測區(qū)，阻尼系數(shù)有助于控制未知數(shù)改正數(shù)的變化幅度，避免因為法方程系數(shù)矩陣的奇異性造成不穩(wěn)定解的情況[17,19,23]。將式(2)簡化之后得到式(3)

(3)

通過對法方程中的地面點未知數(shù)進(jìn)行消元處理，得到改化的法方程，如式(4)所示

(4)

此時，通過解算式(4)可以得到外方位元素未知數(shù)改正數(shù)，再將其回代至式(5)中求解得到地面點坐標(biāo)未知數(shù)

VΔxp=lp-WTΔxc

(5)

但在實際處理中，地面點未知數(shù)往往通過利用新的外方位元素進(jìn)行前方交會得到。

1.2 預(yù)條件共軛梯度法

共軛梯度法是一種迭代求解大型線性方程組的方法，給定一個線性方程組(6)

By=c

(6)

式中，B是線性方程組的系數(shù)矩陣；y是未知數(shù)向量；c是常數(shù)項向量。

共軛梯度法求解的基本思想是，給定線性方程組(6)的初始解x0，然后根據(jù)線性方程組系數(shù)矩陣，常數(shù)項向量，計算共軛梯度法迭代的殘差向量，方向向量，然后計算線性方程組的新解x1，重復(fù)上述過程，直到其方向向量的模小于給定的閾值，也即線性方程組的解的改變量小于給定閾值，即可退出迭代，此時的解xn即為線性方程組的最終解。

理論上共軛梯度法收斂的次數(shù)與線性方程組系數(shù)矩陣的條件數(shù)密切相關(guān)，條件數(shù)越大，需要迭代的次數(shù)越多。因此，為了減少迭代次數(shù)，則需要降低線性方程組系數(shù)矩陣的條件數(shù)?？梢栽诰€性方程組(6)兩端分別左乘一個預(yù)條件矩陣M-1，得到式(7)

M-1By=M-1c

(7)

此時，系數(shù)矩陣變?yōu)镸-1B，其迭代次數(shù)也相應(yīng)變?yōu)榫仃嘙-1B的條件數(shù)，通過選取合適預(yù)條件矩陣M，可以有效降低線性方程組的條件數(shù)，從而減少解算時的迭代次數(shù)。預(yù)條件矩陣M的選取需要滿足一定的原則，即構(gòu)造簡單且易于求逆，同時與系數(shù)矩陣相乘后可以減少系數(shù)矩陣的條件數(shù)。本文采用應(yīng)用范圍較廣，且計算簡單，效果穩(wěn)定的塊狀Jacobi預(yù)條件矩陣。其他的一些預(yù)條件矩陣可能效果更好，但是計算復(fù)雜且穩(wěn)定性差[18]。

在區(qū)域網(wǎng)平差計算中，可以對法方程(4)采用預(yù)條件共軛梯度法迭代求解，塊狀Jacobi預(yù)條件矩陣則是由法方程系數(shù)矩陣中對角線上的塊狀子矩陣構(gòu)成，其結(jié)構(gòu)示意圖如圖1所示。

圖中每一個子塊由3×3個元素組成(本文中每組外方位元素有6個元素，對應(yīng)的塊狀Jacobi預(yù)條件矩陣是由對角線上的6×6大小的子塊組成)。

如前所述，引入預(yù)條件矩陣的目的是為了在共軛梯度法求解法方程的過程中，減少系數(shù)矩陣的條件數(shù)，從而減少共軛梯度法的迭代次數(shù)。塊狀Jacobi預(yù)條件矩陣是構(gòu)造最為簡單，易于求逆，且效果最為穩(wěn)定的預(yù)條件矩陣，滿足預(yù)條件矩陣的選取原則，因此本文選擇塊狀Jacobi預(yù)條件矩陣來進(jìn)行處理。

預(yù)條件共軛梯度法的具體技術(shù)流程如下所示：

針對線性方程組(6)，給定初始參數(shù):y0=0,r0=c-By0,d0=M-1r0,k=0;

Do

2:yk+1=yk+αkdk

3:rk+1=rk-αkBdk

5:dk+1=M-1rk+1+βkdk

6:k=k+1

1.3 不精確牛頓解法

如圖2所示，本文算法包括兩種迭代過程，一種是區(qū)域網(wǎng)平差迭代，一種是預(yù)條件共軛梯度法迭代求解法方程，后者被包含于前者當(dāng)中，預(yù)條件共軛梯度法迭代的目的是求得法方程的解，然后利用這組解繼續(xù)進(jìn)行區(qū)域網(wǎng)平差迭代。由于兩種迭代都是近似逼近的過程，對預(yù)條件共軛梯度法迭代而言，該迭代過程只有在前幾次迭代會對未知數(shù)有較大的改進(jìn)，后面的若干次迭代對未知數(shù)的改進(jìn)并不大，而這后面的若干次迭代仍然會消耗大量的計算資源和時間，因此可以考慮提前結(jié)束預(yù)條件共軛梯度法迭代，此時雖然只得到了一組法方程的不精確解，但如果使用這組不精確的解繼續(xù)進(jìn)行區(qū)域網(wǎng)平差迭代，不影響區(qū)域網(wǎng)平差的最終精度，則可認(rèn)為提前結(jié)束預(yù)條件共軛梯度法是可行的，此時，只需要確定在何時停止迭代即可。

本文引入一種不精確牛頓解法，其基本思想是引入新的迭代收斂條件(判斷強(qiáng)制序列系數(shù)Forcing sequence是否小于指定閾值)，通過設(shè)定強(qiáng)制系數(shù)的閾值，提前終止迭代，用法方程的一種不精確解替代其精確解，從而達(dá)到減少迭代次數(shù)的目的。強(qiáng)制序列系數(shù)ηk的計算方法如式(8)所示

(8)

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

引入預(yù)條件共軛梯度法、不精確牛頓解以及法方程系數(shù)矩陣的壓縮存儲方法后，傳統(tǒng)的區(qū)域網(wǎng)平差流程需要進(jìn)行徹底的改變，傳統(tǒng)的逐點生成法方程后，更新法方程時，必須根據(jù)索引找到該點生成子塊在壓縮后法方程中的位置，然后取出對應(yīng)的子塊，并與當(dāng)前子塊累加，然后將結(jié)果更新至壓縮存儲的法方程中，求解法方程時，不再需要對法方程進(jìn)行求逆運算，而是區(qū)域按照1.2節(jié)中的步驟，通過預(yù)條件共軛梯度法迭代求解法方程。全新的區(qū)域網(wǎng)平差技術(shù)流程如圖2所示。

1.5 法方程系數(shù)矩陣的壓縮存儲與計算

一般來講，在傳統(tǒng)區(qū)域網(wǎng)平差中，如果采用算法對法方程系數(shù)矩陣進(jìn)行壓縮，由于改變了法方程的存儲結(jié)構(gòu)，則在對壓縮的法方程求逆時，其操作非常復(fù)雜，而本文采用預(yù)條件共軛梯度法迭代求解法方程，無須對法方程系數(shù)矩陣進(jìn)行直接求逆運算，求解過程中只需要每次迭代時計算法方程系數(shù)矩陣與某個向量的乘積即可，因此該方法為法方程系數(shù)矩陣的壓縮存儲提供可能，同時，為了適應(yīng)區(qū)域網(wǎng)平差中法方程系數(shù)矩陣的塊狀結(jié)構(gòu)，本文使用一種分塊矩陣壓縮方法[18]。該方法的基本思想是以每一個子塊矩陣為單元(每個子塊對應(yīng)兩組外方位元素的協(xié)方差，即每個子塊大小為6×6)，僅存儲非零子塊，舍去為零的子塊，將法方程系數(shù)矩陣中的上三角部分(由于法方程是對稱矩陣，因此只需要存儲其上三角矩陣：圖3中的黑色填充塊)的非零子塊按照從左到右，從上到下的順序依次存儲至一個一維存儲結(jié)構(gòu)，如圖4所示。同時記錄每個子塊在原始法方程中的位置和大小，即圖4中的第2行，其具體存儲方法如圖3和圖4所示。

圖2 基于預(yù)條件共軛梯度法的區(qū)域網(wǎng)平差技術(shù)流程Fig.2 The whole workflow of the bundle adjustment with PCG algorithm

圖3 待壓縮的系數(shù)塊狀矩陣Fig.3 The uncompressed block matrix

從圖3和表1可以看出，本文使用的法方程系數(shù)矩陣壓縮存儲方法的最小操作單位是一個子塊矩陣，該子矩陣對應(yīng)一組外方位元素未知數(shù)，由于對法方程系數(shù)矩陣的運算均是以子塊矩陣為單位，因此本壓縮算法可以繼承法方程的塊狀結(jié)構(gòu)，保留了法方程系數(shù)矩陣在構(gòu)建，存儲以及運算過程中的塊狀結(jié)構(gòu)特性。同時，本文采用預(yù)條件共軛梯度法求解法方程，無需對法方程系數(shù)矩陣進(jìn)行直接求逆運算，僅需要計算法方程系數(shù)矩陣與某個向量的乘積，因此可以考慮將矩陣與向量的乘積，分解為各個子塊矩陣分別與向量對應(yīng)位置的乘積，然后累加。由于壓縮算法保留了塊狀結(jié)構(gòu)，因此壓縮存儲后，可以方便地進(jìn)行矩陣向量積運算。本算法的壓縮效率與法方程系數(shù)矩陣的稀疏程度密切相關(guān)，法方程系數(shù)矩陣越稀疏，則其壓縮效率越高，反之亦然。

表1 分塊矩陣壓縮的存儲結(jié)構(gòu)

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

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

采用10組數(shù)據(jù)對本文提出的方法進(jìn)行試驗，這10組數(shù)據(jù)中，一部分是由數(shù)碼相機(jī)以及手機(jī)拍攝得到，還有一部分是直接從網(wǎng)絡(luò)上下載的拍攝同一區(qū)域的影像[12]，這些影像從拍攝角度、距離、分辨率以及相機(jī)的焦距等各不相同，影像的分布呈隨機(jī)性，不滿足規(guī)則排列。為了驗證本文提出的方法在內(nèi)存使用效率，計算效率以及平差精度等個方面的性能，并與傳統(tǒng)的算法在上述各個方面進(jìn)行全方位對比。本文使用的數(shù)據(jù)信息如表2所示。

表2 試驗數(shù)據(jù)信息統(tǒng)計

表2中，地面點是指匹配得到的地面連接點數(shù)量，影像點則是這些地面點對應(yīng)于不同影像上的同名像點。另外，除了第3組數(shù)據(jù)中存在少量地面控制點，用于檢查控制點平差精度，其余測區(qū)由于均不是來自于常規(guī)航空攝影，大多為網(wǎng)絡(luò)下載的影像，因此基本上無地面控制點數(shù)據(jù)可用。所有測區(qū)平差均采用無控制點自由網(wǎng)平差方式，第3組數(shù)據(jù)則額外進(jìn)行帶控制點平差，用于檢查控制點平差精度，具體結(jié)果見2.5節(jié)。除了網(wǎng)絡(luò)下載的影像外(網(wǎng)絡(luò)下載影像數(shù)據(jù)由Sameer在其個人網(wǎng)站公開的數(shù)據(jù)，且?guī)в幸呀?jīng)匹配好的連接點數(shù)據(jù))[12]，其余測區(qū)的連接點數(shù)據(jù)均由本課題組研發(fā)的空三匹配軟件匹配得到，所使用的平差軟件是由作者自主開發(fā)的區(qū)域網(wǎng)平差軟件。本文中所有試驗使用的硬件配置為Inter (R) Core(TM) i5-33320M CPU 2.60GHz,8.00GB RAM，軟件環(huán)境為Windows 10操作系統(tǒng)。

2.2 權(quán)策略

在區(qū)域網(wǎng)平差中，觀測值權(quán)的選擇可以在一定程度上削弱粗差對平差系統(tǒng)的影響，本文的給觀測值的初始權(quán)均為1，在每次區(qū)域網(wǎng)平差迭代之后，根據(jù)各觀測值的殘差重新計算各觀測值的權(quán)值，若觀測值的殘差大于給定的閾值，則給該觀測值一個較小的權(quán)，若觀測值的殘差小于給定的閾值，則觀測值的權(quán)保持不變。未知數(shù)的權(quán)為單位1，如1.1節(jié)所述，由于在法方程中加入了阻尼項(damping term)，因此可以約束未知數(shù)的改正數(shù)，使其在特定的范圍變化，以避免平差系統(tǒng)因為沒有控制點而發(fā)生迭代不收斂的情況。

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

對10組試驗數(shù)據(jù)，分別統(tǒng)計傳統(tǒng)算法以及本文提出算法的內(nèi)存使用情況，該內(nèi)存包括了整個算法中各個數(shù)據(jù)項占用內(nèi)存的總和(包括原始數(shù)據(jù)，法方程系數(shù)矩陣，常數(shù)項向量，其他變量等)，也即算法程序運行的實際內(nèi)存，其結(jié)果如表3所示。

表3 傳統(tǒng)算法及本文方法的內(nèi)存使用情況對比

Tab.3 Comparison of memory requirement of conventional method and proposed method

數(shù)據(jù)影像數(shù)傳統(tǒng)方法/MB本文方法/MB1 40 2．0 1．525215．713．238416．314．448817．014．9531457．730．8639463．619．937961313．7134．5819361214．1509．6945855959．51363．2101368252939．23202．6

從表3和圖4中可以看出，本文方法對內(nèi)存的需求明顯低于傳統(tǒng)辦法的內(nèi)存需求，傳統(tǒng)的方法對內(nèi)存的需求隨著影像數(shù)的增加呈近似指數(shù)式增長，而本文方法則僅呈現(xiàn)線性增長。當(dāng)影像數(shù)增加到13 682時，傳統(tǒng)的辦法需要至少53 GB內(nèi)存，才能進(jìn)行運算，而本文方法則只需要約3 GB內(nèi)存即可。

圖4 傳統(tǒng)方法與本文方法的內(nèi)存使用對比圖Fig.4 Comparison of memory requirement of conventional method and proposed method

2.4 平差效率對比

分別使用傳統(tǒng)的方法和本文方法對10組試驗數(shù)據(jù)進(jìn)行區(qū)域網(wǎng)平差處理(圖5)，并記錄每組數(shù)據(jù)采用不同平差方法所使用的時間，如表4所示。

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

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

從表4和圖5可以看出，傳統(tǒng)方法平差的計算耗時也隨著影像數(shù)的增加呈指數(shù)式增加，而本文算法則隨著影像數(shù)的增加呈近似的線性增長趨勢。當(dāng)影像數(shù)較小時，本文算法和傳統(tǒng)的算法計算耗時相當(dāng)，無明顯區(qū)別，而隨著影像數(shù)逐漸增加至200以上時，本文算法逐漸顯現(xiàn)出其優(yōu)越性，繼續(xù)增加影像數(shù)，本文算法的優(yōu)勢更加明顯，計算速度可提升3～10倍。值得指出的是，當(dāng)影像數(shù)達(dá)到4585及以上時，傳統(tǒng)的算法由于受到內(nèi)存大小的限制，無法進(jìn)行平差，而本文算法由于采用了法方程系數(shù)矩陣壓縮存儲算法，仍然可以正常進(jìn)行平差處理，且計算耗時在可接受范圍內(nèi)。此時的計算速度比影像數(shù)為1936時的計算速度還要快，這是由于在本文的預(yù)條件共軛梯度法迭代求解法方程過程中，其迭代次數(shù)根據(jù)不同觀測網(wǎng)型結(jié)構(gòu)、粗差數(shù)量以及內(nèi)外方位元素初始值的準(zhǔn)確度的不同而有明顯區(qū)別，因此出現(xiàn)影像數(shù)較大的測區(qū)比影像數(shù)較小的測區(qū)平差所需時間要少的情況是可能的。在傳統(tǒng)算法中，影響平差系統(tǒng)計算速度的主要是法方程的求逆運算，而隨著影像數(shù)增加，法方程的大小呈指數(shù)級增長，因而其計算耗時也呈指數(shù)級增長。本文算法采用預(yù)條件共軛梯度算法求解法方程，無需對法方程系數(shù)矩陣進(jìn)行直接求逆運算，且由于采用了法方程系數(shù)矩陣壓縮算法，平差系統(tǒng)的耗時大大降低，僅隨著影像數(shù)的增加而呈線性增長。

2.5 平差精度對比

本文算法在內(nèi)存使用和計算效率方面均優(yōu)于傳統(tǒng)的算法，為了全方位驗證本文方法的優(yōu)越性，還需要比較本文算法與傳統(tǒng)算法在平差精度方面的表現(xiàn)。還是對上述10組數(shù)據(jù)分別采用傳統(tǒng)算法和本文算法進(jìn)行處理，并將其精度指標(biāo)進(jìn)行對比分析(本文由于沒有使用控制點和檢查點，因此精度指標(biāo)用相對精度來衡量，即用像點殘差來表示)。具體精度統(tǒng)計信息如表5所示。

表5 本文算法與傳統(tǒng)算法的精度對比

Tab.5 The comparison of accuracy of the conventional method and proposed method 像素

數(shù)據(jù)影像數(shù)傳統(tǒng)算法本文算法xyxy1400．4790．6240．4430．5532520．7510．5640．7500．5633840．4570．1950．4580．1934880．6850．6890．5250．54753140．2820．2470．2820．24763940．6440．6010．6440．60179610．7320．6390．7320．639819360．8170．8090．8180．80994585N/AN/A0．8190．7401013682N/AN/A0．7810．746

如表5和圖6、7所示，本文算法在平差精度上與傳統(tǒng)算法相當(dāng)，二者的精度差異在實際應(yīng)用中可忽略不計。

圖6 本文算法與傳統(tǒng)算法x方向的殘差對比Fig.6 The comparison of reprojection error of the conventional method and proposed method in x direction

圖7 本文算法與傳統(tǒng)算法y方向的殘差對比Fig.7 The comparison of reprojection error of the conventional method and proposed method in y direction

以上試驗均是不帶控制點進(jìn)行的自由網(wǎng)區(qū)域網(wǎng)平差，為了驗證本算法在有控制點情況下的平差精度，本文在第3個測區(qū)采用共27個外業(yè)測量點，其中13個點作為控制點，剩余點作為檢查點。采用帶控制點進(jìn)行平差后，如表6所示，在像點殘差和檢查點殘差方面本文方法與傳統(tǒng)方法仍然具有相當(dāng)?shù)木人健?/p>

表6 第3個測區(qū)帶控制點區(qū)域網(wǎng)平差精度統(tǒng)計

2.6 與同類方法對比

文獻(xiàn)[22]也是用預(yù)條件共軛梯度法進(jìn)行區(qū)域網(wǎng)平差，其使用的數(shù)據(jù)大部分也來自于本文從網(wǎng)上下載的數(shù)據(jù)(由Sameer在其個人網(wǎng)站公開的數(shù)據(jù)),但使用的是GPU并行區(qū)域網(wǎng)平差，在求解的過程中，法方程系數(shù)矩陣是沒有以任何形式存儲在法方程中的，只是在需要的時候，利用GPU實時計算出來的，其計算速度因而比本文算法要快。從該文獻(xiàn)中找到了與本文中相同的3組試驗數(shù)據(jù)，并將其試驗精度結(jié)果與本文的精度結(jié)果進(jìn)行對比，如表7所示。

表7 本文方法與國際上同類方法的對比

文獻(xiàn)[22]中的精度沒有用具體的數(shù)字以表格的形式給出，而是給出了幾張對比圖形，但是可以通過圖形中刻度值得到其精度的大概范圍。該方法由于使用了GPU并行計算，其計算速度也因此比本文的方法要快，但如表7所示，其計算精度普遍較差，可見該方法雖然有很快的計算速度，但是由于GPU的單精度浮點運算，其計算精度明顯低于本文的方法。

綜上所示，本文算法不僅節(jié)省了內(nèi)存使用量，提高了計算速度，同時還保證了平差的最終精度幾乎無損，因此本算法非常適合于除了規(guī)則航空攝影以外的其他攝影測量區(qū)域網(wǎng)平差問題。傳統(tǒng)算法和本文算法的高程精度相當(dāng)，但都不是理想的高程精度，可能是由于沒有對內(nèi)方位元素做自檢校，即沒有對相機(jī)主光軸傾角，鏡頭畸變等系統(tǒng)誤差進(jìn)行有效補(bǔ)償，僅解算了外方位元素。本文方法雖不如GPU并行計算方法效率高，但是在精度上卻有明顯的優(yōu)勢。

3 結(jié) 論

針對各類新型傳感器給攝影測量數(shù)據(jù)源帶來的多源性、隨機(jī)性、大數(shù)據(jù)等方面的挑戰(zhàn)，本文引入預(yù)條件共軛梯度法以及不精確牛頓解法求解區(qū)域網(wǎng)平差中的法方程，使用一種有效的法方程系數(shù)矩陣壓縮格式，并在此基礎(chǔ)之上建立了全新的區(qū)域網(wǎng)平差方法流程，使之在保證計算精度的同時，可以節(jié)省內(nèi)存使用量，并提高計算速度。本文共使用了10組真實數(shù)據(jù)，包括無人機(jī)影像、傾斜影像、手機(jī)影像以及網(wǎng)絡(luò)下載的影像分別進(jìn)行對比試驗。通過對試驗數(shù)據(jù)的對比分析，可以得出以下結(jié)論：①本文算法相較于傳統(tǒng)的算法可以節(jié)省3～15倍內(nèi)存空間，影像數(shù)越大，壓縮效率越高；②本文算法的計算速度比傳統(tǒng)算法大幅提升，且數(shù)據(jù)量越大，計算速度提升倍率越高；③本文算法保證了精度與傳統(tǒng)算法相當(dāng)。

本算法雖然在一定程度上節(jié)省了內(nèi)存空間，提高了計算速度，但是由于本算法使用預(yù)條件共軛梯度法以及不精確牛頓解法迭代求解法方程，因此本算法的穩(wěn)定性仍需要大量真實數(shù)據(jù)進(jìn)行測試驗證。本文僅解算外方位元素，內(nèi)方位元素是在平差之前進(jìn)行檢校的，因此無需解決內(nèi)外方位元素之間的相關(guān)性這一傳統(tǒng)的平差難題，但如果內(nèi)外方位元素同時求解，也可以使用本文的算法提高平差的效率，節(jié)省內(nèi)存空間，只是本文算法對于內(nèi)外方位元素的相關(guān)性問題，也并沒有明確的優(yōu)勢。另外，本算法存儲了法方程，當(dāng)影像數(shù)大于15 000并繼續(xù)增大時，本文算法所需的內(nèi)存仍會繼續(xù)增加，且最終一定會超過普通電腦甚至是專業(yè)工作站的內(nèi)存，因此需要研究新的方法，徹底解決法方程系數(shù)矩陣的內(nèi)存需求問題，即不存儲法方程系數(shù)矩陣，而是在平差迭代過程中實時計算法方程系數(shù)矩陣，此時需要高性能的計算設(shè)備，引入GPU并行計算和大型服務(wù)器并行計算是可能的方法，同時還要解決引入GPU后，僅使用單精度浮點運算帶來的精度損失問題，本文后續(xù)將進(jìn)一步研究此類方法。

[1] 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.

[2] 陳馳, 楊必勝, 彭向陽. 低空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.

[3] 閆利, 費亮, 葉志云, 等. 大范圍傾斜多視影像連接點自動提取的區(qū)域網(wǎng)平差法[J]. 測繪學(xué)報, 2016, 45(3): 310-317. 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. DOI: 10.11947/j.AGCS.2016.20140673.

[4] 季順平, 史云. 車載全景相機(jī)的影像匹配和光束法平差[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.

[5] 李德仁. 展望大數(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.

[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. Bandwidth 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(1): 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(1): 43-45.

[10] FRADSEN P E, JONASSON K, NIELSEN H B, et al. Unconstrained Optimization[M/OL]. 3rd ed. Denmark: Informatics and Mathematical Modelling, Technical University of Denmark, 2004. http:∥www.imm.dtu.dk/courses/02611/uncon.pdf.

[11] HESTENES M R, STIEFEL E. Methods of Conjugate Gradients for Solving Linear Systems[J]. Journal of Research of the National Bureau of Standards, 1952, 49(6): 409-436.

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

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

[14] BYR? D M, ?STR? M K. Conjugate Gradient Bundle Adjustment[M]∥DANIILIDIS K, MARAGOS P, PARAGIOS N. Computer Vision-ECCV 2010. Berlin Heidelberg: Springer, 2010, 6312: 114-127.

[15] JIAN Y D, BALCAN D C, DELLAERT F. Generalized Subgraph Preconditioners for Large-scale Bundle Adjustment[C]∥Proceedings of IEEE International Conference on Computer Vision. Barcelona: IEEE, 2011: 295-302.

[16] BRU R, MARN J, MAS J, et al. Balanced Incomplete Factorization[J]. SIAM Journal on Scientific Computing, 2008, 30(5): 2302-2318.

[17] BYR? D M, ?STR? M K. Bundle Adjustment using Conjugate Gradients with Multiscale Preconditioning[C]∥Proceedings of 2009 British Machine Vision Conference. London: BMVC, 2009.

[18] 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.

[19] WU Changchang, AGARWAL S, CURLESS B, et al. Multicore Bundle Adjustment[C]∥Proceedings of IEEE Conference on Computer Vision and Pattern Recognition. Providence, RI: IEEE, 2011: 3057-3064.

[20] BELL N, GARLAND M. Implementing Sparse Matrix-Vector Multiplication on Throughput-Oriented Processors[C]∥Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis. Portland, OR: IEEE, 2009: 1-11.

[21] SAAD Y. Iterative Methods for Sparse Linear Systems[M]. 2nd ed. Philadelphia, PA: SIAM, 2003.

[22] GRIMES R G, KINCAID D R, YOUNG D M. ITPACK 2.0 User’s Guide[R]. Technical Report CNA-150. Austin, TX: Center for Numerical Analysis, University of Texas, 1980.

[23] NIELSEN H B, Damping Parameter in Marquardt’s Method[R]. Technical Report IMM-REP-1999-05. Denmark: Technical University of Denmark, 1999.

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

A Fast and Effective Block Adjustment Method with Big Data

ZHENG Maoteng1,ZHANG Yongjun2,ZHU Junfeng3,XIONG Xiaodong3,ZHOU Shunping1

1. National Engineering Research Center for Geographic Information System, China University of Geosciences, Wuhan 430074, China； 2. School of Remote Sensing and Information Engineering, Wuhan University, Wuhan 430079, China； 3. Smart Mapping Technology Inc., Beijing 100830, China

To deal with multi-source, complex and massive data in photogrammetry, and solve the high memory requirement and low computation efficiency of irregular normal equation caused by the randomly aligned and large scale datasets, we introduce the preconditioned conjugate gradient combined with inexact Newton method to solve the normal equation which do not have strip characteristics due to the randomly aligned images. We also use an effective sparse matrix compression format to compress the big normal matrix, a brand new workflow of bundle adjustment is developed. Our method can avoid the direct inversion of the big normal matrix, the memory requirement of the normal matrix is also decreased by the proposed sparse matrix compression format. Combining all these techniques, the proposed method can not only decrease the memory requirement of normal matrix, but also largely improve the efficiency of bundle adjustment while maintaining 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 15 minutes while achieving sub-pixel accuracy.

block adjustment; preconditioned conjugate gradient; inexact Newton method; sparse matrix compressing; big data

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

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

鄭茂騰，張永軍，朱俊峰，等.一種快速有效的大數(shù)據(jù)區(qū)域網(wǎng)平差方法[J].測繪學(xué)報，2017,46(2)：188-197.

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.

P237

A

1001-1595(2017)02-0188-10

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

2016-07-26

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

E-mail： tengve@163.com

修回日期： 2016-11-09