袁斌

無人機(jī)影像光束法平差研究

袁斌

（山東理工大學(xué)，山東 淄博 255000）

隨著無人機(jī)技術(shù)的不斷發(fā)展，相機(jī)所獲取的同名像點(diǎn)數(shù)量大大增加，傳統(tǒng)的光束法平差解算方法在面對(duì)大數(shù)據(jù)量所帶來的內(nèi)存消耗大、計(jì)算效率低等問題上遇到了前所未有的挑戰(zhàn)。主要圍繞降低無人機(jī)影像光束法平差內(nèi)存占用以及提高平差效率兩方面展開研究。對(duì)于光束法平差稀疏矩陣特有的塊狀結(jié)構(gòu)，利用一種特殊的分塊存儲(chǔ)計(jì)算方法，大大減少了平差過程中的內(nèi)存占用；利用預(yù)條件共軛梯度法和不精確牛頓解法求解法方程，大大提高了平差效率。

光束法平差；系數(shù)矩陣；預(yù)條件共軛梯度法；不精確牛頓解法

1 引言

無人機(jī)（unmanned aerial vehicle，UAV）依靠高靈活性、低成本且可攜帶高分辨率數(shù)碼相機(jī)等優(yōu)點(diǎn)已成為目前獲取大規(guī)模地理信息的主要方式。因自身能力所限，無人機(jī)所搭載相機(jī)的拍攝區(qū)域面積偏小，在測(cè)區(qū)面積相同的情況下，無人機(jī)影像數(shù)量相較于傳統(tǒng)的航空攝影測(cè)量會(huì)增加很多，同名像點(diǎn)的數(shù)量也會(huì)大大增加，意味著光束法平差需要處理的數(shù)據(jù)量大大增加。這會(huì)導(dǎo)致光束法平差誤差方程系數(shù)矩陣及法方程系數(shù)矩陣規(guī)模十分龐大。隨著大數(shù)據(jù)時(shí)代的到來，如何高效精確地解算大規(guī)模數(shù)據(jù)是如今光束法平差的研究熱點(diǎn)之一[1-2]。

光束法平差系數(shù)矩陣是一種大規(guī)模的稀疏矩陣。目前解決稀疏矩陣的存儲(chǔ)問題已有多種方法，如對(duì)角線存儲(chǔ)法、坐標(biāo)存儲(chǔ)法、三元組存儲(chǔ)法、超矩陣存儲(chǔ)法、CRS存儲(chǔ)法、動(dòng)態(tài)存儲(chǔ)法等[3]。經(jīng)研究發(fā)現(xiàn)，以上方法雖然能夠解決大部分稀疏矩陣的壓縮存儲(chǔ)問題，但這些壓縮存儲(chǔ)方法均無法完整保存光束法平差法方程系數(shù)矩陣特有的塊狀結(jié)構(gòu)（每張影像外方位元素一一對(duì)應(yīng)著一個(gè)大小為6×6的方陣），不能直接適用于光束法平差稀疏矩陣的壓縮與存儲(chǔ)[4]。本文將根據(jù)其特有的塊狀結(jié)構(gòu)，采用一種特殊的稀疏矩陣壓縮存儲(chǔ)方法，僅對(duì)非零元素進(jìn)行存儲(chǔ)，利用較小的內(nèi)存空間存儲(chǔ)計(jì)算大規(guī)模稀疏矩陣。

傳統(tǒng)光束法平差的解算方法主要為列文伯格-馬夸爾特法算法[5]（Levenberg Marquaedt，LM），主要通過對(duì)法方程系數(shù)矩陣的直接求逆得到最終解。但面對(duì)數(shù)據(jù)規(guī)模較大的無人機(jī)影像，LM算法已不再適用。目前，用共軛梯度法（Conjugate Gradient，CG）取代LM算法解算光束法平差的想法引起了人們的關(guān)注[6]。共軛梯度法能充分利用稀疏矩陣的特性，在計(jì)算中僅將非零元素存儲(chǔ)在稀疏矩陣中。研究表明，共軛梯度法存儲(chǔ)量小，計(jì)算簡便，結(jié)果穩(wěn)定，可直接用于光束法區(qū)平差的解算[7]。共軛梯度法的迭代次數(shù)與系數(shù)矩陣的條件數(shù)成正比，系數(shù)矩陣的條件數(shù)越大，迭代次數(shù)越多。因此，可通過減小系數(shù)矩陣的條件數(shù)來達(dá)到提高迭代效率的目的。一般通過在線性方程組左右兩邊同時(shí)左乘一個(gè)預(yù)條件矩陣來減少系數(shù)矩陣的條件數(shù)。常用的預(yù)條件矩陣有雅可比預(yù)條件矩陣（Jacobi）、對(duì)稱逐次超松弛預(yù)條件矩陣（Symmetric Successive Over-relaxation，SSOR）、基于QR分解的預(yù)條件矩陣、平衡不完全分解預(yù)條件矩陣（Balanced Incomplete Factorization）、多尺度預(yù)條件矩陣等[8]。預(yù)條件矩陣避免了對(duì)線性方程組系數(shù)矩陣[9]的直接求逆，大大提高了迭代收斂速度，目前正廣泛應(yīng)用于光束法平差解算[10]。為進(jìn)一步提高平差效率，在預(yù)條件共軛梯度法的基礎(chǔ)上又引入了不精確牛頓解法（Inexact Newton method），其可以在迭代過程中改變迭代收斂條件，提前終止迭代，用不精確解代替精確解。研究證明，在光束法平差過程中，利用不精確牛頓解法提前終止迭代是可行的，且不會(huì)影響最終的平差精度。本文將通過實(shí)驗(yàn)說明上述方法的可行性和高效性。

2 理論與方法

2.1 光束法平差

光束法平差的基本測(cè)量單位是每個(gè)像點(diǎn)所對(duì)應(yīng)的光線束，其基本模型是共線方程，誤差方程如下：

式（1）中：為誤差方程中對(duì)應(yīng)外方位元素的系數(shù)矩陣；為誤差方程中對(duì)應(yīng)加密點(diǎn)地面坐標(biāo)的系數(shù)矩陣；為影像外方位的改正數(shù)向量；為加密點(diǎn)地面坐標(biāo)的改正數(shù)向量；為誤差方程常數(shù)項(xiàng)向量。

將式（1）按照最小二乘原理法化，得到相應(yīng)的法方程：

為表達(dá)簡單，可將式（2）簡寫為：

加密點(diǎn)地面坐標(biāo)改正數(shù)的個(gè)數(shù)遠(yuǎn)遠(yuǎn)大于影像外方位元素改正數(shù)的個(gè)數(shù)，因此，可以利用矩陣Schur補(bǔ)方法消去加密點(diǎn)地面坐標(biāo)改正數(shù)，得到僅含有一類未知數(shù)的改化法方程：

在實(shí)際處理中，通常利用得到的新的外方位元素進(jìn)行多片前方交會(huì)，得到加密點(diǎn)地面坐標(biāo)。

2.2 稀疏矩陣的存儲(chǔ)計(jì)算

分析式（1）和式（3），可以得到法方程系數(shù)矩陣的結(jié)構(gòu)，如圖1所示。

圖1 法方程系數(shù)矩陣結(jié)構(gòu)圖

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

共軛梯度法在求解病態(tài)方程方面具有較大的優(yōu)勢(shì)，同時(shí)對(duì)于稀疏矩陣的處理也十分有效。本文選取共軛梯度法求解大型線性方程組：

=（5）

共軛梯度法收斂速度的快慢很大程度上取決于線性方程組系數(shù)矩陣的條件數(shù)，系數(shù)矩陣的條件數(shù)越小，共軛梯度法所需的迭代次數(shù)就越少，收斂速度越快。在利用共軛梯度法進(jìn)行光束法平差處理時(shí)，由于其未知數(shù)數(shù)量多，改化法方程系數(shù)矩陣的條件數(shù)大，這會(huì)大大降低光束法平差的迭代效率。因此，為了加快迭代收斂速度，有必要尋求一種方法來降低線性方程組系數(shù)矩陣的條件數(shù)。為此，可以在式（5）兩邊同時(shí)乘一個(gè)預(yù)條件矩陣－1，得到：

－1=－1（6）

在式（6）中，系數(shù)矩陣由變成了－1，其相應(yīng)的條件數(shù)也變成了矩陣－1的條件數(shù)。通過選擇合適的預(yù)條件矩陣－1，可以有效減少線性方程組系數(shù)矩陣的條件數(shù)，達(dá)到加快迭代收斂速度的目的。預(yù)條件矩陣－1的選擇一般需要滿足以下兩個(gè)條件：與系數(shù)矩陣相似，便于構(gòu)造和求逆；可以有效降低系數(shù)矩陣的條件數(shù)。

目前常見的預(yù)條件矩陣有Jacobi預(yù)條件矩陣、SSOR預(yù)條件矩陣、基于QR分解的預(yù)條件矩陣、平衡不完全分解預(yù)條件矩陣、多尺度預(yù)條件矩陣等。其中，Jacobi預(yù)條件矩陣最便于構(gòu)造和求逆，且實(shí)踐證明其效果最為穩(wěn)定，滿足預(yù)條件矩陣選取的兩大原則。其他的預(yù)條件矩陣的效果可能會(huì)優(yōu)于Jacobi預(yù)條件，但它們的計(jì)算較為復(fù)雜且穩(wěn)定性差。因此，本文選擇Jacobi預(yù)條件矩陣對(duì)式（5）進(jìn)行預(yù)處理。

Jacobi預(yù)條件矩陣由系數(shù)矩陣A中對(duì)角線上的塊狀子矩陣構(gòu)成，在利用預(yù)條件共軛梯度法求解改化法方程時(shí)，其相應(yīng)的Jacobi預(yù)條件矩陣是由式（4）系數(shù)矩陣中對(duì)角線上的大小為6×6的子矩陣構(gòu)成，其結(jié)構(gòu)如圖2所示。

圖2 Jacobi預(yù)條件矩陣的結(jié)構(gòu)示意圖

2.4 不精確牛頓解法

預(yù)條件矩陣的引入雖然可以在一定程度上加快共軛梯度法的迭代效率，但當(dāng)其方程規(guī)模很大時(shí)，其收斂速度仍不理想。幸運(yùn)的是，經(jīng)研究發(fā)現(xiàn)，在預(yù)條件共軛梯度法的迭代過程中，通過迭代求得的未知數(shù)只有前幾次迭代會(huì)產(chǎn)生較大的改進(jìn)，后續(xù)的迭代過對(duì)未知數(shù)的改進(jìn)很小，而后續(xù)的迭代過程仍會(huì)占用大量的內(nèi)存和時(shí)間，因此，可以考慮提前終止預(yù)條件共軛梯度法迭代，只保留前幾次對(duì)未知數(shù)有較大改進(jìn)效果的迭代。

在光束法平差解算過程中，上一次共軛梯度法得到的改化法方程的解作為下一次光束法平差迭代計(jì)算的初值，因此提前終止迭代雖然只得到一組不精確解，但使用這組不精確解作為初值繼續(xù)進(jìn)行光束法平差迭代并不會(huì)對(duì)最終精度產(chǎn)生影響。所以提前終止預(yù)條件共軛梯度法的迭代，并用其不精確解代替其嚴(yán)密解繼續(xù)進(jìn)行光束法平差迭代是可行的，此時(shí)，只需要確定何時(shí)停止迭代即可。

最常用的共軛梯度法的近似解是不精確牛頓解法。不精確牛頓解法是牛頓法的一種改進(jìn)方法，其基本思想是在迭代求解改化法方程的過程中引入新的迭代收斂條件，將迭代的終止條件從待求參數(shù)殘差向量的閾值更改為強(qiáng)制序列系數(shù)（Forcing sequence）。在預(yù)條件共軛梯度法每次迭代后，需要計(jì)算當(dāng)前的強(qiáng)制序列系數(shù)，并判斷此強(qiáng)制序列系數(shù)k是否小于指定閾值，通過設(shè)置強(qiáng)制序列系數(shù)的閾值，提前終止迭代，用改化法方程的一組不精確解替代其精確解。實(shí)踐表明，不精確牛頓解法既能保證平差解算的最終精度，又能大幅度降低預(yù)條件共軛梯度法的迭代次數(shù)。

強(qiáng)制序列系數(shù)k的計(jì)算方法為：

引入不精確牛頓解后，預(yù)條件共軛梯度法的迭代退出條件為判斷是否k＞t。其中t為強(qiáng)制序列系數(shù)的閾值。關(guān)于如何確定強(qiáng)制序列系數(shù)的閾值t，已有大量的研究，一般按照經(jīng)驗(yàn)法則選擇常數(shù)t=0.1。

引入不精確牛頓解法的預(yù)條件共軛梯度法的迭代求解步驟如下。

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

為了驗(yàn)證本文所提的稀疏矩陣的分塊壓縮與計(jì)算方法能減少計(jì)算機(jī)內(nèi)存占用，以及將預(yù)條件矩陣和不精確牛頓解法引入到共軛梯度法，能夠有效提高大規(guī)模光束法區(qū)域網(wǎng)空中三角測(cè)量平差效率。選取一組包含35張像片、12 063個(gè)加密點(diǎn)，53 075個(gè)影像點(diǎn)的數(shù)據(jù)進(jìn)行實(shí)驗(yàn)。該數(shù)據(jù)共含有3個(gè)地面控制點(diǎn)，控制點(diǎn)在高斯坐標(biāo)系下的高斯平面坐標(biāo)如表1所示。

本次數(shù)據(jù)中加密點(diǎn)的個(gè)數(shù)已達(dá)到上萬個(gè)點(diǎn)位，單純應(yīng)用傳統(tǒng)的光束法區(qū)域網(wǎng)平差方法進(jìn)行計(jì)算涉及部分矩陣以達(dá)到上萬級(jí)別，且矩陣的稀疏化也較為嚴(yán)重，直接進(jìn)行數(shù)據(jù)的存儲(chǔ)無疑會(huì)給計(jì)算機(jī)的存儲(chǔ)帶來巨大的負(fù)荷。若大型稀疏矩陣中存在少數(shù)奇異值時(shí)對(duì)解算結(jié)果會(huì)產(chǎn)生較大的偏差，傳統(tǒng)的最小二乘迭代方法在處理含有奇異值的解算過程，會(huì)受到奇異值的擾動(dòng)，對(duì)解算結(jié)果的精度產(chǎn)生較大的影響，導(dǎo)致解算結(jié)果不可靠的情況發(fā)生，本文采用預(yù)條件共軛梯度法對(duì)法方程進(jìn)行求解，克服了奇異值對(duì)解算結(jié)果的影響，提高了算法的收斂性，詳細(xì)的解算結(jié)果如表2所示。

表1 控制點(diǎn)的高斯平面坐標(biāo)

控制點(diǎn)號(hào)X/mY/mZ/m GCP_0157 8241.361 02 843 930.948 5564.169 9 GCP_02578 365.495 02 843 855.874 4564.260 0 GCP_03578 319.473 12 843 766.751 7582.906 8

表2 光束法區(qū)域網(wǎng)平差解算情況

迭代次數(shù)最大內(nèi)存消耗/MB耗用時(shí)間/s 5107.64444.926

由表2可知，采用本文所述方法，對(duì)上萬地面點(diǎn)的光束法區(qū)域網(wǎng)空三解算僅占用了較小的內(nèi)存空間，迭代次數(shù)為5次，總耗時(shí)44.926 s，保證了較快的收斂速度。此數(shù)據(jù)光束法區(qū)域網(wǎng)平差檢查點(diǎn)殘差方面精度如表3所示。

表3 檢查點(diǎn)坐標(biāo)精度

檢查點(diǎn)號(hào)△X/m△Y/m△Z/m GCP_010.1430.0970.397 GCP_020.1240.1410.352

由表3可以明顯看出，對(duì)于較大規(guī)模數(shù)據(jù)采用本文方法進(jìn)行光束法區(qū)域網(wǎng)平差解算在實(shí)際中是可行的。檢查點(diǎn)的坐標(biāo)精度在，坐標(biāo)分量上的精度明顯高于坐標(biāo)分量上的精度，定位解算中普遍存在坐標(biāo)分量存在較大偏差的情況，光束法區(qū)域網(wǎng)平差中同樣存在此現(xiàn)象。由本實(shí)驗(yàn)的解算結(jié)果的數(shù)值表明，本文所提的方法在較大規(guī)模的空間加密中是可行的，克服了大型稀疏矩陣存儲(chǔ)的負(fù)荷，且精度上滿足現(xiàn)實(shí)的生產(chǎn)實(shí)際精度。

4 結(jié)論

本文針對(duì)無人機(jī)影像大數(shù)據(jù)量給傳統(tǒng)的光束法區(qū)域網(wǎng)平差解算帶來的內(nèi)存消耗大、計(jì)算效率低等問題，提出了一套快速高效的解決方案并加以實(shí)現(xiàn)。針對(duì)稀疏矩陣，利用其系數(shù)矩陣特有的塊狀結(jié)構(gòu)，使用了一種特殊的稀疏矩陣的存儲(chǔ)與計(jì)算方法。在共軛梯度法的基礎(chǔ)上又引入了預(yù)條件矩陣及不精確牛頓解法，進(jìn)一步提高了共軛梯度法求解改化法方程的平差效率，使之在保證解算精度的同時(shí)又能提高運(yùn)算速度。從本文實(shí)驗(yàn)結(jié)果中可以看出，加密點(diǎn)的平面精度良好，但高程精度相對(duì)較低，在今后的研究中，應(yīng)加強(qiáng)對(duì)光束法區(qū)域網(wǎng)平差的點(diǎn)位高程精度的研究。

［1］季順平，史云.車載全景相機(jī)的影像匹配和光束法平差［J］.測(cè)繪學(xué)報(bào)，2013，42（1）：98-104，111.

［2］閆利，費(fèi)亮，葉志云，等.大范圍傾斜多視影像連接點(diǎn)自動(dòng)提取的區(qū)域網(wǎng)平差法［J］.測(cè)繪學(xué)報(bào)，2016（3）：310-317，338.

［3］BELL N，GARLAND M.Implementing sparse matrix- vector multiplication on throughput-oriented processors［C］//Proceedings of the ACM/IEEE Conference on High Performance Computing，2009.

［4］ YI-ZHOU H E，CHEN X，WANG H.Modeling correlated samples via sparse matrix gaussian graphical models［J］.Journal of Zhejiang University Science C，2013（2）：36-46.

［5］NOCEDAL J，WRIGHT S J.Numerical optimization［M］.Berlin：Springer，2000.

［6］AGARWAL S，SNAVELY N，SEITZ S，et al.Bundle adjustment in the large［J］.In ECCV10，2010（2）：29-42.

［7］MARTIN B，KALLE ?.Conjugate gradient bundle adjustment［C］//European Conference on Computer Vision，2010.

［8］BYROD M，?STROM K.Conjugate gradient bundle adjustment［J］.Notes Comput Sci，2010（6312）： 114-127.

［9］鄭茂騰，張永軍，朱俊峰，等.一種快速有效的大數(shù)據(jù)區(qū)域網(wǎng)平差方法［J］.測(cè)繪學(xué)報(bào)，2017，46（2）：58-67.

［10］孫磊，郭海濤，李傳廣，等.大區(qū)域無人機(jī)影像數(shù)據(jù)的快速幾何處理［J］.測(cè)繪工程，2016，25（3）：38-43.

2095－6835（2020）10－0007－03

P237

A

10.15913/j.cnki.kjycx.2020.10.003

袁斌（1995—），女，碩士研究生，研究方向?yàn)閿z影測(cè)量與LiDAR數(shù)據(jù)處理。

〔編輯：嚴(yán)麗琴〕