激光攝像式傳感器標(biāo)定方法研究

劉世望，胡云卿，林 軍

（中車株洲電力機(jī)車研究所有限公司，湖南 株洲 412001）

0 引言

激光攝像式傳感器標(biāo)定是實(shí)現(xiàn)軌道交通視覺檢測的重要前提之一。利用激光攝像技術(shù)進(jìn)行視覺檢測時(shí)，圖像坐標(biāo)系與世界坐標(biāo)系間存在多個(gè)自由度的變換，且攝像機(jī)鏡頭存在非線性畸變，所以在視覺檢測中對激光攝像式傳感器進(jìn)行高精度標(biāo)定的難點(diǎn)在于對激光圖像進(jìn)行多自由度的坐標(biāo)轉(zhuǎn)換及鏡頭畸變補(bǔ)償。

激光攝像式傳感器標(biāo)定常見的方法有拉絲標(biāo)定法[1-2]、鋸齒靶標(biāo)標(biāo)定法[3-4]和二次交比不變標(biāo)定法[5-8]，這些標(biāo)定方法主要存在以下問題：（1）結(jié)構(gòu)光平面的特征點(diǎn)提取困難；（2）結(jié)構(gòu)光平面特征點(diǎn)數(shù)量有限；（3） 計(jì)算步驟復(fù)雜。為解決這些方法的不足，本文采用一種基于針形靶標(biāo)的運(yùn)動(dòng)標(biāo)定方法，其將結(jié)構(gòu)光平面投射在針形靶標(biāo)上，形成近似圓形的光斑，其中心即為標(biāo)定所需提取的特征點(diǎn)；同時(shí)針形靶標(biāo)沿位移平臺(tái)按固定方向等間隔移動(dòng)，以獲取多組光斑特征點(diǎn)。該方法解決了結(jié)構(gòu)光平面的特征點(diǎn)提取困難和標(biāo)定特征點(diǎn)數(shù)量不足的問題。

文獻(xiàn)[9-14]針對視覺檢測實(shí)際應(yīng)用，忽略鏡頭畸變因素，采用攝像機(jī)線性模型進(jìn)行標(biāo)定，但沒有考慮激光攝像式傳感器的非線性畸變；文獻(xiàn)[15]側(cè)重于將激光攝像式傳感器應(yīng)用于圖像特征點(diǎn)提??；文獻(xiàn)[16]考慮鏡頭畸變，建立了非線性標(biāo)定模型，采用最小二乘法和高斯-牛頓迭代法進(jìn)行參數(shù)求解，但未給出詳細(xì)的推導(dǎo)過程。

針對上述現(xiàn)狀，本文綜合考慮激光攝像式傳感器鏡頭徑向、偏心及薄棱畸變問題[17-19]，采用基于針形靶標(biāo)的運(yùn)動(dòng)標(biāo)定方法，建立激光攝像式傳感器非線性標(biāo)定模型；針對模型中待求參數(shù)難求解的問題，采用極大似然估計(jì)法來降低圖像噪聲引起的誤差，然后用Levenberg-Marquardt（簡稱“L-M”）算法求出標(biāo)定參數(shù)最優(yōu)解，并與非線性最小二乘法、高斯-牛頓迭代法進(jìn)行對比，從而得到高精度的激光攝像式傳感器標(biāo)定方法。

1 激光攝像式傳感器標(biāo)定模型

激光攝像式傳感器成像及標(biāo)定模型如圖1 所示。以激光平面作為參考面，建立激光攝像式傳感器的世界坐標(biāo)系OW-XWYW。OW-XWYW平面與激光平面共面，OC-XCYCZC為攝像機(jī)坐標(biāo)系，O-uv為以圖像中心為原點(diǎn)的圖像坐標(biāo)系，同時(shí)以圖像平面左上角OE為坐標(biāo)原點(diǎn)建立坐標(biāo)系OE-XEYE，攝像機(jī)視場范圍內(nèi)任意一點(diǎn)(xw,yw,zw)對應(yīng)圖像坐標(biāo)系中坐標(biāo)(xE,yE)。

圖1 激光攝像式傳感器成像及標(biāo)定示意Fig. 1 Imaging and calibration diagram of laser camera sensor

設(shè)R和T為攝像機(jī)坐標(biāo)系轉(zhuǎn)化為世界坐標(biāo)系的旋轉(zhuǎn)矩陣和平移矩陣，即攝像機(jī)外部參數(shù)，如式（1）。

假設(shè)α,β,μ,ν,γ為攝像機(jī)內(nèi)部參數(shù)，其中α和β分別是u軸和v軸的尺度因子（或稱為有效焦距），（μ,ν）是光學(xué)中心，γ是u軸和v軸的不垂直因子。由于鏡頭畸變的存在，本文將鏡頭徑向畸變系數(shù)(k1和k2)、薄棱畸變系數(shù)(p1和p2)及偏心畸變系數(shù)(s1和s2)考慮進(jìn)去，建立非線性標(biāo)定模型。由攝影投射定理可知，OCXCYCZC與OE-XEYE坐標(biāo)系變換關(guān)系為

由于OW-XWYW平面與激光平面共面，標(biāo)定過程中參考點(diǎn)均在激光平面上，故可令zw=0。根據(jù)式（1）和式（2），可推出激光攝像式傳感器圖像坐標(biāo)系與世界坐標(biāo)系的關(guān)系：

式中：a1～a9為激光攝像式傳感器的攝像機(jī)線性模型參數(shù)。

激光平面上任意點(diǎn)P的圖像像素坐標(biāo)為(xE,yE)?？紤]鏡頭畸變后，實(shí)際成像像素坐標(biāo)為(xD,yD)，則畸變前后關(guān)系如下：

式中：δx與δy分別為XE方向和YE方向的畸變因子。

圖像中心O點(diǎn)像素坐標(biāo)為(x0,y0)，其與考慮畸變后的P點(diǎn)在圖像坐標(biāo)系O-uv下的坐標(biāo)(u,v)關(guān)系如下：

畸變因子δx與δy的表達(dá)式如式（6）所示，其中系數(shù)k1,k2,p1,p2,s1,s2被統(tǒng)稱為鏡頭畸變參數(shù)。

將式（4）和式（6）代入式（3），得到激光攝像式傳感器的非線性標(biāo)定模型，如式（7）所示，其傳感器標(biāo)定參數(shù)（a1～a8,k1,k2,p1,p2,s1及s2）的目標(biāo)函數(shù)如式（8）所示。

式中：n——標(biāo)定特征點(diǎn)數(shù)量。

2 標(biāo)定模型參數(shù)求解方法

視覺標(biāo)定過程中常用的參數(shù)求解方法主要有非線性最小二乘法、高斯-牛頓迭代算法以及L-M 算法。在非線性最小二乘法及高斯-牛頓迭代算法的求解過程中，會(huì)出現(xiàn)矩陣奇異或接近奇異的情形，使方程求解過程變得十分困難；而L-M 算法通過改變矩陣的特征值結(jié)構(gòu)來解決該問題，同時(shí)L-M 算法在求解過程中可以通過參數(shù)設(shè)置靈活調(diào)節(jié)收斂速度。

2.1 非線性最小二乘法[20]

由于式（8）是非線性方程組，無法直接通過線性最小二乘法求解，故通過非線性最小二乘法求解激光攝像式傳感器標(biāo)定參數(shù)，求解步驟如下：

（1）將傳感器標(biāo)定參數(shù)分成兩部分，即線性模型參數(shù)(a1,a2, …,a8)和畸變參數(shù)(k1,k2,p1,p2,s1,s2)。

（2）將畸變參數(shù)作為已知量代入式（8）并將非線性方程組轉(zhuǎn)化為線性方程組，采用線性最小二乘法求解線性模型參數(shù)。

（3）將步驟（2）中計(jì)算得到的線性模型參數(shù)作為已知量代入式（8），將式（8）轉(zhuǎn)化為關(guān)于畸變參數(shù)的線性方程組并采用線性最小二乘法求解畸變參數(shù)。

（4）重復(fù)步驟（2）和步驟（3），通過多次迭代，最終接近最優(yōu)解。

2.2 高斯-牛頓迭代法

高斯-牛頓迭代法是一種非線性優(yōu)化算法，用來求解激光攝像式傳感器的標(biāo)定參數(shù)，具體步驟如下：

（1）采用2.1 節(jié)中非線性最小二乘法迭代5 次，將求解得到的激光攝像式傳感器的標(biāo)定參數(shù)作為高斯-牛頓迭代的初值并根據(jù)式（8）建立高斯-牛頓優(yōu)化計(jì)算模型[21]：

利用式（10）可求得Guass-Newton 方向參數(shù)Di，其中F=[fx1,fy1, …,fxi,fyi]T。

（2）沿著Guass-Newton 方向Di作一維搜索，求得步長m使式（11）成立。

2.3 基于極大似然估計(jì)的L-M 算法

在高斯-牛頓迭代計(jì)算過程中，求解矩陣(ATA)-1時(shí)，容易出現(xiàn)矩陣ATA奇異或接近奇異的情形，這時(shí)求解式（10）得到的Guass-Newton 方向Di誤差較大，極易導(dǎo)致最終迭代計(jì)算結(jié)果無法收斂。

針對高斯-牛頓迭代法中矩陣ATA奇異或接近奇異的情形，L-M 算法通過增加正定對角矩陣的方式改變矩陣ATA的特征值結(jié)構(gòu)，使其變成條件數(shù)較好的對稱正定矩陣，彌補(bǔ)了高斯-牛頓迭代法的不足。

L-M 算法以高斯-牛頓迭代法為基礎(chǔ)，令

式中：I為i階單位矩陣；α為正實(shí)數(shù)。

當(dāng)α=0 時(shí)，Di即為Guass-Newton 方向；當(dāng)α足夠大時(shí)，逆矩陣(ATA+αI)-1主要取決于αI，Di接近最速下降方向。一般而言，α初值為0.001，變換因子β=10，允許誤差ε=0.000 01。當(dāng)F(X i+1)>F(X i)且|ATA|≤ε時(shí)，停止計(jì)算；否則，置α=αβ，繼續(xù)迭代。當(dāng)F(X i+1)＜F(X i)且|ATA|≤ε時(shí)，停止計(jì)算；否則，置α=α/β，繼續(xù)迭代。

本算法提出先采用極大似然估計(jì)法求得由高斯噪聲和散粒噪聲引起的圖像坐標(biāo)誤差，進(jìn)而在標(biāo)定模型參數(shù)求解過程中剔除高斯噪聲和散粒噪聲的干擾，將實(shí)際成像像素坐標(biāo)減去由高斯噪聲和散粒噪聲引起的像素坐標(biāo)誤差，然后用L-M 法求出參數(shù)(a1, …,a8,k1, …,s2)的最優(yōu)解。

設(shè)激光平面上任意點(diǎn)P的圖像像素坐標(biāo)為(xE,yE)，考慮鏡頭畸變后，實(shí)際成像像素坐標(biāo)為(xD,yD)。實(shí)際上，圖像坐標(biāo)的提取過程會(huì)受到高斯噪聲、散粒噪聲的干擾[22-23]。一般而言，高斯噪聲引起的像素坐標(biāo)誤差服從正態(tài)分布，而散粒噪聲引起的像素坐標(biāo)誤差服從伽馬分布，誤差(εx,εy)如式（14）所示。

任意點(diǎn)P對應(yīng)的由高斯噪聲引起的像素坐標(biāo)誤差概率密度函數(shù)f(p1)如式（15）所示。

式中：σ——標(biāo)準(zhǔn)差。

任意點(diǎn)P對應(yīng)的由散粒噪聲引起的像素坐標(biāo)誤差概率密度函數(shù)f(p2)如下：

構(gòu)造極大似然函數(shù)L(p1)和L(p2)，其等于所有點(diǎn)(p1,p2, …,pb)對應(yīng)的像素坐標(biāo)誤差概率密度函數(shù)乘積，具體如式（17）和式（18）所示。

根據(jù)極大似然原理[23-25]，當(dāng)L(p1)和L(p2)取最大值時(shí)，所得的估計(jì)值即為像素坐標(biāo)誤差值；根據(jù)式（17）和式（18）可知，當(dāng)L(p1)和L(p2)取最大值時(shí)，式（19）取值最小。

由激光攝像式傳感器模型可知，式（19）（目標(biāo)函數(shù)）為非線性函數(shù)，采用L-M 算法對該目標(biāo)函數(shù)進(jìn)行求解，容易得到目標(biāo)函數(shù)的極大似然估計(jì)值，也就可以得到像素坐標(biāo)誤差值。

3 激光攝像式傳感器標(biāo)定試驗(yàn)

3.1 標(biāo)定試驗(yàn)

針形靶標(biāo)標(biāo)定原理如圖2 所示，采用針形靶標(biāo)對標(biāo)定平臺(tái)、刻度尺、滑動(dòng)絲桿、計(jì)算機(jī)進(jìn)行標(biāo)定，再通過數(shù)據(jù)采集軟件和Matlab 實(shí)現(xiàn)對標(biāo)定數(shù)據(jù)的采集和處理。

圖2 針形靶標(biāo)標(biāo)定平臺(tái)示意圖Fig. 2 Schematic diagram of the calibration platform for needle target diagram

激光攝像式傳感器針形靶標(biāo)標(biāo)定步驟如下：

步驟一，安裝激光攝像式傳感器。調(diào)整傳感器在工裝固定臺(tái)的位置，確保結(jié)構(gòu)光始終垂直投射于針形靶標(biāo)上，再將激光攝像式傳感器固定在平臺(tái)上。

步驟二，攝像機(jī)捕捉結(jié)構(gòu)光平面投射于針形靶標(biāo)形成的光斑圖像，并提取光斑中心作為特征點(diǎn)的圖像坐標(biāo)，同時(shí)記錄光斑在世界坐標(biāo)系下的世界坐標(biāo)。

步驟三，以相同間距移動(dòng)針形靶標(biāo)，每次移動(dòng)后重復(fù)步驟二，盡量使所得特征點(diǎn)覆蓋攝像機(jī)全部視場范圍。

步驟四，根據(jù)上述的標(biāo)定模型參數(shù)求解方法，將獲取的特征點(diǎn)對代入計(jì)算模型，獲得最優(yōu)解。

結(jié)構(gòu)光平面投射在針形靶標(biāo)上形成的光斑圖像如圖3 所示，圖中僅列出了針形靶標(biāo)在位移平臺(tái)5 個(gè)不同位置的光斑圖像。選取光斑圖像的中心作為標(biāo)定點(diǎn)對的特征點(diǎn)圖像坐標(biāo)，特征點(diǎn)世界坐標(biāo)則為結(jié)構(gòu)光與針形靶針相交位置距離固定平臺(tái)的水平距離和垂直距離。等間隔地移動(dòng)靶標(biāo)，共獲取800 個(gè)標(biāo)定特征點(diǎn)對，特征點(diǎn)的圖像坐標(biāo)如圖4（a）所示，世界坐標(biāo)如圖4（b）所示。

圖3 針形靶標(biāo)光斑圖像Fig. 3 Spot images of needle target

圖4 特征點(diǎn)坐標(biāo)標(biāo)定Fig. 4 Coordinate calibration of feature points

針對實(shí)驗(yàn)獲得的800 組標(biāo)定數(shù)據(jù)，采用非線性最小二乘法迭代5 次，發(fā)現(xiàn)傳感器標(biāo)定模型參數(shù)(a1, …,a8,k1,k2,p1,p2,s1,s2)的值比較穩(wěn)定，迭代5 次得到的傳感器標(biāo)定模型參數(shù)值如表1、表2 所示。

表1 標(biāo)定線性模型參數(shù)Tab. 1 Calibration linear model parameters

表2 標(biāo)定鏡頭畸變參數(shù)Tab. 2 Calibration lens distortion parameters

將傳感器標(biāo)定模型參數(shù)作為基于極大似然估計(jì)的L-M 優(yōu)化算法的初值，通過L-M 算法尋找傳感器標(biāo)定模型參數(shù)最優(yōu)值，標(biāo)定線性模型參數(shù)收斂曲線如圖5 所示，標(biāo)定鏡頭畸變參數(shù)收斂曲線如圖6 所示。

圖5 標(biāo)定線性模型參數(shù)收斂曲線Fig. 5 Convergence curves of calibration linear model parameters

圖6 標(biāo)定鏡頭畸變參數(shù)收斂曲線Fig. 6 Convergence curves of calibration lens distortion parameters

迭代10 步之后，參數(shù)基本收斂，此時(shí)所得到的標(biāo)定參數(shù)為最優(yōu)解，如表3 和表4 所示。

表3 標(biāo)定鏡頭畸變參數(shù)Tab. 3 Calibration lens distortion parameters

表4 標(biāo)定線性模型參數(shù)Tab. 4 Calibration linear model parameters

3.2 標(biāo)定試驗(yàn)結(jié)果誤差分析

針對標(biāo)定試驗(yàn)中得到的標(biāo)定數(shù)據(jù)，分別采用非線性最小二乘法、高斯-牛頓迭代法和基于極大似然估計(jì)的L-M 優(yōu)化法這3 種方法進(jìn)行處理，得到模型參數(shù)(a1, …,a8,k1, …,s2)的最優(yōu)值，然后將參數(shù)代回標(biāo)定模型中，進(jìn)行逆向運(yùn)算。將800 組特征點(diǎn)的像素坐標(biāo)分別代入標(biāo)定模型，分別求得特征點(diǎn)在世界坐標(biāo)系下的坐標(biāo)值(xw1,yw1), (xw2,yw2), (xw3,yw3)。以原始世界坐標(biāo)(xw,yw)為標(biāo)準(zhǔn)，(Δx, Δy)為世界坐標(biāo)系下的誤差，采用非線性最小二乘法得到的誤差為(Δx1, Δy1)=(xw-xw1,yw-yw1)，分別如圖7（a）和圖7（b）所示；采用高斯-牛頓迭代法得到的誤差為(Δx2, Δy2)=(xw-xw2,yw-yw2)，分別如圖7（c）和圖7（d）所示；采用基于極大似然估計(jì)的L-M 優(yōu)化法得到的誤差為(Δx3, Δy3)=(xw-xw3,yw-yw3)，分別如圖7（e）和圖7（f）所示。

圖7 800 組特征點(diǎn)世界坐標(biāo)系下的誤差曲線Fig. 7 Error curves of 800 groups of characteristic points in world coordinate system

表5 800 組特征點(diǎn)坐標(biāo)算法誤差 Tab. 5 Errors of 800 groups of feature points（單位：mm）

為確保實(shí)驗(yàn)的可靠性和算法的有效性，另獲取1 000 組未參與標(biāo)定模型計(jì)算的數(shù)據(jù)進(jìn)行誤差分析（圖8），采用3 種方法分別計(jì)算得到該組數(shù)據(jù)的世界坐標(biāo)為(xw4,yw4)，(xw5,yw5)和(xw6,yw6)。

圖8 1 000 組未參與標(biāo)定模型計(jì)算的數(shù)據(jù)誤差分析Fig. 8 Error analysis of 1 000 groups of data not involved in calibration model calculation

非線性最小二乘法得到的誤差為(Δx4, Δy4)，如圖8（a）和圖8（b）所示；高斯-牛頓迭代法得到的誤差為(Δx5, Δy5)，分別如圖8（c）和圖8（d）所示；基于極大似然估計(jì)的L-M 優(yōu)化法得到的誤差(Δx6, Δy6)，分別如圖8（e）和圖8（f）所示。采用這3 種方法得到的誤差最大值分別為εmax4，εmax5，εmax6；誤差平均值分別為εave4，εave5，εave6；重投影誤差方均根值分別εw4，εw5，εw6，具體如表6 所示。

表6 1 000 組特征點(diǎn)坐標(biāo)算法誤差 Tab. 6 Errors of 1 000 groups of feature points（單位：mm）

綜合表5 和表6，對兩組數(shù)據(jù)進(jìn)行平均計(jì)算，得到非線性最小二乘法、高斯-牛頓法、基于極大似然估計(jì)的L-M 算法誤差最大值分別是（0.493, 0.584），（0.436, 0.586），（0.081, 0.112）；誤差平均值分別是（0.156, 0.121），（0.147, 0.120），（0.011, 0.012）；重投影誤差方均根值分別是0.245 7, 0.236 6, 0.024 3。可見，基于極大似然估計(jì)L-M 算法的標(biāo)定精度誤差最大值較非線性最小二乘法的和高斯-牛頓法的分別減小了（0.412 mm, 0.472 mm）和（0.355 mm, 0.474 mm）； 誤差平均值分別減小了（0.145 mm, 0.109 mm）和（0.136 mm, 0.108 mm）；重投影誤差方均根值分別減小了0.221 4 mm 和0.212 3 mm。

4 結(jié)語

本文建立了基于非線性最小二乘法、高斯-牛頓迭代法及基于極大似然估計(jì)的L-M 優(yōu)化算法的激光攝像式傳感器標(biāo)定模型；并采用3 種不同方法進(jìn)行了標(biāo)定模型參數(shù)求解。其中非線性最小二乘法迭代300 次后模型參數(shù)基本收斂，高斯-牛頓迭代法迭代20 次后模型參數(shù)基本收斂，基于極大似然估計(jì)的L-M 算法迭代10 次后模型參數(shù)基本收斂。對于基于極大似然估計(jì)的L-M算法，因考慮高斯噪聲與散粒噪聲影響，對已有標(biāo)定算法進(jìn)行了優(yōu)化，使得標(biāo)定精度有效提升；但該算法的計(jì)算復(fù)雜度相對較高，后續(xù)研究中需繼續(xù)優(yōu)化。