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

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

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

0 引 言

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

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

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

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

1 基本理論

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

1.1 ECAH模型

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

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

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

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

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

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

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

1.2 超聲衰減實(shí)驗(yàn)

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

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

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

圖2 超聲測(cè)量裝置Fig.2 Ultrasonic measuring device

2 多輸出回歸算法

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

2.1 梯度提升決策樹(shù)

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

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

2.2 GBDT回歸組合

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

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

3 特征選取

3.1 實(shí)驗(yàn)信號(hào)

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

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

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

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

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

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

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

3.2 特征選取

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

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

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

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

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

為了避免特征數(shù)據(jù)量綱不同導(dǎo)致在回歸樹(shù)預(yù)測(cè)時(shí)帶來(lái)的偏差，需要對(duì)數(shù)據(jù)進(jìn)行預(yù)處理。對(duì)特征為種類的頻率進(jìn)行 onehot編碼，將 660、750、830 kHz三種頻率轉(zhuǎn)化為數(shù)值1、2、3；對(duì)其他連續(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ù)集組成

本文的樣品通過(guò)篩分法進(jìn)行配置，即先通過(guò)孔篩篩選出在各個(gè)粒徑區(qū)間范圍內(nèi)的顆粒粒徑，然后由精密天平稱重組合成近似各種分布類型的懸移質(zhì)樣品。以對(duì)照組的三種懸移質(zhì)樣本為例進(jìn)行說(shuō)明，三個(gè)樣本分別服從近似單峰分布、均勻分布和不規(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í)驗(yàn)共配置了36組實(shí)驗(yàn)樣品，其中33組用于制作訓(xùn)練數(shù)據(jù)集，其余3組用于制作驗(yàn)證數(shù)據(jù)集。對(duì)所有樣品進(jìn)行三種頻率下的超聲衰減實(shí)驗(yàn)，然后針對(duì)每一次實(shí)驗(yàn)分別提取 10項(xiàng)特征作為數(shù)據(jù)集的輸入，并且將 17種粒徑分布的值作為數(shù)據(jù)集的輸出。

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

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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

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