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

（上海理工大學(xué) 能源與動(dòng)力工程學(xué)院，上海 200093）

定量而精確的顆粒速度測(cè)量是獲取顆粒流動(dòng)特性、高效組織工質(zhì)流動(dòng)的前提，是研究與優(yōu)化系統(tǒng)結(jié)構(gòu)的設(shè)計(jì)基礎(chǔ)，在能源、動(dòng)力、化工與環(huán)境等領(lǐng)域具有廣泛而重要的工程應(yīng)用價(jià)值［1-2］。

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

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

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

1.1 單幀長曝光成像原理

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

圖1 背光式測(cè)量系統(tǒng)示意圖Fig.1 Schematic diagram of measurement system with backward illumination

圖2 運(yùn)動(dòng)模糊軌跡示意圖Fig.2 Schematic diagram of blurred motion trajectory

本文對(duì)處于運(yùn)動(dòng)狀態(tài)下的顆粒進(jìn)行如下三點(diǎn)假設(shè)：①運(yùn)動(dòng)顆粒為圓形，直徑為D；②由于曝光時(shí)間很短，故可以假設(shè)顆粒在曝光時(shí)間內(nèi)為近似的勻速直線運(yùn)動(dòng)，速度為v；③顆粒在垂直于成像平面方向上速度為0，即為二維運(yùn)動(dòng)。

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

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

1.2 軌跡圖像灰度分布

根據(jù)運(yùn)動(dòng)軌跡成像原理，對(duì)圓點(diǎn)軌跡圖像的灰度分布進(jìn)行分析。假設(shè)有一不透光圓點(diǎn)半徑為R，以速度v由左向右水平運(yùn)動(dòng)，成像系統(tǒng)放大倍率為M，圖像傳感器的像元大小為P，圓點(diǎn)在曝光時(shí)間t內(nèi)形成的理想軌跡如圖3（a）所示。假設(shè)運(yùn)動(dòng)軌跡長度大于直徑（實(shí)驗(yàn)時(shí)可通過調(diào)節(jié)曝光時(shí)間以滿足該條件），可將其分成Ⅰ和Ⅱ兩個(gè)區(qū)域分別進(jìn)行研究，其中區(qū)域Ⅰ為軌跡的兩端[見圖3（b）]，區(qū)域Ⅱ?yàn)檐壽E的中間部分[見圖 3（c）]。

圖3 圓點(diǎn)運(yùn)動(dòng)模糊軌跡區(qū)域劃分及參數(shù)示意Fig.3 Division of blurred motion dot trajectory and the relative parameters

由于運(yùn)動(dòng)軌跡具有對(duì)稱性，以軌跡幾何中心為原點(diǎn)，以顆粒運(yùn)動(dòng)方向?yàn)閤軸建立像素坐標(biāo)系，取軌跡的左上部分區(qū)域進(jìn)行分析。以區(qū)域Ⅰ為例，在背光成像方式下，區(qū)域Ⅰ中坐標(biāo)點(diǎn)（x，y）在成像過程中被圓點(diǎn)遮光的時(shí)間tz為

式中，L為該像素被圓點(diǎn)遮光時(shí)間下，圓點(diǎn)運(yùn)動(dòng)的軌跡長度。由于像素點(diǎn)相對(duì)灰度（該像素點(diǎn)灰度與背景灰度的差值） ΔGI與該點(diǎn)被遮光的時(shí)間tz成線性關(guān)系，即

式中，k為系統(tǒng)常數(shù)，與光源強(qiáng)度和系統(tǒng)響應(yīng)有關(guān)。同理推導(dǎo)得到區(qū)域Ⅱ中的圖像灰度值分布為

可見，區(qū)域Ⅰ的灰度在同一y值下與x成線性關(guān)系，且線性系數(shù)為v的函數(shù)；而區(qū)域Ⅱ的灰度與x無關(guān)，只與y有關(guān)。本文基于區(qū)域Ⅰ的灰度分布進(jìn)行速度測(cè)量。當(dāng)y= 0時(shí)，即軌跡區(qū)域Ⅰ長軸線上的像素點(diǎn)灰度分布為

1.3 圖像處理流程

針對(duì)基于灰度分布擬合的圓點(diǎn)速度測(cè)量過程，對(duì)圓點(diǎn)運(yùn)動(dòng)軌跡圖像的信息提取算法進(jìn)行研究，形成的圖像處理流程如圖4所示。

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

具體步驟為：

（1）圖像預(yù)處理：讀取圓點(diǎn)軌跡圖像IP[如圖5（a）所示]和背景圖像Ib，對(duì)Ib-IP的圖像經(jīng)過維納去噪后得到圖像I，如圖5（b）所示；

圖5 圓點(diǎn)運(yùn)動(dòng)模糊軌跡Fig.5 Blurred motion dot trajectory

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

（3）針對(duì)單個(gè)軌跡圖像，利用Radon算法確定軌跡中心點(diǎn)并建立坐標(biāo)系；

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

2 仿真圖像分析

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

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

圖6 半徑81個(gè)像素的圓點(diǎn)水平運(yùn)動(dòng)301個(gè)像素的運(yùn)動(dòng)軌跡仿真圖片F(xiàn)ig.6 Simulated motion trajectory of a dot with radius 81 pixels and moving length 301 pixels

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

圖7 仿真圖片的 k 值標(biāo)定及測(cè)速結(jié)果Fig.7 Calibration of k values and measured velocities for synthetic images

3 實(shí)驗(yàn)裝置及參數(shù)標(biāo)定

3.1 背光式圖像法圓點(diǎn)測(cè)速實(shí)驗(yàn)裝置

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

圖8 背光式圖像法圓點(diǎn)速度測(cè)量系統(tǒng)Fig.8 Measurement system of the dot velocity with backward illumination imaging

3.2 參數(shù) k 的標(biāo)定

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

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

圖9 不同圓點(diǎn)直徑和速度下對(duì)應(yīng)的 kFig.9 k values for different dot diameters and velocities

3.3 速度測(cè)量結(jié)果與誤差分析

根據(jù)圖像處理算法原理結(jié)合實(shí)驗(yàn)圖像，采用灰度分布擬合法、二值化結(jié)合Regionprops 函數(shù)法以及二值化結(jié)合Radon變換函數(shù)法等三種測(cè)速方法進(jìn)行圖像處理并獲得速度測(cè)量誤差，結(jié)果如圖10所示。

由圖10中可知，當(dāng)圓點(diǎn)直徑較小時(shí)，三種圖像處理方法的相對(duì)誤差相差無幾，但隨著圓點(diǎn)直徑增加，二值化結(jié)合Regionprops函數(shù)法和Radon變換函數(shù)法相對(duì)誤差大大增加，而灰度分布擬合法的測(cè)速相對(duì)誤差則無明顯增加，甚至有所降低。究其原因，如圖11所示，黑色區(qū)域?yàn)槎祷蟮能壽E圖像（Radaon函數(shù)處理基于該區(qū)域），橢圓為與二值化區(qū)域具有相同標(biāo)準(zhǔn)二階中心矩的橢圓（Regionprops函數(shù)處理基于該區(qū)域）；利用Radon函數(shù)法處理運(yùn)動(dòng)軌跡時(shí)，其徑向線方向的投影無法解決二值化過程中的長軸損失問題，所以，其對(duì)測(cè)速結(jié)果帶來的誤差最大；而Regionprops函數(shù)法測(cè)速結(jié)果雖優(yōu)于Radon變換函數(shù)法，究其原因是由于其等效橢圓長軸拉伸了部分條件下的二值化軌跡圖像，但是不同情況下長軸損失的程度各不相同。相比較而言，本文提出的利用軌跡成像原理速度法由于避免了直接基于二值化圖進(jìn)行測(cè)量，測(cè)速誤差可降低5%～25%。

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

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

4 結(jié) 論

本文從圖像灰度分布識(shí)別角度，對(duì)單幀長曝光軌跡法測(cè)速誤差進(jìn)行理論和實(shí)驗(yàn)研究，提出一種基于軌跡圖像灰度分布的測(cè)速方法。具體結(jié)論為：

（1）單幀長曝光圓點(diǎn)軌跡的長軸線兩端區(qū)域（區(qū)域I），像素點(diǎn)相對(duì)灰度與該點(diǎn)距軌跡質(zhì)心的距離在理論上呈較為嚴(yán)格的線性關(guān)系，且標(biāo)定系數(shù)k與圓點(diǎn)大小R和曝光時(shí)間t均無關(guān)；

（2）在對(duì)仿真圖像的分析中發(fā)現(xiàn)，圓點(diǎn)直徑越大，本文提出的基于軌跡圖像灰度分布的測(cè)速方法計(jì)算結(jié)果越準(zhǔn)確；在測(cè)速范圍內(nèi)軌跡長度即顆粒速度對(duì)其影響不大；圖像噪聲使得誤差整體增大，但不影響上述規(guī)律；

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