大角度立體像對相對定向的混合共軛梯度算法

李佳田，王聰聰，賈成林，牛一如，王 瑜，張文靖，吳華靜，李 鍵

1. 昆明理工大學(xué)國土資源工程學(xué)院，云南 昆明 650093； 2. 昆明理工大學(xué)云南省高校高原山區(qū)空間信息測繪技術(shù)應(yīng)用工程研究中心，云南 昆明 650093

立體像對相對定向是恢復(fù)攝影時刻相鄰兩像片攝影光束的相互關(guān)系，從而使同名光線對對相交[1]。相對定向算法研究具有典型的實(shí)際意義。

在航空攝影測量中，立體像對相對定向時像片間的旋偏角一般小于6°，可對相對定向元素中的基線分量近似，以簡化相對定向模型，在相對定向元素初值為0的前提下，通常采用最速下降法迭代求解[1-2]。隨著傾斜攝影測量和低空攝影測量的發(fā)展，立體像對相對定向時像片間的旋偏角可能達(dá)到甚至超過10°，最速下降法常因相對定向元素初值無法獲取而失敗[3]，一般采用相對定向線性變換法求解[4]。上述兩種算法的不足之處是：最速下降法具有線性收斂速度、求解精度高，但過分依賴相對定向元素初值[5-6]，使得大角度相對定向時收斂速度減慢、收斂至局部極值甚至不收斂[7]；相對定向線性變換是一種直接解法，雖然無須預(yù)知相對定向元素的初值，但求解精度不高，一般作為光束法平差的初值[8]。因此，兩種算法均較難兼顧收斂性與求解精度兩個方面的要求。

關(guān)于初值依賴問題，國內(nèi)外學(xué)者已取得豐碩的成果，主要分為代數(shù)法和矩陣分解法兩類[9]。代數(shù)法是根據(jù)共面方程的約束關(guān)系，將相對定向元素引入約束方程，構(gòu)造含有未知數(shù)的約束方程組[7，10-16]。其中，文獻(xiàn)[7]將計(jì)算機(jī)視覺中的8點(diǎn)算法引入相對定向中，無須初值，僅需8個同名點(diǎn)即可解算相對定向元素；較之8點(diǎn)算法，5點(diǎn)算法因有較高精度而被關(guān)注，因其采用多項(xiàng)式求解，存在多個解及誤解排除的問題[13-16]。矩陣分解法是通過分解基礎(chǔ)矩陣或本質(zhì)矩陣解算相對定向元素[17-19]。文獻(xiàn)[18]在顧及基礎(chǔ)矩陣元素非獨(dú)立條件下實(shí)現(xiàn)無須初值的相對定向元素求解。文獻(xiàn)[19]從基礎(chǔ)矩陣中直接導(dǎo)出立體像對像片間的平移和旋轉(zhuǎn)量，較難得到精確的相對定向元素。近年來，針對最速下降法具有較高計(jì)算精度，依賴于初值的特點(diǎn)，有學(xué)者提出混合算法以克服初值依賴并提高相對定向的解算精度。文獻(xiàn)[20]通過分解本質(zhì)矩陣得到相對定向元素初值，然后采用最速下降法迭代求解。文獻(xiàn)[21]結(jié)合遺傳算法的全局收斂性與最速下降法的局部收斂性，利用兩者的混合算法求解大角度立體像對相對定向。

初值依賴性、數(shù)值收斂性及求解精度是解算相對定向需考慮的3個問題。文獻(xiàn)[22]對常規(guī)模擬退火法產(chǎn)生新的待估參數(shù)及接受概率進(jìn)行改進(jìn)，提出了隨機(jī)爬山(stochastic hill climbing，SHC)算法。它具有全局收斂性，且較之遺傳算法與模擬退火法，其具有收斂速度快的優(yōu)點(diǎn)。因此，為解決大角度相對定向時，最速下降法對相對定向元素初值依賴的問題，并加快算法的收斂速度，本文提出一種基于SHC算法與共軛梯度的混合共軛梯度算法，采用超線性收斂的共軛梯度法替代最速下降法，加快局部收斂速度。當(dāng)相對定向元素初值未知時，為避免共軛梯度法陷入局部極值并加快全局收斂速度，利用SHC算法為其提供迭代初值。

1 混合共軛梯度法解算大角度相對定向

1.1 相對定向模型

連續(xù)像對相對定向是以左像片為基準(zhǔn)，確定右像片相對于左像片的相對定向元素(BX,BY,BZ,φ,ω,κ)，如圖1所示。

圖1 相對定向示意圖Fig.1 Schematic diagram of relative orientation

圖1中以左像片的像空間坐標(biāo)系為像對的像空間輔助坐標(biāo)系，記為S1-X1Y1Z1，以右攝影中心S2為原點(diǎn)建立另一個像空間輔助坐標(biāo)系，記為S2-X2Y2Z2。兩坐標(biāo)系互相平行，且左右像點(diǎn)a1,a2在各自像空間輔助坐標(biāo)系中的坐標(biāo)分別為(X1,Y1,Z1)和(X2,Y2,Z2)，像點(diǎn)坐標(biāo)分別為(x1,y1)和(x2,y2)，S2在S1-X1Y1Z1中的坐標(biāo)為(BX,BY,BZ)，則每對同名像點(diǎn)的共面條件方程為[1]

(1)

式中

(2)

式中，f為攝像機(jī)焦距；R為相對定向φ、ω、κ角元素方向余弦構(gòu)成的矩陣。

相對定向中，由于立體像對像片間的旋偏角較小，一般采用簡化模型求解，即基線分量BX已知，BY和BZ用小角度μ、ν近似表示

(3)

對式(3)中5個相對定向元素BY、BZ、φ、ω、κ分別求導(dǎo)，采用最速下降法迭代求解。在大角度相對定向情況之下，存在如下問題：

(1) 立體像對像片間的旋偏角較大，式(3)的簡化模型不再適用。

(2) 最速下降法具有線性收斂速度，因此，收斂速率方面仍有改進(jìn)的空間。

(3) 最速下降法是局部收斂算法，在相對定向元素初值難以獲取的情況下，計(jì)算會陷入局部極值，甚至不收斂。

1.2 混合共軛梯度法

SHC算法是一種改進(jìn)后的模擬退火法，同樣具有全局收斂的性質(zhì)，它通過對模擬退火算法中新的待估參數(shù)產(chǎn)生過程和概率重新修改，獲得比模擬退火法更快的收斂速度。設(shè)待估參數(shù)當(dāng)前值和新值分別為Xk和Xk+1，f為目標(biāo)函數(shù)，則其只接受f(Xk+1)

Xk+1=Xk+sign·random·dx

(4)

式中，random為[0，1]均勻分布的隨機(jī)數(shù)；當(dāng)x≥0時，sign=1，否則sign=-1；dx為給定的擾動步長，控制待估參數(shù)的修改幅度。

為進(jìn)一步加快算法收斂速度，在局部搜索過程中，考慮采用共軛梯度法進(jìn)行計(jì)算，即混合共軛梯度算法，計(jì)算過程如圖2所示。

(1) 按式(4)擾動當(dāng)前待估參數(shù)值來產(chǎn)生新值。

(2) 判斷產(chǎn)生的新值是否符合要求。

(3) 若符合要求，則接受新值，并以新值為初值進(jìn)行共軛梯度法迭代，否則，返回步驟(1)繼續(xù)擾動產(chǎn)生新值。

(4) 判斷是否滿足全局收斂條件，滿足則結(jié)束計(jì)算，否則，則返回(1)。

結(jié)合混合共軛梯度算法的特點(diǎn)及1.1節(jié)大角度相對定向存在的3個問題，本文具體改進(jìn)是：

(1) 采用具有全局收斂性的SHC算法進(jìn)行全局搜索，為局部搜索提供迭代初值。

(2) 不采用小角度μ、ν近似，BX已知的簡化模型，而直接對BX、BY、BZ、φ、ω、κ這6個相對定向元素求導(dǎo)。

(3) 采用共軛梯度法代替最速下降法，加快局部搜索的收斂速度。

圖2 混合共軛梯度法Fig.2 Hybrid conjugate gradient

1.3 共軛梯度法

由上文可知，求解式(1)中6個相對定向元素BX、BY、BZ、φ、ω、κ，至少需要6對同名像點(diǎn)。為提高解算精度，常有多余觀測。設(shè)有n對同名像點(diǎn)(n>6)，則有由n個式(1)構(gòu)成的方程組

H(p)=[h1(p),h2(p),…,hn(p)]T=0

(5)

式中，p=[BX,BY,BZ,φ,ω,κ]T為相對定向元素。式(5)是典型的非線性優(yōu)化過程[23]。共軛梯度法的求解思路是：給定一個相對定向元素初值p0，假定初值p0在真值p*鄰域內(nèi)，則相對定向元素可通過式(6)迭代計(jì)算

pk+1=pk+αkPk

(6)

式中

P0=-g0

(7)

Pk+1=-gk+1+βkPk

(8)

式中，αk為1維搜索步長，通過線性搜索獲得；Pk為搜索方向；βk為參數(shù)。不同的βk對應(yīng)不同算法，文中采用Polak-Ribiere-Polyak算法[24]

(9)

(10)

進(jìn)一步由式(1)可得

(11)

因此，在給定1維搜索步長αk的情況下，若相對定向元素初值充分接近真值，則相對定向問題可利用共軛梯度法多次迭代式(6)解算。

1.4 SHC算法[25]

在大角度相對定向中，若相對定向元素初值充分接近真值，則大角度相對定向問題可由超線性收斂的共軛梯度法解算。但其相對定向元素初值一般較難獲取，若將相對定向?qū)⒊踔翟O(shè)為0，則會使共軛梯度法收斂于局部極小值，甚至不收斂。為此，在大角度相對定向元素初值未知的情況下，采用SHC算法計(jì)算相對定向元素初值。具體過程如下：

(1) 已知左右像片n對同名像點(diǎn)坐標(biāo)和攝影機(jī)焦距f，給定相對定向元素初值pk。

(2) 由式(2)計(jì)算左右像對同名像點(diǎn)在各自像空間輔助坐標(biāo)系的坐標(biāo)(X1,Y1,Z1)和(X2,Y2,Z2)。

(3) 由式(1)、式(5)計(jì)算得到H(pk)。

(4) 按式(4)修改當(dāng)前相對定向元素pk以產(chǎn)生新的相對定向元素pnew，計(jì)算得到ΔH=H(pnew)-H(pk)。

2 試驗(yàn)與分析

為驗(yàn)證本文算法的有效性，分別進(jìn)行模擬試驗(yàn)和實(shí)測試驗(yàn)。

2.1 模擬試驗(yàn)

已知攝像機(jī)焦距f及立體像對左像片像點(diǎn)坐標(biāo)，以左像片為基準(zhǔn)，采用交向攝影模擬產(chǎn)生3組右像片分別與左像片構(gòu)成虛擬立體像對，相對定向元素模擬值見表1，其中僅考慮了大角度相對定向。然后分別計(jì)算不同相對定向元素下右像片的像點(diǎn)坐標(biāo)，并在其中加入均值為0，方差為0.05的高斯噪聲。為模擬大角度立體像對，假設(shè)成像不受圖幅限制。以左右像片的像點(diǎn)坐標(biāo)、攝像機(jī)焦距f為已知條件，采用本文算法、最速下降法及文獻(xiàn)[21]中基于最速下降法的混合遺傳算法(簡稱為算法3)分別對相對定向元素求解，并將相對定向元素計(jì)算值與模擬值進(jìn)行比較。試驗(yàn)中，相對定向線元素和角元素的模型擾動步長dx分別為0.05 mm和0.02 rad。

由表1可知，由于相對定向元素模擬值傾角較大，為保證最速下降法正確收斂，試驗(yàn)中采用相對定向線性變換法[4]為其提供初值。而本文算法和算法3中相對定向元素初值均設(shè)為0，即BX=BY=BZ=φ=ω=κ=0。3種算法計(jì)算結(jié)果見表2，本文算法和算法3的迭代次數(shù)見表3。

表1 相對定向元素模擬值

表2 相對定向元素計(jì)算結(jié)果

表3 迭代次數(shù)

由表2、表3可知，本文算法的計(jì)算精度高于算法3，與最速下降法相當(dāng)，但最速下降法依賴于相對定向元素初值，為保證其正確收斂，需預(yù)知相對定向元素初值；而本文算法無須提供初值，且較之算法3，收斂速度較快。原因在于：混合共軛梯度算法中采用了共軛梯度法，具有超線性收斂速度，雖然可能陷入局部極值，但由于采用SHC算法通過擾動產(chǎn)生新的隨機(jī)相對定向元素，并按目標(biāo)函數(shù)值對產(chǎn)生的相對定向元素進(jìn)行取舍，保證收斂方向和速度。因此，在多次迭代中，能夠跳出共軛梯度法陷入的局部極值。試驗(yàn)中，當(dāng)相對定向元素初值距離模擬值較遠(yuǎn)時，混合共軛梯度算法在不同極值處進(jìn)行多次局部優(yōu)化，并在每次局部搜索過程中，因SHC算法對所產(chǎn)生的相對定向元素進(jìn)行判斷，越來越接近模擬值，經(jīng)多次迭代后，最終收斂。

2.2 實(shí)測試驗(yàn)

如圖3所示，為云南省某地居民樓房兩條基線下交向攝影獲得的3幅像片，構(gòu)成2對連續(xù)立體像對。航拍相機(jī)為Canon 5D Mark Ⅱ，鏡頭焦距為35 mm，相對航高445 m，圖像大小為449×300像素，像素大小為6.4 μm。為減少同名點(diǎn)定位帶來的外部誤差，采用人工量測的8個同名像點(diǎn)來進(jìn)行相對定向元素的解算。采用最速下降法、算法3和本文算法分別對立體像對進(jìn)行相對定向元素解算。由于立體像對像片間的旋偏角較大，試驗(yàn)中采用相對定向線性變換法[4]為最速下降法提供初值，而本文算法和算法3中，相對定向元素初值設(shè)為0，即BX=BY=BZ=φ=ω=κ=0，3種算法的計(jì)算結(jié)果見表4，本文算法和算法3的迭代次數(shù)見表5。試驗(yàn)中相對定向線元素和角元素的模型擾動步長dx分別為0.05 mm和0.02 rad。

圖3 像片序列Fig.3 Sequence of images

立體像對算法BX/mmBY/mmBZ/mmφ/radω/radκ/rad精度/μm像片2/像片1最速下降法-1.3675-0.01620.0524-1.06470.0698-1.57082.1算法3-1.3674-0.01620.0521-1.06450.0698-1.57062.4本文算法-1.3675-0.01630.0525-1.06480.0699-1.57072.2像片3/像片1最速下降法-1.39110.02450.0411-1.02970.02540.04361.8算法3-1.39100.02450.0413-1.02960.02540.04382.1本文算法-1.39110.02440.0410-1.02970.02550.04371.8

表5 迭代次數(shù)

從表4、表5可知，本文算法計(jì)算精度與最速下降法相當(dāng)，相對定向達(dá)到約±0.5像素的精度，而最速下降法因依賴于相對定向元素初值，需要相對定向線性變換法提供初值；較之算法3，本文算法具有計(jì)算精度高且收斂快的優(yōu)點(diǎn)。

3 結(jié) 論

提出了一種全局收斂的混合共軛梯度算法。試驗(yàn)表明，在無須預(yù)知相對定向元素初值的前提下，算法能夠正確解算大角度立體像對相對定向，并具有解算精度高和收斂速度快的特點(diǎn)。此外，合理選取SHC算法中的擾動步長將會進(jìn)一步加快收斂速度。針對大角度相對定向解算，擾動步長的定量分析與選取有待進(jìn)一步研究。