□ 李 楠

長(zhǎng)安大學(xué) 汽車學(xué)院 西安 710064

1 噴霧錐角定義

汽車發(fā)動(dòng)機(jī)的燃油霧化效果直接影響發(fā)動(dòng)機(jī)的經(jīng)濟(jì)性、動(dòng)力性能及排放性能[1]。

目前,發(fā)動(dòng)機(jī)噴霧錐角的定義比較多,研究人員普遍采用Naber提出的噴霧錐角定義:在油束貫穿距離的1/2處作水平直線,與油束寬度邊緣交于B、C兩點(diǎn),與噴嘴A處形成的∠BAC即為噴霧錐角,如圖1所示。

油束貫穿距離定義為在開(kāi)始噴油后指定時(shí)刻的噴霧邊緣與噴油器頂端之間沿噴油器軸線方向的最大距離,即圖1中的L。

圖1 噴霧錐角示意圖

2 現(xiàn)有測(cè)量方法

文獻(xiàn)[2]通過(guò)圖1中B、C兩點(diǎn)坐標(biāo)與噴嘴A處坐標(biāo)構(gòu)成一個(gè)三角形,利用余弦定理計(jì)算出噴霧錐角。由于只考慮三個(gè)點(diǎn),并沒(méi)有考慮油束整個(gè)邊緣,因此隨著油束長(zhǎng)度增加,這一測(cè)量方法的誤差較大。文獻(xiàn)[3]通過(guò)對(duì)噴霧圖像進(jìn)行直方圖均衡化,增強(qiáng)油束與背景的對(duì)比度,用自適應(yīng)濾波函數(shù)對(duì)圖像進(jìn)行平滑處理,從而去除噪聲,再用Otsu法二值化處理,用坎尼算子進(jìn)行邊緣檢測(cè),結(jié)合霍夫變換擬合油束邊界線,計(jì)算兩條直線的斜率,從而得出噴霧錐角。這一方法雖然考慮了油束的邊界整體特性,但是由于前期預(yù)處理會(huì)損失油束邊緣,而且霍夫變換參數(shù)難以控制,因此導(dǎo)致擬合出的直線不止一條,測(cè)量精度不高。

針對(duì)上述問(wèn)題,筆者提出一種改進(jìn)的汽車發(fā)動(dòng)機(jī)噴霧錐角測(cè)量算法。這一算法根據(jù)原噴霧圖像的灰度直方圖特點(diǎn),將像素值不為0的點(diǎn)都提取出來(lái),全部置為1,再通過(guò)連通域、孔填充等方法完整分割出油束區(qū)域,并結(jié)合卡爾曼濾波獲取邊界數(shù)據(jù),用最小二乘法計(jì)算出噴霧錐角。

3 噴霧圖像獲取

噴霧圖像的獲取方法主要有機(jī)械測(cè)量法、電子測(cè)量法、光學(xué)測(cè)量法。機(jī)械測(cè)量法與電子測(cè)量法在測(cè)量精度與穩(wěn)定性方面都不如光學(xué)測(cè)量法[1],因此筆者采用光學(xué)測(cè)量法,利用高速攝像機(jī)拍攝噴霧場(chǎng)。噴霧圖像拍攝系統(tǒng)如圖2所示。

圖2 噴霧圖像拍攝系統(tǒng)

噴霧圖像拍攝系統(tǒng)主要設(shè)備有高壓油泵、電磁閥、高速攝像機(jī)、可控燃油噴射系統(tǒng)等,在一個(gè)標(biāo)準(zhǔn)大氣壓下進(jìn)行噴射。高速攝像機(jī)采用MotionXtra HG_100K型高速電荷耦合器件(CCD)數(shù)字?jǐn)z像機(jī),如圖3所示。

圖3 MotionXtra HG_100K型高速CCD數(shù)字?jǐn)z像機(jī)

獲取噴霧圖像時(shí),打開(kāi)氙燈向系統(tǒng)提供充分穩(wěn)定的光照。將攝像機(jī)的分辨率調(diào)整為1 600像素×1 200像素,拍攝頻率為1 000 幀/s。燃料為70%柴油與30%乙醇混合燃料,高壓油泵調(diào)節(jié)噴射壓力。試驗(yàn)所用噴嘴的孔徑為0.366 mm,噴射壓力為25 MPa。

在噴霧圖像中,噴霧錐角、油束貫穿距離等信息以像素為單位,筆者設(shè)計(jì)的測(cè)量算法同樣基于像素進(jìn)行編寫,所以需要對(duì)尺寸進(jìn)行標(biāo)定,進(jìn)而確定每個(gè)像素所代表的實(shí)際尺寸。在獲取噴霧圖像的同時(shí),將刻度尺固定在噴嘴中軸線正上方,確?？潭瘸吆陀褪休S線相對(duì)于攝像機(jī)鏡頭在同一水平面上。通過(guò)單位刻度統(tǒng)計(jì)像素?cái)?shù)量,進(jìn)而確定單位像素代表的實(shí)際距離。噴霧圖像與刻度尺如圖4所示,刻度尺局部放大如圖5所示。

圖4 噴霧圖像與刻度尺

圖5 刻度尺局部放大

通過(guò)Photoshop軟件將圖片放大,1 cm所占像素為31個(gè),即單位像素代表的實(shí)際距離為1/31 cm。

4 噴霧圖像處理

4.1 圖像預(yù)處理

由于高速攝像機(jī)拍攝到的噴霧圖像為.rgb格式,因此需要將其轉(zhuǎn)換為灰度圖像以減小儲(chǔ)存量,然后提取圖像中的感興趣區(qū)域。用MATLAB軟件中的Rgb2gray函數(shù)將拍攝到的原始.rgb格式圖像轉(zhuǎn)換為灰度圖像,并用Imcrop函數(shù)分割出感興趣區(qū)域,同時(shí)顯示該區(qū)域的灰度直方圖。原始.rgb格式圖像如圖6所示,感興趣區(qū)域圖像如圖7所示,灰度直方圖如圖8所示。

4.2 圖像二值化

觀察灰度直方圖可以發(fā)現(xiàn),像素值基本落在[0,50]之間,此時(shí)目標(biāo)與背景對(duì)比度不大,如果用Otsu法進(jìn)行全閾值分割,將會(huì)把部分油束當(dāng)作背景,造成油束邊緣損失。Otsu法閾值分割圖像如圖9所示。用筆者的測(cè)量算法,首先提取灰度圖像中灰度值不為0的像素點(diǎn),將所有不為0的像素點(diǎn)全部置為1,得到如圖10所示的二值化圖像。然后用連通域去除小于一定面積的區(qū)域,如圖11所示。經(jīng)過(guò)多次試驗(yàn),筆者選擇去除面積小于100像素×100像素的區(qū)域。最后進(jìn)行孔填充,得到如圖12所示孔填充后的圖像。

圖6 原始.rgb格式圖像

圖7 感興趣區(qū)域圖像

圖8 灰度直方圖

圖9 Otsu法閾值分割圖像

圖10 二值化圖像

圖11 去除小于一定面積區(qū)域

圖12 孔填充后圖像

4.3 獲取油束輪廓邊界

利用MATLAB軟件中的Bwboundaries函數(shù)獲取油束的邊界輪廓,如圖13所示。再根據(jù)油束貫穿距離的1/2上下對(duì)稱分割油束,得到油束分割圖,如圖14所示。

圖13 油束邊界輪廓

4.4 卡爾曼濾波處理

卡爾曼濾波是一種基于線性無(wú)偏最小方差估計(jì)的遞推濾波方法,將信號(hào)看作一個(gè)線性系統(tǒng)在白噪聲作用下的輸出,利用狀態(tài)方程來(lái)描述這一系統(tǒng)的輸入輸出關(guān)系[4]?？柭鼮V波原理是利用之前一個(gè)狀態(tài)的估計(jì)值和當(dāng)前狀態(tài)的測(cè)量值,來(lái)估算出當(dāng)前狀態(tài)的估計(jì)值。由于卡爾曼濾波主要針對(duì)離散變量,因此需要一個(gè)離散型控制過(guò)程[5]。這一系統(tǒng)可以用式(1)線性隨機(jī)方程表示,并用式(2)計(jì)算當(dāng)前時(shí)刻系統(tǒng)的測(cè)量值。

圖14 油束分割圖

式中:X(k)為k時(shí)刻的系統(tǒng)狀態(tài);X(k-1)為k-1時(shí)刻的系統(tǒng)狀態(tài);U(k)為k時(shí)刻對(duì)系統(tǒng)的控制量;A、B為狀態(tài)轉(zhuǎn)移矩陣;Z(k)為k時(shí)刻的測(cè)量值;H為觀測(cè)矩陣;W(k)為過(guò)程噪聲;V(k)為測(cè)量噪聲。

筆者假設(shè)W(k)和V(k)為高斯白噪聲,在式中以方差Q、R代替,其中Q為系統(tǒng)過(guò)程噪聲的方差,R為系統(tǒng)測(cè)量噪聲的方差。測(cè)量噪聲方差R與儀器測(cè)量精度息息相關(guān),一般可以從統(tǒng)計(jì)學(xué)方面去理解,即對(duì)試驗(yàn)數(shù)據(jù)進(jìn)行大量統(tǒng)計(jì),然后得出測(cè)量方差。過(guò)程噪聲方差Q可以通過(guò)對(duì)比試驗(yàn)獲得。通過(guò)大量試驗(yàn),取R為0.35、Q為0.01時(shí),濾波效果較好。

用k-1時(shí)刻的系統(tǒng)狀態(tài)預(yù)測(cè)k時(shí)刻的系統(tǒng)狀態(tài)X(k|k-1):

X(k|k-1)=AX(k-1|k-1)+BU(k)

(3)

式中:X(k-1|k-1)為k-1時(shí)刻系統(tǒng)狀態(tài)最優(yōu)值。

計(jì)算對(duì)應(yīng)X(k|k-1)的協(xié)方差矩陣P(k|k-1):

P(k|k-1)=AP(k-1|k-1)A′+Q

(4)

式中:P(k-1|k-1)為X(k-1|k-1)狀態(tài)對(duì)應(yīng)的協(xié)方差矩陣;A′為A的轉(zhuǎn)置矩陣;Q為系統(tǒng)過(guò)程的協(xié)方差矩陣。

有了當(dāng)前k時(shí)刻預(yù)測(cè)值X(k|k-1)與k時(shí)刻系統(tǒng)測(cè)量值Z(k),就可以得到k時(shí)刻系統(tǒng)的估計(jì)值X(k|k):

(5)

式中:Kg為卡爾曼增益。

(6)

式中:R為系統(tǒng)測(cè)量過(guò)程的協(xié)方差矩陣。

由于離散控制過(guò)程結(jié)束后卡爾曼濾波器才結(jié)束運(yùn)行,因此要不斷更新k時(shí)刻系統(tǒng)狀態(tài)X(k|k)的協(xié)方差矩陣P(k|k):

P(k|k)=[1-Kg(k)H]P(k|k-1)

(7)

隨著狀態(tài)的不斷更新,卡爾曼濾波器進(jìn)行回歸運(yùn)算,直到離散控制過(guò)程結(jié)束[5]。之所以引入?yún)f(xié)方差,是因?yàn)橄到y(tǒng)更新存在一定的誤差,考慮協(xié)方差可以提高當(dāng)前狀態(tài)估計(jì)值的精度[6]。

由于卡爾曼濾波在處理噴霧油束邊界數(shù)據(jù)時(shí)回歸性和誤差都好于均值濾波與中值濾波,因此筆者選擇卡爾曼濾波進(jìn)行處理。對(duì)上半部分油束邊界數(shù)據(jù)進(jìn)行分析,由于油束的對(duì)稱性,下半部分的分析情況與上半部分相同。采用不同濾波方法對(duì)油束邊界數(shù)據(jù)的處理結(jié)果如圖15所示。

圖15 油束邊界數(shù)據(jù)濾波處理結(jié)果

4.5 最小二乘法擬合

提取出噴霧油束邊界數(shù)據(jù)后,需要得到油束兩側(cè)邊界的直線,然后根據(jù)兩條直線的斜率計(jì)算出夾角,這一夾角即為噴霧錐角。目前求油束邊界擬合直線的方法主要有霍夫變換法及基于霍夫變換的改進(jìn)方法,如當(dāng)一側(cè)的擬合直線不止一條時(shí),選取最長(zhǎng)的一條直線作為最佳直線[7]?；舴蜃儞Q擬合油束邊界如圖16所示。雖然霍夫變換法的抗干擾能力較強(qiáng),但是擬合直線不止一條時(shí),下半部分油束邊界檢測(cè)效果較差,噴霧錐角測(cè)量誤差較大,而且通過(guò)調(diào)整霍夫變換法的參數(shù)來(lái)控制擬合直線的條數(shù)往往效果不佳[8]。因此,考慮到油束經(jīng)過(guò)卡爾曼濾波后噪聲點(diǎn)明顯減少,加上對(duì)擬合直線的精度要求較高,筆者采用誤差更小的最小二乘法來(lái)擬合油束邊界[9]。最小二乘法的定義為MIN∑[yi-Φ(xi)]2,(xi,yi)為已知數(shù)據(jù)點(diǎn),Φ(xi)為擬合函數(shù)。

圖16 霍夫變換擬合油束邊界

擬合直線為一次函數(shù),最小二乘法考慮已知數(shù)據(jù)點(diǎn)到擬合直線距離的二次方和最小[10]。最小二乘法擬合油束邊界如圖17所示。

圖17 最小二乘法擬合油束邊界

通過(guò)計(jì)算最小二乘法得到的兩條直線的斜率,可以求得對(duì)應(yīng)時(shí)刻的噴霧錐角,試驗(yàn)計(jì)算結(jié)果為15.231 7°。

5 測(cè)量結(jié)果及分析

對(duì)所有噴霧圖像進(jìn)行批處理,分別應(yīng)用文獻(xiàn)[2]余弦定理測(cè)錐角方法、文獻(xiàn)[3]霍夫變換測(cè)錐角方法,以及筆者所提出的測(cè)量算法進(jìn)行噴霧錐角測(cè)量,結(jié)果如圖18所示。

圖18 噴霧錐角測(cè)量結(jié)果

測(cè)量結(jié)果表明,三種測(cè)量方法都能在總體上反映噴霧錐角的變化趨勢(shì),即先迅速增大后逐漸趨于平穩(wěn)。在啟噴階段,由于噴射壓力較大,噴霧錐角快速增大。在增大到一定階段后,由于噴射壓力減小及環(huán)境背壓的存在,導(dǎo)致噴霧錐角逐漸穩(wěn)定。應(yīng)用筆者測(cè)量算法得到的噴霧錐角曲線穩(wěn)定性要優(yōu)于文獻(xiàn)[2]與文獻(xiàn)[3]方法。文獻(xiàn)[2]方法在前幾幀測(cè)得的噴霧錐角大于筆者測(cè)量算法,這是因?yàn)閱婋A段,油束速度快,此時(shí)利用油束貫穿距離1/2處的最小行與最大行像素點(diǎn)作為計(jì)算噴霧錐角的參照點(diǎn)是不合理的。文獻(xiàn)[3]在啟噴階段的噴霧錐角與筆者測(cè)量算法所得結(jié)果相差不大,但是之后明顯小于筆者測(cè)量算法,且波動(dòng)較大,次要原因是噴霧圖像不規(guī)則,主要原因是預(yù)處理階段對(duì)噴霧圖像進(jìn)行二值化時(shí)丟失了油束邊界信息,所以后期總體噴霧錐角變小。

6 結(jié)束語(yǔ)

針對(duì)測(cè)量汽車發(fā)動(dòng)機(jī)噴霧錐角過(guò)程中直接使用Otsu法二值化會(huì)損失油束邊界信息,以及霍夫變換檢測(cè)油束邊界直線精度不高這兩個(gè)問(wèn)題,基于MATLAB軟件數(shù)字圖像處理技術(shù)提出了一種新的測(cè)量算法:首先考慮噴霧圖像實(shí)際灰度像素值比較小的特征,將像素值不為0的都設(shè)為1;然后用孔填充、連通域去除小于一定面積的區(qū)域,完整分割出噴霧油束,提取油束邊界數(shù)據(jù),選擇回歸性較好的卡爾曼濾波進(jìn)行濾波去噪;最后基于最小二乘法擬合油束邊界數(shù)據(jù),根據(jù)直線斜率計(jì)算出噴霧錐角。為了充分驗(yàn)證這一測(cè)量算法的可行性及準(zhǔn)確性,在同一工況下拍攝完整噴霧過(guò)程圖像進(jìn)行噴霧錐角測(cè)量,確認(rèn)筆者測(cè)量算法能夠有效保護(hù)油束邊界,測(cè)量精度較高,具有一定的實(shí)用價(jià)值。