無人機(jī)光束法空三稀疏矩陣處理與實現(xiàn)

亓晨 ，趙西安，董友強(qiáng)，陳柯

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

1 引言

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

2 無人機(jī)光束法空三平差

2.1 光束法空三平差模型

光束法空三平差的基本數(shù)學(xué)模型是共線方程:

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

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

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

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

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

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

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

3 大規(guī)模稀疏矩陣處理

3.1 大規(guī)模稀疏矩陣內(nèi)存管理

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

現(xiàn)在有一個由2條航帶，每航帶3張影像組成的區(qū)域，如圖1所示。

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

圖1 區(qū)域網(wǎng)點位分布圖

圖2 誤差方程系數(shù)矩陣結(jié)構(gòu)

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

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

對于A的轉(zhuǎn)置AT的存儲:

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

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

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

對于系數(shù)矩陣BT的存儲:

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

法方程式系數(shù)矩陣結(jié)構(gòu)如圖3所示:

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

(3)對于N11的存儲:

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

(4)對于N12的存儲:

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

等同于對N12的處理，只是基本單位是轉(zhuǎn)置關(guān)系，記為。

(6)對于N22的存儲:

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

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

3.2 大規(guī)模稀疏矩陣算法

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

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

(1)N11的解求方法:

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

(2)N12的解求方法:

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

在求得N12以后，將N12的基本單位轉(zhuǎn)置便可得到。

(4)N22的解求方法:

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

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

(5)NN的解求方法:

易知，NN是一個36×36的方陣，由36個6×6的方陣組成，將每個方陣看作一個元素，NN便簡化為一個6×6的二維數(shù)組，解求NN的任務(wù)即簡化為計算NN中的36個元素，記NN中第i行第j列的元素為NNij。首先，找到第i像片與第j像片的所有公共點，記作p1，p2，…，pn(假設(shè)有 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取反，便得到改化法方程式系數(shù)矩陣N。

3.3 矩陣運算穩(wěn)定性分析

對N22進(jìn)行結(jié)構(gòu)剖析，由圖3可知，N22是由11個三階方陣對角排列組成的稀疏矩陣，因而對其求逆等同于對各個三階方陣求逆。在數(shù)據(jù)不好的情況下，個別三階方陣的對應(yīng)的行列式值幾乎為零，接近奇異，導(dǎo)致無法求逆或求逆結(jié)果非常不準(zhǔn)確［10］，N22就會呈現(xiàn)病態(tài)，導(dǎo)致平差精度下降，結(jié)果不穩(wěn)定。

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

4 無人機(jī)光束法空三平差實現(xiàn)

4.1 內(nèi)存管理實現(xiàn)

內(nèi)存管理的實現(xiàn)，主要體現(xiàn)在方程系數(shù)矩陣的VC++內(nèi)存設(shè)置管理。3.1中所述各矩陣VC++內(nèi)存管理方法相同，其實現(xiàn)程序如下。

稀疏矩陣內(nèi)存管理實現(xiàn)程序:

4.2 矩陣運算穩(wěn)定性實現(xiàn)

以下為矩陣運算穩(wěn)定性實現(xiàn)方法:

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

5 實驗與分析

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

數(shù)據(jù)采集基本信息 表1

本實驗數(shù)據(jù)含有12個地面控制點，控制點在高斯坐標(biāo)系下的高斯平面坐標(biāo)如表2所示:

控制點的高斯平面坐標(biāo) 表2

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

光束法空三平差數(shù)據(jù)及解算 表3

由表3可知，實驗數(shù)據(jù)量巨大，如果在光束法平差過程中不對稀疏矩陣進(jìn)行處理，是不可能完成解算的。而采用本文所述算法，不僅使解算變?yōu)榭赡?，還保證了較快的收斂速度及較高的穩(wěn)定性，節(jié)省了計算機(jī)內(nèi)存。

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

控制點精度 表4

檢查點精度 表5

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

6 結(jié)論與展望

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

［1］魯恒，李永樹，何敬等．無人機(jī)低空遙感影像數(shù)據(jù)的獲取與處理［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］吳杰，李明峰，余騰．測量數(shù)據(jù)處理中病態(tài)矩陣和正則化方法［J］．大地測量與地球動力學(xué)，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］ 李德仁，鄭肇葆．解析攝影測量學(xué)［M］．北京:測繪出版社，1992:161～166．

［10］ 同濟(jì)大學(xué)應(yīng)用數(shù)學(xué)系．線性代數(shù)［M］．北京:高等教育出版社，2003:43．