多輸出回歸算法在超聲粒徑分布反演中的研究

應(yīng)啟帆，謝代梁，徐志鵬，徐 雅，劉鐵軍，黃震威

(中國計量大學(xué)浙江省流量計量技術(shù)研究重點(diǎn)實(shí)驗室，浙江 杭州 310018)

0 引 言

水體懸移質(zhì)粒徑分布測量在河流、湖泊、海洋和生態(tài)系統(tǒng)等領(lǐng)域中有著重要價值，其準(zhǔn)確測量有助于認(rèn)識懸移質(zhì)沉降、擴(kuò)散規(guī)律[1]，對水土流失規(guī)律研究[2]和水污染治理[3]等領(lǐng)域的發(fā)展有重大意義。

粒徑測量方法按照工作原理[4]，可分為篩分法、顯微鏡法、沉降法、光散射法和超聲法等。其中，超聲法測量具有非浸入式、穿透性強(qiáng)、操作簡單等優(yōu)點(diǎn)，受到了廣泛關(guān)注。Shukla等[5]使用超聲波譜在線監(jiān)測乙酰氨基酚結(jié)晶過程中的顆粒生長；章維等使用超聲波衰減譜和相速度譜及Epstein-Carhart-Allegra-Hawley (ECAH)模型反演計算懸濁液顆粒粒徑分布[6]；Ding等[7]使用超聲頻譜測量了水分散體系中重油的粒徑；李燁明等采用人工蜂群算法優(yōu)化超聲粒徑分布計算結(jié)果[8]。

對粒徑分布的測量，傳統(tǒng)學(xué)者多聚焦在選用更精準(zhǔn)的模型和更加優(yōu)化的算法來提升粒徑分布結(jié)果的準(zhǔn)確度。蘇明旭等[9]通過構(gòu)造與貝塞爾(Bessel)函數(shù)有關(guān)的變量改造系數(shù)矩陣，降低了矩陣求解的條件數(shù)，拓展了模型適用的頻率和顆粒粒徑。Silva等[10]使用 6～14 MHz頻率光譜法和擴(kuò)展多重散射模型測量體積分?jǐn)?shù)為10%～50%的葵花籽油包水型乳液的液滴尺寸分布。由于機(jī)器學(xué)習(xí)的發(fā)展和測量數(shù)據(jù)的提升，有學(xué)者嘗試用數(shù)據(jù)來尋找顆粒粒徑分布和輸入特征之間的關(guān)系，通過訓(xùn)練機(jī)器學(xué)習(xí)模型來進(jìn)行預(yù)測。Thompson等[11]利用基因表達(dá)編程和人工神經(jīng)網(wǎng)絡(luò)開發(fā)了一種新模型，用于預(yù)測活躍施工現(xiàn)場和露天采礦作業(yè)中裸露的土壤表面產(chǎn)生的雨水徑流中的總懸浮泥沙粒徑分布；Sch?fer等[12]使用卷積神經(jīng)網(wǎng)絡(luò)和人工圖像進(jìn)行訓(xùn)練，通過分割重疊顆粒來預(yù)測多相流中液滴的粒徑分布；Manee等[13]提出了一種深度學(xué)習(xí)傳感器，可對晶體生長過程中的顆粒進(jìn)行在線監(jiān)測。

在本文中，我們提出了一種結(jié)合多輸出回歸算法和超聲衰減實(shí)驗的粒徑分布測量方法。首先根據(jù)超聲衰減實(shí)驗和 ECAH模型等先驗信息獲取實(shí)驗信號并提取特征，然后利用梯度提升決策樹算法預(yù)測單種粒徑，并組合全部粒徑預(yù)測模型構(gòu)建多輸出回歸模型預(yù)測樣品的粒徑分布。最后將近似單峰分布、均勻分布和不規(guī)則分布的三種樣品粒徑分布預(yù)測結(jié)果與篩分法確定的粒徑分布進(jìn)行對比，驗證多輸出回歸模型的準(zhǔn)確性，可為粒徑分布測量提供一種新的參考方式。

1 基本理論

超聲法測量顆粒粒徑主要根據(jù)理論模型和優(yōu)化算法來獲得結(jié)果。理論模型介紹了超聲波在不同粒徑顆粒液體中的衰減情況，也成為了在機(jī)器學(xué)習(xí)算法中的特征基礎(chǔ)。接下來主要介紹 ECAH模型和超聲衰減實(shí)驗，了解粒徑測量的過程。

1.1 ECAH模型

ECAH模型說明了超聲波在水體中受到的粘滯力、熱損失和聲散射等因素的影響，可以描述超聲波具體的衰減過程。它忽略了顆粒之間的相互作用，適用于5%體積濃度以下、粒徑在1 000 um以下的懸移質(zhì)溶液環(huán)境，具體說明如式(1)[14-15]所示：

式中：α為超聲波在水中的衰減系數(shù)(單位：Np·m-1)；kc為連續(xù)介質(zhì)中的波數(shù)；φ為懸移質(zhì)溶液體積濃度，(單位：g·L-1)；R為懸移質(zhì)顆粒半徑(單位：μm)；An為懸移質(zhì)溶液中的超聲衰減系數(shù)，不同的下標(biāo)n代表不同種衰減情況。

顆粒粒徑的求解過程是第一類 Fredholm 方程計算問題，可利用超聲衰減系數(shù)、濃度等參數(shù)之間的相關(guān)關(guān)系來求解，是典型的反問題，參數(shù)關(guān)系如式(2)[16]所示：

式中：q是顆粒粒度的頻度分?jǐn)?shù)，是要求解的結(jié)果；K是消聲系數(shù)，是關(guān)于頻率和粒徑區(qū)間的函數(shù)。

通過對式(2)的離散化處理，可以轉(zhuǎn)化成矩陣形式來求解[16]：

式中： fi為不同測量頻率(單位：kHz)；Rj為不同粒徑區(qū)間(單位：Np·m-1)；α( fi, Rj)為在特定頻率和特定粒徑區(qū)間時的衰減系數(shù)(單位：Np·m-1)。

等式左邊是聲衰減過程，是超聲波在不同頻率和粒徑里的衰減過程，通過實(shí)驗獲得。等式右邊是理論部分，可以結(jié)合ECAH模型和各項物性參數(shù)來進(jìn)行求解。式(1)～(3)描述了采用 ECAH模型求解顆粒粒徑分布的過程，是先驗信息。后續(xù)的特征部分就是根據(jù)這部分理論來選取與超聲衰減過程密切的相關(guān)參數(shù)作為數(shù)據(jù)集的輸入。

1.2 超聲衰減實(shí)驗

超聲衰減實(shí)驗采用兩個聚焦換能器進(jìn)行聲電信號轉(zhuǎn)換，通過一發(fā)一收的工作方式，對電壓幅值進(jìn)行測量，可探測到超聲波在懸移質(zhì)溶液中的衰減過程。聚焦換能器采用的材料為壓電材料PZT82；球面內(nèi)徑為 100 mm；投影直徑為 90 mm，換能器中心頻率為 750 kHz，單個聚焦換能器具體的結(jié)構(gòu)如圖1所示。

圖1 超聲換能器結(jié)構(gòu)Fig.1 Ultrasonic transducer structure

發(fā)射換能器接收信號發(fā)生器的電信號，將其轉(zhuǎn)化成聲信號發(fā)射，聲信號在穿過懸移質(zhì)溶液后到達(dá)接收換能器。在這個過程中超聲波會發(fā)生反射，從而多次到達(dá)接收換能器，采用示波器測量接收換能器提供的電壓信號，就可看出超聲波衰減過程，具體的測量裝置如圖2所示。

圖2 超聲測量裝置Fig.2 Ultrasonic measuring device

2 多輸出回歸算法

粒徑分布反演問題在機(jī)器學(xué)習(xí)中的表現(xiàn)是多輸入多輸出回歸問題，根據(jù)多個輸入找到對應(yīng)的多個目標(biāo)預(yù)測結(jié)果。本文采用組合梯度提升決策樹(Gradient Boosting Decision Tree, GBDT)算法來解決這一問題。首先利用 GBDT算法來找到單種粒徑的回歸預(yù)測結(jié)果，然后將所有GBDT算法組合起來，對整體粒徑分布進(jìn)行預(yù)測。

2.1 梯度提升決策樹

梯度提升決策樹是一種集成學(xué)習(xí)算法，算法的基礎(chǔ)是分類回歸決策樹，利用損失函數(shù)的梯度方向決定后續(xù)回歸樹的優(yōu)化方向，通過貪心算法每一次找到下一步更優(yōu)的回歸樹，最后對所有回歸樹進(jìn)行集成，輸出最后結(jié)果，就構(gòu)成了梯度提升決策樹模型。算法訓(xùn)練[17]的步驟如下所示：

算法每次對負(fù)梯度和輸入特征進(jìn)行擬合，構(gòu)建新的回歸樹來修正模型，其問題類型是多輸入單輸出的回歸預(yù)測模型，因此不適用于粒徑分布問題，只適用于單種粒徑分布的預(yù)測，需要將多個GBDT模型進(jìn)行組合來解決問題。

2.2 GBDT回歸組合

對每種粒徑，單獨(dú)訓(xùn)練一個GBDT模型，可以根據(jù)輸入特征和粒徑的數(shù)值進(jìn)行回歸預(yù)測。為了預(yù)測樣品的粒徑分布，需要所有模型組合，每個輸出對應(yīng)粒徑的結(jié)果，最后按照對整體粒徑的權(quán)重確定整體粒徑的頻度分布，多輸出回歸模型示意圖如圖3所示。

圖3 多輸出回歸模型示意圖Fig.3 Diagram of the multi-output regression model

3 特征選取

3.1 實(shí)驗信號

為了獲得實(shí)驗信號提取特征來制作數(shù)據(jù)集，搭建了如圖4所示的實(shí)驗測試系統(tǒng)。首先接收端聚焦換能器接收功率放大器放大的猝發(fā)波信號。超聲信號到達(dá)接收換能器時部分轉(zhuǎn)換成電壓幅值，部分反射回接收聚焦換能器。最后呈現(xiàn)在示波器上的是超聲波到達(dá)接收端聚焦換能器時的電壓幅值，最終將數(shù)據(jù)讀取并存儲在電腦中。實(shí)驗器材：發(fā)射和接收聚焦換能器、功率放大器、信號發(fā)生器、電磁攪拌器和示波器。

圖4 實(shí)驗系統(tǒng)Fig.4 Experimental system

實(shí)驗獲得的電壓信號經(jīng)去噪后的波形如圖5所示。同時考慮到換能器中心頻率的影響，提取其在頻域下的幅值，如圖6所示。

圖5 去噪后接收信號波形圖Fig.5 Waveforms of the received signals after denoising

圖6 去噪后接收信號頻譜圖Fig.6 Spectrums of the received signals after denoising

圖5中三個回波表示超聲波三次到達(dá)接收聚焦換能器時的電壓幅值，聲衰減系數(shù)就可通過這三個峰峰值計算而出，如式(7)[8]所示：

式中：αw為清水中的衰減系數(shù)(單位：Np·m-1)；V0是發(fā)射換能器電壓幅值(單位：V)，V1為接收換能器電壓幅值(單位：V)；L為兩個聚焦換能器之間的距離(單位：m)。

3.2 特征選取

多輸出回歸模型的特征來源于先驗信息，即顆粒的 ECAH模型和超聲衰減實(shí)驗的信號特征。接下來以各個參數(shù)與粒徑之間的關(guān)系來說明選擇它作為特征的理由。

根據(jù)ECAH模型反演來看，式(1)說明了粒徑與超聲衰減系數(shù)和體積濃度之間的關(guān)系，由圖5中的三次回波可計算出兩個衰減系數(shù)，因此將兩個聲衰減系數(shù)和體積濃度作為特征。式(3)說明了頻率可以影響聲衰減系數(shù)的計算，同時考慮到超聲換能器中心頻率對信號的影響，在選擇特征時加入頻率這一參數(shù)。

圖5說明了聲衰減變化的電壓衰減過程。圖6說明了圖5中三次回波在各頻率下的幅值，結(jié)果表明在中心頻率附近的傅里葉變換后幅值最大。三次回波的電壓峰峰值和三次回波的傅里葉變換最大幅值，是超聲衰減的關(guān)鍵因素，因此也選擇其作為特征。

綜上所述，選擇了懸移質(zhì)溶液體積濃度、兩次聲衰減系數(shù)、頻率、三個電壓峰峰值幅值和三次回波傅里葉變換后的最大幅值共 10個屬性作為數(shù)據(jù)集的特征，具體說明如表1所示。

表1 數(shù)據(jù)集輸入特征Table 1 Input features of data

為了避免特征數(shù)據(jù)量綱不同導(dǎo)致在回歸樹預(yù)測時帶來的偏差，需要對數(shù)據(jù)進(jìn)行預(yù)處理。對特征為種類的頻率進(jìn)行 onehot編碼，將 660、750、830 kHz三種頻率轉(zhuǎn)化為數(shù)值1、2、3；對其他連續(xù)數(shù)據(jù)進(jìn)行標(biāo)準(zhǔn)化處理，將數(shù)據(jù)映射到[0,1]，轉(zhuǎn)換公式如式(8)[18]所示：

式中：x*是數(shù)據(jù)組中數(shù)據(jù)的轉(zhuǎn)換結(jié)果，x是數(shù)據(jù)組中數(shù)據(jù)初始值，xmin是數(shù)據(jù)組中數(shù)據(jù)最小值，xmax是數(shù)據(jù)組中數(shù)據(jù)最大值。

4 結(jié)果分析

4.1 數(shù)據(jù)集組成

本文的樣品通過篩分法進(jìn)行配置，即先通過孔篩篩選出在各個粒徑區(qū)間范圍內(nèi)的顆粒粒徑，然后由精密天平稱重組合成近似各種分布類型的懸移質(zhì)樣品。以對照組的三種懸移質(zhì)樣本為例進(jìn)行說明，三個樣本分別服從近似單峰分布、均勻分布和不規(guī)則分布，粒徑區(qū)間有 17種。單種粒徑區(qū)間質(zhì)量占樣本總質(zhì)量的百分比為該種粒徑區(qū)間的頻度，將 17種粒徑區(qū)間頻度進(jìn)行組合，就構(gòu)成了樣本的粒徑分布，具體的粒徑組成如表2所示。

表2 懸移質(zhì)顆粒樣品粒徑分布組成Table 2 Composition of the particle size distribution of suspended particles

本次實(shí)驗共配置了36組實(shí)驗樣品，其中33組用于制作訓(xùn)練數(shù)據(jù)集，其余3組用于制作驗證數(shù)據(jù)集。對所有樣品進(jìn)行三種頻率下的超聲衰減實(shí)驗，然后針對每一次實(shí)驗分別提取 10項特征作為數(shù)據(jù)集的輸入，并且將 17種粒徑分布的值作為數(shù)據(jù)集的輸出。

本文將訓(xùn)練數(shù)據(jù)集用于訓(xùn)練多輸出回歸模型，同時將驗證數(shù)據(jù)集輸入模型，即可得到粒徑分布預(yù)測結(jié)果。因為考慮到有三種頻率下的實(shí)驗，所以為了得到更加精確的結(jié)果，對三次預(yù)測結(jié)果進(jìn)行均值化處理，作為真正的粒徑分布結(jié)果。

4.2 預(yù)測結(jié)果和誤差分析

中位徑D50代表一個樣品的累計粒度分布百分?jǐn)?shù)達(dá)到50%時所對應(yīng)的粒徑，中位徑誤差ε是常用的粒徑誤差計算公式，如式(9)[19]所示：

式中：D50,act為實(shí)際的中位徑(單位：μm)；D50,pre是模型預(yù)測確定的中位徑(單位：μm)。

首先利用訓(xùn)練數(shù)據(jù)集訓(xùn)練多輸出回歸模型，然后利用驗證數(shù)據(jù)集和模型預(yù)測三種樣品的粒徑分布結(jié)果。預(yù)測得到的單峰分布、均勻分布和不規(guī)則分布樣品篩分結(jié)果和多輸出回歸模型預(yù)測結(jié)果對比圖分別如圖7～9所示。

圖 7～9的橫坐標(biāo)是粒徑區(qū)間的組合，柱狀圖分別為篩分法確定的該種粒徑頻率分布和模型預(yù)測的粒徑頻率分布，縱坐標(biāo)為粒徑的頻率分布，即單種粒徑區(qū)間質(zhì)量占樣品質(zhì)量的百分比。所有粒徑的結(jié)果進(jìn)行組合就構(gòu)成了樣本粒徑分布的結(jié)果。接下來統(tǒng)計了每種樣品的單個粒徑的最大相對誤差、決定系數(shù)和中位徑誤差，如表3所示。

表3 不同樣本中位徑和測量誤差Table 3 Median diameters and measurement errors of different samples

圖7 單峰分布樣品兩種方法測量結(jié)果對比Fig.7 Comparison of measurement results of two methods for unimodal samples

三種樣品的單種顆粒相對誤差范圍都在±10%以內(nèi)，只有單峰分布的個別樣本在粒徑權(quán)重較小時出現(xiàn)了較大的偏差，為 12.74%。這表明，單種粒徑占總體粒徑分布的權(quán)重質(zhì)量較小時會導(dǎo)致相對較大的偏差。中位徑的誤差都較小，表明預(yù)測分布與實(shí)際分布較為一致，可以很好地反應(yīng)懸移質(zhì)顆粒粒徑分布組成。

圖8 均勻分布樣品兩種方法測量結(jié)果對比Fig.8 Comparison of measurement results of two methods for uniformly distributed samples

圖9 不規(guī)則分布樣品兩種方法測量結(jié)果對比Fig.9 Comparison of measurement results of two methods for irregularly distributed samples

5 結(jié) 論

本文通過 ECAH模型和超聲衰減實(shí)驗等先驗信息，了解傳統(tǒng)粒徑分布測量的過程，并說明了顆粒粒徑與體積濃度、聲衰減系數(shù)之間存在的耦合關(guān)系，使用機(jī)器學(xué)習(xí)算法對粒徑分布進(jìn)行預(yù)測來為粒徑分布測量提供一種新的方法。

粒徑分布預(yù)測作為一種多輸入多輸出回歸問題，采用組合梯度提升決策樹來進(jìn)行預(yù)測。首先使用 GBDT對單種粒徑進(jìn)行預(yù)測，然后組合所有粒徑預(yù)測過程，就構(gòu)成了對粒徑分布結(jié)果的預(yù)測。然后選擇了與粒徑關(guān)系密切的體積濃度、二次聲衰減系數(shù)和超聲衰減實(shí)驗中最明顯的三次回波電壓峰峰值和傅里葉變換最大幅值共 10項參數(shù)作為數(shù)據(jù)集的特征。最后，采用近似單峰分布、均勻分布和不規(guī)則分布的三種樣品作為預(yù)測結(jié)果，三種樣品中單種粒徑相對誤差范圍在±10%以內(nèi)，中位徑誤差分別為0.07%、-0.1%和-2.2%，表明預(yù)測樣品結(jié)果與實(shí)際樣品結(jié)果較為一致，可以很好地測量樣品的粒徑分布。

本文利用的特征根據(jù)先驗信息進(jìn)行選取，也可以進(jìn)一步通過對超聲信號的分解來選取與粒徑關(guān)系密切的特征來進(jìn)行訓(xùn)練。同時，使用更大規(guī)模的數(shù)據(jù)集來進(jìn)行訓(xùn)練，可得到某個區(qū)域更準(zhǔn)確的粒徑分布結(jié)果。本文的嘗試，可以為粒徑分布測量問題提供一種新的思路。