基于軌跡圖像灰度分布的顆粒速度測量

周 騖，邵星翔，陳本珽，蔡小舒

（上海理工大學 能源與動力工程學院，上海 200093）

定量而精確的顆粒速度測量是獲取顆粒流動特性、高效組織工質流動的前提，是研究與優化系統結構的設計基礎，在能源、動力、化工與環境等領域具有廣泛而重要的工程應用價值［1-2］。

基于數字圖像處理的顆粒檢測方法是圖像處理技術的重要分支，相比于其他測量方法圖像法具有直觀、可靠、簡便等優點［3-4］，尤其適用于稀疏顆粒相多參數在線測量［5-6］。單幀長曝光軌跡圖像法［7-8］通過適當延長工業相機的曝光時間以獲得離散顆粒的單幀“拖影”圖像，即運動軌跡，通過軌跡寬度和長短同時獲取顆粒粒徑和運動速度信息，從而可以降低對測量系統幀率的要求。已經有許多學者對圖像法運動場測量的方法及其應用做了大量研究，例如：吳學成等［9］對煤粉顆粒進行了基于軌跡圖像的速度和粒徑測量；王若琳等［10］采用單幀運動模糊圖像處理測量了顆粒粒度信息；張超等［11］利用軌跡圖像實現了噴嘴液滴粒徑和速度測量系統構建；馮明亮等［12］從軌跡識別的角度，對單幀長曝光軌跡法測速上限的影響因素進行了理論分析和實驗研究。但是，軌跡法的測量精度還有進一步提升的空間，尤其需要從圖像處理方面做進一步的研究以獲得更準確的結果。

本文基于顆粒運動中的成像過程，即軌跡成像法，首先，通過理論分析，提出基于軌跡圖像灰度分布的顆粒測速方法，推導顆粒速度測量的理論公式；其次，基于仿真圖像的處理分析，證明上述方法的可行性和可靠性；最后，通過實驗研究，對基于該方法的顆粒速度測量誤差進行分析，并與前期工作中采用的直接二值化結合Regionprops函數或Radon變換測速方法的誤差進行對比研究，為利用背光照明式成像系統測量運動顆粒提供參考和依據。

1 基于軌跡圖像灰度分布的顆粒測速原理

1.1 單幀長曝光成像原理

圖1為典型的背光式單幀長曝光圖像法測量系統示意圖。光源直接朝向成像系統照明，物體（如顆粒）對光線形成遮擋從而在相機中投影成像，即形成背景亮物體暗的“剪影”圖像。當在曝光時間t內顆粒平行于成像平面有運動時，像場的光能量（或能量的遮擋效應）于曝光時間t內在圖像傳感器上的累積，產生電荷信號，經過數模轉換，最終由計算機量化為數字圖像。因此，顆粒圖像實際上是顆粒在曝光時間內經過的所有位置的像的疊加，如圖2所示，運動模糊的圖像正蘊含了物體的速度信息。

圖1 背光式測量系統示意圖Fig.1 Schematic diagram of measurement system with backward illumination

圖2 運動模糊軌跡示意圖Fig.2 Schematic diagram of blurred motion trajectory

本文對處于運動狀態下的顆粒進行如下三點假設：①運動顆粒為圓形，直徑為D；②由于曝光時間很短，故可以假設顆粒在曝光時間內為近似的勻速直線運動，速度為v；③顆粒在垂直于成像平面方向上速度為0，即為二維運動。

基于上述假設，顆粒粒徑等于軌跡圖像的寬度；顆粒在曝光時間t內通過的位移L等于S和D的差，則顆粒在曝光時間內的平均速度為

常規速度測量的方法即通過圖像二值化后識別軌跡的長軸長度和短軸長度，已知相機的曝光時間，基于式（1）獲得速度。獲取S和D的圖像處理算法主要有兩種：①Matlab軟件自帶的圖像處理工具箱中的Regionprops函數；②圖像在一個特定角度下的徑向線方向投影得到軌跡圖像的Radon變換函數。Regionprops函數是Matlab軟件中一個重要的圖像分析函數，它可以快速度量圖像區域的像素數、重心、等效圓直徑等參數［13］。學者們利用函數中的等效短軸和等效長軸這兩個參數分別來近似計算顆粒軌跡圖像的寬度和長度，從而進一步提取出顆粒的粒徑和速度信息。Radon變換函數由奧地利數學家約翰·雷登于1917年提出，其數學意義是平面內函數f（x，y）沿特定直線的線積分，在二維圖像中，一幅圖像的Radon變換結果就是這幅圖像在一個特定角度下的徑向線方向的投影［14］。但由成像原理導致軌跡邊緣灰度與背景灰度接近，且存在系統噪聲，難以合理選取閾值大小，二值化操作常常帶來較大誤差。因此，本文擬從圖像灰度分布的角度，探索顆粒速度測量的方法，避免或減小二值化過程帶來的誤差。

1.2 軌跡圖像灰度分布

根據運動軌跡成像原理，對圓點軌跡圖像的灰度分布進行分析。假設有一不透光圓點半徑為R，以速度v由左向右水平運動，成像系統放大倍率為M，圖像傳感器的像元大小為P，圓點在曝光時間t內形成的理想軌跡如圖3（a）所示。假設運動軌跡長度大于直徑（實驗時可通過調節曝光時間以滿足該條件），可將其分成Ⅰ和Ⅱ兩個區域分別進行研究，其中區域Ⅰ為軌跡的兩端[見圖3（b）]，區域Ⅱ為軌跡的中間部分[見圖 3（c）]。

圖3 圓點運動模糊軌跡區域劃分及參數示意Fig.3 Division of blurred motion dot trajectory and the relative parameters

由于運動軌跡具有對稱性，以軌跡幾何中心為原點，以顆粒運動方向為x軸建立像素坐標系，取軌跡的左上部分區域進行分析。以區域Ⅰ為例，在背光成像方式下，區域Ⅰ中坐標點（x，y）在成像過程中被圓點遮光的時間tz為

式中，L為該像素被圓點遮光時間下，圓點運動的軌跡長度。由于像素點相對灰度（該像素點灰度與背景灰度的差值） ΔGI與該點被遮光的時間tz成線性關系，即

式中，k為系統常數，與光源強度和系統響應有關。同理推導得到區域Ⅱ中的圖像灰度值分布為

可見，區域Ⅰ的灰度在同一y值下與x成線性關系，且線性系數為v的函數；而區域Ⅱ的灰度與x無關，只與y有關。本文基于區域Ⅰ的灰度分布進行速度測量。當y= 0時，即軌跡區域Ⅰ長軸線上的像素點灰度分布為

1.3 圖像處理流程

針對基于灰度分布擬合的圓點速度測量過程，對圓點運動軌跡圖像的信息提取算法進行研究，形成的圖像處理流程如圖4所示。

圖4 基于灰度分布擬合的圓點速度測量流程Fig.4 Flow chart of velocity measurement for the dot based on gray distribution fitting

具體步驟為：

（1）圖像預處理：讀取圓點軌跡圖像IP[如圖5（a）所示]和背景圖像Ib，對Ib-IP的圖像經過維納去噪后得到圖像I，如圖5（b）所示；

圖5 圓點運動模糊軌跡Fig.5 Blurred motion dot trajectory

（2）采用Otsu算法計算整幅圖像閾值，并進行圖像二值化，二值化圖像邊界位置如圖5（b）中黑色閉合曲線所示，若一幅圖像中有多條軌跡則進行軌跡圖的分割；

（3）針對單個軌跡圖像，利用Radon算法確定軌跡中心點并建立坐標系；

（4）如圖5（b）所示，提取顆粒軌跡的長軸長度L1和短軸長度L2，以長軸一端點xA作為選取像素點的起始位置（xA=L1/2）；從該點起向坐標中心方向共選取L2/2個像素位置，作為待擬合區域，即 (xB,xA)；并讀取圖像I中該區域的灰度值ΔG（由于第一步預處理，此處已為相對灰度值）。

2 仿真圖像分析

為了對基于軌跡圖像灰度分布的測速方法可行性進行驗證，仿真生成不同圓點半徑和不同速度下的軌跡圖片，按照1.3節的圖像處理流程進行處理獲得速度，并與理論值進行比較，分析噪聲的影響。

基于Matlab軟件平臺，采用Motion算子卷積清晰圓點圖片形成運動模糊軌跡。以圖6（a）為例，它是半徑為81個像素的圓點水平運動301個像素時形成的圖片，其中運動模糊前的圖片背景灰度為0，顆粒靜止成像的前景灰度為0.8；添加高斯噪聲（噪聲方差3.25）后如圖6（b）所示。假設成像系統曝光時間為20 μs，放大倍率為1倍，像元大小為5.3 μm，根據圖6（a）獲得的標定系數k為1.0268×107。

圖6 半徑81個像素的圓點水平運動301個像素的運動軌跡仿真圖片Fig.6 Simulated motion trajectory of a dot with radius 81 pixels and moving length 301 pixels

假設圓點半徑由11變化到101個像素，運動模糊軌跡長度由101變化至1001個像素，獲得不同圓點半徑和不同速度下標定的k值，如圖7（a）所示?？梢姡涸撓禂蹬c所選圓點大小及運動速度基本無關；在圓點較小時，由于擬合數據點數量過少（約為圓點半徑對應的像素個數），標定的k值變化范圍較大；在圓點較大時，若運動長度接近或小于圓點直徑，由于所擬合區域將包含圖3中的區域Ⅱ， ΔGI0～x數據將偏離線性關系，上述方法將不再適用，即該方法的前提是運動軌跡足夠長（這一點在實際測量中可以通過調節曝光時間來保證）。取標定系數的平均k值（1.0306×107），根據1.3節所述的處理流程，得到速度計算值，并與理論值進行比較，結果如圖 7（b）和（c）所示。圖 7（b）為不含噪聲圖像的測速結果誤差，其值等于各k值與平均k值的相對誤差；隨著圓點直徑的增大，測量相對誤差有整體減小的趨勢，這一前提是運動軌跡長度大于圓點直徑；但隨著軌跡長度的進一步增加，由于運動模糊導致前景圖像灰度逐漸接近背景，即所擬合曲線的斜率逐漸降低，誤差又逐漸增大，即測速范圍存在上限，而該上限與粒徑有關。圖7（c）為仿真圖片添加方差為3.25的高斯噪聲后的測速結果。在軌跡長度較長時，擬合曲線的灰度變化可能僅有幾個灰度級，噪聲造成的影響非常大，極個別到100%以上。為更好地顯示可測量區域的誤差，此處僅顯示501個運動長度以下的測速誤差。可見，誤差水平整體提高，但圓點直徑和軌跡長度的影響趨勢與圖7（b）類似。

圖7 仿真圖片的 k 值標定及測速結果Fig.7 Calibration of k values and measured velocities for synthetic images

3 實驗裝置及參數標定

3.1 背光式圖像法圓點測速實驗裝置

本文搭建了一套背光式圖像法圓點速度測量系統，如圖8所示。該系統包括工業相機（FL3-U3-13Y3M）、遠心鏡頭（2 ×）、直流電機、光度計、轉速儀和平行光源等。拍攝對象為固定在鋁板扇葉上標定板中的黑色圓點。實驗時標準圓點板隨著鋁板扇葉的轉動獲得轉速，并被調至鏡頭焦平面，以保證相機鏡頭軸線、光源軸線及標定圓點板所拍攝區域中心位于同一水平線上，且與標定板平面垂直，避免遠離光軸位置的離焦模糊現象。相機采集圖片位深選為10位。值得注意的是，測量的軌跡長度一般為圓點直徑的3～8倍，在此區間內運動軌跡的平均曲率僅為0.0035，因此，圓點可認為是沿著水平方向運動，此時適用1.2節中的運動軌跡成像原理。

圖8 背光式圖像法圓點速度測量系統Fig.8 Measurement system of the dot velocity with backward illumination imaging

3.2 參數 k 的標定

基于背光式圖像法實驗裝置，首先對該實驗系統進行標定，實驗選取直徑分別為80、100、150、200、250和400 μm共6種圓點，每種圓點運動速度分別取 5.0、6.3、7.6、8.9、10.1、11.4和12.6 m·s-1，拍攝并篩選出含有運動軌跡的517張實驗圖片進行處理，獲取相應的x與為橫坐標，以 ΔGx為縱坐標繪制散點圖并擬合出斜率k，結果如圖9（a）所示。該圖中整體趨勢與仿真結果一致。當圓點直徑較小時，用于擬合的數據較少，且在高速條件下相對灰度值變化較小[如圖9（b）所示]，擬合誤差增加；即測速范圍上限隨著粒徑的降低而降低，如本文實驗條件下，直徑100 μm圓點的運動速度上限約為8.9 m·s-1。粒徑較大或速度較低時，擬合曲線如圖9（c）所示，符合基于灰度分布擬合速度測量方法的測速范圍。

上述仿真和實驗結果均表明：本文提出的基于軌跡圖像灰度分布的測速方法更適合使用在圓點直徑較大的速度場檢測中；測速范圍存在上、下限，下限要保證運動長度大于顆粒直徑，上限需保證擬合區域的圖片灰度有較明顯的變化；還需要深入研究確定上、下限的方法；此外，噪聲的存在明顯影響圖像的質量，進而影響測速結果的準確性。圖像去噪方法值得進一步深入研究。

圖9 不同圓點直徑和速度下對應的 kFig.9 k values for different dot diameters and velocities

3.3 速度測量結果與誤差分析

根據圖像處理算法原理結合實驗圖像，采用灰度分布擬合法、二值化結合Regionprops 函數法以及二值化結合Radon變換函數法等三種測速方法進行圖像處理并獲得速度測量誤差，結果如圖10所示。

由圖10中可知，當圓點直徑較小時，三種圖像處理方法的相對誤差相差無幾，但隨著圓點直徑增加，二值化結合Regionprops函數法和Radon變換函數法相對誤差大大增加，而灰度分布擬合法的測速相對誤差則無明顯增加，甚至有所降低。究其原因，如圖11所示，黑色區域為二值化之后的軌跡圖像（Radaon函數處理基于該區域），橢圓為與二值化區域具有相同標準二階中心矩的橢圓（Regionprops函數處理基于該區域）；利用Radon函數法處理運動軌跡時，其徑向線方向的投影無法解決二值化過程中的長軸損失問題，所以，其對測速結果帶來的誤差最大；而Regionprops函數法測速結果雖優于Radon變換函數法，究其原因是由于其等效橢圓長軸拉伸了部分條件下的二值化軌跡圖像，但是不同情況下長軸損失的程度各不相同。相比較而言，本文提出的利用軌跡成像原理速度法由于避免了直接基于二值化圖進行測量，測速誤差可降低5%～25%。

圖10 三種圖像處理方法誤差結果對比Fig.10 Comparison of errors among three image processing methods

圖11 測量軌跡長短軸誤差示意圖Fig.11 Measuring errors of major axis length and minor axis length of the trajectory

4 結 論

本文從圖像灰度分布識別角度，對單幀長曝光軌跡法測速誤差進行理論和實驗研究，提出一種基于軌跡圖像灰度分布的測速方法。具體結論為：

（1）單幀長曝光圓點軌跡的長軸線兩端區域（區域I），像素點相對灰度與該點距軌跡質心的距離在理論上呈較為嚴格的線性關系，且標定系數k與圓點大小R和曝光時間t均無關；

（2）在對仿真圖像的分析中發現，圓點直徑越大，本文提出的基于軌跡圖像灰度分布的測速方法計算結果越準確；在測速范圍內軌跡長度即顆粒速度對其影響不大；圖像噪聲使得誤差整體增大，但不影響上述規律；

（3）搭建了圓點測速實驗裝置并通過實驗比較了三種測速方法；在實驗條件下，對比Otsu二值化結合Regionprops函數法以及二值化結合Radon函數法。本文提出的基于灰度分布擬合速度法測速誤差可降低5%～25%，可靠性更高。