王文顥， 劉振丙， 雒志超

(1.桂林電子科技大學 電子工程與自動化學院,廣西 桂林 541004；2.桂林電子科技大學 計算機與信息安全學院,廣西 桂林 541004)

由于近紅外技術(shù)的非接觸、速度快等特點，被廣泛使用在農(nóng)產(chǎn)品[1]、食品[2]和藥物[3]的質(zhì)量檢測。近紅外光譜定量分析法[4]利用已知的樣本及性質(zhì)通過化學計量學方法可建立定量分析模型，用于未知樣本的性質(zhì)預測。

在應用中，由于源機與目標機的測量光譜存在偏差，所以不能直接使用源機上建立的定量分析模型。若在目標機上再次測量多組樣本數(shù)據(jù)，建立新的定量分析模型又費時費力。為了解決這一問題，模型傳遞技術(shù)[5]應運而生。模型傳遞技術(shù)首先建立源機的測量光譜與目標機的測量光譜間關(guān)系的數(shù)學模型，再通過該模型對目標機上的光譜進行轉(zhuǎn)換。用在源機建立的定量分析模型，分析傳遞后的光譜便可獲得更準確的預測結(jié)果。

模型傳遞分為有標樣法和無標樣法。有標樣法需要使用標準樣本在多機上分別檢測，而無標樣法不需要標樣集。有標樣方法包括Shenk’s算法[6]、直接校正法(DS)[7]、分段校正法(PDS)[8]等；無標樣方法包括有限脈沖響應算法(FIR)[9]。其中直接校正法和分段校正法是主流的模型傳遞方法。直接校正法直接使用源機和目標機上測量光譜整體來構(gòu)建校正模型；分段校正法法是直接校正法的改進，增加了源機光譜波長點和目標機光譜波長點近鄰的限制窗口，窗口的大小直接影響校正的結(jié)果。

動態(tài)時間規(guī)整算法(dynamic time warping,簡稱 DTW)[10-11]在語音分析中經(jīng)常使用。因為人聲的隨機性，同一單詞的2次發(fā)音并不相同，動態(tài)時間規(guī)整算法可以尋找兩組語音序列中幀的關(guān)系，計算序列的相似度。

提出的模型傳遞算法基于動態(tài)時間規(guī)整算法。該算法首先計算雙機光譜上波長點的相關(guān)距離。再使用動態(tài)時間規(guī)整算法尋找使整體的相關(guān)距離最小的源機光譜與目標機光譜的關(guān)聯(lián)關(guān)系。最后根據(jù)關(guān)聯(lián)關(guān)系，構(gòu)建回歸模型。

1 算法描述

設Sm(i,j)為源機的光譜矩陣，Ss(i,j)為目標機的光譜矩陣。其中S(i,:)是第i個光譜樣本的波長數(shù)據(jù)，S(S,j)是每個光譜樣本的第j個波長點的吸光度。

1.1 計算波長點的相關(guān)距離

近紅外光譜設備受到自身老化、環(huán)境溫濕度等條件的影響，源機與目標機測得的光譜之間不但會產(chǎn)生基線漂移，而且波長點也會產(chǎn)生偏移，故需要建立源機光譜Sm與目標機光譜Ss的對應關(guān)系。當源機光譜第p個波長點Sm(:,p)與目標機光譜第q個波長點Ss(:,q)對應時，向量間的相關(guān)距離小，則相關(guān)系數(shù)大。相關(guān)系數(shù)使用皮爾遜相關(guān)系數(shù)計算,

(1)

根據(jù)相關(guān)系數(shù)計算相關(guān)距離,

DSm(:,p),Ss(:,q)=1-ρSm(:,p),Ss(:,q)。

(2)

1.2 動態(tài)時間規(guī)整算法

對于給定的2個序列X=(x1,x2,,xN)和Y=(y1,y2,,yM)，使用動態(tài)時間規(guī)整得到2個序列的匹配關(guān)系。當2個序列間的距離最小時，匹配關(guān)系是最佳的。兩個序列的匹配關(guān)系如圖1所示。

圖1 兩個序列的匹配關(guān)系

光譜數(shù)據(jù)Sm(:,i)可以看作序列X中xi，Ss(:,j)可以看作序列Y中yi。

1.2.1 計算代價矩陣

首先構(gòu)造一個代價矩陣C∈RN×M,它由序列中的元素的距離構(gòu)成。將其元素定義為：

ci,j=Dxi,yj,xi∈X,yj∈Y。

(3)

樸素的動態(tài)時間規(guī)整算法使用Dxi,yj=‖xi-yj‖作為距離。這里更關(guān)注元素的相關(guān)性，故

Dxi,yj=1-ρxi,yj。

(4)

1.2.2 計算規(guī)整路徑

動態(tài)時間規(guī)整方法可以依據(jù)代價矩陣得到一條規(guī)整路徑。圖2為一條從序列X到序列Y的規(guī)整路徑。

圖2 最佳規(guī)整路徑

設P=(p1,p2,,pL)為規(guī)整路徑，pl=(xi,yj)為路徑中的第l個元素，表示序列X的元素xi和序列Y的元素yj相關(guān)。規(guī)整路徑必須滿足3個條件：

1)邊界條件。起點只能是(x1,y1)點，終點只能是(xN,yM)點。

2)順序限制。序列中的點必須依次逐個匹配，禁止跳過序列中的點。

3)步長限制。一次只能移動一步進行匹配，方向為向右、向上、向右上，例如當前點為(xi,yj)，那么下一步只能為(xi+1,yj)、(xi,yj+1)或(xi+1,yj+1)。

使用規(guī)整路徑的代價和作為算法的損失函數(shù)，

(5)

當pl=(xi,yj)時，c(pl)=c(xi,yi)=cxi,yj。動態(tài)時間規(guī)整算法尋找代價和最小的規(guī)則路徑：

DTW(X,Y)=min{CP(X,Y)}。

(6)

構(gòu)造累計代價矩陣Dacc，其計算式為：

當i=1,j=1時，

Dacc(i,j)=c(xi,yj)；

(7)

當i=1,2≤j≤M時，

Dacc(i,j)=Dacc(i,j-1)+c(xi,yj)；

(8)

當2≤i≤N,j=1時，

Dacc(i,j)=Dacc(i-1,j)+c(xi,yj)；

(9)

當2≤i≤N,2≤j≤M時，

Dacc(i,j)=min{Dacc(i-1,j-1),Dacc(i-1,j),

Dacc(i,j-1)}+c(xi,yj)。

(10)

當前點累計的代價是前一步的最小累計代價與當前點代價的和。累計代價矩陣構(gòu)造完成后，從終點pL=(xN,yM)按最小累計代價回溯到起始點p1=(x1,y1)即可得到規(guī)整路徑。

1.2.3 根據(jù)規(guī)整路徑建立校正模型

源機與目標機測量光譜的波長點的關(guān)系由規(guī)整路徑確定。根據(jù)此關(guān)系來構(gòu)建校正模型。

當Sm(:,i)與Ss(:,j)單一對應時，構(gòu)造一元回歸方程：

Sm(:,i)=bi,M+1+bi,jSs(:,j)。

(11)

當Sm(:,i)與Ss(:,j-n)到Ss(:,j+m)相匹配時，構(gòu)建多元回歸方程：

Sm(:,i)=bi,M+1+bi,j-nSs(:,j-n)++

bi,j+mSs(:,j+m)。

(12)

式(11)、(12)通過最小二乘法求得系數(shù)。完成每個波長點的計算后，構(gòu)造轉(zhuǎn)移矩陣F∈R(M+1)×N，其元素為回歸方程的系數(shù)，定義為：

(13)

(14)

2 實驗部分

2.1 數(shù)據(jù)集

為驗證基于動態(tài)時間規(guī)整算法的模型傳遞方法，實驗數(shù)據(jù)使用的藥物數(shù)據(jù)集和玉米數(shù)據(jù)集來自網(wǎng)頁(http://www.eigenvector.com/data/tablets/index.html)和網(wǎng)頁(http://www.eigenvector.com/data/Corn/index.html)。

藥物數(shù)據(jù)集來自2臺近紅外光譜儀對654片藥片樣本采樣獲得的數(shù)據(jù)。每個樣本的光譜波長從600～1898 nm，采樣間隔2 nm，共650維。采集到的654對光譜數(shù)據(jù)被分為訓練155對、測試460對和驗證40對。此外，還包含藥物質(zhì)量、藥物硬度、活性物質(zhì)含量3種性質(zhì)供分析使用。

玉米數(shù)據(jù)集是3臺近紅外光譜儀對80個玉米樣本采樣獲得的數(shù)據(jù)。每個樣本的光譜波長從1100～2498 nm，采樣間隔2 nm，共700維。此外，還包含油、水、淀粉和蛋白質(zhì)含量4種性質(zhì)供分析使用。

2.2 建立定量模型

對于藥物數(shù)據(jù)集，先對光譜數(shù)據(jù)做零均值化處理。選擇儀器1作為源機，在訓練集上使用PLS算法對藥物數(shù)據(jù)創(chuàng)建定量分析模型，對藥物的3種性質(zhì)進行預測。

對于玉米數(shù)據(jù)集，先對光譜數(shù)據(jù)做零均值化處理。選擇m5作為源機，使用K-S算法劃分訓練集和測試集。其中訓練集樣本50個，測試機樣本30個。使用PLS算法對訓練集創(chuàng)建定量分析模型，對玉米的4種性質(zhì)進行預測。

使用10折交叉驗證創(chuàng)建定量分析模型，最優(yōu)的PLS模型主因子數(shù)由平均預測均方差確定，結(jié)果如表1所示。

表1 使用源機光譜建立的PLS模型

2.3 使用動態(tài)時間規(guī)整算法建立傳遞模型

使用動態(tài)時間規(guī)整算法得到規(guī)整路徑，以此得到波長點的對應關(guān)系。

藥物數(shù)據(jù)集以波長600～920 nm為例，對應關(guān)系和規(guī)整路徑如圖3、4所示。圖3中圓形點為源機波長點，星形點為目標機波長點，二者的細線是其對應關(guān)系?？梢钥闯?，在700 nm附近的峰值存在變化，動態(tài)時間規(guī)整算法能找到準確的對應關(guān)系和規(guī)整路徑，并建立回歸方程。

圖3 波長點的匹配關(guān)系

圖4 兩臺近紅外光譜儀器的波長的規(guī)整路徑

玉米數(shù)據(jù)集里的源機與目標機偏差較小，其光譜波長點單一對應。動態(tài)時間規(guī)整算法也能找到準確的規(guī)整路徑，得出準確的關(guān)聯(lián)關(guān)系，回歸方程部分退化為一元線性回歸。

求解上述回歸方程，得到轉(zhuǎn)移矩陣F。

3 結(jié)果與討論

3.1 模型傳遞平均光譜差異對比

使用動態(tài)時間規(guī)整算法得到轉(zhuǎn)移矩陣F，通過式(14)可計算轉(zhuǎn)移后的目標機光譜。藥物數(shù)據(jù)集上的源機光譜、目標機光譜、傳遞后的目標機光譜如圖5所示。其中目標機光譜和源機光譜差異較大，目標機光譜在傳遞后和源機光譜的差異則顯著減小。

圖5 模型傳遞后的平均光譜

表2為直接校正算法、分段校正算法、動態(tài)時間規(guī)整算法得到的傳遞光譜與未傳遞光譜的平均差異(ARMS)。由表2可看到，光譜的ARMS在應用了模型傳遞算法后較未傳遞的光譜均有下降。其中動態(tài)時間規(guī)整算法傳遞后的光譜的ARMS整體較小，意味著算法的適用性更好，傳遞的結(jié)果和原機光譜匹配度高、有效。

表2 各算法傳遞光譜的平均差異(ARMS)

3.2 與傳統(tǒng)模型傳遞算法對比

將傳遞前的光譜和使用直接校正算法、分段校正算法、動態(tài)時間規(guī)整算法傳遞后的光譜輸入如表1所示的PLS回歸模型中預測樣本的性質(zhì)，其預測結(jié)果的預測均方偏差(RMSEP)如表3、4所示。其中分段校正算法的窗口大小通過交叉驗證確定。

表3 在藥物數(shù)據(jù)集上使用各算法傳遞光譜的預測均方差(RMESP)

結(jié)果表明，直接使用目標機光譜進行預測會由于其和源機光譜的差異產(chǎn)生不客觀的預測結(jié)果。在使用模型傳遞算法后，預測結(jié)果均有提升。在不同的預測項目中有不同的算法取得優(yōu)勢，使用動態(tài)時間規(guī)整算法得到的傳遞后光譜在8項預測中的4項取得領先，故準確度優(yōu)于直接校正算法和分段校正算法，但由于其運算量大于直接校正算和分段校正法算法，故需要更多的計算時間。

表4 在玉米數(shù)據(jù)集上使用各算法傳遞光譜的預測均方差(RMESP)

4 結(jié)束語

提出一種新的近紅外光譜模型傳遞算法。該算法依次得到規(guī)整路徑、匹配關(guān)系、回歸模型、傳遞矩陣。使用傳遞矩陣便可對目標機光譜進行傳遞，以獲得更準確的定量分析的預測結(jié)果。在2個數(shù)據(jù)集上的實驗證實了算法的有效性，但計算時間略長。