區(qū)域網(wǎng)空中三角測量數(shù)據(jù)平差計(jì)算效率提升策略

胡海彥，余旭初，陳　虹，江振治，方　勇

(1. 信息工程大學(xué)地理空間信息學(xué)院,河南 鄭州 450052; 2. 地理信息工程國家重點(diǎn)實(shí)驗(yàn)室,

陜西 西安 710054; 3. 西安測繪研究所,陜西 西安 710054)

HU Haiyan,YU Xuchu,CHEN Hong,JIANG Zhenzhi,FANG Yong

區(qū)域網(wǎng)空中三角測量數(shù)據(jù)平差計(jì)算效率提升策略

胡海彥1,2,3，余旭初1，陳虹2,3，江振治2,3，方勇2,3

(1. 信息工程大學(xué)地理空間信息學(xué)院,河南 鄭州 450052; 2. 地理信息工程國家重點(diǎn)實(shí)驗(yàn)室,

陜西 西安 710054; 3. 西安測繪研究所,陜西 西安 710054)

An Efficient Computational Strategy for Aero-triangulation of Block Bundle Adjustment

HU Haiyan,YU Xuchu,CHEN Hong,JIANG Zhenzhi,FANG Yong

摘要：針對光束法空中三角測量區(qū)域網(wǎng)平差數(shù)值計(jì)算中計(jì)算量占比很大的法化方程組系數(shù)矩陣，提出了一種攝影區(qū)域影像片號編排下的法化方程組系數(shù)矩陣降維精簡策略，可以有效降低數(shù)據(jù)處理的時(shí)間復(fù)雜度和空間復(fù)雜度，從而明顯提升平差計(jì)算效率。利用實(shí)際作業(yè)數(shù)據(jù)進(jìn)行了計(jì)算效率試驗(yàn)，與常規(guī)空中三角測量區(qū)域網(wǎng)數(shù)據(jù)平差相比，計(jì)算效率提高了近11倍。

引文格式： 胡海彥，余旭初，陳虹,等. 區(qū)域網(wǎng)空中三角測量數(shù)據(jù)平差計(jì)算效率提升策略[J].測繪通報(bào)，2015(9)：75-78.DOI:10.13474/j.cnki.11-2246.2015.0284

關(guān)鍵詞：光束法平差；空三；線性代數(shù)；計(jì)算效率

中圖分類號：P231

文獻(xiàn)標(biāo)識碼：B

文章編號：0494-0911(2015)09-0075-04

收稿日期：2014-09-09

作者簡介：胡海彥(1977—)，男，博士生，研究方向?yàn)楹娇蘸教旄呔榷ㄎ焕碚撆c方法。E-mail:173412586@qq.com

一、引言

光束法空中三角測量(簡稱空三)區(qū)域網(wǎng)平差計(jì)算屬于計(jì)算密集型數(shù)據(jù)處理，如果不配合良好的算法設(shè)計(jì)和計(jì)算策略，僅靠硬件性能的提升難以滿足日益增長的計(jì)算需求，在大范圍區(qū)域攝影和高冗余(多基線)攝影條件下尤其如此[1]。雖然硬件性能的提升基本遵循摩爾定律，但即使性能再好的硬件，如果運(yùn)行了指數(shù)級的算法(算法的時(shí)間復(fù)雜度按級別可大致分為對數(shù)級、線性級、冪數(shù)級及指數(shù)級)，當(dāng)問題規(guī)模較大時(shí)，運(yùn)行效率也是十分低下的，甚至出現(xiàn)計(jì)算崩潰的現(xiàn)象[2-3]。

從解析攝影測量誕生時(shí)起，圍繞區(qū)域網(wǎng)空三數(shù)據(jù)平差處理中的計(jì)算策略或軟件算法優(yōu)化等方面的研究就從未停止過，有大量文獻(xiàn)和教材可供查閱[4-5]。這些研究有的圍繞獨(dú)立模型法或航帶法展開具體研究，其本質(zhì)上是攝影區(qū)域的分區(qū)分塊處理，一定程度上影響了區(qū)域平差的整體性，導(dǎo)致了一定的誤差傳遞與累積[6-7]；有的就是進(jìn)行“平-高”交替分求，并未有效降低問題規(guī)模，因此并不是解決計(jì)算效率提升的關(guān)鍵所在[8-10]。本文針對光束法整體平差求解中影響計(jì)算效率的關(guān)鍵點(diǎn)——法化方程系數(shù)矩陣，提出一種法化系數(shù)矩陣分塊降維精簡及攝影區(qū)域影像片號編排影響下的矩陣再精簡策略，希望在解決影響光束法區(qū)域網(wǎng)平差計(jì)算效率的關(guān)鍵問題上有所突破。

二、平差計(jì)算策略

1. 分塊精簡法化系數(shù)矩陣

在量取同名像點(diǎn)坐標(biāo)及控制點(diǎn)對應(yīng)像點(diǎn)坐標(biāo)后，空三數(shù)據(jù)平差的主要技術(shù)環(huán)節(jié)是組建觀測方程組，實(shí)現(xiàn)途徑是對共線條件方程線性化后依據(jù)最小二乘原理得到法化方程組，它是關(guān)于未知數(shù)改正量的線性方程組

NΔ=K

(1)

參考文獻(xiàn)式中，N為法化后方程組系數(shù)矩陣；Δ為未知數(shù)改正向量；K為常數(shù)向量，矩陣中各元素的具體描述可見[11]。高效求解該線性系統(tǒng)的關(guān)鍵是法化系數(shù)矩陣N的結(jié)構(gòu)分析與研究[12-13]。對于區(qū)域網(wǎng)平差而言，法化系數(shù)矩陣N往往特別巨大，基本策略是進(jìn)行矩陣分塊，以達(dá)到矩陣降維，從而減少計(jì)算時(shí)間及內(nèi)存消耗的目的。

式(1)可根據(jù)外方位元素和物空間坐標(biāo)兩類未知數(shù)進(jìn)行分塊

(2)

塊對角矩陣在對角線方向上包含非零子矩陣，其他元素均為0，這種矩陣的最大特點(diǎn)是其逆矩陣也為塊對角矩陣，且子矩陣為原矩陣對應(yīng)子矩陣的逆[14]。因此，塊對角矩陣的求逆相對其他常規(guī)非零矩陣就容易得多。基于此，可對式(2)進(jìn)行重組，以便提高求解效率。首先將式(2)拆解為兩組矩陣式

(3)

(4)

(5)

再將式(5)代入式(3)中可得

(6)

即

(7)

類似矩陣分塊的方法同樣可用于計(jì)算協(xié)方差矩陣。可將N-1分塊，記為

(8)

由于CN=NC=I，因此式(8)中關(guān)于矩陣C的各分塊矩陣可按式(9)—式(11)求得

(9)

(10)

(11)

2. 像片編排次序優(yōu)化

圖1為一個(gè)攝影區(qū)域內(nèi)3條航線(每條航線9張像片)規(guī)則攝影下的影像分布情況，航向、旁向重疊率分別為60%和30%。圖中僅示意了1、2條航線的前3張像片地面覆蓋狀況，虛線部分為像對所對應(yīng)規(guī)則立體模型排列示意。具有代表性的連接點(diǎn)A映射在像片1—1、1—2、1—3、2—1、2—2、2—3 6張像片上，這6張像片能夠組成的可能立體像對都被A點(diǎn)“鏈接”了。整個(gè)攝影區(qū)域的影像鏈接(重疊)關(guān)系如圖2所示，示意了整個(gè)攝影區(qū)域所有連接點(diǎn)形成的拓?fù)溥B接關(guān)系“圖”(直線段或弧線段表示)。

這些影像鏈接(重疊)關(guān)系即確定了精簡法化方程系數(shù)矩陣中相應(yīng)位置的非零子矩陣。至于非零子矩陣在精簡法化方程系數(shù)矩陣中的具體矩陣位置則是由像片號的排序方式所決定的。一般有兩種排序策略—沿航線排序、跨航線排序，見表1。沿航線像片排序指先依航線序，再依像片序進(jìn)行整個(gè)攝影區(qū)域的像片排序；而跨航線像片排序指先依像片序，再依航線序進(jìn)行整個(gè)攝影區(qū)域的像片排序。

圖1　每航線9張像片的3條航線分布

圖2　關(guān)于連接點(diǎn)的像片間重疊拓?fù)潢P(guān)系

表1　沿/跨航線圖1攝影區(qū)域影像排序

圖3(a)為沿航線像片排序策略下的精簡法化方程系數(shù)矩陣結(jié)構(gòu)，非零元素趨向?qū)蔷€形成帶狀排列。帶寬為對角線到最遠(yuǎn)非零元素距離，所示矩陣的帶寬為6×12=72。圖3(b)為跨航線像片排序策略下的精簡法化方程系數(shù)矩陣結(jié)構(gòu)，其帶寬為6×8=48，相比圖3(a)，其帶寬大大地減小了，而這對矩陣存儲(chǔ)及解算耗時(shí)都大有益處。

根據(jù)矩陣運(yùn)算操作次數(shù)分析可知，對于非帶狀的精簡法化方程系數(shù)矩陣，線性系統(tǒng)的解算時(shí)間正比于未知數(shù)個(gè)數(shù)(6m)的3次方。例如，此處的27張航片的情況為解算時(shí)間正比于(6×27)3=4.2×106。對于帶狀矩陣，解算時(shí)間正比于帶寬平方乘以未知數(shù)個(gè)數(shù)。例如，圖3(a)對應(yīng)的解算時(shí)間正比于722×(6×27)=8.4×105，其比非帶狀矩陣情況解算時(shí)間快了5倍。而圖3(b)對應(yīng)的解算時(shí)間正比于482×(6×27)=3.7×105，其比非帶狀矩陣情況解算時(shí)間快11倍多。

三、試驗(yàn)及分析

有關(guān)試驗(yàn)數(shù)據(jù)的測區(qū)情況與區(qū)域網(wǎng)基本參數(shù)見表2，每條航線18張影像，野外測量點(diǎn)69個(gè)(29個(gè)用于控制，40個(gè)用于精度檢查)，采用Geolord-AT商用平差軟件進(jìn)行作業(yè)計(jì)算，結(jié)果見表3。

表2　區(qū)域網(wǎng)基本參數(shù)

圖3　沿航線、跨航線像片排序策略下的精簡法化方程系數(shù)矩陣結(jié)構(gòu)

耗時(shí)/sσ0/μmX/mY/mZ/m1328.770.1090.0970.4248.927.800.0810.0920.383

為了驗(yàn)證本文提出的跨航線影像片號編排下的精簡法化系數(shù)矩陣空三計(jì)算效率提升效果，采用以下4種試驗(yàn)方案進(jìn)行計(jì)算效率比對，見表4。

方案1：原始法化方程組(式(1))的直接解算(片號隨機(jī)編排下的暴力解)。

方案2：未考慮影像片號編排策略的非帶狀對角精簡法化方程組解算(片號編排隨機(jī))，這是對式(7)、式(5)的直接應(yīng)用。

方案3：沿航線影像片號編排下的非帶狀對角精簡法化方程組解算，這是依據(jù)圖3(a)對式(7)、式(5)的間接應(yīng)用。

方案4：跨航線影像片號編排下的非帶狀對角精簡法化方程組解算，這是依據(jù)圖3(b)對式(7)、式(5)的間接應(yīng)用。

根據(jù)上述4種方案，采用CLAPACK線性代數(shù)數(shù)值計(jì)算軟件包，分別實(shí)現(xiàn)了各方案對應(yīng)的4個(gè)軟件模塊用于數(shù)值計(jì)算效率評估，該軟件包提供了可靠的矩陣基本操作和運(yùn)算，并針對稀疏矩陣和帶狀塊對角矩陣提供了專門的API函數(shù)調(diào)用。基礎(chǔ)計(jì)算平臺軟硬件為：32位微軟XP操作系統(tǒng)、2.4 GHz Pentium(R)4CPU、1GB RAM。

實(shí)際作業(yè)數(shù)據(jù)的規(guī)模為6m+3n=6×108+3×2140=7068。為了分析計(jì)算效率隨問題規(guī)模大小變化的規(guī)律，在此基礎(chǔ)上抽稀加密點(diǎn)個(gè)數(shù)為1070個(gè)，則問題規(guī)模減小為6×108+3×1070=3858；再利用ERDAS遙感軟件以同樣的匹配精度(1/3像元)再匹配加密，使得加密點(diǎn)個(gè)數(shù)增加到3985個(gè)，則問題規(guī)模增加到6×108+3×3985=12 603。試驗(yàn)過程中，同一問題規(guī)模下不同試驗(yàn)方案的初始值等輸入條件是完全相同的。

表4　4種試驗(yàn)方案下的計(jì)算耗時(shí)統(tǒng)計(jì)　s

從表4可以看出，方案1無論是時(shí)間消耗還是空間消耗都是所有方案中最大的，計(jì)算效率最低。理論上隨問題規(guī)模的增長呈冪數(shù)級時(shí)間復(fù)雜度規(guī)律，而當(dāng)問題規(guī)模為12 603時(shí)，甚至出現(xiàn)內(nèi)存溢出、程序崩潰現(xiàn)象。相比方案1，方案2計(jì)算耗時(shí)陡降。而且更進(jìn)一步，方案3和方案4分別又比方案2計(jì)算效率提高了約5倍和11倍，與理論分析基本一致，達(dá)到了預(yù)期大幅度提高平差計(jì)算效率的目的，而且與作業(yè)單位提供的成果相比，平差精度相當(dāng)，甚至還略高些(見表3)，表明本文所提出策略實(shí)現(xiàn)的平差軟件不會(huì)額外損失平差精度。

四、結(jié)束語

法化系數(shù)矩陣的結(jié)構(gòu)和規(guī)模是影響光束法區(qū)域網(wǎng)平差空三計(jì)算效率的最主要因素之一，針對這種情況，本文提出了一種攝影區(qū)域跨航線影像片號編排下的法化系數(shù)矩陣分塊降維精簡策略，通過理論分析和實(shí)際作業(yè)航攝影像數(shù)據(jù)平差試驗(yàn)，表明基于該策略的數(shù)據(jù)平差處理算法相比常規(guī)平差算法，在時(shí)間效率和存儲(chǔ)效率上都有大幅度提升。這對空三數(shù)據(jù)處理相關(guān)程序算法設(shè)計(jì)與軟件實(shí)現(xiàn)都有很好的借鑒意義，對航攝規(guī)劃、影像獲取及數(shù)據(jù)處理過程中的片號編排也有一定的指導(dǎo)作用。當(dāng)然，這里的兩種像片排序策略比對分析只是針對規(guī)則區(qū)域攝影進(jìn)行了討論。對于非規(guī)則攝影區(qū)域，應(yīng)該研究更為適用的方法以求得最小帶寬，從而節(jié)省解算時(shí)間，后續(xù)還存在更進(jìn)一步的研究空間。

[1]MCGLONE J C，MIKHAIL E M，BETHEL J S，et al. Manual of Photogrammetry [M]. Maryland：American Society for Photogrammetry and Remote Sensing, 2004.

[2]SHAFFER C A. 數(shù)據(jù)結(jié)構(gòu)與算法分析: C++ 版[M].北京：電子工業(yè)出版社, 2002.

[3]嚴(yán)蔚敏，吳偉民. 數(shù)據(jù)結(jié)構(gòu)[M].北京：清華大學(xué)出版社, 1992.

[4]WOLF P，DEWITT B，WILKINSON B. Elements of Photogrammetry with Application in GIS[M]. New York: McGraw-Hill Education, 2013.

[5]EGELS Y，KASSER M. Digital Photogrammetry[M]. London: Taylor & Francis, 2010.

[6]MORGAN D，F(xiàn)ALKNER E. Aerial Mapping: Methods and Applications[M]. London: Taylor & Francis, 2010.

[7]LUHMANN T，ROBSON S，KYLE S. Close Range Photogrammetry: Principles, Techniques and Applications[M]. Caithness:Whittles Publishing, 2006.

[8]江延川. 解析攝影測量學(xué)[M]. 鄭州: 信息工程大學(xué),2001.

[9]王樹根. 攝影測量原理與應(yīng)用[M]. 武漢: 武漢大學(xué)出版社, 2009.

[10]王之卓. 攝影測量原理[M]. 武漢: 武漢大學(xué)出版社, 2007.

[11]BEHAN A. Digital Photogrammetry: A Practical Course[J]. The Photogrammetric Record, 2010, 25(132): 476-477.

[12]STRANG G. Introduction to Linear Algebra[M]. Wellesley: Wellesley-Cambridge Press, 2009.

[13]DEMMEL J W. Applied Numerical Linear Algebra[M]. Philadelphia: Siam, 1997.

[14]MARCUS M. Introduction to Linear Algebra[M]. Mineola: Courier Dover Publications, 1988.

[15]FOX L. An Introduction to Numerical Linear Algebra[M]. Oxford: Clarendon Press, 1964.

[16]SEARLE S R. Matrix Algebra Useful for Statistics[M]. New York: Wiley-Interscience, 2006.

[17]SCHWERDTFEGER H. Introduction to Linear Algebra and the Theory of Matrices[M]. Groningen: P. Noordhoff, 1950.引文格式： 桂新，戴微. MASK在遙感影像水系提取中的應(yīng)用[J].測繪通報(bào)，2015(9)：79-82.DOI:10.13474/j.cnki.11-2246.2015.0285