基于最小二乘配置的光學(xué)鏡頭畸變校正方法

何子清, 葛 超,王春陽

(中國白城兵器試驗中心，吉林 白城 137001)

1 引 言

短焦大視場光學(xué)測試設(shè)備在靶場應(yīng)用中具有舉足輕重的地位。但是由于電荷耦合器件制造誤差、鏡片曲面誤差、各鏡片間軸向安裝誤差以及對中誤差等，破壞了攝影中心、待檢測點及其像點間的共線關(guān)系，產(chǎn)生了畸變[1]；并隨著視場的增大，畸變值隨之增大，雖然不影響圖像的清晰度，卻嚴重影響設(shè)備的測角精度[2]，而測角精度正是檢核設(shè)備測試能力的重要技術(shù)指標，因此，不能忽視畸變對其產(chǎn)生的影響。為了提高光學(xué)設(shè)備的測角精度，必須采用正確合理的方法對光學(xué)鏡頭進行畸變校正，近些年，眾多學(xué)者對此做了大量的研究工作[3-10]。

鏡頭畸變主要有3種:徑向畸變、偏心畸變和薄棱鏡畸變，靶場光學(xué)鏡頭受徑向畸變影響較大，即離鏡頭中心越遠，變形越大。文獻[3]提出了一種無需昂貴的器材，便可完成高精度畸變校正的的方法，目前已在靶場使用，該方法以校正點附近若干已知點的畸變值作為觀測值，以二次多項式的系數(shù)作為估計參數(shù)，采用最小二乘估計得到修正模型[3]。該方法雖然能一定程度上抑制畸變產(chǎn)生的影響，但是仍然存在一些問題。比如，校正模型僅是針對局部有效，每個校正點的修正模型均不同，沒有整體連續(xù)性，無法充分利用已知數(shù)據(jù)的相關(guān)性；另外，這種模型只顧及了畸變的趨勢項，而沒有考慮隨機項。事實上，這種隨機項與概率統(tǒng)計意義上的隨機誤差不同，雖然在數(shù)值上呈隨機特性，但在光學(xué)鏡頭制作完成后就已經(jīng)固定了，它對測角精度的影響與在視場中的位置有關(guān)，因此，亦可稱為似系統(tǒng)誤差。這樣稱呼它是因為，它對已知校正點附近的點有系統(tǒng)性影響，對距離較遠的點則呈隨機性影響[11]。而文獻[3]中的方法將畸變誤差作為偶然誤差處理，以距離待校正點越近的網(wǎng)格點對其修正值影響越大為依據(jù)，雖符合一定的邏輯，但并沒有將誤差系統(tǒng)性和隨機性因素體現(xiàn)在模型中，這顯然是不合理的。因此，需要建立一種既能顧及畸變趨勢性影響，又能顧及隨機性影響的修正模型，另外，由于畸變值形成的曲面較復(fù)雜，尤其視場邊緣，相對起伏較大，二次多項式難以描述整個視場內(nèi)畸變值的變化趨勢，甚至可能出現(xiàn)“龍格現(xiàn)象”[12]。

為了提高光學(xué)設(shè)備的測角精度，本文提出一種基于多面函數(shù)的最小二乘配置光學(xué)鏡頭畸變校正方法，與其他校正方法不同，該方法既顧及畸變值的趨勢特性，又顧及其隨機特性，以多面函數(shù)代替二次多項式，能更好的反應(yīng)畸變值曲面的整體連續(xù)變化情況；根據(jù)擬合殘差求出協(xié)方差函數(shù)的估計參數(shù)。最后，利用本文方法和傳統(tǒng)方法分別對靶場使用的光學(xué)鏡頭進行校正處理，驗證了本文方法的有效性與可行性。

2 原理及數(shù)學(xué)模型

2.1 靶場畸變校正方法

目前，正在使用的修正模型如下：

Δα=ax2+by2+cxy+dx+ey+f，

(1)

其中，Δα為畸變值，x、y為編碼器讀數(shù)，a、b、c、d、e為系數(shù)。

將整個視場以網(wǎng)格的形式進行劃分，網(wǎng)格間距根據(jù)不同的設(shè)備取不同值，選取待校正點附近4×4網(wǎng)格點作為擬合節(jié)點，網(wǎng)格點對校正點修正值的影響程度被認為是不一樣的，依據(jù)到待校正點距離越近影響越大為原則，進行定權(quán)處理。由于網(wǎng)格點之間相互獨立，因此，權(quán)陣上的元素，除對角線外，其他均為0，對角線元素值為網(wǎng)格點與測量點距離的倒數(shù)。

根據(jù)式(1)建立誤差方程：

(2)

其中，A為系數(shù)矩陣。

在滿足VTPV=min條件下，建立法方程：

(3)

式中,P為權(quán)矩陣,求解上述法方程，得到函數(shù)的系數(shù)估計值，模型確定以后，代入編碼器讀數(shù)，即可得到對應(yīng)位置點的角度修正值。

2.2 基于最小二乘配置的畸變校正模型

根據(jù)光學(xué)鏡頭畸變的特點，建立校正模型：

(4)

式中，L為觀測向量，即網(wǎng)格點處的畸變值；AX為畸變值曲面連續(xù)的整體變化趨勢；A為趨勢項的系數(shù)矩陣；X為非隨機參數(shù)；如果選二次多項式作為變化趨勢項，AX與式(2)中一致；BY表示畸變值曲面不規(guī)則的隨機部分；S為已測點信號向量，即網(wǎng)格點處畸變值；S′為未測點信號向量，即待校正點處畸變值；Δ為隨機誤差。

Δ和Y的隨機特性為：

(5)

將式(4)改成誤差方程的形式：

(6)

按廣義最小二乘原理：

(7)

可得到：

(8)

其中：∑S為網(wǎng)格點的協(xié)方差矩陣，∑SS′為網(wǎng)格點與待校正點的互協(xié)方差矩陣，∑S′為待校正點的協(xié)方差矩陣。協(xié)方差矩陣中的元素是根據(jù)一個隨位置變化的協(xié)方差函數(shù)計算得到，精確確定協(xié)方差函數(shù)是最小二乘配置的關(guān)鍵。

在已有的經(jīng)驗協(xié)方差函數(shù)中，高斯函數(shù)應(yīng)用最廣，其表達式如式(9)所示:

σ2(d)=σ2(0)e-K2d2，

(9)

式中：σ2(0)表示距離為0的方差，K為待定的未知參數(shù)，d為視場內(nèi)兩點之間的距離。

根據(jù)如下公式計算σ2(0)和σ(di)：

(10)

多面函數(shù)法在解決不規(guī)則的連續(xù)曲面函數(shù)逼近問題有很好的效果，因此，本文選擇用多面函數(shù)代替?zhèn)鹘y(tǒng)的二次多項式對畸變值曲面進行擬合。

假設(shè)畸變值曲面是一個連續(xù)不規(guī)則變化的表面，利用多面函數(shù)對其進行逼近[12-13]：

(11)

式中：γi為待求系數(shù)，R為核函數(shù)，(xi,yi)為核函數(shù)的中心點。核函數(shù)表達式為：

R(x,y,xi,yi)=[(x-xi)2+(y-yi)2+δ2]ω，

(12)

式中：δ為平滑因子，它的數(shù)值大小影響核函數(shù)形狀；ω為非零實數(shù)。在本文中，核函數(shù)參數(shù)取值為：δ=0.1,ω=0.5。

將式(11)轉(zhuǎn)成矩陣形式，將多面函數(shù)融入最小二乘配置中，模型為：

L=Rγ+BY+Δ，

(13)

式中，Rγ表示畸變值曲面的趨勢項部分；R為n×m維核函數(shù)矩陣；γ為1×m維的待定系數(shù)向量；BY為畸變值曲面的隨機部分。

根據(jù)廣義最小二乘原理，解得最小二乘配置解為：

(14)

待校正點的修正值為：

2.3 協(xié)方差函數(shù)擬合

協(xié)方差函數(shù)模型的準確性是最小二乘配置的關(guān)鍵。為避免線性化帶來的模型誤差,文章對高斯函數(shù)做一些處理。

對式(9)兩邊取對數(shù)，即可得到如下形式：

(15)

建立平差模型：

(16)

通過最小二乘求得σ2(0)和K。

3 實驗數(shù)據(jù)分析與驗證

選用靶場某型炸點坐標測量系統(tǒng)的光學(xué)經(jīng)緯儀作為實驗對象，對其進行畸變校正。具體實驗步驟如下:圖1為在實驗室布設(shè)的場景圖，在平行光管一端放置12W聚光照明燈，作為無窮遠目標，盡量使燈光亮點在圖像上僅占1到2個像素。調(diào)整經(jīng)緯儀位置，使光軸盡量與平行光管主軸重疊，并將經(jīng)緯儀整平，水平角置零。將整個視場劃分為若干網(wǎng)格，轉(zhuǎn)動經(jīng)緯儀使光源位于網(wǎng)格點處，單幀采集圖像，取光源的像面坐標，并代入式(18)中，得到目標相對經(jīng)緯儀的方位角和垂直角[14]:

(18)

其中：α0為方位角，λ0為垂直角，f為焦距，α為目標相對經(jīng)緯儀的方位角，λ為目標相對經(jīng)緯儀的垂直角。

圖1 實驗場景圖Fig.1 Map of experimental scene

由于經(jīng)緯儀測角精度為0.5"，遠遠高于設(shè)備測角精度40"，因此，可假設(shè)經(jīng)緯儀測角誤差為零，通過式(18)和(19)聯(lián)合求得畸變值:

(19)

其中：Δαi為水平方向的畸變值，Δλi為垂直方向的畸變值，αi為目標相對經(jīng)緯儀的方位角，λi為目標相對經(jīng)緯儀的垂直角。

(a)水平方向(a) Horizontal direction

(b)垂直方向(b) Vertical direction圖2 畸變值曲面的Delaunay三角形網(wǎng)圖Fig.2 3D mesh surface figure of distortion

以水平角0°、垂直角90°處為原點，建立二維坐標系，水平方向為x軸，垂直方向為y軸，視場范圍為-21°

圖3 待校正點在視場內(nèi)的分布圖Fig.3 Distribution map of the points to be corrected in the field

分別用二次多項式局部最小二乘法和本文方法對光學(xué)經(jīng)緯儀設(shè)備進行畸變校正。在進行最小二乘配置前，需要確定經(jīng)驗協(xié)方差函數(shù)，以得到準確的協(xié)方差矩陣∑S、∑SS′、∑S′。通過1.3小節(jié)介紹的協(xié)方差擬合公式(16)，求得水平方向的σ2(0)和K為107、0.8，垂直方向的σ2(0)和K為109、0.78。

在鏡頭視場內(nèi)隨機選取300個點，如圖3所示，圖中紅色“×”代表待校正點，藍色“+”代表用于多面函數(shù)擬合的核函數(shù)中心點。

采用3種方案進行畸變校正，具體步驟如下：

方案1：基于二次多項式的局部最小二乘法，該法選取待校正點周圍4×4網(wǎng)格點數(shù)據(jù)作為插值節(jié)點擬合趨勢項，計算所有待校正點的修正值，將結(jié)果與實測畸變值進行求差，得到校正精度；

方案2：多面函數(shù)擬合法，該法選取357個網(wǎng)格點作為擬合點，用多面函數(shù)擬合畸變值曲面，計算所有待校正點的擬合值，將結(jié)果與實測畸變值進行求差，得到擬合精度。

方案3：基于多面函數(shù)的最小二乘配置法，該法選取357個網(wǎng)格點作為核函數(shù)中心點，用多面函數(shù)擬合趨勢項，計算所有待校正點的修正值，將結(jié)果與實測畸變值進行求差，得到校正精度。

通過以上3種方案對設(shè)備進行校正，3種方案的各項指標如表1所示，誤差統(tǒng)計如圖4所示。

表1 3種方案處理結(jié)果統(tǒng)計表

從表1可以得出，方案3的各項指標均優(yōu)于方案1和方案2。水平方向上，兩者的RMS值分別為31.3"、23.1"和11.7"，垂直方向上，分別為18.6"、11.2"和7.7"。相比于傳統(tǒng)的方法，采用多面函數(shù)擬合畸變值曲面，精度在兩個方向上分別提高8.2"和7.4"，說明多面函數(shù)作為曲面的趨勢項效果更好；相比于多面函數(shù)擬合，基于多面函數(shù)的最小二乘配置法的校正精度在兩個方向上分別提高11.4"和3.5"，說明去除趨勢項之后，最小二乘配置會對剩余的不規(guī)則隨機變化項進行一定的修正，使校正精度更高。以上數(shù)據(jù)證明了以多面函數(shù)擬合趨勢項更能反映畸變值的整體連續(xù)變化趨勢，去除趨勢項后得到的隨機變化部分也更為平穩(wěn)。另外，從誤差統(tǒng)計圖中可以看出，基于多面函數(shù)的最小二乘配置法的誤差明顯小于傳統(tǒng)方法和多面函數(shù)擬合法，進一步驗證了本文方法的有效性。

(a)水平方向(a) Horizontal direction

(b)垂直方向(b) Vertical direction圖4 3種方案待校正點誤差統(tǒng)計圖Fig.4 Error statistical chart of the points to be corrected in three method

從實驗數(shù)據(jù)的分析中可以總結(jié)出如下結(jié)論：

(1)采用二次多項式局部最小二乘法，對每個待校正點在小區(qū)域內(nèi)進行校正，存在無法充分利用數(shù)據(jù)之間的關(guān)聯(lián)性的問題。而采用基于多面函數(shù)的最小二乘配置法可以解決上述問題，它是以多面函數(shù)的形式，從整體的角度，反映畸變值的整體連續(xù)變化趨勢。

(2)傳統(tǒng)方法只顧及觀測值中的偶然誤差，而基于多面函數(shù)的最小二乘配置法不但能顧及畸變值曲面的趨勢項，而且也能顧及不規(guī)則變化的隨機項。因此，相比于傳統(tǒng)方法，以多面函數(shù)擬合趨勢項更能反映畸變值的整體連續(xù)變化趨勢，去除趨勢項后得到的隨機變化部分也更為平穩(wěn)。

(3)與傳統(tǒng)方法相比，基于多面函數(shù)的最小二乘配置法不但能較好的反映畸變值的整體連續(xù)變化趨勢，校正精度也明顯優(yōu)于傳統(tǒng)方法。

4 結(jié) 論

通過對光學(xué)測試設(shè)備進行畸變校正，驗證了本文方法的可行性與有效性。實驗結(jié)果證明：本文提出以多面函數(shù)代替?zhèn)鹘y(tǒng)的方法擬合畸變值曲面的變化趨勢，能達到更好的效果，另外，最小二乘配置法不但能顧及畸變值曲面的趨勢項，而且也能顧及不規(guī)則變化的隨機項，使待校正點的推估精度更高。