無人機光束法空三稀疏矩陣處理與實現

亓晨 ，趙西安，董友強，陳柯

(1.北京建筑大學測繪與城市空間信息學院，北京 100044;2.北京靈圖軟件技術有限公司，北京 100193)

1 引言

無人機航測具有高機動性、靈活性、使用成本低等特點，在航測領域得到越來越廣泛的應用［1～3］。然而無人機航測與傳統航測相比具有以下不足:無人機影像像幅小，數量多，信息采集量大，后期在光束法空三平差中需要處理大量的數據，致使誤差方程系數矩陣規模龐大，可達幾萬階或幾十萬階;無人機飛行姿態不穩定，傾角過大，影像重疊部分不規則，畸變明顯，嚴重影響數據質量，致使光束法空三平差中出現接近奇異的病態方陣，此類矩陣無法求逆或求逆結果非常不準確，影響平差結果的穩定性。因此，如何高效利用有限的計算機內存完成如此大規模稀疏矩陣的存儲和運算以及對病態矩陣的避免，保證運算的穩定性，值得我們探討和研究。對稀疏矩陣的壓縮存儲有多種方案，較常用的有對角線存儲法，對稱矩陣的變帶寬存儲法，CSR存儲法，坐標存儲法，Ellpack－Itpack存儲法等［4～6］。經研究發現，上述技術雖然能夠解決多數大規模稀疏矩陣的存儲和解算問題，卻不能針對光束法平差系數矩陣的獨特結構直接使用。提高光束法空三平差穩定性的措施，多數集中在對病態方程的解算方面，如嶺估計法、截斷奇異值法和正則化方法等［7～8］，卻忽略了對平差前期對病態矩陣的避免。對于無人機光束法空三平差，內存管理、矩陣解算技術以及對病態矩陣的避免，對于實現其解的高效性和穩定性都極為重要。

2 無人機光束法空三平差

2.1 光束法空三平差模型

光束法空三平差的基本數學模型是共線方程:

對上式進行線性化，可得到誤差方程式，寫成矩陣的形式為:

在式(3)這類誤差方程式中含有兩類未知數:外方位元素t和加密點的地面攝影測量坐標X。相應的法方程式為:

消去未知數X后，可得改化法方程式:

對改化法方程式(6)采用循環分塊法求出全區每片的外方位元素，然后采用前方交會法求出加密點的地面坐標［9］。

2.2 光束法空三平差存在問題

光束法空三平差具有嚴密的理論基礎，是目前空三解算最常用的方法。但無人機航測的大數據量會導致誤差方程系數矩陣規模龐大，且矩陣規模隨著測區大小和像點數量的增加而增加。如果一個區域網由7航帶，每航帶50張影像組成，每像對匹配得到300對同名像點，地面點數量為 40 000左右，那么誤差方程式的數量至少為300×50×7×2=210 000，相應的誤差方程系數矩陣則是一個數十萬階的大規模稀疏矩陣。對這樣一個矩陣的存儲，計算機內存顯然是不夠的，需要我們針對矩陣的特點，進行有序分塊存儲。

光束法空三平差需要通過參數消去法將法方程式(5)轉換為改化法方程式(6)，這個過程中需要對N22求逆，得到。而在無人機航測數據質量不好的情況下，方陣N22的行列式值幾乎為零，接近奇異，呈現病態，無法求逆或求逆結果非常不準確［10］，影響解算的精度和穩定性。因此需要對N22的結構深入剖析，從而采取措施避免病態。

3 大規模稀疏矩陣處理

3.1 大規模稀疏矩陣內存管理

本文根據無人機光束法空三平差中遇到的大規模稀疏矩陣的特點，設計了一種內存管理方法。下面用一個實例來詳細介紹這種內存管理方法。

現在有一個由2條航帶，每航帶3張影像組成的區域，如圖1所示。

首先，將所有像片進行編號，除第1航帶，每航帶接由上一航帶像片編號依次增加，就本例而言，所有像片依次編號為 1，2，3，4，5，6。依次按照像片順序，每像片所包含的像點順序，列出誤差方程式，該誤差方程式系數矩陣結構如圖2所示:

圖1 區域網點位分布圖

圖2 誤差方程系數矩陣結構

(1)對于系數矩陣A的存儲:

將A定義為一個4維數組，首先以像片為單位將A 分為 A1，A2，A3，A4，A5，A6。假設第 i像片中有 n 個像點，則將Ai按像點點號大小依次分為Ai1，…，Ain，則Aik(1≤k≤m)表示第i像片上的第k個像點所列誤差方程系數矩陣中的A矩陣。

對于A的轉置AT的存儲:

此外，對于AT的存儲，還要考慮光束法空三平差的特點以及法方程式(5)中N12的求解。光束法空三平差中的未知參數分為像片的外方位元素與加密點坐標兩類，但不包含控制點坐標，對控制點所列誤差方程式系數矩陣中的B為零矩陣，致使N12=AT·B的求解過程中與控制點對應的N12全為零矩陣，所以可以認為與控制點對應的AT并未參與N12的解算。因此為了N12的解算方便，需對AT進行另外一種方式的存儲，使得AT中不包含控制點對應的AT。將按該方式存儲的AT記為(A_E)T，下面給出(A_E)T的存儲方式:

(2)對于系數矩陣B的存儲:

針對上文提到的光束法空三平差的特點，為了解算方便，對B進行分塊存儲時，也不再考慮控制點。將B定義為一個四維數組，首先以像片為單位將B分為B1，B2，B3，B4，B5，B6。假設第 i像片有 m 個像點(非控制點)，則將Bi按像點點號大小依次分為Bi1，…，Bim，Bik(1≤k≤m)則代表第i像片中第k個像點(非控制點)所列誤差方程系數矩陣中的B矩陣。

對于系數矩陣BT的存儲:

將BT定義為一個四維數組，等同于對B的處理，假設第i像片有m個像點(非控制點)，則將按像點點號大小依次分為(1≤k≤m)則代表第i像片中第k個像點所列誤差方程系數矩陣中B矩陣的轉置BT。

法方程式系數矩陣結構如圖3所示:

圖3 法方程式系數矩陣結構

(3)對于N11的存儲:

由圖可知，N11是由6個6×6的方陣沿對角線排列組成的對角方陣(圖中左上角6個黑色方塊)。將N11定義為一個三維數組，依次按像片順序將N11分為N(11)1，N(11)2，N(11)3，N(11)4，N(11)5，N(11)6，分別對應6 個方陣。

(4)對于N12的存儲:

由圖可知，N12的基本單位是6×3的小矩陣。將N12定義為一個四維數組，按像片順序將N12依次分為N(12)1，N(12)2，N(12)3，N(12)4，N(12)5，N(12)6。假設第 i像片有m個像點(非控制點)，則按照該像片像點序號依次分為 N(12)i1，…，N(12)im，N(12)ik(1≤k≤m)表示第 i像片中第k個像點(非控制點)對應的6×3的矩陣ATB。

等同于對N12的處理，只是基本單位是轉置關系，記為。

(6)對于N22的存儲:

由圖可知，N22是由若干3×3的方陣對角排列組成的對角方陣，每個小方陣與區域網中每個加密點是一一對應的關系。將N22定義為一個三維數組，數組中每個元素表示區域網中每個加密點對應的3×3方陣，假設某點是區域網所有待定點中第k個點，則該點對應的小方陣記為N(22)k。

由法方程式(5)可得改化法方程式(6)，改化法方程式的系數矩陣規模有限，且通常情況下，非零元素較少，故不屬于大規模稀疏矩陣，所以這里不需討論。

3.2 大規模稀疏矩陣算法

對于大規模稀疏矩陣已經實現分塊存儲，對于方程式的解算，普通的矩陣運算方法已不適合，需要設計一套適用于上述矩陣分塊存儲方式的矩陣運算方法。這里重點討論兩個過程:由誤差方程系數矩陣得到法方程系數矩陣的過程以及由法方程系數矩陣得到改化法方程式系數矩陣的過程。

首先給出由誤差方程系數矩陣A和B得到法方程系數矩陣 N11，N12，，N22的矩陣運算方法。

(1)N11的解求方法:

由N11的存儲方法可知，N11的基本單位是對應第i像片的6×6的方陣N(11)i，所以這里只需討論對于N(11)i的解算方法。假設第i像片有n個像點，則:

(2)N12的解求方法:

由N12的存儲方法可知，N12的基本單位是對應第i像片中第k個像點(非控制點)的6×3的矩陣N(12)ik，所以這里只需討論對于N(12)ik的解算方法。假設第i像片中某點是第k個像點(非控制點)，則:

在求得N12以后，將N12的基本單位轉置便可得到。

(4)N22的解求方法:

由N22的存儲方法可知，N22的基本單位是與區域網中第k個加密點(非控制點)對應的3×3的方陣N(22)k，所以這里只需討論N(22)k的解算方法。首先，需要通過遍歷像點數據得到包含該像點(非控制點)的像片數n以及每張像片中該像點對應的誤差方程的系數矩陣B，依次記為B1，B2，…，Bn，則

在法方程系數矩陣全部解求完畢之后，需要由法方程系數矩陣解求改化法方程式的系數矩陣N11－N12，為表達簡便，暫且將 N12簡記為 NN，將 N11－N12簡記為 N。

(5)NN的解求方法:

易知，NN是一個36×36的方陣，由36個6×6的方陣組成，將每個方陣看作一個元素，NN便簡化為一個6×6的二維數組，解求NN的任務即簡化為計算NN中的36個元素，記NN中第i行第j列的元素為NNij。首先，找到第i像片與第j像片的所有公共點，記作p1，p2，…，pn(假設有 n 個公共點)，pk(1≤k≤n)在第 i像片所有非控制點中的序號記作pki，在第j像片所有非控制點中的序號記作pkj，則:

(6)N的解求方法:

將NN中對角線上6×6的方陣NN11，…，NN66依次作 NNii－N(11)i(i=1，2，…，6)的處理，再將此時的 NN取反，便得到改化法方程式系數矩陣N。

3.3 矩陣運算穩定性分析

對N22進行結構剖析，由圖3可知，N22是由11個三階方陣對角排列組成的稀疏矩陣，因而對其求逆等同于對各個三階方陣求逆。在數據不好的情況下，個別三階方陣的對應的行列式值幾乎為零，接近奇異，導致無法求逆或求逆結果非常不準確［10］，N22就會呈現病態，導致平差精度下降，結果不穩定。

由3.1中所提N22的存儲方式可知，N22中每個三階方陣與區域網中每個加密點是一一對應的關系。為了避免N22呈現病態，只需將其中病態的三階方陣所對應的加密點剔除即可。采用的方法是，在列出誤差方程式之前，求出每個加密點在N22中對應三階方陣的行列式值，將行列式值的絕對值小于0.1的三階方陣對應的待定點直接由原始像點數據中剔除。如此一來，避免了方陣N22會出現病態的情況，保證了平差解算的精度和穩定性。

4 無人機光束法空三平差實現

4.1 內存管理實現

內存管理的實現，主要體現在方程系數矩陣的VC++內存設置管理。3.1中所述各矩陣VC++內存管理方法相同，其實現程序如下。

稀疏矩陣內存管理實現程序:

4.2 矩陣運算穩定性實現

以下為矩陣運算穩定性實現方法:

說明:①PhotoN為所有航帶像片的總數;②PointN是一維數組，PointN［i］表示第i張像片包含的像點數;③ExPointN是一維數組，ExPointN［i］表示第i張像片包含的非控制點像點數;④SPointN是整個區域網中加密點的個數。⑤ClmofB是一個二維數組，ClmofB［i］［j］表示第i像片中第j個非控制點像點對應的加密點在區域網所有加密點中的序號。⑥Calculating()是一個用于求解N22的函數。⑦valueHLS()是一個用于求解N22行列式值的函數。⑧SickPoint是一個一維數組，存儲需要剔除的加密點的序號。⑨PointData是一個vector數組，用于存儲每個像點的像平面坐標x、y和對應的地面攝影測量坐標 X、Y、Z，如 PointData［i］［j］．x表示第i像片第j個像點的像平面坐標x。

5 實驗與分析

本次無人機實驗數據采集區域為新疆沙爾拜客水電站壩區，地形平坦，以農田為主，包含河流及少量民房。無人機所攜帶數碼相機型號為JC4，焦距 24 mm，影像分辨率為 3888 xs×2592 xs，CCD尺寸 22.2 mm×14.8 mm。測區含有7個航帶，每航帶像片數在45～50之間，數據采集基本信息如表1所示:

數據采集基本信息 表1

本實驗數據含有12個地面控制點，控制點在高斯坐標系下的高斯平面坐標如表2所示:

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

為了檢查解算精度，將其中CP_03，CP_06，CP_09，CP_11這4個控制點作為檢查點。光束法空三平差解算的數據及解算情況如表3所示:

光束法空三平差數據及解算 表3

由表3可知，實驗數據量巨大，如果在光束法平差過程中不對稀疏矩陣進行處理，是不可能完成解算的。而采用本文所述算法，不僅使解算變為可能，還保證了較快的收斂速度及較高的穩定性，節省了計算機內存。

表4是本次光束法空三平差實驗的控制點精度報告，表5是檢查點的精度報告:

控制點精度 表4

檢查點精度 表5

由以上兩表明顯可以看出，采用本文所述算法進行無人機光束法空三平差，精度遠遠高于 1∶2 000的測圖要求。

6 結論與展望

本文提出了針對無人機光束法空三大規模稀疏矩陣的一種內存管理與解算方法并加以實現。通過對稀疏矩陣進行分塊存儲與運算以及部分加密點的剔除，達到了只對稀疏矩陣的非零元素進行存儲的目的。有效降低了無人機光束法空三解算過程中對計算機內存的消耗，提高了運算效率，減少了運算時間，同時避免了病態矩陣的出現，保證了解算結果的穩定性。最后，通過對無人機實驗數據進行測試，驗證了該方法的可行性與高效性，平差的精度也符合要求。但由精度報告也可看到，高程精度總體要低于平面精度，同時穩定性也較差。今后將會就如何改善無人機光束法空三的高程精度及穩定性作相應研究并加以實現。

［1］魯恒，李永樹，何敬等．無人機低空遙感影像數據的獲取與處理［J］．測繪工程，2011(1):51～54．

［2］Laliberte A S，Rango A．Texture and Scale in Object－based Analysis of Subdecimeter Resolution Unmanned Aerial Vehicle(UAV)Imagery［J］．Geoscience and Remote Sensing，2009，761 ～769．

［3］Laliberte A S，Herrick J E，Rango A，et al．Acquisition，or Thorectification，and Object－ based Classification of Unmanned Aerial Vehicle(UAV)Imagery for Rangeland Monitoring［J］．Photogrammetric Engineering ＆ Remote Sensing，2010，661 ～672．

［4］Ichitaro Yamazaki，Richard，Scalettar．A High－Quality Preconditioning Technique for Multi－Length－Scale Symmetric Positive Definite Linear Systems［J］．Numerical Mathematics:Theory，Methods and Applications，2009，469 ～484．

［5］A sparse matrix model－based optical proximity correction algorithm with model－based mapping between segments and control sites［J］．Journal of Zhejiang University－Science C(Computers＆ Electronics)，2011，436～442．

［6］Yi－zhou HE，Xi CHEN，Hao Wang．Modeling correlated samples via sparse matrix Gaussian graphical models［J］．Journal of Zhejiang University－Science C(Computers＆ E-lectronics)，2013，107 ～117．

［7］吳杰，李明峰，余騰．測量數據處理中病態矩陣和正則化方法［J］．大地測量與地球動力學，2010(04):102～105．

［8］N．S．Mera，L．Elliott，D．B．Ingham．On the use of genetic algorithms for solving ill－posed problems［J］．Inverse Problems in Engineering，2003，105 ～121．

［9］ 李德仁，鄭肇葆．解析攝影測量學［M］．北京:測繪出版社，1992:161～166．

［10］ 同濟大學應用數學系．線性代數［M］．北京:高等教育出版社，2003:43．