無人機影像光束法平差研究

袁斌

無人機影像光束法平差研究

袁斌

（山東理工大學，山東 淄博 255000）

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

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

1 引言

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

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

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

2 理論與方法

2.1 光束法平差

光束法平差的基本測量單位是每個像點所對應的光線束，其基本模型是共線方程，誤差方程如下：

式（1）中：為誤差方程中對應外方位元素的系數矩陣；為誤差方程中對應加密點地面坐標的系數矩陣；為影像外方位的改正數向量；為加密點地面坐標的改正數向量；為誤差方程常數項向量。

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

為表達簡單，可將式（2）簡寫為：

加密點地面坐標改正數的個數遠遠大于影像外方位元素改正數的個數，因此，可以利用矩陣Schur補方法消去加密點地面坐標改正數，得到僅含有一類未知數的改化法方程：

在實際處理中，通常利用得到的新的外方位元素進行多片前方交會，得到加密點地面坐標。

2.2 稀疏矩陣的存儲計算

分析式（1）和式（3），可以得到法方程系數矩陣的結構，如圖1所示。

圖1 法方程系數矩陣結構圖

2.3 預條件共軛梯度法

共軛梯度法在求解病態方程方面具有較大的優勢，同時對于稀疏矩陣的處理也十分有效。本文選取共軛梯度法求解大型線性方程組：

=（5）

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

－1=－1（6）

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

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

Jacobi預條件矩陣由系數矩陣A中對角線上的塊狀子矩陣構成，在利用預條件共軛梯度法求解改化法方程時，其相應的Jacobi預條件矩陣是由式（4）系數矩陣中對角線上的大小為6×6的子矩陣構成，其結構如圖2所示。

圖2 Jacobi預條件矩陣的結構示意圖

2.4 不精確牛頓解法

預條件矩陣的引入雖然可以在一定程度上加快共軛梯度法的迭代效率，但當其方程規模很大時，其收斂速度仍不理想。幸運的是，經研究發現，在預條件共軛梯度法的迭代過程中，通過迭代求得的未知數只有前幾次迭代會產生較大的改進，后續的迭代過對未知數的改進很小，而后續的迭代過程仍會占用大量的內存和時間，因此，可以考慮提前終止預條件共軛梯度法迭代，只保留前幾次對未知數有較大改進效果的迭代。

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

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

強制序列系數k的計算方法為：

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

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

3 實驗與分析

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

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

表1 控制點的高斯平面坐標

控制點號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 光束法區域網平差解算情況

迭代次數最大內存消耗/MB耗用時間/s 5107.64444.926

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

表3 檢查點坐標精度

檢查點號△X/m△Y/m△Z/m GCP_010.1430.0970.397 GCP_020.1240.1410.352

由表3可以明顯看出，對于較大規模數據采用本文方法進行光束法區域網平差解算在實際中是可行的。檢查點的坐標精度在，坐標分量上的精度明顯高于坐標分量上的精度，定位解算中普遍存在坐標分量存在較大偏差的情況，光束法區域網平差中同樣存在此現象。由本實驗的解算結果的數值表明，本文所提的方法在較大規模的空間加密中是可行的，克服了大型稀疏矩陣存儲的負荷，且精度上滿足現實的生產實際精度。

4 結論

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

［1］季順平，史云.車載全景相機的影像匹配和光束法平差［J］.測繪學報，2013，42（1）：98-104，111.

［2］閆利，費亮，葉志云，等.大范圍傾斜多視影像連接點自動提取的區域網平差法［J］.測繪學報，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］鄭茂騰，張永軍，朱俊峰，等.一種快速有效的大數據區域網平差方法［J］.測繪學報，2017，46（2）：58-67.

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

2095－6835（2020）10－0007－03

P237

A

10.15913/j.cnki.kjycx.2020.10.003

袁斌（1995—），女，碩士研究生，研究方向為攝影測量與LiDAR數據處理。

〔編輯：嚴麗琴〕