基于改進風驅(qū)動算法的攝像機標定優(yōu)化方法

任久斌，曹中清

(西南交通大學 機械工程學院，四川 成都 610031)

0 引 言

現(xiàn)有的攝像機標定方法主要有直接線性變換法、Tsai兩步標定法和張正友標定法[1-4]。由于張正友標定法使用的是平面標定板，操作簡單，成本低廉，所以它在實際攝像機標定中被廣泛應(yīng)用。但因其采用的非線性優(yōu)化算法具有局限性，常常陷入局部最優(yōu)解，所以得到的攝像機標定結(jié)果精度較差。近年來，許多研究者提出了基于遺傳算法(GA)、模擬退火算法(SA)、粒子群算法(PSO)等的改進方法，盡管這些算法在一定程度上提升了標定精度，但是其改進效果并不明顯[5-8]。

風驅(qū)動優(yōu)化算法是近年來提出的一種全局優(yōu)化算法，它具有開發(fā)能力強、收斂速度快、尋優(yōu)效率高的特點，非常適合于求解多維多模態(tài)的問題[9-11]。然而基本風驅(qū)動算法中固有參數(shù)的人為設(shè)置非常困難，不同的設(shè)置值通常會帶來差異性的結(jié)果?；诖耍疚奶岢龌趨f(xié)方差矩陣自適應(yīng)進化策略的改進算法，通過利用協(xié)方差矩陣中不同維度的約束關(guān)系，實現(xiàn)參數(shù)采樣空間的選擇和重組，確保參數(shù)的采樣分布朝著適應(yīng)度值減小的方向移動。在獲得自適應(yīng)選擇的固有參數(shù)后，再利用基本風驅(qū)動算法優(yōu)化攝像機內(nèi)參。實驗結(jié)果表明，改進算法的收斂速度、收斂精度以及穩(wěn)定性都優(yōu)于其它算法。

1 攝像機成像模型

為了表達場景中的目標點與其投影點之間的關(guān)系，需要首先建立攝像機的成像模型[3]。攝像機成像模型中最常用的是線性模型(即針孔模型)，如圖1所示。圖中Ow-XwYwZw表示世界坐標系，OC-XCYCZC表示攝像機坐標系，O0-uv表示圖像像素坐標系，O1-xy表示圖像物理坐標系。

圖1 攝像機線性模型

假設(shè)空間點P在世界坐標系中的坐標為(Xw,Yw,Zw)，其在圖像像素坐標系中的投影點p的坐標為(u,v)，它們之間的關(guān)系可以表示為

(1)

式中：αx=f/dx，αy=f/dy，它們分別表示攝像機在x方向和y方向上以像素為單位的焦距，γ為坐標軸的不垂直度系數(shù)；M1為3×3的矩陣，由攝像機的內(nèi)部參數(shù)αx、αy、u0、v0、γ構(gòu)成，稱為攝像機的內(nèi)參矩陣；M2為3×4的矩陣，由旋轉(zhuǎn)矩陣R和平移向量t構(gòu)成，描述了攝像機坐標系相對于世界坐標系的關(guān)系，稱為攝像機的外參矩陣；M為3×4的矩陣，即內(nèi)參矩陣和外參矩陣的乘積，稱為攝像機的投影矩陣。

(2)

其中，δx和δy分別表示圖像中某一點在x方向和y方向的畸變值。它們的形式如式所示

(3)

式(3)中3項分別是徑向畸變、切向畸變和薄棱鏡畸變，p1、p2、k1、k2、s1、s2為非線性畸變系數(shù)。本文只考慮徑向畸變和切向畸變的前兩項，即k1、k2、p1、p2，因為這兩項足以描述高精度測量中的透鏡畸變。

2 攝像機模型參數(shù)初值的求解

根據(jù)張正友的攝像機標定方法，為了得到高精度的標定結(jié)果，需要首先為優(yōu)化算法計算提供內(nèi)參初值。通過將世界坐標系固定在平面標定板上，令標定板平面為Zw=0的平面，則上述攝像機的成像模型可以簡化為

(4)

其中，r1和r2為旋轉(zhuǎn)矩陣R的第一列和第二列。再利用圖像平面和標定板平面間的單應(yīng)性關(guān)系可得

(5)

(6)

其中，λ為尺度因子。由于旋轉(zhuǎn)矩陣的兩個列向量r1和r2正交，于是有

(7)

然后再將式(6)和式(7)結(jié)合，整理變換后可得如下兩個約束

(8)

由于攝像機內(nèi)參矩陣有5個參數(shù)，而每個單應(yīng)性矩陣只能夠提供兩個約束，因此至少需要3個不同的單應(yīng)性矩陣才能求解，這里可以通過攝像機采集至少3幅不同方位的標定板圖像實現(xiàn)[8]。

為了便于計算，令

(9)

由式(9)可以看到B是一個對稱矩陣，它可以用一個6維的向量表示為

(10)

hiTBhj=vijTb

(11)

其中，vij為

(12)

由式(8)的約束關(guān)系可得

(13)

在采集到標定板的n(n≥3)幅不同圖像后，利用最小二乘法可以求解出b，進而求解得到攝像機內(nèi)參的初始值。在得到單應(yīng)性矩陣和內(nèi)參矩陣的基礎(chǔ)上，還可以求得攝像機的外部參數(shù)

(14)

3 基于改進風驅(qū)動算法的攝像機參數(shù)優(yōu)化

3.1 基本風驅(qū)動優(yōu)化算法

風驅(qū)動優(yōu)化算法(wind driven optimization，WDO)是由Bayraktar提出的一種基于群體的啟發(fā)式迭代搜索算法，在求解全局優(yōu)化問題時具有簡單易行和高效的特點。Bayraktar首先通過對地球表面對流層空氣運動狀況的分析，抽象出簡化的大氣模型，然后結(jié)合牛頓第二定律和理想氣體狀態(tài)方程，推導出WDO算法中空氣微團的速度更新公式和位置更新公式[11]。

對某一單位質(zhì)量的空氣微團受力分析，由牛頓第二定律可得

ρa=∑Fi

(15)

式中：ρ為空氣的密度；a為空氣微團運動的加速度，其大小為a=Δu/Δt，Δu表示Δt時間內(nèi)的速度變化；Fi為作用在空氣微團上的力。

由空氣動力學知識，式(15)中主要包含了4種主要的作用力，使得空氣微團以一定速度向某個方向移動或者使它偏離現(xiàn)有路徑，這4種力分別是：氣壓梯度力(FPG)、摩擦力(FF)、重力(FG)和科里奧利力(FC)。利用它們的簡化表示，式(15)可以變換為

(16)

氣壓梯度的變化是因為空氣溫度差異引起的，為了將空氣微團的密度和溫度聯(lián)系起來，可以利用如下所示的理想氣體狀態(tài)方程

P=ρRT

(17)

其中，P為大氣壓強，R為理想氣體常數(shù)，T為氣體的溫度。

為了簡化式(16)，令δV=1， Δt=1，并將式(17)代入式(16)后得

(18)

式中：Pcur表示當前位置空氣微團的壓力值；Δu=unew-ucur，ucur表示空氣微團當前迭代的速度，unew表示空氣微團下一次迭代的速度。

氣壓梯度是使空氣微團從當前位置xcur移向最優(yōu)位置xopt的原因，其大小為|Popt-Pcur|，因此氣壓梯度可表示為

(19)

其中，Popt表示最優(yōu)位置的壓力值。將重力加速度矢量g改寫為g=|g|(0-xcur)，代入到式(18)中得

(20)

(21)

其中，c=-2|Ω|RT。式(21)中使用了實際的壓力值，當壓力值非常大的時候，空氣微團更新的速度也會變得不切實際的大，從而削弱了WDO的搜索性能。取而代之的是使用基于壓力值排序的方法，即通過對所有空氣微團的壓力值按升序排序，用排序后當前位置空氣微團的序號代替上述公式中的壓力值，這樣式(21)可以改寫為

(22)

式中：i代表當前位置xcur空氣微團壓力值Pcur排序后的序號，最優(yōu)位置xopt處的最優(yōu)壓力值Popt排序序號始終為1。

在得到空氣微團的速度更新方程以后，根據(jù)式(23)可以得出空氣微團的位置更新方程

xnew=xcur+(unewΔt)

(23)

式中：通常假設(shè)Δt=1。

式(22)和式(23)即為最終得到的基本W(wǎng)DO算法的速度更新方程和位置更新方程。

3.2 協(xié)方差矩陣自適應(yīng)進化策略

基本W(wǎng)DO算法中包含了4個參數(shù)α、g、RT和c，不同的優(yōu)化問題通常需要將它們設(shè)置為不同的值，而且在整個優(yōu)化迭代過程中它們固定不變[11]。為了減輕手動設(shè)置的負擔，避免人為選擇參數(shù)帶來的不確定性，本文提出基于協(xié)方差矩陣自適應(yīng)進化策略(CMAES)自動地選擇WDO算法的4個參數(shù)。

CMAES將α、g、RT和c作為待優(yōu)化參數(shù)，群體的采樣分布用協(xié)方差矩陣描述，通過自適應(yīng)地改變采樣分布的位置和形狀，引導搜索群體朝著適應(yīng)度值減小的方向進化，直到達到最大迭代次數(shù)為止[13-15]。CMAES算法的具體操作步驟如下：

(1)參數(shù)初始化。設(shè)置群體的大小λ，重組操作的子群大小μ=λ/2，重組權(quán)值wi(i=1,2,…,μ)以及自適應(yīng)調(diào)整所需的常量cσ、μeff、dσ、cc、c1見表1，設(shè)置算法的最大迭代次數(shù)G和問題的維度N(N=4，即待優(yōu)化的4個參數(shù))。在搜索空間中隨機選擇采樣分布的初始均值m(0)，采樣初始步長σ(0)設(shè)置為最大參數(shù)范圍的1/3，設(shè)置初始步長進化路徑Pσ(0)和初始協(xié)方差矩陣進化路徑Pc(0)為0，將初始協(xié)方差矩陣C(0)設(shè)置為單位矩陣I，進化代數(shù)g的初值設(shè)置為0。

表1 CMAES內(nèi)部常量參數(shù)的計算公式

(2)采樣操作。按照下式從采樣分布空間中抽取λ個個體形成搜索種體

(24)

(25)

(3)選擇和重組操作。計算群體中所有個體的適應(yīng)度值(即空氣微團的壓力值)，并對其進行升序排序，選擇適應(yīng)度較小的前μ個個體形成子群，然后對子群個體加權(quán)求和更新采樣分布的均值

(26)

(4)采樣步長調(diào)整

(27)

(28)

其中，協(xié)方差矩陣的逆(C(g))-1/2=B(g)(D(g))-1(B(g))T。

(5)協(xié)方差矩陣調(diào)整，引導群體搜索全局最優(yōu)解。先對協(xié)方差矩陣進化路徑Pc進行更新

(29)

Pc的進化方向保證了采樣分布沿均值更新方向移動。協(xié)方差矩陣自適應(yīng)更新公式為

(30)

其中，第三項利用了適應(yīng)度較好的μ個個體，因此采樣分布朝著適應(yīng)度減小的方向進化[15]。

(6)若g達到最大迭代次數(shù)G，則退出循環(huán)；否則，g=g+1，并跳轉(zhuǎn)到步驟(2)繼續(xù)執(zhí)行。

3.3 攝像機內(nèi)參優(yōu)化過程

為了避免過多的參數(shù)優(yōu)化導致算法尋優(yōu)能力下降，本文采用外參固定的優(yōu)化模型，在第2節(jié)計算的初值基礎(chǔ)之上，優(yōu)化攝像機的內(nèi)部參數(shù)。假設(shè)采集了n幅標定板圖像，每幅圖像有m個特征點，定義優(yōu)化的目標函數(shù)(即重投影誤差)如下所示

(31)

基于改進風驅(qū)動算法的攝像機內(nèi)參優(yōu)化過程如下：

(1)算法初始化。設(shè)置空氣微團種群的大小M，優(yōu)化問題的維度D(即待優(yōu)化的9個攝像機內(nèi)參：αx、αy、u0、v0、γ、k1、k2、p1、p2)，最大進化代數(shù)Gen，將當前進化代數(shù)gen初始化為0。在攝像機內(nèi)參初值的基礎(chǔ)上產(chǎn)生空氣微團的位置種群和速度種群，并為每個空氣微團隨機初始化一個參數(shù)(α、g、RT和c)，所有參數(shù)組成參數(shù)種群；

(2)利用目標函數(shù)(31)評估初始空氣微團種群的適應(yīng)度值，并對它們進行升序排序，記錄下適應(yīng)度值最小的個體gbest；

(3)根據(jù)參數(shù)種群提供的WDO算法參數(shù)，由式(22)和式(23)更新空氣微團種群中每個個體的速度和位置，再次利用式(31)計算種群更新后的適應(yīng)度值并排序，更新全局最優(yōu)個體gbest；

(4)利用CMAES自適應(yīng)更新參數(shù)種群，具體實現(xiàn)如3.2節(jié)所示，其中參數(shù)種群的大小λ和空氣微團種群大小M相同，維度N=4；

(5)令gen=gen+1，若gen

4 實驗與結(jié)果分析

為了驗證本文改進算法的性能，下面進行了攝像機標定實驗。首先，用攝像機從不同方位拍攝棋盤格標定板(9×12，每個方格大小為10 mm)的20幅圖像，由于世界坐標系固定在標定板上(Zw=0)，因此可以很方便地計算出標定板上每個角點的世界坐標，而對于標定板圖像中角點的像素坐標，本文利用Harris角點檢測算法提取。然后再根據(jù)第2節(jié)所描述的方法計算出攝像機內(nèi)參的初值，所得結(jié)果見表2。

表2 攝像機內(nèi)參初值

在得到各個參數(shù)的初值后，利用下面的方式對它們進行簡單地擴展，以此創(chuàng)建每個參數(shù)的搜索邊界：[αxαyu0v0]±100，γ∈[-1,1]， [k1k2p1p2]∈[-2,2]。實驗中，本文設(shè)置的種群大小為50，最大迭代次數(shù)為100。

接下來分別利用張正友標定法(Zhang)、粒子群算法(PSO)、基本風驅(qū)動算法(WDO)和本文改進的風驅(qū)動算法(CMAES-WDO)進行實驗，各算法分別運行50次，將各次所得的參數(shù)優(yōu)化結(jié)果和重投影誤差值取平均，結(jié)果見表 3。

從表3中可以看出，張正友標定法的重投影誤差均值為0.217 69，經(jīng)過PSO算法優(yōu)化的標定重投影誤差均值為0.179 32，經(jīng)過基本W(wǎng)DO算法優(yōu)化的標定重投影誤差均值為0.158 12，而經(jīng)過本文改進的風驅(qū)動算法的標定重投影誤差均值僅為0.079 61。多次實驗的結(jié)果表明，盡管經(jīng)過基本W(wǎng)DO算法優(yōu)化標定的攝像機內(nèi)參精度有所改善，但因為是人為選擇的算法參數(shù)，所以改進的程度并不高。

表3 各算法的優(yōu)化結(jié)果

圖2～圖5中展示了4種方法的標定重投影誤差分布。從圖2和圖3的對比可以看出，PSO優(yōu)化標定的重投影誤差分布相較于Zhang標定法的重投影誤差分布更加密集，說明PSO算法在一定程度上能夠提高攝像機的標定精度。從圖3和圖4的對比可以看出，基本W(wǎng)DO算法相較于PSO算法的性能更優(yōu)，但標定精度的改善程度并不大。從圖5中可以看出，經(jīng)過CMAES-WDO算法優(yōu)化標定后的重投影誤差在x方向和y方向均集中在[-0.15,0.15]之間，和其它3種方法相比其分布密度最小，所以它最終得到的攝像機標定參數(shù)精度最高。圖6中展示了PSO算法、基本W(wǎng)DO算法和CMAES-WDO算法用于攝像機標定優(yōu)化時目標值的變化過程，可以發(fā)現(xiàn)CMAES-WDO算法迭代約10次就可以完全收斂，而且收斂精度約為0.08，相較于其它兩種算法，CMAES-WDO算法的收斂速度和收斂精度都更優(yōu)。

圖2 Zhang標定法的重投影誤差分布

圖3 PSO算法的重投影誤差分布

圖4 WDO算法的重投影誤差分布

圖5 CMAES-WDO算法的重投影誤差分布

圖6 目標值優(yōu)化過程

5 結(jié)束語

本文利用協(xié)方差矩陣自適應(yīng)進化策略自動地確定了基本風驅(qū)動算法中的固有參數(shù)，避免了因為人為參數(shù)選擇帶來的差異性結(jié)果，并將改進后的算法用于攝像機參數(shù)的標定優(yōu)化中。實驗結(jié)果表明，改進風驅(qū)動算法相較于基本風驅(qū)動算法的性能更好，攝像機參數(shù)標定優(yōu)化后的精度更高，得到的攝像機內(nèi)參結(jié)果更準確。通過與粒子群算法和基本風驅(qū)動算法對比，改進算法的收斂速度更快，收斂精度更高。此外，多次重復實驗的結(jié)果表明，改進算法的穩(wěn)定性也更好。