劉世望，胡云卿，林 軍

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

0 引言

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

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

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

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

1 激光攝像式傳感器標定模型

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

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

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

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

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

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

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

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

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

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

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

式中：n——標定特征點數(shù)量。

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

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

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

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

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

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

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

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

2.2 高斯-牛頓迭代法

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

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

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

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

2.3 基于極大似然估計的L-M 算法

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

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

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

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

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

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

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

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

式中：σ——標準差。

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

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

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

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

3 激光攝像式傳感器標定試驗

3.1 標定試驗

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

圖2 針形靶標標定平臺示意圖Fig. 2 Schematic diagram of the calibration platform for needle target diagram

激光攝像式傳感器針形靶標標定步驟如下：

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

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

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

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

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

圖3 針形靶標光斑圖像Fig. 3 Spot images of needle target

圖4 特征點坐標標定Fig. 4 Coordinate calibration of feature points

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

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

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

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

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

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

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

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

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

3.2 標定試驗結(jié)果誤差分析

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

圖7 800 組特征點世界坐標系下的誤差曲線Fig. 7 Error curves of 800 groups of characteristic points in world coordinate system

表5 800 組特征點坐標算法誤差 Tab. 5 Errors of 800 groups of feature points（單位：mm）

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

圖8 1 000 組未參與標定模型計算的數(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）所示；基于極大似然估計的L-M 優(yōu)化法得到的誤差(Δx6, Δy6)，分別如圖8（e）和圖8（f）所示。采用這3 種方法得到的誤差最大值分別為εmax4，εmax5，εmax6；誤差平均值分別為εave4，εave5，εave6；重投影誤差方均根值分別εw4，εw5，εw6，具體如表6 所示。

表6 1 000 組特征點坐標算法誤差 Tab. 6 Errors of 1 000 groups of feature points（單位：mm）

綜合表5 和表6，對兩組數(shù)據(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?？梢?，基于極大似然估計L-M 算法的標定精度誤差最大值較非線性最小二乘法的和高斯-牛頓法的分別減小了（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é)語

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