基于灰度圖差分梯度的霧化角測量方法

魏永杰，李明凱，王浩然，張樹日

（河北工業大學機械工程學院，天津300130）

引言

噴霧技術在燃油、消防、冶金、環保、航空航天、化工等領域有了越來越廣泛的應用。由噴嘴噴出的液滴在周圍空氣流的作用下形成復雜的氣液兩相流，霧場特性分析是提高霧化效果的重要途徑，其中霧化角是一項重要指標。劉祺等分析了在不同環境條件下對航空渦輪風扇發動機中所使用的離心式噴嘴噴霧特性，并開展了實驗研究[1]。梁博健等以噴射角和射流流量作為評價指標，考慮多項參考因素的影響，研究高壓水除鱗噴嘴的射流性能[2]。蔣仲安等采用一次霧化和二次多級霧化理論得出影響噴嘴霧化性能的主要因素，并通過實驗分析各因素對噴嘴霧化性能的影響[3]。楊國華等基于高速攝影技術和圖像處理方法，以不同螺旋升角作為變量，對氨法脫硫用螺旋形實心錐噴嘴的霧化特性進行了分析和預測[4]。郭鵬宇等分析了變壓器滅火領域噴嘴的安裝角度、噴嘴直徑和噴霧壓力等對噴嘴霧化特性的影響規律，得到了最佳的噴嘴安裝角度[5]。賈衛東等使用激光相位多普勒粒子測速儀直接測定噴霧霧化角，精度較高，直觀地顯示了霧滴數分布狀況[6]。

根據流體力學原理，液體以一定壓力經噴嘴噴射后，與空氣流相互作用，形成一定范圍內的霧場。由于噴射出的是微米級的小顆粒，而噴霧流的速度較高，會擾動周圍空氣形成差壓氣流，尤其是多相流形成的細小霧滴隨空氣湍流飄灑，造成噴霧場的邊緣不清晰，測量比較困難。

噴嘴霧化角即噴嘴形成的霧場邊界所形成的夾角。除目視測量外，測量噴嘴霧化角主要是采用圖像法[7]。除普通可見光圖像法外，Wei Yijie 等采用高速相機對燃油噴霧使用粒子追蹤測速法進行圖像分析的同時，獲得了噴霧的瞬態輪廓和霧化角[8]；Chen Run等根據激光散射吸收法(laser absorption-scattering,LAS)，采用紫外、可見光雙光束成像技術分析了燃油噴霧的特性和形貌特征[9]；Gr?ger Karsten 等采用X 射線圖像方法對噴霧的近場形貌進行了分析[10]。與上述方法相比，采用普通可見光圖像法分析霧場邊界和霧化角是最經濟快速的方法，但由于噴霧場是由液滴顆粒組成的不連續流體場，采用通常的圖像邊緣算子，如Sobel算子、Canny 算子、Laplace 算子等直接在灰度圖上提取到的霧場邊界并不是很理想。

基于上述原因，一般采用的處理過程是通過對噴霧圖像降噪濾波后進行二值化處理，然后根據圖像形態學針對二值化的圖像進行邊緣提取和計算[11]。由于噴霧圖像對比度較差且霧滴粒度、濃度分布不同，造成灰度圖像不均勻，圖像二值化過程比較困難。目前有多種算法進行復雜灰度圖的二值化方法，如全局閾值二值化、自適應二值化、OTSU 算法等。但這些方法主要是針對圖像進行后期處理，計算噴霧場邊界缺少客觀的評價依據。

根據圖像形態學梯度原理，將灰度圖像分別膨脹、腐蝕后與原圖像相減得到梯度圖像，理論上也可以得到圖像邊界，但由于霧場圖像灰度值分布復雜，難以得到有效邊界。

圖像的差分梯度算法即針對有某種關聯特性的兩幅圖像或經過處理的兩幅圖像，將對應像素或相鄰某個區域的圖像做差分，從而得到圖像的特征[12]。

本文提出根據噴霧灰度圖得到霧場范圍擴大的二值化掩模板，并將此掩模板作用于灰度圖像，針對原始灰度圖進行迭代，并對噴霧范圍邊緣的平均灰度值進行梯度計算，認為梯度最大時的圖像為實際噴霧圖像，從而得到霧場邊界，并可以進行噴霧角度計算。

1 算法原理

1.1 霧場邊界特征分析

圖1 是拍攝的噴霧圖像和采用邊緣算子直接處理灰度圖的結果。圖1（a）是拍攝的噴霧原圖。該噴嘴是一個小型的細水霧噴嘴，噴霧由左到右水平噴出。由于噴霧濃度較小、氣流影響較大，因此在噴射一定距離后霧場和背景的灰度值相差不大。為了更清晰地看清楚噴霧圖像，采用直方圖均衡化對圖1（a）進行了處理，如圖1（b）所示，圖中可以看到上、下2 個灰度值較大的邊緣。由于射流卷吸問題[13]，受氣液兩相流和空氣剪應力的共同作用[14]，造成空氣流在2 個邊緣附近形成負壓且細水霧中的較小顆粒受阻力大，飄散慢，停留時間長，故在2 個灰度值較大的邊緣外各有一條灰度值較小且分布不連續的細水霧飄散帶[15]。因此，認為灰度值較大的邊緣為實際噴霧的有效邊界范圍。

圖1（c）、圖1（d）、圖1（e）和圖1（f）分別是采用4 種邊緣算子對圖1（b）灰度圖進行處理的結果。但這4 種算子只得到霧場邊緣的1 個范圍或1 個局部，無法得到清晰完整的邊緣，更無法判斷實際霧場的有效邊界。

圖1 噴霧灰度圖和邊緣檢測結果Fig.1 Grayscale image of spray and boundary detection results

圖2 是將圖1（b）膨脹后與原圖差分得到的梯度圖像。圖2 中可以看到模糊的邊界，與圖1 中Canny 算子、Robert 算子得到的邊緣相似，但仍不能判斷實際邊界。

圖2 形態學梯度處理結果Fig.2 Gradient processing results by morphology

通過上述分析可知，這些方法是在一定灰度值變化范圍內尋找圖像灰度差異點，并認為這些點是圖像的邊界。因此，會找到多個邊界點，從而得到邊界是1 個范圍，而不是某條確定的邊界線。

1.2 灰度圖差分梯度算法

為了得到有效的邊界線，采用在灰度圖上直接統計一定噴霧邊緣范圍內的灰度值，取灰度差的最大值作為判據提取邊界的方法。

霧場和背景相比，霧場是高灰度值，而背景是低灰度值。為此，首先在原始灰度圖像上，通過圖像形態學的膨脹操作獲得比實際霧場區域大的假設霧場區，然后采用全閾值二值化得到掩模板圖像，其中假設的霧場區域為“白”，背景區域為“黑”，并將此二值化圖作為迭代算法的初始掩模板。

直方圖均衡化并不改變實際霧場的分布，為觀察方便，從獲得初始掩模板的步驟開始，均采用圖1（b）直方圖均衡化的圖像作為原始圖像。將該掩模板迭代進行腐蝕操作，依次得到邊界范圍逐漸減小的掩模板；采用相鄰的2 個掩模板分別作用于圖1（c）的圖像，將得到的2 個圖像做差分后，計算差分灰度圖像的累計灰度值；將2 個相鄰差分圖像的平均灰度之差作為灰度梯度。

由于飄散的細水霧從有效霧場邊界向外濃度逐漸降低，因此該灰度梯度從迭代開始逐漸增大；而在有效霧場范圍內，灰度梯度依次減小，在飄散的細水霧顆粒和有效霧場邊界處有最大值。因此，通過迭代可以得到該有效邊界。

圖3所示為迭代原理。由圖3 可知，在圖像處理整個過程中，只有模板是二值化的圖像。將原始圖像多次膨脹后取全閾值二值化，得到模板1，將此模板“罩”在原始灰度圖上，得到噴霧區的灰度圖1，同時將噴霧區外的灰度值置為0，并計算噴霧區的平均灰度值1；將二值化模板1 再腐蝕一次，作為二值化模板2，重復上述過程，得到噴霧區的灰度圖2 和平均灰度值2。與灰度值1 相減，得到灰度圖的梯度1；同理，得到第n次迭代的梯度n。由上面的分析可知，當梯度n取得最大值時停止迭代，得到噴霧有效邊界。

圖3 灰度圖梯度迭代原理框圖Fig.3 Block diagram of gradient iteration on grayscale images

2 實驗及分析

采用上述提出的方法，針對圖1（a）實際采集的霧場圖像，經直方圖均衡化后進行迭代處理，得到霧場有效邊界，如圖4所示。

圖4（a）是迭代算法采用的掩膜板，是將圖1（b）的灰度圖通過膨脹操作后采用全局二值化獲得。比較圖4（a）和圖1（b），噴霧范圍擴大，尤其從左側噴嘴位置觀察更加明顯。噴霧是由滿足一定分布的液滴粒子構成，在拍攝的原始圖像和直方圖均衡化圖像中，在噴霧區會有一定的孔洞，這些孔洞是由于該處不存在液滴粒子或液滴粒子數量較少形成的。圖像形態學的膨脹操作恰好具有填補孔洞的作用，因此，圖4（a）是連續性較好的圖像。一般來說，對于霧場范圍的圖像測量方法需要采用平滑濾波，本文采用膨脹法提高了連通性，且后續步驟針對灰度圖進行迭代，所以未采用平滑濾波方法。

圖4 噴霧邊界迭代結果Fig.4 Iterative results of spray boundary

圖4（b）是針對灰度圖進行迭代的結果。比較圖4（b）和圖1（b），處理后的結果與圖1（b）中較亮的邊界一致，因此認為此范圍是霧場噴射的邊界，而不是細水霧飄散形成的邊界。比較圖4（a）和圖4（b），由于在獲得初始掩模板時，在圖4（a）右側中部有1 小塊未經膨脹去除的黑斑。因此，在迭代過程中，將掩膜板逐次腐蝕后在圖4（b）右側對應實際噴霧位置生成了1 塊錯誤的小黑塊。同時，在圖4（a）掩膜板下端中間位置，代表擴大的假設噴霧區出現了較大的直線條，經多次腐蝕操作后，對應的圖4（b）代表實際噴霧邊界的下端中間位置出現了較大誤差。迭代過程中的梯度值如表1所示。

表1 灰度梯度迭代數據Table 1 Results of gradient iteration on grayscale images

由表1 可知，在迭代第8 次時，灰度梯度有最大值，此時即霧場有效邊界。通過比較第8 次迭代以后的數據，尤其是第8 次和第9 次迭代的數據可知，由于霧場的連續性特性，在邊界附近的值較大，而霧場內部液滴粒子分布較為一致，因此灰度梯度逐漸減??；在霧場外部，由于細水霧的飄灑，使霧場有效邊界外的附近區域中，細水霧液滴逐漸減小，且細水霧液滴少、與噴射的水霧液滴不存在連續性，因此第7 次迭代數值為2.503，遠小于第8 次迭代數值100.094，且由第7 次到第1 次依次減小。因此，將第8 次迭代的灰度圖作為識別得到的霧場區域圖，即圖4（b）對應的灰度圖。

表1 數據中，第2、4 次迭代數據略小于前一次迭代數據，第14 次迭代數據略小于后一次迭代數據，與霧場實際分布有波動有關。

為了與傳統的針對二值化圖像處理方法進行比較，采用二值化方法對圖1（a）的灰度原始圖先進行二值化邊界識別，如圖5所示。

圖5 傳統二值化霧場邊界處理結果Fig.5 Results of traditional binary spray field boundary

圖5（a）是直接觀察和對比原始圖設定閾值的結果；圖5（b）是選取比圖5（a）更大閾值的結果；圖5（c）是采用OTSU 算法自動選取閾值的結果。從圖5 的3 個結果可知，難以采用數值判據來分析霧場邊界。

為驗證上述方法的有效性，對另外3 種不同濃度的噴嘴霧場進行迭代運算。圖6 是噴霧邊界處理結果。其中圖6（a）、圖6（b）、圖6（c）是拍攝的原始圖，對應的霧場濃度依次減小。如圖6（a）所示，在產生水噴霧場時，施加的壓力越大，霧滴能完全破碎，霧場越濃；如圖6（b）和圖6（c）所示，施加的壓力越小，霧滴不能破碎，且噴出的水量少、霧場稀。這兩種霧場濃度依次減小，拍攝圖片的灰度值也減小，同時霧滴因未完全破碎呈近細流狀，霧場邊界更模糊，直接用邊界提取也更加困難。

圖6（d）、圖6（e）、圖6（f）分別是采用本文方法的處理結果，可見隨壓力減小，霧場區域逐漸減??；由于邊界附近霧滴少，因此圖6（e）、圖6（f）中得到的霧場邊界更偏離直線。

圖6 不同濃度霧場對比結果Fig.6 Comparison results of spray field with different density

為了擬合噴霧角，將圖4（b）進行二值化處理，如圖4（c）所示。由于處理過程中，已將圖4（b）非霧場區域置0。因此，二值化后的圖4（c）與圖4（b）霧場區域邊界是相同的。

提取圖4（c）的邊界點S，進行分段曲線擬合，再分析曲線斜率，分別找到上、下兩條曲線斜率變化的最大值?？紤]到噴霧角一般為錐形，因此取斜率變化最大值左邊的上、下兩段曲線，分別進行直線擬合，取兩條擬合直線的夾角作為噴霧角。經計算，得到該噴霧角為82.65°，如圖7所示。

圖7 噴霧角計算結果Fig.7 Calculation results of spray angle

圖7 中，取斜率變化最大值并截取由噴霧開始位置之后的平直段，得到兩條由邊界點構成的曲線。圖中的兩條直線分別是對上、下邊界曲線擬合得到的直線，圖中坐標代表的均是像素數。

由于霧場本身的動態特性且受周圍氣流擾動影響較大，因此對霧場邊界和角度進行準確測量和標定還比較困難。本文提出用灰度梯度的最大值作為判據和圖像處理后進行分界角度計算，提供了一種數值判斷依據。

上述方法是針對深色背景下的水霧場圖像進行處理，當背景較復雜時，用本方法與光譜圖像的主成分分析方法等結合使用，也可以實現霧場邊界定位。

3 結論

由于噴霧場液滴的非連續性以及兩相流的作用，造成霧場空間分布動態波動特性明顯，且霧場邊界不清晰，灰度值差別較大。將灰度圖像進行二值化再提取邊界或計算噴霧角，缺少定量評判依據，得到的結果與實際霧場有偏離。

根據灰度圖梯度算法，除掩模板處理的二值化圖像外，其他處理步驟均是針對灰度圖進行處理，并通過構造灰度圖梯度，當灰度圖梯度為最大值時，認為噴霧圖像邊界是實際噴霧場的有效邊界，在此基礎上進行噴霧角擬合計算。該方法提出了一種采用數值判斷噴霧場有效邊界的計算方法，通過灰度梯度最大值確定霧場邊界，從而通過程序自動實現霧場邊界提取和噴霧角計算，避免傳統方法中需要通過直接觀察缺乏評判依據的問題。