基于近紅外高光譜成像技術(shù)的 寧夏羊肉產(chǎn)地鑒別

王 靖,丁佳興,郭中華,*,何鳳杰,梁曉燕

(1.寧夏大學(xué)物理與電子電氣工程學(xué)院,寧夏銀川 750021;2.寧夏大學(xué)農(nóng)學(xué)院,寧夏銀川 750021)

羊肉肌纖維較細(xì),肉質(zhì)細(xì)膩,口感俱佳;并且含有礦物質(zhì)、維生素以及多種人體必需的氨基酸等營養(yǎng)成分,具有很高的營養(yǎng)價值,也因其高蛋白質(zhì)、低膽固醇、低脂肪的特點(diǎn)深受人們喜愛[1]。羊肉的品質(zhì)跟其產(chǎn)地來源密切相關(guān),同一品種不同產(chǎn)地羊肉品質(zhì)大不相同,如寧夏鹽池灘羊肉因其低脂低膽固醇、鮮香細(xì)嫩和口感釀?wù)奶攸c(diǎn),產(chǎn)品價格比其他產(chǎn)地羊肉每公斤高出10元以上[2]。近年來市場上利用其他產(chǎn)地的羊肉冒充鹽池灘羊獲取暴利的事件時常發(fā)生,因此研究羊肉產(chǎn)地快速鑒別方法對打擊非法商販,保護(hù)消費(fèi)者利益,維護(hù)產(chǎn)品口碑意義重大。

目前,國內(nèi)外對羊肉產(chǎn)地鑒別的方法主要包括:水解氨基酸[3]、同位素鑒別方法[4]、指紋圖譜溯源技術(shù)[5]。目前,應(yīng)用光譜技術(shù)主要對肉類品種鑒別較多。其中,劉瑋等[6-7]利用傅里葉變換近紅外光譜分析方法分別建立了新鮮豬肉和新鮮雞肉的判別模型;王文秀[8]等利用近紅外光譜技術(shù)對冷鮮豬肉和解凍肉進(jìn)行了鑒別;牛曉穎[9]等人利用近紅外光譜技術(shù)對驢肉、牛肉、豬肉和羊肉建立了判別模型;Cozzolino等[10]分別對豬肉、羊肉和雞肉的樣本進(jìn)行品種判別分析。然而,利用光譜技術(shù)對肉類產(chǎn)地判別較少。高光譜的定性分析研究主要應(yīng)用于物質(zhì)判別分析方面,其融合了電子學(xué)、光學(xué)、計算機(jī)科學(xué)和信息處理,將光譜技術(shù)和二維成像技術(shù)結(jié)合在一起,因此,具有快速、無損、連續(xù)多波段、高分辨率和光譜圖像合一的特點(diǎn)[11-14]。

K最近鄰分類算法(K-nearest neighbor,KNN)是一種典型的較流行、非參數(shù)、有效的分類方法,屬于傳統(tǒng)統(tǒng)計模式識別算法[15]。偏最小二乘判別分析(Partial Least Squares-Discriminant Analysis,PLS-DA)基于 PLS回歸方程,其作為判別分析的常用方法,具有簡單高效的特點(diǎn)[16-18]。不同生長環(huán)境的羊肉從外觀雖然很難辨別,但其所處的生長環(huán)境(氣候、土壤、水質(zhì))有所不同,會導(dǎo)致主要化學(xué)成分(蛋白質(zhì)、脂肪、水分等)的結(jié)構(gòu)和含量存在一定差異,這種差異可在高光譜上得到反映。故本文研究利用近紅外高光譜成像技術(shù)對寧夏同一品種不同產(chǎn)地的羊肉進(jìn)行鑒別,采集不同產(chǎn)地的3類羊肉樣品原始光譜,利用SG卷積平滑(Savitzky-Golay smoothing,SGS)、標(biāo)準(zhǔn)正態(tài)變化(standardized normal variate,SNV)、面積歸一化法(Area-Normalization)、多元散射校正(Multiplicative scatter correction,MSC)預(yù)處理方法處理,結(jié)合立偏最小二乘判別分析PLS-DA模型優(yōu)選出最優(yōu)預(yù)處理方法;然后利用連續(xù)投影算法(Successive Projection Algorithm,SPA)、競爭性正自適應(yīng)加權(quán)算法(Competitive Adaptive Reweighted Sampling,CARS)和無信息變量消除法(Uninformative Variable Elimination,UVE)選取的特征波長,分別建立PLS-DA、KNN判別分析模型,優(yōu)選出最優(yōu)模型,從而實(shí)現(xiàn)固原、鹽池、銀川羊肉產(chǎn)地的鑒別。

1 材料與方法

1.1 材料與儀器

實(shí)驗(yàn)羊 選自寧夏鹽池、銀川、固原,羊被屠宰后,取羊后腿肌肉放于低溫保鮮盒運(yùn)至實(shí)驗(yàn)室,用手術(shù)刀片去除樣本表面的脂肪和肌膜,整形切塊,肉樣大小約為20 mm×30 mm×10 mm,分別取各產(chǎn)地樣本75個,總樣本共225個,用密封袋密封、編號。

近紅外高光譜成像系統(tǒng) 900～1700 nm,光譜分辨率3 nm,256個波段,如圖1所示,該系統(tǒng)由高光譜成像儀(Imspector N17E)、CCD相機(jī)(Zelos-285GV)、鹵鎢燈(SIA-LS-TDIF)、電控位移平臺(PSA200-11-X)組成。

圖1 近紅外高光譜成像系統(tǒng)Fig.1 Near infrared hyperspectral imaging system

1.2 實(shí)驗(yàn)方法

1.2.1 高光譜圖像的采集 為保證采集到的圖像清晰不失真,經(jīng)過預(yù)實(shí)驗(yàn)嘗試,最終將成像光譜儀的曝光時間設(shè)為10 ms,物鏡高度為385 mm,電控位移平臺速度為15 mm/s。

由于光源在各波段下強(qiáng)度分布不均勻及箱體中暗電流存在會導(dǎo)致圖像光譜中含有較大的噪聲。因此要對高光譜圖像進(jìn)行黑白校正[19],以消除噪聲的影響。黑白校正公式為:

式(1)

其中:I校正后的漫反射光譜圖像;R樣本原始的漫反射光譜圖像;D暗圖像;W白板的漫反射圖像。獲得樣本高光譜圖像后,利用ENVI 4.6軟件選取整塊肉作為感興趣區(qū)域(Region of Interest,ROI),計算出每張ROI的平均反射光譜,并將其作為該樣本的反射光譜。

1.2.2 光譜預(yù)處理 由于儀器噪音、樣本表面凹凸不平等因素對光譜的影響,需要將原始光譜數(shù)據(jù)進(jìn)行預(yù)處理。選取SGS、SNV、Area-Normalization、MSC預(yù)處理方法處理,結(jié)合偏最小二乘判別分析PLS-DA模型選出最優(yōu)預(yù)處理方法。

1.2.3 樣本劃分方法 樣本集需要分成校正集和預(yù)測集,用于建立模型及驗(yàn)證預(yù)測模型的預(yù)測能力。采用Galvao等提出的SPXY(Sample Set Partitioning Based on Joint X-Y distance)法劃分樣本。其是在KS(Kennard Stone)法的基礎(chǔ)上發(fā)展起來的,其優(yōu)點(diǎn)是將變量X和Y均考慮在內(nèi),劃分得到的校正集代表性強(qiáng),模型性能好。

1.2.4 特征波長提取 由于全光譜是包含所有光譜變量,信息冗余,選用適當(dāng)?shù)奶卣鞑ㄩL提取方法剔除不相關(guān)或者非線性變量,實(shí)現(xiàn)用少數(shù)關(guān)鍵變量代替全光譜,達(dá)到降低模型運(yùn)算量和復(fù)雜度、提高模型穩(wěn)定性和預(yù)測準(zhǔn)確性的目的。選用SPA、CARS和UVE選取特征波長,分別建立PLS-DA、KNN判別分析模型。

1.3 數(shù)據(jù)處理

光譜預(yù)處理在The Unscrambler X 10.4上實(shí)現(xiàn),其余算法在Matlab R2014a上完成。

2 結(jié)果與分析

2.1 預(yù)處理方法選擇

本文采用SGS(SG-smoothing)、Area-Normalization、SNV和MSC預(yù)處理方法對原始光譜處理,建立PLS-DA模型,每次建立PLS-DA模型之前,首先找到最佳主成分?jǐn)?shù),方法如下:設(shè)定最大主成分?jǐn)?shù)20、數(shù)據(jù)不縮放、交叉驗(yàn)證組數(shù)為10,原始光譜的交互驗(yàn)證的錯誤率隨主成分?jǐn)?shù)的變化如圖2,根據(jù)交叉驗(yàn)證中的錯誤率對應(yīng)的主成分?jǐn)?shù)(16)確定為最佳。同樣的方法找到不同預(yù)處理數(shù)據(jù)的最佳主成分?jǐn)?shù)并建立PLS-DA模型,結(jié)果統(tǒng)計如表1所示。

圖2 原始光譜的交叉驗(yàn)證中不同主成分?jǐn)?shù)下的錯誤率Fig.2 Error rate of different principal components in cross validation of original spectrum

預(yù)處理方法主成分?jǐn)?shù)交互驗(yàn)證中最小錯誤率準(zhǔn)確率OriginalSpectrum160.16440.9378SG-smoothing140.12440.9244Area-Normalization120.15560.9333SNV160.14670.9244MSC150.15110.9289

從表1可以看出,原始光譜的PLS-DA模型主成分?jǐn)?shù)為16時,交互驗(yàn)證錯誤率最低為0.1644,經(jīng)過上述4種預(yù)處理方法建立的PLS-DA模型交互驗(yàn)證的最小的錯誤率均低于原始光譜建立的PLS-DA模型,SG-smoothing的PLS-DA模型錯誤率最低為0.1244;但SG-smoothing、SNV和MSC的PLS-DA模型準(zhǔn)確率均小于0.93，明顯低于原始光譜PLS-DA模型;經(jīng)SG-smoothing、Area-Normalization和MSC預(yù)處理后的PLS-DA模型最優(yōu)主成分?jǐn)?shù)低于原始光譜PLS-DA模型,一般來說,主成分?jǐn)?shù)越少,模型越穩(wěn)定,Area-Normalization的PLS-DA模型的主成分?jǐn)?shù)最少(12)。綜合主成分?jǐn)?shù)、交互驗(yàn)證錯誤率和模型準(zhǔn)確率判定Area-Normalization為最優(yōu)預(yù)處理方法。

2.2 樣本劃分

采用SPXY方法對225個樣本按照3∶1比例劃分校正集和預(yù)測集,劃分結(jié)果如表2所示。

表2 利用SPXY法劃分樣本結(jié)果Table 2 Results of sample division by SPXY method

2.3 特征波長提取

由于全光譜包含所有光譜信息,數(shù)據(jù)量大、信息冗余且存在共線性變量,以全光譜建模會增加建模的復(fù)雜度,降低計算速度,影響建模效果。所以,本文選擇3種方法提取特征波長,從全光譜中去除共線性變量,挑選有用變量,最終實(shí)現(xiàn)用關(guān)鍵波長代替全光譜建模,從而達(dá)到對全光譜降維、提高模型穩(wěn)定性和建模速度的目的。

2.3.1 利用SPA提取特征波長 應(yīng)用SPA選取波長數(shù)時,設(shè)置波長范圍是1～25,數(shù)據(jù)不縮放,計算每個變量數(shù)下的RMSECV,根據(jù)RMSECV最小,確定最佳特征波長數(shù)。不同波長數(shù)下的RMSECV值如圖3所示,可以看出當(dāng)波長數(shù)為17時,RMSECV最小,故確定經(jīng)SPA提取出17個最佳特征波長，依次為:924.6、927.5、930.5、933.5、936.5、939.5、945.4、951.4、960.3、987.2、996.1、1008.0、1046.8、1124.3、1151.1、1297.1、1675.6 nm。

圖3 SPA中不同有效波長數(shù)下的RMSECVFig.3 Variation of RMSECV with number of effective wavelengths in SPA

2.3.2 利用CARS提取特征波長 在應(yīng)用CARS提取之前,首先確定PLS模型中最優(yōu)主成分?jǐn)?shù),設(shè)置最大主成分?jǐn)?shù)15,蒙特卡洛采樣次數(shù)1000,從總樣本提取3/4作為校正集,得到不同主成分?jǐn)?shù)下的RMSECV值,結(jié)果如圖4所示。從圖4中可知,當(dāng)主成分?jǐn)?shù)為7時,RMSECV值最小(0.5447),故最優(yōu)主成分?jǐn)?shù)為7。然后設(shè)定CARS參數(shù):蒙特卡洛采樣次數(shù)為1000,主成分?jǐn)?shù)為7,交叉驗(yàn)證組數(shù)為10。對全波段樣本光譜進(jìn)行篩選,羊肉產(chǎn)地鑒別篩選過程見圖5a。結(jié)果表明,隨著變量篩選過程的進(jìn)行被挑選的波長數(shù)逐漸下降,下降趨勢由快變慢,體現(xiàn)了波長變量篩選的粗選與精選。圖5b 為變量篩選過程中交互驗(yàn)證錯誤率的變化趨勢。結(jié)果顯示,交互驗(yàn)證錯誤率隨著篩選過程的進(jìn)行先下降后上升,蒙特卡洛(Monte Carlo,MC)采樣次數(shù)為36～61時,交互驗(yàn)證錯誤率均為0.5171且最小,MC采樣次數(shù)繼續(xù)增加后,交互驗(yàn)證錯誤率隨挑選變量數(shù)減少而增大,再結(jié)合篩選過程中波長變量回歸系數(shù)變化趨勢(見圖5c),“*”所對應(yīng)的位置為36次MC采樣,所以應(yīng)用CARS方法在MC采樣為36次時,交互驗(yàn)證錯誤率最小,挑選出40個特征波長,分別為:924.6、927.5、966.3、984.2、1002、1005、1008、1011、1049.7、1052.7、1055.7、1067.6、1070.6、1079.5、1154.1、1160、1163、1177.9、1183.9、1198.8、1213.7、1216.6、1219.6、1222.6、1225.6、1228.6、1231.5、1288.2、1297.1、1303.1、1312、1315、1320.9、1323.9、1392.5、1401.4、1505.7、1663.7、1675.6、1678.6 nm。

圖4 RMSECV隨PLS主成分?jǐn)?shù)的變化規(guī)律Fig.4 Variation of RMSECV with number of principal components of PLS

圖5 羊肉產(chǎn)地鑒別的CARS特征波長篩選過程Fig.5 Process of CARS characteristic wavelength selection for identification of mutton origin注:a為羊肉產(chǎn)地鑒別篩選過程； b為變量篩選過程中交互驗(yàn)證錯誤率的變化趨勢； c為篩選過程中波長變量回歸系數(shù)變化趨勢。

2.3.3 利用UVE提取特征波長 首先根據(jù)PLS交互驗(yàn)證模型中RMSECV最小確定PLS的最佳主成分?jǐn)?shù),本研究中當(dāng)主成分?jǐn)?shù)為7時,RMSECV最小,因此將主成分確定為7。設(shè)定交互驗(yàn)證組數(shù)為25,隨機(jī)波長數(shù)為256,運(yùn)行UVE計算256個輸入變量的穩(wěn)定性結(jié)果如圖6所示。圖中豎線左側(cè)為256個光譜變量,右邊為256個隨機(jī)變量。兩條水平虛線為變量選擇閾值(10.20),閾值的選擇標(biāo)準(zhǔn)為隨機(jī)變量穩(wěn)定性最大絕對值的99%。兩條虛線之間的信息被認(rèn)定是無用信息,兩條虛線外的被認(rèn)定是有用信息,其對應(yīng)的波長被挑選出來。最終,應(yīng)用UVE共選出121個特征波長,這些波長主要分布在999.1～1401.5,1502.8～1562.4,1633.9～1678.6 nm區(qū)間。

圖6 PLS主成分?jǐn)?shù)為7時UVE的穩(wěn)定性分布曲線Fig.6 Stability distribution of variables selected by UVE at seven principal components in PLS

2.4 建模算法的比較與分析

根據(jù)原始光譜和預(yù)處理基于PLS-DA羊肉產(chǎn)地鑒別模型建模對比發(fā)現(xiàn)SNV法為最優(yōu)預(yù)處理方法,分別建立全波段光譜(Full Spectrum,FS)和3種特征波長提取方法基于PLS-DA和KNN的羊肉產(chǎn)地鑒別模型。

2.4.1 KNN建模結(jié)果 在運(yùn)行KNN算法時,首先選擇合適的K值,K值的大小影響模型的穩(wěn)定性和預(yù)測能力。圖7為全光譜交互驗(yàn)證的錯誤率與K值大小的分布圖,選擇交互驗(yàn)證最低錯誤率對應(yīng)的K值,所以選擇K=4建立KNN模型。同樣的方法確定CARS、SPA、UVE特征波長的K值分別為5、6、5,分別建立對應(yīng)的KNN模型,結(jié)果統(tǒng)計如表3。

圖7 KNN算法中K值的選擇Fig.7 K values selection for KNN algorithm

表3 基于不同特征波長挑選方法的KNN羊肉產(chǎn)地鑒別模型Table 3 KNN model based on different characteristic wavelength selection method for identification of geographical origins of mutton

表4 基于不同特征波長挑選方法的PLS-DA羊肉產(chǎn)地鑒別模型Table 4 PLS-DA model based on different characteristic wavelength selection method for identification of geographical origins of mutton

從表3可以看出,基于全光譜和特征波長提取方法建立的KNN模型中校正集的正確率均小于80%,預(yù)測集正確率均小于70%,說明模型穩(wěn)定性和預(yù)測能力較差。對比全光譜和特征波長建立的KNN模型發(fā)現(xiàn),SPA-KNN和UVE-KNN的模型效果明顯遜于FS-KNN模型,而CARS-KNN模型校正集正確率為70.83%略小于FS-KNN模型,且預(yù)測集正確率為68.42%大于FS-KNN模型,說明CARS-KNN模型校正性能略差于FS-KNN模型,但模型預(yù)測能力優(yōu)于FS-KNN模型。所以,CARS特征波長提取方法優(yōu)于SPA和UVE,得到的40個特征波長包含了大量的有用信息,可以代替全光譜建模,達(dá)到減少模型的復(fù)雜性,提高建模速度的目的。

2.4.2 PLS-DA建模結(jié)果 建立PLS-DA判別模型時,首先需要確定最優(yōu)主成分?jǐn)?shù),選擇RMSECV最小時對應(yīng)的主成分?jǐn)?shù)確定為最優(yōu)。本研究建立PLS-DA模型時,確定參數(shù):數(shù)據(jù)不縮放、交叉驗(yàn)證組數(shù)為10。應(yīng)用CARS、SPA、UVE方法挑選特征波長建立PLS-DA羊肉產(chǎn)地鑒別模型并與全光譜對比,模型效果統(tǒng)計如表4。

從表4可以看出,在3種特征波長提取方法中,SPA提取的特征波長17個為最少,SPA-PLS-DA模型主成分?jǐn)?shù)為8低于FS-PLS-DA模型,說明通過SPA確實(shí)降低了模型的復(fù)雜性,但其校正集、預(yù)測集準(zhǔn)確率分別為77.98%、66.67%,明顯低于FS-PLS-DA模型,說明SPA方法剔除了過多變量信息,其中包括部分有用信息;經(jīng)過UVE提取出121個特征波長,保留的變量較多,UVE-PLS-DA模型主成分?jǐn)?shù)為12略低于FS-PLS-DA模型,但其校正集、預(yù)測集準(zhǔn)確率明顯低于FS-PLS-DA模型,說明該模型保留了部分無用信息,篩選能力不強(qiáng);經(jīng)過CARS方法挑選出40個特征波長,就數(shù)量上而言對全光譜起到了降維的作用,雖然CARS-PLS-DA模型的校正集準(zhǔn)確率為90.48%，低于FS-PLS-DA模型的94.05%,預(yù)測集準(zhǔn)確率為84.21%，略低于FS-PLS-DA模型的85.96%,說明CARS-PLS-DA模型的準(zhǔn)確性稍差,但是其最佳主成分?jǐn)?shù)為10小于FS-PLS-DA模型(13),通常情況下模型主成分?jǐn)?shù)越小,模型越穩(wěn)定,所以綜合考慮模型的穩(wěn)定性和準(zhǔn)確性,說明CARS特征波長提取方法較好,提取得到的特征波長能代表全光譜建立PLS-DA模型。

3 結(jié)論

對比預(yù)處理光譜的PLS-DA模型建模效果可知,經(jīng)過面積歸一化法預(yù)處理方法最好;對經(jīng)過面積歸一化法的光譜數(shù)據(jù)挑選特征波長,SPA、CARS、UVE算法提取特征波長分別為17、40、121個;分別建立基于特征波長的KNN、PLD-DA判別模型,結(jié)果表明所建立的KNN模型效果較差,3種特征波長中利用CARS提取的特征波長建模效果最佳,代替全光譜建立PLS-DA判別模型是可行的;對比模型效果,CARS-PLS-DA為最優(yōu)模型,校正集正確率90.48%,預(yù)測集正確率84.21%。

[1]烏仁張嘎. 羊肉的營養(yǎng)價值分析及飼養(yǎng)管理水平對羊肉品質(zhì)的影響[J]. 當(dāng)代畜禽養(yǎng)殖業(yè),2017(2):3-5.

[2]王家云,王松磊,賀曉光,等. 基于NIR高光譜成像技術(shù)的灘羊肉內(nèi)部品質(zhì)無損檢測[J]. 現(xiàn)代食品科技,2014(6):257-262,249.

[3]劉興勇,林濤,劉宏程,等.基于水解氨基酸分析山羊肉的產(chǎn)地溯源[J].現(xiàn)代食品科技,2013,29(11):2788-2792.

[4]孫淑敏,郭波莉,魏益民,等.多種穩(wěn)定性同位素(C、N、H)分析在羊肉產(chǎn)地溯源中的應(yīng)用[J].2011年第四屆國際食品安全高峰論壇論文集.

[5]孫淑敏. 羊肉產(chǎn)地指紋圖譜溯源技術(shù)研究[D].楊凌:西北農(nóng)林科技大學(xué),2012.

[6]Dubois A J,Lewisb E N,Jr F S F,et al. Bacterial identification by near-infrared chemical imaging of food-specific cards[J]. Food Microbiology,2005,22(6)577-583.

[7]Swatland HJ.Observations on rheological,electrical,and optical changes during rigor development in pork and beef[J].J Anim Sci,1997,75:975-985.

[8]王文秀,彭彥昆. 基于近紅外光譜的冷鮮肉-解凍肉的判別研究[J]. 食品安全質(zhì)量檢測學(xué)報,2014(3):754-760.

[9]牛曉穎,邵利敏,董芳,等. 基于近紅外光譜和化學(xué)計量學(xué)的驢肉鑒別方法研究[J]. 光譜學(xué)與光譜分析,2014(10):2737-2742.

[10]Cozzolino D,Murray I.Identiflcation of animal meat muscles by visle and near infrared reflectance spectroscopy[J].Swks Society of Food Science and Technology,2004,37(4):447-452.

[11]王雷,喬曉艷,董有爾,等.高光譜圖像技術(shù)在農(nóng)產(chǎn)品檢測中的應(yīng)用進(jìn)展[J].應(yīng)用光學(xué),2009,30(4):639-645.

[12]逛全勝,陳蕾,王平,等.島光譜遙感技術(shù)在海洋研究的應(yīng)用及展望機(jī)[J].海洋湖沼通報,2008(3):168-173.

[13]田有文,王曉娟.基于高光譜圖像技術(shù)的農(nóng)產(chǎn)品品質(zhì)無損檢測[J].農(nóng)機(jī)化研究,2009(10):220-222.

[14]楊國鵬,余旭初,馮伍法,高光譜遙感技術(shù)的發(fā)展與應(yīng)用現(xiàn)狀[J].測繪通報,2008(10):1-4.

[15]劉應(yīng)東,?；菝?基于k-最近鄰圖的小樣本KNN分類算法[J].計算機(jī)工程,2011,37(9):198-200.

[16]唐果,田曠達(dá),李祖紅,等. 近紅外光譜結(jié)合PLS-DA劃分煙葉等級[J].農(nóng)業(yè)工程學(xué)報,2013(4):60-62.

[17]王逸之,董文淵,李永和,等.基于近紅外光譜結(jié)合PLS-DA法的野外竹種識別技術(shù)研究[J].竹子研究匯刊,2014,33(4):16-20.

[18]楊忠,任海青,江澤慧.PLS-DA法判別分析木材生物腐朽的研究[J].光譜學(xué)與光譜分析,2008,28(4):793-796.

[19]ElMasry G,Wang N,ElSayed A. Hyperspectral imaging for nondestructive determination of some quality attributes for strawberry[J]. Journal of Food Engineering,2007,81:98-107.