丁 聰，邸南思

(上?？睖y設計研究院有限公司，上海 200434)

1 概述

水泵、水輪機特性曲線是泵站、電站設計的重要基礎資料，廣泛用于原型機組參數(shù)計算、水力過渡過程分析[1-2]、經(jīng)濟效益評價[3]等。日常設計需要從期刊、書籍等資料中收集水泵、水輪機特性曲線。但資料多為紙質、圖片格式，需進行人工讀點、數(shù)值化。而人工數(shù)值化費時費力，且主觀誤差較大。為減少圖像曲線數(shù)據(jù)提取工作量，提高曲線數(shù)據(jù)提取準確性，已有圖像曲線數(shù)據(jù)提取的軟件開發(fā)，如Getdata Dragh Digitizer、FindGraph、Digitize。但此類軟件一般為手動選點或無法較好地適用于含縱橫坐標軸網(wǎng)格線的圖像。

Matlab因其強大的矩陣處理能力，在圖像處理方面具有較大優(yōu)勢。MATLAB圖像處理研究已經(jīng)涉及到醫(yī)學[4]、食品[5]、交通[6]等各行各業(yè)。并且已有研究利用Matlab進行圖像曲線數(shù)據(jù)提取，但多數(shù)仍停留在手動選點、曲線擬合或僅適用含單條曲線的圖像，或未考慮曲線走向[7-9]。為此提出利用已識別曲線像素，尋求周圍曲線像素并結合曲線斜率不可突變作為判斷條件，甄別取舍尋求到的周圍像素點進行曲線自動追蹤，并根據(jù)像素坐標與真實坐標一一對應關系進行曲線數(shù)據(jù)換算，以避免手工操作的繁瑣和曲線追蹤的偏離。

2 圖像處理難點及思路

若曲線與背景絕對分開，可降低曲線追蹤難度。收集的曲線圖像多為彩色圖像，因此，需先將圖像灰度化、二值化處理。

若曲線很細，僅為曲線中心線，一可提高曲線數(shù)據(jù)提取準確性；二可方便曲線追蹤。因此，需先對曲線細化處理。

收集的曲線圖像可能有褪色、變質、偏轉等問題，致使影響曲線提取的因素較多(噪聲大)，需先進行圖像預處理(如：去噪、旋轉、裁剪)。

多數(shù)收集的曲線圖像內含有坐標軸、網(wǎng)格線，且曲線、網(wǎng)格線交叉較多。較多的交叉點增加了曲線追蹤難度，因此，需將圖像中的坐標軸、網(wǎng)格線先行刪除。

處理好的曲線圖像便可利用遍歷法及斜率預測法進行曲線的追蹤識別，并可根據(jù)坐標、像素點位置對應關系進行曲線數(shù)據(jù)的換算導出。

水泵特性曲線一般比較復雜、且有多條，智能分條識別難度較大，可采用人為指定需識別曲線減小識別難度。

依據(jù)上述圖像處理難點及解決方法，曲線處理可采用圖1流程。

圖1 曲線提取處理流程示意

3 圖像處理過程

3.1 圖像讀入

圖像處理前，需將圖像讀入Matlab工作空間。Matlab提供有圖像文件的讀入和顯示函數(shù)imread、imshow，用來讀取與顯示如bmp、jpg、tiff、hdf等格式文件。

3.2 圖像灰度化及灰度增強

圖像灰度化就是由RGB三通道數(shù)據(jù)的彩色圖像變?yōu)閱瓮ǖ罃?shù)據(jù)圖像。灰度圖像上每個像素的顏色值稱為灰度，范圍為0～255。圖像灰度化常用方法有3種：① 分量法，將彩色圖像的R、G、B三分量的亮度之一作為灰度圖像的灰度值；② 最大值法，將彩色圖像的R、G、B三分量亮度最大值作為灰度圖像的灰度值；③ 加權平均法，將彩色圖像的R、G、B三分量以不同的權重進行加權平均。Matlab提供的灰度化函數(shù)為rgb2gray，使用的是加權平均法。由于人眼對綠色敏感最高，對藍色敏感最低，故采用心理學灰度公式：

Gray= 0.114B+ 0.587G+ 0.299R

(1)

灰度增強是指按一定變換關系逐點改變原圖中每一個像素的灰度值，以改善圖像視覺效果?；叶仍鰪姵Ｓ玫姆椒ㄓ芯€性灰度變化、非線性灰度變化、直方圖均衡化等。直方圖是圖像的重要統(tǒng)計特征，是表示數(shù)字圖像中每一灰度級與該灰度級出現(xiàn)的頻數(shù)間的統(tǒng)計關系。直方圖均衡化是通過灰度變化將一幅圖像轉換為另一幅均衡直方圖，即在每個灰度級上都具有相同的像素點數(shù)的過程。采用直方圖均衡化可使原始圖像的灰度直方圖從比較集中的某個灰度區(qū)間變成在全部灰度范圍內的均勻分布。Matlab直方圖均衡化函數(shù)為histeq。

3.3 圖像去噪

圖像去噪是減少數(shù)字圖像中噪聲的過程。圖像去噪的方法有均值濾波、自適應維納濾波、中值濾波、形態(tài)學噪聲濾波、小波變換等方法。小波變換在時域、頻域同時具有良好的局部化性質和多分辨率，可有效地把信號和噪聲區(qū)別開來。

小波分析進行圖像去噪復原主要有3個步驟[10]：① 對圖像信號進行小波分解。選擇合適的小波和恰當?shù)姆纸鈱哟?，然后對含有噪聲的圖形進行N層小波分解。② 對經(jīng)過層次分解后的高頻系數(shù)進行閾值量化。對分解的每一層，選擇一個恰當?shù)拈撝?，并對該層高頻系數(shù)進行閾值量化處理。③ 利用二維小波重構圖像信號。同樣，根據(jù)小波分解后的第N層近似(低頻系數(shù))和經(jīng)過閉值量化處理后的各層細節(jié)(高頻系數(shù))，來計算二維信號的小波重構。

Matlab小波去噪函數(shù)為ddencmp(閾值選取函數(shù))、wdencmp(去噪函數(shù))。

3.4 圖像二值化

灰度閾值法是最古老的圖像二值化方法[11]。選取適當閾值，像素灰度值超過閾值則分配到最大灰度，低于閾值分配到最小灰度。

閾值的選取方法可分為全局閾值法和局部閾值法。全局閾值法指對整幅圖像選取單一的閾值進行二值化，在物體和背景的灰度差別明顯時效果較好。曲線圖像中曲線與背景的灰度差別較大，全局閾值法較為合適。典型的全局閾值法包括Otsu 法、全局迭代法、最大墑方法等。Otsu法主要思想是從圖像灰度直方圖中把灰度值的集合用閾值T分為兩類，然后根據(jù)兩個類別的平均值方差(類間方差)和各類的方差(類內方差)的比為最大來確定閾值T。

假定給定的圖像有L個灰度級，設閾值為T，把具有T以下灰度值的像素和具有比T大的灰度值的像素分為兩類，設為類1和類2，類1中的像素數(shù)量為W1(T)，灰度值的均值為M1(T)，方差為σ1(T)，同樣類2中的像素量為W2(T)，灰度值的均值為M2(T)，方差為σ2(T)。全體像素的平均值為MT，則類內方差由下式計算：

(2)

類間方差由下式計算：

(3)

Matlab最大類間差法全局閾值求解函數(shù)為graythesh。

3.5 曲線細化

圖像的細化算法較多，數(shù)學形態(tài)學細化算法提取的圖像骨架較連續(xù)，斷點少，對噪聲不敏感，且提取的邊緣較光滑[12-13]。細化過程根據(jù)擊中或擊不中變換定義，圖像二值文件A除去被結構元素序列B擊中的元素。

A?B=A-(A#B)=A∩(A#B)C

(4)

圖像細化是對圖像文件進行基于結構元素序列{B}的細化

A?{B}=((…((A?B1)?B2)…)?Bn

(5)

即連續(xù)使用{B}對A細化。

Matlab圖像數(shù)學形態(tài)學函數(shù)為bwmorph，其中細化操作為thin(細化)和skel(骨架化)。

3.6 圖像旋轉、裁剪、坐標軸網(wǎng)格線刪除

由于曲線圖像中的坐標軸、網(wǎng)格線一般均為直線?？赏ㄟ^尋找計算坐標軸、網(wǎng)格線所在直線與水平、豎直存在的夾角進行圖像旋轉。

直線尋找可采用Hough變換，其原理是利用圖像空間和Hough參數(shù)空間的點、線對偶性，把圖像空間中的檢測問題轉換到參數(shù)空間，通過在參數(shù)空間里進行簡單的累加統(tǒng)計，然后在Hough參數(shù)空間尋找累加器峰值的方法檢測直線。

Hough變換的參數(shù)空間為極坐標空間，極坐標中的方程為：

ρ=xcosθ+ysinθ

(6)

其中ρ是直線離原法線距離，θ是法線對x軸的角度(如圖2所示)。

圖2 直角坐標、極坐標對應關系示意

可見直線的Hough變換在極坐標中的一點，而點的Hough變化是正弦函數(shù)。由此可知在圖像空間中共線的點對應在參數(shù)空間中相交的線，參數(shù)空間中相交于同一點的所有曲線在圖像空間里都有共線的點與之對應。利用相交點的累加統(tǒng)計便可檢測直線。

依次利用Matlab圖像函數(shù)hough(Hough變換)、houghpeaks(提取峰值點)、houghlines(提取峰值點對應的直線)提取圖像文件中的直線。然后根據(jù)夾角采用imrotate函數(shù)進行圖像旋轉。

圖像中若含有縱橫網(wǎng)格線，對曲線追蹤干擾較大?？上葘⒕W(wǎng)格線刪除，具體方法為：首先利用Hough變化識別圖像中的直線，然后對直線長度進行排序，最后根據(jù)給定的縱、橫網(wǎng)格線條數(shù)刪除識別出的長度較長的直線。

3.7 曲線區(qū)域特征點獲得

圖像文件中曲線數(shù)量可能為多條，若采用遍歷算法尋找圖像中曲線，其計算量較大，并且曲線易混淆。因此，提出在手動選擇需要提取數(shù)據(jù)曲線的基礎上進行曲線追蹤。

MATALB函數(shù)getpts可獲得用戶在當前圖像上使用鼠標點擊點的坐標，利用此函數(shù)便可獲得鼠標點擊曲線處的像素點位置。為避免點擊出現(xiàn)偏差，可以點擊像素點位置為圓心逐漸向外尋找曲線像素點。

3.8 曲線追蹤識別

獲得區(qū)域特征點后可利用遍歷法尋找曲線像素點，并確定為曲線追蹤的起始點。然后利用遍歷法在已識別曲線像素周圍尋求下一曲線像素并結合曲線斜率不可突變作為判斷條件，甄別取舍尋求到的周圍像素點進行曲線追蹤。曲線斜率可通過綜合比較統(tǒng)計方向碼、已得曲線擬合計算端點斜率獲得。

方向碼是邊界跟蹤中描述行走方向規(guī)定的概念，可規(guī)定如圖3所示：東、東北、北、西北、西、西南、南、東南，分別使用數(shù)字0、1、2、3、4、5、6、7表示。統(tǒng)計已得曲線最近幾個像素點的方向碼，并計算其平均值，便可預測下一像素點的方向碼。

圖3 方向碼示意

曲線擬合方法可利用Matlab多項式擬合函數(shù)polyfit對已得曲線像素點進行擬合，然后利用擬合曲線得到端點斜率或下一像素點位置。

3.9 數(shù)據(jù)換算導出

圖像是由像素點構成，像素點位置與實際坐標值(x軸、y軸)存在一一對應關系。因此，將識別出曲線的像素點位置與坐標范圍進行比例換算便可得到曲線坐標值。假設曲線圖像像素尺寸(長度×寬度)為(A，B)；曲線的像素點位置，行、列值為(a、b)；曲線X軸坐標范圍為(Xmin，Xmax)，Y軸坐標范圍為(Ymin，Ymax)。則數(shù)據(jù)換算公式如下。

(7)

最后將換算得到的曲線數(shù)據(jù)導出。Matlab導出的數(shù)據(jù)文件可為mat(Matla數(shù)據(jù))、Excel等形式。可利用uiputfile函數(shù)保存為mat文件，xlswrite函數(shù)保存為Excel文件。

4 水泵特性曲線提取實例

根據(jù)上述圖像處理過程，利用MATLAB的GUI編輯界面，編寫曲線數(shù)據(jù)提取軟件包。利用此軟件包進行水泵特性曲線提取。編制的軟件界面如圖4所示。

圖4 曲線提取軟件GUI界面示意

假設需處理的水泵特性曲線如圖5所示，二值化后圖像如圖6所示，細化并刪除坐標網(wǎng)格線的圖像如圖7所示。圖像文件中有多條曲線，可依次選擇需提取曲線，得到曲線數(shù)據(jù)繪制圖形如圖8所示，其中紅色曲線是依據(jù)提取出的數(shù)據(jù)點重新繪制而成。

圖5 水泵特性曲線原圖像示意

圖6 二值化處理后圖像示意

圖7 細化并刪除坐標網(wǎng)格線后示意

圖8 讀取所有曲線后圖像示意

比較原始圖片和提取繪制的圖像可以得出，提取繪制圖像與原圖像無偏差。

5 結語

1)利用已識別曲線像素尋求周圍曲線像素并結合曲線斜率不可突變作為判斷條件，甄別取舍尋求到的周圍像素點可保證曲線追蹤的連續(xù)性。

2)利用MATLAB軟件進行圖像曲線數(shù)據(jù)提取，能克服肉眼判斷帶來的誤差，提取數(shù)據(jù)與原圖像一致，大大提高了準確性及提取效率。

3)MATLAB語句簡單，GUI軟件包可實現(xiàn)軟件的迅速調用及推廣應用。