區域網空中三角測量數據平差計算效率提升策略

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

(1. 信息工程大學地理空間信息學院,河南 鄭州 450052; 2. 地理信息工程國家重點實驗室,

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

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

區域網空中三角測量數據平差計算效率提升策略

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

(1. 信息工程大學地理空間信息學院,河南 鄭州 450052; 2. 地理信息工程國家重點實驗室,

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

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

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

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

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

關鍵詞：光束法平差；空三；線性代數；計算效率

中圖分類號：P231

文獻標識碼：B

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

收稿日期：2014-09-09

作者簡介：胡海彥(1977—)，男，博士生，研究方向為航空航天高精度定位理論與方法。E-mail:173412586@qq.com

一、引言

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

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

二、平差計算策略

1. 分塊精簡法化系數矩陣

在量取同名像點坐標及控制點對應像點坐標后，空三數據平差的主要技術環節是組建觀測方程組，實現途徑是對共線條件方程線性化后依據最小二乘原理得到法化方程組，它是關于未知數改正量的線性方程組

NΔ=K

(1)

參考文獻式中，N為法化后方程組系數矩陣；Δ為未知數改正向量；K為常數向量，矩陣中各元素的具體描述可見[11]。高效求解該線性系統的關鍵是法化系數矩陣N的結構分析與研究[12-13]。對于區域網平差而言，法化系數矩陣N往往特別巨大，基本策略是進行矩陣分塊，以達到矩陣降維，從而減少計算時間及內存消耗的目的。

式(1)可根據外方位元素和物空間坐標兩類未知數進行分塊

(2)

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

(3)

(4)

(5)

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

(6)

即

(7)

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

(8)

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

(9)

(10)

(11)

2. 像片編排次序優化

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

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

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

圖2　關于連接點的像片間重疊拓撲關系

表1　沿/跨航線圖1攝影區域影像排序

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

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

三、試驗及分析

有關試驗數據的測區情況與區域網基本參數見表2，每條航線18張影像，野外測量點69個(29個用于控制，40個用于精度檢查)，采用Geolord-AT商用平差軟件進行作業計算，結果見表3。

表2　區域網基本參數

圖3　沿航線、跨航線像片排序策略下的精簡法化方程系數矩陣結構

耗時/sσ0/μmX/mY/mZ/m1328.770.1090.0970.4248.927.800.0810.0920.383

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

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

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

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

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

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

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

表4　4種試驗方案下的計算耗時統計　s

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

四、結束語

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

[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. 數據結構與算法分析: C++ 版[M].北京：電子工業出版社, 2002.

[3]嚴蔚敏，吳偉民. 數據結構[M].北京：清華大學出版社, 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，FALKNER 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]江延川. 解析攝影測量學[M]. 鄭州: 信息工程大學,2001.

[9]王樹根. 攝影測量原理與應用[M]. 武漢: 武漢大學出版社, 2009.

[10]王之卓. 攝影測量原理[M]. 武漢: 武漢大學出版社, 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在遙感影像水系提取中的應用[J].測繪通報，2015(9)：79-82.DOI:10.13474/j.cnki.11-2246.2015.0285