基于OpenCV-Python的土壤顆粒動態休止角試驗*

周鵬飛，蒙賀偉，2，3，梁榮慶，張炳成，坎雜，2，3

(1. 石河子大學機械電氣工程學院,新疆石河子,832000; 2. 農業農村部西北農業裝備重點實驗室,新疆石河子,832000; 3. 綠洲特色經濟作物生產機械化教育部工程研究中心,新疆石河子,832000)

0 引言

新疆是中國的棉花主產區,受地理環境條件影響,棉花普遍采用覆膜種植[1-2]。棉花收獲期后田間鋪設的殘膜經機械化回收,其中夾雜裹挾大量土壤,致使殘膜初清理及資源化利用難度大[3-5]。針對該類土壤的物料特性進行研究,是設計開發相關篩分、輸送機械的基礎,對解決殘膜雜質分離問題具有重要意義。

在農業機械實際作業中,土壤顆粒大多處于快速運動的動態變化中,如篩分作業等,靜態休止角不足以充分表征土壤顆粒的流變特性,因此需針對動態休止角進行研究。土壤動態休止角的傳統測定方法是使用測量儀器對土壤運動時形成的休止角進行測量,這種方法誤差較大,難以滿足實際研究需求。近年來,隨著計算機技術的不斷發展,應用于研究的手段和途徑也在不斷擴充和豐富,相關學者積極探索利用計算機圖像處理技術輔助研究的途徑[6]。向偉等[7]采用Matlab軟件對南方黏壤土靜態休止角單側圖像進行處理,測得了該類土壤的靜態休止角。王涵斌等[8]利用基于Sobel算子的計算機圖像處理技術獲取了沙壤土靜態休止角邊界輪廓,通過線性擬合實現對沙壤土靜態休止角的測量。戴飛等[9]通過將采集的土壤靜態休止角圖像導入CAD軟件,利用內置標注工具實現對土壤靜態休止角的測量。田辛亮等[10]利用Matlab軟件對玉米秸稈—土壤混料堆積圖像進行處理,提取邊界輪廓并進行線性擬合,獲得了該混料的靜態休止角。鹿芳媛等[11]將采集圖像導入In-Sight Explorer軟件,利用該軟件的幾何測量模塊實現對水稻芽種靜態休止角的測量。劉凡一等[12]依據k均值聚類算法,獲取了麥粒靜態休止角和動態休止角的邊界,經線性擬合完成麥粒靜、動態休止角的測量。由上可知,將計算機圖像處理技術與基礎試驗研究相結合從而提升研究的準確性和便捷性,已然成為一種發展趨勢。

本文以新疆棉田機收膜雜中的土壤為主要試驗對象,搭建了一套土壤動態休止角測量裝置,采用基于OpenCV-Python的計算機圖像處理技術對圖像進行有效處理,獲取不同轉速下的土壤動態休止角,對結果進行對比分析,匹配適宜土壤動態休止角測量的最佳轉速,為同類散粒物料動態休止角測定以及流動特性研究提供參考。

1 動態休止角測量裝置搭建與圖像采集

1.1 動態休止角測量裝置搭建

本文采用轉鼓法[13]對土壤動態休止角進行測量,搭建的土壤動態休止角測量裝置如圖1所示。土壤物料裝填于空心圓筒內部,空心圓筒前端設有透明蓋板,且通過螺栓進行緊固連接,以便于土壤裝填;空心圓筒后端與電機輸出軸相連,并保持空心圓筒與電機輸出軸在同一軸線。工作時,空心圓筒水平放置于4個對稱布置的托輪上,在電機的驅動下繞自身軸線旋轉。土壤在自身重力、顆粒間摩擦力和離心力的共同作用下,沿空心圓筒內壁進行傾落運動。

1.2 圖像采集與攝像頭標定

數字圖像采集環節采用CMOS攝像頭(分辨率:1 280 像素×720像素;幀率:60 fps),攝像頭水平固定于三腳架,并與空心圓筒保持在同一軸線上,通過USB接口實現與計算機間的通信,并由OpenCV攝像頭驅動程序啟動。所采集的數字圖像在計算機端進行后續處理和分析,處理流程如圖2所示。

圖2 土壤動態休止角測量流程圖

由于光線進入攝像鏡頭時會產生畸變,從而導致進行圖像處理獲取特征信息時出現偏差,對成像造成直接影響。因此,在對采集的圖像進行處理前,需對攝像頭進行標定。本文采用Zhang[14]提出的基于黑白棋盤格的攝像頭標定方法(簡稱張氏標定法),該方法具有精度較高、便于操作等優點。根據張氏標定法,選取平面二維12×10黑白相間的棋盤格為標定板,單個方格尺寸為20 mm×20 mm,其內角點數為11×9,通過OpenCV-Python庫函數驅動攝像頭獲取不同角度的圖像共計9張,并對圖片進行角點檢測并進行標識,從而確定角點坐標,其檢測結果如圖3所示。利用OpenCV-Python庫函數對攝像頭的內參數矩陣進行求解,利用非線性最小二乘法估計畸變系數,使用極大似然估計法優化參數,解得攝像頭的內參數矩陣Q如式(1)所示,畸變系數K如式(2)所示。

圖3 角點提取圖像

(1)

K=[-0.0431.032-0.009-0.011-6.541]

(2)

2 基于OpenCV-Python的圖像處理

2.1 圖像預處理

本文主要針對土壤在空心圓筒內運動時的動態休止角,而在圖像采集過程中難免會產生冗余信息,對所采集圖像進行預處理可以有效剔除不需要的圖像信息、保留某些對于后續處理有重要作用的圖像特征。

1) RGB圖像是3通道圖像,灰度圖是單通道圖像,通過將RGB圖像向灰度圖轉換,可在保留圖像梯度信息的前提下,提高運算速度,縮短處理時間?；叶葓D像與RGB圖像之間的映射關系如式(3)所示。

Gray(x,y)=T(R(x,y),G(x,y),B(x,y))

(3)

式中:Gray(x,y)——圖像灰度值;

R(x,y)、G(x,y)、B(x,y)——各通道分量;

T——圖像灰度值與各通道分量間的函數關系。

本文采用OpenCV庫函數對圖像進行灰度化處理,處理后的圖像如圖4所示。

圖4 灰度圖

2) 土壤動態休止角形成于空心圓筒內部,空心圓筒外部的圖像信息對本文無實際意義。將無用信息剔除,保留所需圖像特征,可提高圖像處理效率。以空心圓筒內壁為邊界,將圖像進行圓形裁剪,保留空心圓筒內圖像特征信息。具體實施方式是:(1)利用OpenCV內置cv2.HoughCircles()函數通過霍夫變換確定空心圓筒的半徑與圓心所在位置;(2)根據空心圓筒半徑與圓心位置剪切去除其余圖像,保留圓內圖像信息,處理后的圖像如圖5所示。

圖5 圓形裁剪后的圖像

2.2 圖像降噪

圖像噪聲是指存在于圖像數據中的不必要的或多余的干擾信息。由于在圖像采集或者傳輸過程中受到電磁環境、光照條件及各種隨機信號等的影響,采集的圖像將存在一定噪聲,直接對后續處理結果的準確性產生影響。對圖像進行降噪處理,可在一定程度上還原圖像信息,降低噪聲帶來的不良影響,滿足圖像后續處理需要。根據統計分布特征,圖像噪聲主要可以分為高斯噪聲、泊松噪聲和脈沖噪聲等。針對不同噪聲,常用高斯濾波和中值濾波等方法對圖像進行降噪處理。利用該類方法進行圖像降噪時,會以損失部分圖像特征信息為代價。為避免降噪處理對圖像特征信息的影響,本文采用三維塊匹配濾波(Block-matching and 3D filtering,BM3D)算法對圖像進行降噪處理。BM3D算法是一種基于圖像塊間相似性進行降噪的三維濾波算法,實現過程可分為基礎估計階段和最終估計階段兩個階段,而每個階段可以分為塊匹配、協同濾波和聚集三個部分,其具體實現過程如下[15-16]。

1) 在整幅圖像中搜尋符合參考塊的相似度要求的匹配塊,將所有匹配塊組合形成一個三維矩陣,其結果可由式(4)求解得出。

SxR={x∈X:d(FxR,Fx)≤τmatch}

(4)

式中:SxR——與參考塊類似匹配塊的集合;

d(FxR,Fx)——參考塊和匹配塊間的距離;

τmatch——距離閾值;

FxR、Fx——三維濾波算法的參考塊、匹配塊。

2) 對三維矩陣進行線性變換,利用閾值收縮變換域的系數,做逆線性變換,從而達到降低圖像噪聲的目的,其過程如式(5)所示。

(5)

T3D——三維域變換;

T3D-1——T3D逆變換;

γ3D——硬閾值;

λth3D——三維閾值參數;

σ——圖像噪聲的標準差;

N1——塊邊長。

3) 由于同一參考塊可能會存在若干估計值,導致各像素點同樣會產生若干估計值。因此,需對估計值進行加權平均計算,如式(6)所示,通過該方式可得到降噪圖像。

(6)

(7)

式中:ωxR——局部估計的權重;

Nhar——降噪后所有非零元素的個數。

由于無法通過主觀對圖像降噪效果進行定量評價,因此需借助客觀手段對原始圖像和降噪處理后圖像進行對比分析,根據結果確定降噪效果[17]。峰值信噪比(Peak signal-to-noise ratio,PSNR)和結構相似指數(Structural Similarity index,SSIM)是目前常見的評價圖像質量的參數。峰值信噪比是評價圖像噪聲水平的參數,其計算公式如式(8)所示,結果取值范圍為[0,100]。當使用峰值信噪比衡量圖像的質量時,其值越大,說明所使用的降噪算法對噪聲的抑制效果越好。結構相似指數是評價降噪后圖像與原圖像相似度的參數,其計算公式如式(9)所示,結果取值范圍為[0,1]。結構相似指數越大,說明降噪算法對圖像特征的保持情況越好,越有利于后續圖像處理對圖像特征的提取。

(8)

式中:u——N×N的數字圖像;

v——降噪后的圖像。

(9)

式中:μx、μy——圖像u和v的均值;

σx、σy——圖像u和v的方差;

σxy——圖像u和v的協方差;

c1、c2——為避免分母為0而引入的常數。

利用動態休止角測量裝置對土壤動態休止角進行圖像采集,隨機選取圖像10幅,經圖像預處理后,采用BM3D算法對其進行降噪處理,通過對比降噪處理前后的圖像,峰值信噪比和結構相似指數結果如表1所示。

表1 BM3D算法對圖像降噪結果的客觀評價Tab. 1 Objective evaluation of BM3D algorithm on the denoising results of experimental images

由表1可知,峰值信噪比和結構相似指數的均值分別為36.48 dB、0.88,可滿足后續圖像處理需求。

2.3 土壤邊界特征提取

圖像邊界指圖像中灰度或結構發生突變的位置,是圖像最基本的特征,也是紋理特征和形狀特征提取等圖像分析的基礎。本文所需提取的邊界為土壤顆粒外緣輪廓,邊界特征提取的準確性將會直接影響動態休止角的測量結果。目前,常見的邊緣檢測方法主要有Roberts、Sobel、Prewitt和Canny邊緣檢測[18-20]。由于Roberts算法易受噪聲影響;Sobel算法獲取的邊緣位置偏差較大;Prewitt算法在檢測過程中容易出現偽邊緣,造成過分割現象。因此,本文采用OpenCV中cv2.Canny()函數對土壤動態休止角圖像進行邊界檢測并獲取土壤邊界,其結果如圖6所示。

圖6 土壤邊界提取圖

3 基于改進的最小二乘法的線性擬合

針對OpenCV圖像處理后所得的土壤顆粒邊界,以空心圓筒圓心處為原點,水平方向為x軸,豎直方向為y軸建立坐標系,確定土壤邊界各像素點坐標數據,分別利用一般最小二乘法和改進的最小二乘法對坐標數據進行線性擬合,通過對求解所得線性擬合模型進行對比分析,獲取擬合度較優的土壤顆粒邊界方程。

3.1 改進的最小二乘法線性擬合

最小二乘法是一種通過對誤差平方進行最小化處理,根據實測數據估計數學模型中未知參數的數學統計方法[21-22]。該方法便于分析與計算,在數據處理方面應用廣泛,其具體實現過程是:(1)假設存在n個樣本,樣本坐標分別為(pi,qi)(i=1,2,3,…,n),式(10)為該樣本的擬合函數,則樣本與目標函數的誤差平方和為式(11);(2)對誤差平方和函數中每個未知系數求偏導,令偏導數為0,求解方程;(3)將方程的解代入擬合函數,該方程即為此樣本通過最小二乘法所得擬合方程。

(10)

式中:φ(q)——擬合函數;

m——擬合函數的擬合次數;

ak——擬合系數。

(11)

式中:δ(p)——樣本與擬合函數的誤差平方和函數。

在實際應用中,由于隨機噪點和干擾點等因素的影響,使部分實測數據存在較大誤差,導致求解獲取的最小二乘法擬合函數偏離原始數據。為提高擬合效果,本文對最小二乘法算法進行優化改進,具體思路是通過最小二乘法對土壤邊界像素點進行線性擬合,獲取擬合函數后分別計算各像素點與擬合點間距離,根據拉依達準則剔除異常點,對剩余像素點再次進行最小二乘法線性擬合,獲取剔除異常點后的擬合函數,盡可能提高函數擬合度,確保結果能充分反映原始數據。

改進的最小二乘法擬合步驟如下。

1) 以空心圓筒圓心處為原點,以水平方向為x軸、豎直方向為y軸,以相鄰像素點間距為單位距離建立坐標系,其中,土壤邊界各像素點坐標為(xi,yi)(i=1,2,3,…,n)。由于需對土壤邊界進行一階線性擬合,故在式(10)中令m=1,得擬合方程如式(12)所示,進而確定誤差平方和函數如式(13)所示。

φ1(x)=a0+a1x

(12)

(13)

對誤差平方和函數中每個未知系數分別求偏導,令偏導數為0,即

(14)

求解方程得

(15)

此時,誤差平方和函數值最小。

2) 分別計算各像素點與擬合點間距離li(i=1,2,3,…,n),如式(16)所示。根據拉依達準則對各像素點進行判定,將各異常點剔除,獲得新的土壤邊界像素坐標點(xj,yj)(j=1,2,3,…,n′,且n′≤n),再次利用最小二乘法進行線性擬合,求解得剔除異常點后的土壤邊界擬合函數為式(17),其中一次項系數為a1′。

li=|yi-(a0+a1xi)|

(16)

φ1′(x)=a0′+a1′x

(17)

3) 根據擬合函數一次項系數,對土壤動態休止角θ求解,得

θ=arctana1′

(18)

3.2 擬合結果分析

通過自制的土壤動態休止角測量裝置采集獲取土壤動態休止角數字圖像,隨機選取10幀OpenCV處理后圖像,分別利用一般最小二乘法和改進的最小二乘法進行線性擬合,根據概率學理論對線性擬合模型分析,結果如表2所示。

表2 擬合結果比較Tab. 2 Comparison of fitting results

通常情況下,利用一般最小二乘法對數據進行擬合時,易受極端異常點的影響,從而降低所構建模型的準確性。由表2可知,相對于一般最小二乘法擬合模型,改進的最小二乘法擬合模型的決定系數R2平均提升了0.29%,均方根誤差平均下降了1.74,擬合效果和精度優于一般最小二乘法擬合模型。該結果表明,改進的最小二乘法可以有效提升線性方程擬合度。

4 試驗結果與分析

4.1 試驗材料及儀器

試驗材料取自新疆瑪納斯縣北五岔鎮棉田機收膜雜中的土壤,該區域地處天山北麓、準噶爾盆地西南緣,土壤以棉田深度為0～50 mm的表層沙壤土為主[23]。在試驗樣本采集過程中按土壤在機收膜雜中所處的不同位置分別進行隨機取樣,通過人工對棉田機收膜雜進行篩揀,剔除機收膜雜中殘膜、棉秸稈等成分。

主要儀器設備:標準土壤篩(孔徑范圍:2～5 mm)、JMB-5003型電子天平(測量范圍:0～500 g,測量精度:0.001 g)、Sartorius MA-45型快速含水率測定儀(質量精度:0.001 g,含水率精度:0.01%)、UT372型激光轉速儀(測量范圍:10～99 999 r/min,測量精度:0.04%)、自封袋。試驗裝置為自制的動態休止角測量裝置。

針對所采集的土壤樣本物料,利用Sartorius MA-45 型快速含水率測定儀測量土壤的含水率,并隨機取樣3次,每次取樣500 g,采用標準土壤篩(孔徑范圍:2～5 mm)對其進行篩分,利用JMB5003型電子天平(測量精度為0.001 g)對不同粒徑下土壤顆粒進行稱重。測定的土壤基本參數如表3所示。

表3 土壤粒徑分布及含水率統計Tab. 3 Soil particle size distribution and moisture content statistics

4.2 試驗方法

試驗采用自制的動態休止角測量裝置進行土壤動態休止角的測定。根據土壤粒徑分布測定結果,選取粒徑在(0,2)區間內的土壤顆粒(圖7)進行動態休止角試驗。填充率過大或過小均不利于對動態休止角的觀測,參照文獻[24]確定本試驗填充率為50%。測量時,通過步進電機調節圓筒旋轉速度,并利用UT372型激光轉速儀進行測量驗證。將攝像頭固定于三腳架上,調至水平,且攝像頭與圓筒處于同一軸線。待裝置運轉平穩后,計算機通過USB接口獲取攝像頭所采集圖像,利用基于OpenCV的計算機圖像處理技術對土壤動態休止角進行測定。

圖7 土壤顆粒

4.3 結果與分析

通過調整裝填土壤后空心圓筒的旋轉速度,使土壤在透明端蓋的內壁上形成平穩流動的斜坡,獲取斜坡和水平方向的夾角,即為對應轉速下的土壤動態休止角,不同轉速下的土壤動態休止角結果如表4所示。

表4 不同轉速下土壤動態休止角測量值Tab. 4 Measured value of soil dynamic angle of repose at different rotational speed

由表4可知,土壤動態休止角隨著轉速的增加呈先減后增趨勢。通過對該現象產生的原因進行分析可知,當轉速較低時,土壤在空心圓筒內呈間歇性滑移運動,土壤顆粒在重力、摩擦力等的共同作用下,與空心圓筒保持相對穩定,并隨空心圓筒共同繞其軸線旋轉,當到達某一臨界高度后,土壤邊界顆粒迅速滑落至顆粒間再次形成穩態,土壤顆粒的間歇運動致使所測得的動態休止角較大;隨著轉速的增大,土壤顆粒由間歇運動逐步轉變為連續運動,故測量所得動態休止角減小;當轉速進一步增大時,部分土壤邊界顆粒的運動狀態由滾落開始向拋落轉變,拋落過程中部分土壤顆粒存在短暫的懸空,該現象致使土壤邊界向外發生偏移,從而導致所測得的動態休止角增大。

目前,相關學者采用線性關系來描述轉速與動態休止角間的關系,但根據轉速與實測所得土壤動態休止角的變化規律發現,單一的線性關系不能很好地進行描述。本文依據不同轉速所對應動態休止角分別進行線性擬合與多項式擬合,擬合所得數學模型分別為式(19)和式(20)。

y1=38.309+0.255x

(19)

y2=39.709-0.548x+0.111x2-

(4.28×10-3)x3

(20)

通過對擬合模型進行方差分析可知,線性擬合模型的R2為89.35%,均方差為0.145;多項式擬合模型的R2為95.72%,均方差為0.061,多項式擬合模型相較于線性擬合模型,擬合度更優。通過所得線性擬合模型和多項式擬合模型繪制轉速與土壤動態休止角關系如圖8所示。

圖8 土壤動態休止角與轉速關系圖

根據文獻[13]可知,不同轉速條件下圓筒內顆粒群具有不同運動特征,而圓筒中顆粒群運動特征可通過式(21)進行預測。

(21)

式中:Fr——Froude數;

ω——轉鼓轉動角速度,rad/s;

R——轉鼓半徑,R=7.5×10-2m;

g——重力加速度,m/s2。

由梅爾曼準則可知,當10-4

(22)

式中:n*——最優試驗轉速;

N(σmin)——標準差最小值時所對應的轉速。

由上可知,本文土壤動態休止角最優試驗轉速為7 r/min。

5 結論

1) 針對棉田機收膜雜中土壤動態休止角研究需求,結合基于OpenCV-Python的計算機圖像處理技術,搭建了一種土壤動態休止角測量裝置,通過該裝置可快速準確獲取土壤動態休止角,滿足試驗測量需求。

2) 利用OpenCV-Python進行圖像處理,通過將采集圖像灰度化,根據cv2.HoughCircles()函數通過霍夫變換裁切圖像,去除多余圖像,保留所需圖像特征信息;利用三維塊匹配濾波算法對獲取的圖像進行降噪處理,相較于原圖,降噪處理后的圖像峰值信噪比為36.48 dB,結構相似指數為0.88;采用Canndy算法檢測并提取土壤邊界像素,利用改進的最小二乘法對土壤邊界進行線性擬合,從而獲取土壤動態休止角。

3) 采用土壤動態休止角測量裝置對不同轉速下的機收膜雜中土壤動態休止角進行測定,研究分析了轉速對動態休止角的影響規律并分別構建了線性擬合模型和多項式擬合模型,通過對比分析得到了擬合度較高的擬合模型,即多項式擬合模型(R2為95.72%,均方差為0.061),在此基礎上,獲取了土壤動態休止角最優測量轉速,即7 r/min。該結果可為相關散粒物料動態休止角的測量提供了參考依據。