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

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

1 引言

無人機(jī)航測(cè)具有高機(jī)動(dòng)性、靈活性、使用成本低等特點(diǎn)，在航測(cè)領(lǐng)域得到越來越廣泛的應(yīng)用［1～3］。然而無人機(jī)航測(cè)與傳統(tǒng)航測(cè)相比具有以下不足:無人機(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ì)算機(jī)內(nèi)存完成如此大規(guī)模稀疏矩陣的存儲(chǔ)和運(yùn)算以及對(duì)病態(tài)矩陣的避免，保證運(yùn)算的穩(wěn)定性，值得我們探討和研究。對(duì)稀疏矩陣的壓縮存儲(chǔ)有多種方案，較常用的有對(duì)角線存儲(chǔ)法，對(duì)稱矩陣的變帶寬存儲(chǔ)法，CSR存儲(chǔ)法，坐標(biāo)存儲(chǔ)法，Ellpack－Itpack存儲(chǔ)法等［4～6］。經(jīng)研究發(fā)現(xiàn)，上述技術(shù)雖然能夠解決多數(shù)大規(guī)模稀疏矩陣的存儲(chǔ)和解算問題，卻不能針對(duì)光束法平差系數(shù)矩陣的獨(dú)特結(jié)構(gòu)直接使用。提高光束法空三平差穩(wěn)定性的措施，多數(shù)集中在對(duì)病態(tài)方程的解算方面，如嶺估計(jì)法、截?cái)嗥娈愔捣ê驼齽t化方法等［7～8］，卻忽略了對(duì)平差前期對(duì)病態(tài)矩陣的避免。對(duì)于無人機(jī)光束法空三平差，內(nèi)存管理、矩陣解算技術(shù)以及對(duì)病態(tài)矩陣的避免，對(duì)于實(shí)現(xiàn)其解的高效性和穩(wěn)定性都極為重要。

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

2.1 光束法空三平差模型

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

對(duì)于A的轉(zhuǎn)置AT的存儲(chǔ):

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

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

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

對(duì)于系數(shù)矩陣BT的存儲(chǔ):

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

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

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

(3)對(duì)于N11的存儲(chǔ):

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

(4)對(duì)于N12的存儲(chǔ):

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

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

(6)對(duì)于N22的存儲(chǔ):

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

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

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

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

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

(1)N11的解求方法:

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

(2)N12的解求方法:

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

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

(4)N22的解求方法:

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

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

(5)NN的解求方法:

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

(6)N的解求方法:

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

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

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

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

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

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

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

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

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

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

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

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

本次無人機(jī)實(shí)驗(yàn)數(shù)據(jù)采集區(qū)域?yàn)樾陆碃柊菘退娬緣螀^(qū)，地形平坦，以農(nóng)田為主，包含河流及少量民房。無人機(jī)所攜帶數(shù)碼相機(jī)型號(hào)為JC4，焦距 24 mm，影像分辨率為 3888 xs×2592 xs，CCD尺寸 22.2 mm×14.8 mm。測(cè)區(qū)含有7個(gè)航帶，每航帶像片數(shù)在45～50之間，數(shù)據(jù)采集基本信息如表1所示:

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

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

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

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

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

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

表4是本次光束法空三平差實(shí)驗(yàn)的控制點(diǎn)精度報(bào)告，表5是檢查點(diǎn)的精度報(bào)告:

控制點(diǎn)精度 表4

檢查點(diǎn)精度 表5

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

6 結(jié)論與展望

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

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

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