基于二維dl的雙片相機(jī)主縱線方程

基于二維dl的雙片相機(jī)主縱線方程

與攝像頭的內(nèi)部相位元素和光學(xué)變形參數(shù)的校正過(guò)程稱(chēng)為攝像頭校正。對(duì)于正在使用的非量規(guī)相機(jī)，校正過(guò)程是從二維圖像中接收三維信息是必要的。隨著攝影測(cè)量和計(jì)算機(jī)視覺(jué)理論的發(fā)展,許多學(xué)者對(duì)攝像機(jī)標(biāo)定技術(shù)進(jìn)行了深入的研究。計(jì)算機(jī)視覺(jué)界最近提出了利用主動(dòng)視覺(jué)技術(shù)進(jìn)行攝像機(jī)自標(biāo)定的一些算法,不過(guò)在少數(shù)情況下,自標(biāo)定技術(shù)得到的解不是惟一的,即使在圖像噪聲很小時(shí),求解的值與實(shí)際值也有較大的差別。BillTriggs提出利用絕對(duì)二次曲線原理進(jìn)行攝像機(jī)自標(biāo)定,但該方法的初始化比較困難。張正友提出了利用旋轉(zhuǎn)矩陣的正交性條件和非線性最優(yōu)化進(jìn)行攝像機(jī)標(biāo)定,其標(biāo)定精度為0.35像素左右。直接線性變換(directlineartransformation,DLT)是建立像點(diǎn)坐標(biāo)和物點(diǎn)坐標(biāo)直接線性關(guān)系的算法,處理時(shí)不需要攝像機(jī)內(nèi)外方位元素的初值,因而在近景攝影測(cè)量中被廣泛使用。國(guó)際上利用三維DLT進(jìn)行攝像機(jī)標(biāo)定已經(jīng)有較多的文獻(xiàn),但到現(xiàn)在為止還沒(méi)有見(jiàn)到利用二維DLT進(jìn)行攝像機(jī)標(biāo)定的算法。1網(wǎng)點(diǎn)可轉(zhuǎn)化的形式二維DLT可表示為:x=h1X+h2Y+h3h7X+h8Y+1y=h4X+h5Y+h6h7X+h8Y+1(1)式中,h1、h2、h3、h4、h5、h6、h7、h8為二維DLT的8個(gè)變換參數(shù)H=(h1,h2,h3,h4,h5,h6,h7,h8)T;X、Y為平面控制(格網(wǎng))點(diǎn)空間坐標(biāo)(Z坐標(biāo)為零);x,y為相應(yīng)的像坐標(biāo),當(dāng)像片點(diǎn)數(shù)大于4個(gè)時(shí),可將式(1)進(jìn)行適當(dāng)變換,則二維DLT參數(shù)可通過(guò)解超定方程AH=0求得。攝影測(cè)量中最常用的共線方程為:x-x0=-f[a1(X-Xs)+b1(Y-Ys)+c1(Ζ-Ζs)]/[a3(X-Xs)+b3(Y-Ys)+c3(Ζ-Ζs)]y-y0=-f[a2(X-Xs)+b2(Y-Ys)+c2(Ζ-Ζs)]/[a3(X-Xs)+b3(Y-Ys)+c3(Ζ-Ζs)](2)式中,x0、y0、f為攝像機(jī)的內(nèi)方位元素;(XS,YS,ZS)為攝站坐標(biāo);(X,Y,Z)為物方空間坐標(biāo);(x,y)為相應(yīng)的像坐標(biāo);R={ai,bi,ci,i=1,2,3}為攝影測(cè)量中常用的旋轉(zhuǎn)角φ、ω、κ(Y為主軸)構(gòu)成的旋轉(zhuǎn)矩陣。對(duì)于平面格網(wǎng)而言,一般將坐標(biāo)系建立在格網(wǎng)中心,X、Y軸在平面內(nèi),此時(shí)所有格網(wǎng)點(diǎn)的Z坐標(biāo)為零。式(2)可轉(zhuǎn)化與式(1)類(lèi)似的形式:x=(fa1λ-a3λx0)X+(fb1λ-b3λx0)Y+(x0-fλ(a1Xs+b1Ys+c1Ζs))-a3λX-b3λY+1y=(fa2λ-a3λy0)X+(fb2λ-b3λy0)Y+(y0-fλ(a2Xs+b2Ys+c2Ζs))-a3λX-b3λY+1(3)式中,λ=(a3XS+b3YS+c3ZS)。比較式(1)和式(3)可知:{h1=fa1/λ-a3x0/λh2=fb1/λ-b3y0/λ(4){h4=fa2/λ-a3y0/λh5=fb2/λ-b3y0/λ(5){h3=x0-f(a1Xs+b1Ys+c1Ζs)/λh6=y0-f(a2Xs+b2Ys+c2Ζs)/λ(6){h7=-a3/λh8=-b3/λ(7)由式(4)、式(5)、式(7)得:{(h1-h7x0)/f=a1/λ(h2-h8x0)/f=b1/λ(8){(h4-h7y0)/f=a2/λ(h5-h8y0)/f=b2/λ(9){-h7=a3/λ-h8=b3/λ(10)將式(8)、式(9)、式(10)上下兩式分別相乘并相加,顧及a1b1+a0b2+a3b3=0,得:(h1-h7x0)?(h2-h8x0)/f2+(h4-h7y0)?(h5-h8y0)/f2+h7h8=0(11)則有∶f=√(-(h1-h7x0)?(h2-h8x0)-(h4-h7y0)?(h5-h8y0))/h7h8(12)當(dāng)主點(diǎn)(x0,y0)已知或者通過(guò)某種方法求得后,即可通過(guò)式(12)求出焦距f。將式(8)、式(9)、式(10)上下兩式分別自乘并相加,顧及a21+a22+a23=1和b21+b22+b23=1并消去λ得:(h1-h7x0)2-(h2-h8x0)2+(h4-h7y0)2-(h5-h8y0)2f2+(h27-h28)=0(13)利用式(11)和式(13)消去焦距即可得:Fh=(h1h8-h2h7)(h1h7-h27x0+h2h8-h28x0)+(h4h8-h5h7)(h4h7-h27y0+h5h8-h28y0)=0(14)不考慮鏡頭畸變影響時(shí),攝像機(jī)的實(shí)際未知數(shù)為9個(gè)(f,x0,y0,φ,ω,κ,XS,YS,ZS),而二維DLT共有8個(gè)參數(shù),則必然無(wú)法惟一分解出攝像機(jī)的9個(gè)未知參數(shù)。事實(shí)上,在給定二維DLT的8個(gè)參數(shù)時(shí),主點(diǎn)(x0,y0)可在主縱線上自由移動(dòng),從而造成外方位元素分解的不惟一性。證明:現(xiàn)取用以Z軸為主軸的A、υ、κ轉(zhuǎn)角系統(tǒng),有tanκ=c1c2,其中κ角為像片上主縱線與y軸的夾角。將下文的式(22)代入可得tanκ=c1c2=a2b3-a3b2a3b1-a1b3=a2/a3-b2/b3b1/b3-a1/a3,將式(21)的相應(yīng)項(xiàng)代入可得tanκ=h4h8-h5h7h1h8-h2h7。顯然,式(14)可表達(dá)為直線y0=Ax0+C的形式,其中A=tanα=h1h8-h2h7h4h8-h5h7為直線的斜率,α為直線與x軸的夾角。不難看出,tanκ=1/tanα,即κ=90°-α,可見(jiàn)式(14)實(shí)質(zhì)上就是主縱線在圖像坐標(biāo)系下的表達(dá)式。因而只要主點(diǎn)(x0,y0)的位置在主縱線y0=Ax0+C上,所有分解出的攝像機(jī)參數(shù)都滿(mǎn)足投影關(guān)系,這也說(shuō)明單張像片無(wú)法利用平面控制網(wǎng)進(jìn)行攝像機(jī)標(biāo)定。主縱線方程式(14)也可表示為(Lx,Ly)(x0y0)=Lc的形式,下文簡(jiǎn)寫(xiě)為L(zhǎng)X=c,其中的各個(gè)系數(shù)可由式(14)得出。并且攝像機(jī)的主點(diǎn)必定位于主縱線上,也即主點(diǎn)(x0,y0)必然滿(mǎn)足方程LX=c。因此,如果有兩張或兩張以上像片,就可以利用多條主縱線通過(guò)解超定方程LX=c求得主點(diǎn)(x0,y0)的坐標(biāo)。值得注意的是,在利用超定方程LX=c求解主點(diǎn)(x0,y0)時(shí)應(yīng)該避免所謂的臨界運(yùn)動(dòng)序列(criticalmotionsequences)。當(dāng)攝像機(jī)固定而標(biāo)定格網(wǎng)只繞其Z軸作旋轉(zhuǎn)時(shí),各像片的二維DLT參數(shù)之間線性相關(guān),此時(shí)各主縱線互相重合,因而無(wú)法求出主點(diǎn)的位置。證明如下。采用齊次坐標(biāo),針孔(pinhole)模型下的平面投影關(guān)系可表示為:s(xy1)=A(r1,r2,t)(XY1)=Η(XY1)(15)式中,A為內(nèi)方位元素矩陣;r1、r2和t分別為旋轉(zhuǎn)矩陣的前兩列和平移向量;H=A(r1,r2,t)稱(chēng)為同形矩陣(homography)。不難看出,利用像點(diǎn)齊次坐標(biāo)的第三行消去比例因子s后的式(15)就是二維DLT方程,因此,二維DLT參數(shù)實(shí)際上與計(jì)算機(jī)視覺(jué)中常用的同形矩陣(第9個(gè)元素取1時(shí))等價(jià)?，F(xiàn)取用同形矩陣的表達(dá)形式,當(dāng)格網(wǎng)只繞Z軸旋轉(zhuǎn)θ角時(shí),兩張像片的旋轉(zhuǎn)矩陣間存在如下關(guān)系:R2=R1(cosθ-sinθ0sinθcosθ0001)(16)用上標(biāo)(1)、(2)分別表示對(duì)應(yīng)于第一和第二張像片的同形矩陣元素,將式(16)代入H=A(r1,r2,t)得:(h(2)1h(2)4h(2)7)=(h(1)1h(1)4h(1)7)cosθ+(h(1)2h(1)5h(1)8)sinθ(h(2)2h(2)5h(2)8)=-(h(1)1h(1)4h(1)7)sinθ+(h(1)2h(1)5h(1)8)cosθ(17)主縱線可表示為y0=Ax0+C的形式,若設(shè)α為直線與x軸的夾角,則直線的斜率為:tanα=(h1h8-h2h7)/(h4h8-h5h7)(18)將式(17)中的相應(yīng)項(xiàng)代入式(18)中可得:h(2)1h(2)8-h(2)2h(2)7=(h(1)1cosθ+h(1)2sinθ)(-h(1)7sinθ+h(1)8cosθ)-(-h(1)1sinθ+h(1)2cosθ)(h(1)7cosθ+h(1)8sinθ)=h(1)1h(1)8-h(1)2h(1)7(19)同理可得:h(2)4h8(2)-h5(2)h7(2)=h4(1)h(1)8-h5(1)h7(1)。由此可知兩張像片的主縱線斜率相等,也即主縱線互相平行。而且,由式(14)知,常數(shù)項(xiàng)C可表示為:C=[(h1h8-h2h7)(h1h7+h2h8)+(h4h8-h5h7)(h4h7-h5h8)]/(h72+h82)(20)由式(17)不難得出,C(2)=C(1)。由此可知兩張像片的主縱線完全重合,當(dāng)然無(wú)法求交點(diǎn)(即主點(diǎn))。因而在實(shí)踐中應(yīng)設(shè)法消除二維DLT參數(shù)之間的相關(guān)性,最簡(jiǎn)單的方法是手持?jǐn)z像機(jī),每次拍攝時(shí)都采用3個(gè)不同的旋轉(zhuǎn)角進(jìn)行像片拍攝,使旋轉(zhuǎn)矩陣之間失去相關(guān)性,從而使主縱線的斜率間有較大差異,以便穩(wěn)定地計(jì)算主點(diǎn)位置。以上詳細(xì)討論了攝像機(jī)內(nèi)方位元素初值(x0,y0)及f的求解方法及實(shí)踐中應(yīng)該注意的臨界序列問(wèn)題,下文將討論外方位元素的分解。將式(10)分別代入式(8)、式(9)可得:a1/a3=-(h1-h7x0)/fh7b1/b3=-(h2-h8x0)/fh8a2/a3=-(h4-h7y0)/fh7b2/b3=-(h5-h8y0)/fh8(21)由b12+b22+b32=1可得:b32=11+(h2-h8x0)2f2h82+(h5-h8y0)2f2h82在以Y為主軸的轉(zhuǎn)角系統(tǒng)下,tanκ=b1b2,由式(21)知,tanκ=b1b2=h2-h8x0h5-h8y0,此式可惟一確定κ角。在求解ω角時(shí),b3的值在開(kāi)平方后首先取正號(hào),將已確定的κ角與通過(guò)b3求得的b1、b2算出的κ相比較,若κ!=κ則說(shuō)明b3應(yīng)取負(fù)號(hào),然后重新計(jì)算b1、b2的值。通過(guò)sinω=-b3即可計(jì)算出ω角。據(jù)旋轉(zhuǎn)矩陣的正交性可得:(c1c2c3)=(a1a2a3)×(b1b2b3)=(a2b3-a3b2a3b1-a1b3a1b2-a2b1)(22)因?yàn)閠anφ=-a3c3=1a1a3b2-a2a3b1,b1、b2已在求ω角時(shí)確定,而a1a3和a2a3可由式(21)確定,因而?角也可惟一確定?？梢钥闯?在求解φ,ω,κ角時(shí),并沒(méi)有計(jì)算整個(gè)旋轉(zhuǎn)矩陣中的所有9個(gè)元素,因而在計(jì)算XS、YS、ZS的初值時(shí),需利用如上計(jì)算出的角φ,ω,κ,重新計(jì)算旋轉(zhuǎn)矩陣的各元素。由式(8)、式(9)、式(10)求平均可得到λ,則XS、YS、ZS的初值可通過(guò)解如下線性方程組獲得:{h3=x0-fλ(a1XS+b1YS+c1ΖS)h6=y0-fλ(a2XS+b2YS+c2ΖS)λ=a3XS+b3YS+c3ΖS(23)2誤差方程的解析運(yùn)用解算出攝像機(jī)的內(nèi)外方位元素初值后,即可用攝影測(cè)量中常用的光束法平差進(jìn)行攝像機(jī)標(biāo)定。由于非量測(cè)數(shù)字?jǐn)z像機(jī)鏡頭一般存在畸變,標(biāo)定時(shí)需求定畸變參數(shù)。引入畸變差后的共線方程為:x-x0-Δx=-fxa1(X-XS)+b1(Y-YS)+c1(Ζ-ΖS)a3(X-XS)+b3(Y-YS)+c3(Ζ-ΖS)=-fxXˉΖˉy-y0-Δy=-fya2(X-XS)+b2(Y-YS)+c2(Ζ-ΖS)a3(X-XS)+b3(Y-YS)+c3(Ζ-ΖS)=-fyYˉΖˉ(24)式中,fx、fy分別為x、y兩個(gè)方向的焦距,此處假設(shè)其有效值可能不相等。Δx=(x-x0)(Κ1r2+Κ2r4)+Ρ1(r2+2(x-x0)2)+2Ρ2(x-x0)(y-y0)Δy=(y-y0)(Κ1r2+Κ2r4)+Ρ2(r2+2(y-y0)2)+2Ρ2(x-x0)(y-y0)(25)式中,r2=(x-x0)2+(y-y0)2;K1,K2稱(chēng)為徑向畸變差(radialdistortion);P1、P2稱(chēng)為偏心畸變差(decenteringdistortion)。用泰勒級(jí)數(shù)將式(24)線性化,即可得到用于標(biāo)定的誤差方程式:vx=?x?XSΔXS+?x?YSΔYS+?x?ΖSΔΖS+?x?φΔφ+?x?ωΔω+?x?κΔκ+?x?XΔX+?x?YΔY+?x?ΖΔΖ+?x?fxΔfx+?x?fyΔfy+?x?x0Δx0+?x?y0Δy0+?x?Κ1ΔΚ1+?x?Κ2ΔΚ2+?x?Ρ1ΔΡ1+?x?Ρ2ΔΡ2-lxvy=?y?XSΔXS+?y?YSΔYS+?y?ΖSΔΖS+?y?φΔφ+?y?ωΔω+?y?κΔκ+?y?XΔX+?y?YΔY+?y?ΖΔΖ+?y?fxΔfx+?y?fyΔfy+?y?x0Δx0+?y?y0Δy0+?y?Κ1ΔΚ1+?y?Κ2ΔΚ2+?y?Ρ1ΔΡ1+?y?Ρ2ΔΡ2-ly(26)誤差方程式各項(xiàng)的系數(shù)均由共線方程求偏導(dǎo)數(shù)得到。若控制點(diǎn)可認(rèn)為無(wú)誤差時(shí),則相應(yīng)的誤差方程式中消去ΔX、ΔY、ΔZ項(xiàng)。在進(jìn)行光束法攝像機(jī)標(biāo)定時(shí),大多情況下法方程的狀態(tài)不良,因而需要對(duì)內(nèi)方位元素加以適當(dāng)?shù)臋?quán)進(jìn)行平差,以保證解算的穩(wěn)定性。3內(nèi)方位元素及光柵成像的實(shí)際標(biāo)定對(duì)于上文提出的數(shù)字?jǐn)z像機(jī)標(biāo)定算法,筆者用一組實(shí)際圖像進(jìn)行了實(shí)驗(yàn)。平面標(biāo)定格網(wǎng)大小約為45cm×45cm,共有900個(gè)設(shè)計(jì)格網(wǎng)點(diǎn),每個(gè)格網(wǎng)點(diǎn)的設(shè)計(jì)坐標(biāo)均用坐標(biāo)量測(cè)儀量測(cè),精度為0.1mm左右。實(shí)際上,對(duì)于自檢校光束法平差而言,格網(wǎng)點(diǎn)完全不需要如此高的精度,精確到0.5mm已經(jīng)足夠,而這對(duì)于工業(yè)制造等領(lǐng)域來(lái)說(shuō)很容易做到。CCD攝像機(jī)的像幅為1300像素×1030像素,標(biāo)定格網(wǎng)置于水平轉(zhuǎn)盤(pán)上,每次都采用不同的旋轉(zhuǎn)角度拍攝,共拍攝了8張像片,每張像片的可見(jiàn)格網(wǎng)點(diǎn)數(shù)約為500左右。每個(gè)可見(jiàn)格網(wǎng)點(diǎn)的像坐標(biāo)都通過(guò)直線特征提取后求交點(diǎn)得到,精度高于0.1像素。如果將平面格網(wǎng)點(diǎn)的量測(cè)坐標(biāo)作為已知值進(jìn)行標(biāo)定,考慮格網(wǎng)點(diǎn)像坐標(biāo)的誤差及空間坐標(biāo)量測(cè)誤差的綜合影響,標(biāo)定后精度必然低于0.3像素,無(wú)法滿(mǎn)足應(yīng)用要求。嚴(yán)密的自檢校光束法平差模型應(yīng)該將格網(wǎng)點(diǎn)的設(shè)計(jì)坐標(biāo)作為未知數(shù),并給以適當(dāng)?shù)臋?quán)值進(jìn)行平差,以精確求解攝像機(jī)參數(shù)及格網(wǎng)點(diǎn)空間坐標(biāo)。利用嚴(yán)密的自檢校光束法平差模型進(jìn)行標(biāo)定后的單位權(quán)中誤差為0.08像素,內(nèi)方位元素的標(biāo)定值及其實(shí)際精度(平差結(jié)果中給出的未知數(shù)中誤差)見(jiàn)表1,其中fx、fy分別為x、y兩個(gè)方向的焦距,(x0、y0)為主點(diǎn)