表征超聲衰減譜粒度的改進(jìn)和聲搜索算法?

蔣 瑜 曲佩玙 賈 楠 蘇明旭

(上海理工大學(xué) 動(dòng)力工程多相流與傳熱重點(diǎn)實(shí)驗(yàn)室 上海 200093)

0 引言

顆粒物廣泛存在于各個(gè)領(lǐng)域，能源領(lǐng)域水煤漿高效燃燒，食品行業(yè)食品結(jié)晶、均質(zhì)化，醫(yī)藥行業(yè)藥物粒徑大小及其在生物體內(nèi)運(yùn)輸?shù)榷忌婕暗筋w粒表征技術(shù)[1?4]。典型顆粒測(cè)量技術(shù)有沉降法[5]、顯微鏡法[6]、動(dòng)態(tài)光散射法[7]、超聲衰減譜法等，其中超聲衰減譜法具有非侵入性、樣品無(wú)需稀釋、裝置成本低等優(yōu)點(diǎn)。超聲波在介質(zhì)中傳播產(chǎn)生衰減，攜帶粒徑信息，因此通過(guò)實(shí)驗(yàn)獲得超聲衰減譜，結(jié)合McClements 模型[8]、Epstein-Carhart-Allegra-Hawley(ECAH)模型[9?10]，可進(jìn)一步反演得到顆粒粒徑分布。

反演算法可分為獨(dú)立模式和非獨(dú)立模式兩類，獨(dú)立模式無(wú)需預(yù)設(shè)分布，求解離散方程組得到粒徑分布；非獨(dú)立模式通常假定顆粒服從高斯分布、羅辛-拉姆勒(Rosin-Rammler,R-R)分布或?qū)?shù)正態(tài)分布，然后定義理論預(yù)測(cè)與實(shí)測(cè)衰減譜的誤差作為目標(biāo)函數(shù)來(lái)構(gòu)建優(yōu)化問(wèn)題，求取最佳粒徑分布。人工智能算法的迅速發(fā)展，為推動(dòng)此類方法提供了巨大的發(fā)展?jié)摿ΑＤ壳胺仟?dú)立模式下粒徑反演中人工魚群算法(Artificial fish school algorithm,AFSA)[11]、蟻群算法(Ant colony optimization algorithm,ACO)[12]和遺傳算法(Genetic algorithm,GA)[13]等得到了較好的研究。與前述算法相比，和聲搜索算法(Harmony search algorithm,HS)作為一種啟發(fā)式全局搜索算法[14]，具有結(jié)構(gòu)簡(jiǎn)單、全局搜索能力強(qiáng)等優(yōu)點(diǎn)。歐陽(yáng)海濱等[15]研究和聲搜索算法和反向?qū)W習(xí)策略，提出反向?qū)W習(xí)產(chǎn)生新和聲與之前的和聲進(jìn)行末位淘汰競(jìng)爭(zhēng)機(jī)制，并應(yīng)用于優(yōu)化熱交換器和減速器設(shè)計(jì)問(wèn)題；Alireza 等[16]利用和聲搜索算法優(yōu)化配電網(wǎng)絡(luò)中配電機(jī)組配置，以降低系統(tǒng)功耗、改善系統(tǒng)電壓分布、提高電壓穩(wěn)定性；Hussein 等[17]利用和聲搜索算法和主成分分析相結(jié)合(Hybrid principal component analysis,HPCA)，利用HPCA 加快對(duì)圖像處理；Shams等[18]將和聲搜索算法應(yīng)用油藏工程歷史匹配問(wèn)題；陳濤[19]針對(duì)特征基因選擇問(wèn)題提出用概率決定對(duì)最差和聲或最好和聲進(jìn)行變異來(lái)優(yōu)化和聲搜索算法。上述在一定程度上提升了算法的性能，但是未能解決算法容易陷入局部最優(yōu)解、求解精度不高的問(wèn)題。

為此，提出一種超聲衰減非獨(dú)立模式反演的改進(jìn)，其核心為通過(guò)和聲搜索算法與擬牛頓法(Broyden-Fletcher-Goldfarb-Shanno,BFGS)[20]并行計(jì)算，合理調(diào)節(jié)和取舍全局搜索與局部搜索能力，使算法兼顧全局搜索效率和逼近精度。針對(duì)R-R分布、正態(tài)分布、對(duì)數(shù)正態(tài)分布函數(shù)的多種粒徑分布情況進(jìn)行模擬，并通過(guò)實(shí)驗(yàn)超聲衰減譜法的反演和圖像法結(jié)果對(duì)比分析，充分驗(yàn)證了改進(jìn)算法的有效性、抗噪性。

1 目標(biāo)函數(shù)構(gòu)建

當(dāng)超聲波通過(guò)含顆粒物介質(zhì)時(shí)，由于顆粒與聲波的作用而產(chǎn)生衰減，可從理論上較為全面地考慮聲波和顆粒相互作用的多種機(jī)制，多分散顆粒系中聲衰減系數(shù)α表達(dá)式為

其中，Kext是模型定義的消聲系數(shù)，dV/dD是顆粒體積頻度分布，為待求值。式(1)在數(shù)學(xué)上屬第一類Fredholm 積分方程，該方程目前沒有理論解析解，將式(1)離散，可寫成矩陣形式：

其中，G為聲衰減分布列向量，A為衰減系數(shù)矩陣，向量W為被測(cè)顆粒的粒徑分布。式(2)是典型的病態(tài)方程，常采用最優(yōu)化算法求解。

在非獨(dú)立求解模式下，通常顆粒系被描述為粒徑頻度分布服從一定的函數(shù)形式，例如正態(tài)分布、對(duì)數(shù)正態(tài)分布或R-R分布，后者由式(3)給出：

為便于采用最優(yōu)化理論求解，往往依據(jù)超聲頻率、粒徑分布以及模型參數(shù)結(jié)合ECAH[21?22]模型預(yù)測(cè)理論聲衰減譜，同時(shí)將其和實(shí)測(cè)聲衰減譜比較，定義誤差函數(shù)如下：

其中，αsim為理論聲衰減預(yù)測(cè)譜，αmeas為實(shí)驗(yàn)衰減譜，通過(guò)誤差函數(shù)最小化，即可求解得最優(yōu)解。

2 改進(jìn)和聲搜索算法

2.1 標(biāo)準(zhǔn)和聲搜索算法

和聲搜索算法源于對(duì)音樂(lè)演奏的模擬，設(shè)樂(lè)器演奏的和聲i的狀態(tài)表示為Xi=(xi1,xi2,···,xiD)，美學(xué)評(píng)估設(shè)為Yi=f(Xi)，即待求目標(biāo)函數(shù)。算法流程如下：

(1)初始化和聲記憶庫(kù)HM：{xi(0)|xLi,j≤xi,j(0)≤xUi,j,i=1,2,···,HMS;j=1,2,···,D}，按照如下方式隨機(jī)產(chǎn)生：

其中，HMS 為記憶庫(kù)大小，xLi,j和xUi,j分別為xi,j取值的下限和上限，unifrnd為在上下限生成均勻分布的隨機(jī)數(shù)。

(2)基于差分變異創(chuàng)作一個(gè)新和聲：

其中，HMCR 為和聲記憶庫(kù)取值概率，PAR 為微調(diào)概率，bw為微調(diào)帶寬，r ∈{1,2,···,HMS}。

(3)更新和聲記憶庫(kù)HM：判斷新和聲目標(biāo)函數(shù)值f(Xnew)是否優(yōu)于當(dāng)前和聲庫(kù)中的f(Xworst)，若是，則用新和聲代替當(dāng)前和聲庫(kù)內(nèi)的最差和聲。

(4)是否滿足終止條件：若滿足停止條件，算法結(jié)束；否則重復(fù)步驟(2)和步驟(3)，直到滿足終止條件。

2.2 擬牛頓法

擬牛頓法是一種較成熟算法，擁有牛頓法收斂速度快的優(yōu)點(diǎn)，同時(shí)用每步梯度信息去近似Hesse 矩陣，減少計(jì)算量[23]，基本迭代公式為dk=?Hk?f(xk)，dk為優(yōu)化方向，由式(7)BFGS修正公式得Hk：

2.3 改進(jìn)和聲搜索算法

作為全局優(yōu)化算法，和聲搜索算法進(jìn)化到后期，記憶庫(kù)多樣性降低，貪婪選擇容易錯(cuò)過(guò)一些有潛力的解向量，陷入極小值，致求解精度不高；而擬牛頓BFGS 算法局部開發(fā)能力強(qiáng)，搜索結(jié)果更依賴于給定初值。由于搜索過(guò)程，局部和全局搜索未設(shè)也不易設(shè)清晰分界線，且前期與后期算法需求不同，為了提高反演準(zhǔn)確性，提出將和聲搜索算法和擬牛頓BFGS 進(jìn)行并行計(jì)算，HS 全局搜索能力強(qiáng)，將其近似解給BFGS 賦為初值利用其局部搜索能力，提高得到精確解的概率，兩者的融合有利于同步提高全局和局部效率。

圖1為改進(jìn)和聲算法(Improved harmony search algorithm,IHS)流程，在迭代前期僅使用和聲搜索算法，基于差分策略生成新和聲全局尋優(yōu)，待一定迭代次數(shù)時(shí)，擬牛頓BFGS算法介入，與和聲搜索算法并行計(jì)算。具體的，將和聲搜索算法尋得當(dāng)前最優(yōu)解賦給擬牛頓BFGS 作為初值，若BFGS 反演目標(biāo)函數(shù)值小于和聲搜索算法所得，則以BFGS 的最優(yōu)解替代和聲搜索種群庫(kù)里最優(yōu)解，更新記憶庫(kù)，直到滿足停止條件。

圖1 和聲搜索算法和BFGS 擬牛頓算法并行計(jì)算流程圖Fig.1 Harmony search and BFGS quasi-Newton algorithm parallel computing flow chart

3 仿真模擬

3.1 迭代次數(shù)

為使IHS 算法搜索解向量過(guò)程優(yōu)化，確定BFGS 介入時(shí)迭代次數(shù)非常關(guān)鍵。其需滿足兩個(gè)條件：反演的解向量誤差小且較為穩(wěn)定；反演所需時(shí)間不能過(guò)長(zhǎng)。為此，對(duì)多個(gè)算例進(jìn)行迭代次數(shù)100～900 間的反演，以50 次反演數(shù)據(jù)均值、方差與反演時(shí)長(zhǎng)加以比較。表1為R-R 分布(= 20 μm,K= 10)的反演結(jié)果，各組特征尺寸均值與設(shè)定值偏差為0.05%。

表1 不同迭代次數(shù)反演結(jié)果Table 1 Inversion results with different iterations

圖2給同一R-R分布不同迭代次數(shù)的反演結(jié)果方差和計(jì)算時(shí)長(zhǎng)，可看出當(dāng)?shù)螖?shù)為200 的時(shí)候，方差最小，50 次反演結(jié)果較為穩(wěn)定。反演時(shí)間隨BFGS 介入的迭代次數(shù)不同而不同，顯然BFGS 在迭代進(jìn)行到900次時(shí)再介入，所需時(shí)間最短，但反演結(jié)果波動(dòng)性相較于其他情況更大。介于9 種情況下所需時(shí)間相差不大，在一個(gè)量級(jí)上，迭代次數(shù)選擇應(yīng)先滿足結(jié)果的精確性與穩(wěn)定性，其次是所需時(shí)間，因此可選擇迭代次數(shù)200次時(shí)介入BFGS算法。

圖2 不同迭代次數(shù)的方差和時(shí)長(zhǎng)Fig.2 Variance and computation time required under different iterations

3.2 單峰分布

表2 三種分布函數(shù)下的單峰參數(shù)反演結(jié)果Table 2 Inversion results of three kinds of distribution functions by HS and IHS

進(jìn)一步分析，圖3比較了真值σ為3、正態(tài)分布不同(10～90 μm)粒徑反演，從圖3(a)可知，僅粒徑= 90 μm 時(shí)呈現(xiàn)隨機(jī)性波動(dòng)，但最大誤差小于5%，其他情況反演結(jié)果穩(wěn)定。對(duì)于σ，圖3(b)中為70 μm 時(shí)，σ值會(huì)有少許波動(dòng)，50次計(jì)算中6 次出現(xiàn)大幅度偏差，為90 μm時(shí)，增至13次。

圖3 不同粒徑下50 次反演變化Fig.3 The change of inversion value of 50 times under different particle size

3.3 雙峰分布

反演雙峰分布時(shí)，需求解參數(shù)增至5個(gè)(分別用以描述兩個(gè)分布函數(shù)及其配比)。表3數(shù)據(jù)為雙峰R-R 分布反演結(jié)果，HS 反演無(wú)法取得合理近似值，其求解與真值誤差極大，最大誤差可至170%。相比之下，IHS 反演峰值位置()與真值偏離較小，最大誤差為10%，分布參數(shù)K值亦有較明顯偏差。

表3 R-R 分布函數(shù)下的5 個(gè)參數(shù)反演值Table 3 Inversion results of R-R distribution functions by HS and IHS

3.4 抗噪性分析

由于實(shí)驗(yàn)超聲譜測(cè)量中往往會(huì)引入誤差(如實(shí)驗(yàn)噪聲或信號(hào)譜分析誤差)，為模擬真實(shí)環(huán)境中隨機(jī)噪聲或其他干擾，考核算法實(shí)用性，在計(jì)算衰減譜時(shí)人為加入了1%～5%的隨機(jī)誤差。

按照表4設(shè)定參數(shù)，將50 次分析結(jié)果和初始設(shè)定進(jìn)行對(duì)比，結(jié)果如圖4所示。當(dāng)設(shè)定R-R 分布= 70 μm、K= 10，在1%、3%隨機(jī)噪聲下反演值為69.83 μm 和69.77 μm，在95%置信水平下的置信區(qū)間分別為(69.825～69.838 μm)和(69.765～69.777 μm)，可見反演結(jié)果基本穩(wěn)定在噪聲水平內(nèi)，沒有放大。當(dāng)隨機(jī)噪聲繼續(xù)增至5%，尋得平均值為69.36 μm，誤差從0.2%增大到0.9%，粒徑參數(shù)保持穩(wěn)定。對(duì)于分布參數(shù)K，隨機(jī)信號(hào)1%、3%時(shí)誤差為2.15%、15.7%，隨機(jī)信號(hào)5% 時(shí)，誤差可增至19.66%，可見對(duì)R-R 分布而言，分布參數(shù)的反演難度更大。

圖4 不同噪聲水平下50 次反演變化Fig.4 The change of inverse value of 50 times under different SNR

表4 加入不同隨機(jī)噪聲IHS 反演結(jié)果Table 4 Inversion results of IHS with different random noise

由于理論上超聲有效帶寬越寬，所能獲得顆粒信息越充分，更利于提高反演精度。因此，針對(duì)同樣算例，仍在上述3 種噪聲情況下將頻譜范圍拓寬50%(1～15 MHz)加以分析。反演結(jié)果的誤差減小為0.11%、0.28%、0.36%；分布寬度K誤差為則依次2.1%、4.49%，7.36%?？梢娫黾宇l譜范圍，反演誤差變小，提升了抗噪性。同時(shí)在噪聲情況下，進(jìn)行了迭代次數(shù)為100～900的9 種BFGS算法介入情況反演分析，仍是迭代次數(shù)為200時(shí)，反演結(jié)果方差最小，最為穩(wěn)定?？紤]到實(shí)驗(yàn)室條件下超聲譜的測(cè)量噪聲水平可控制在1%～2%，結(jié)果表明IHS 對(duì)顆粒系的反演是有效的。

4 實(shí)測(cè)顆粒反演計(jì)算

4.1 實(shí)驗(yàn)原理和裝置

圖5為測(cè)量超聲反演信號(hào)的裝置簡(jiǎn)圖，測(cè)量區(qū)為6.5 mm。為避免溫度波動(dòng)影響， 實(shí)驗(yàn)中將樣品池置于循環(huán)低溫恒溫水槽內(nèi)，溫度設(shè)定為(20±0.5)?C。使用了兩塊等厚度的有機(jī)玻璃緩沖塊，安裝發(fā)射/接收信號(hào)超聲探頭(Olympus V310、中心頻率為5 MHz)，在激勵(lì)源(汕頭超聲，型號(hào)CTS-8077PR)作用下自發(fā)自收模式工作。雙通道數(shù)據(jù)采集卡(NI PCI-5133，100 MS/s)分別采集樣品池內(nèi)介質(zhì)為空氣或被測(cè)樣品時(shí)的緩沖塊-樣品界面反射信號(hào)A1g、A1s，以及樣品池內(nèi)介質(zhì)分別為水或被測(cè)樣品時(shí)的樣品-緩沖塊界面反射信號(hào)A2w、A2s，衰減系數(shù)表達(dá)式為

圖5 超聲測(cè)量原理示意圖Fig.5 Schematic diagram of ultrasonic measurement principle

其中，αw是水的衰減系數(shù)，Rw是緩沖塊與水之間反射系數(shù)，L為樣品池厚度。

4.2 樣品反演計(jì)算

實(shí)驗(yàn)樣品為兩種不同粒徑的硅-水懸濁液，顆粒密度為2.34 g/cm3，體積濃度均為10%，圖6為測(cè)得的實(shí)驗(yàn)超聲衰減譜。圖7是分別用HS、IHS 法對(duì)衰減譜的反演結(jié)果，均較好地區(qū)分出兩種顆粒。對(duì)于標(biāo)稱尺寸為(10±0.6)μm、(20±1.2)μm 硅顆粒由IHS 反演得體積中位徑(Dv50)分別為10.61 μm和21.64 μm，反演結(jié)果與標(biāo)稱值誤差分別為6.0%和8.2%，HS 算法反演結(jié)果則較明顯偏離真值。相同樣品的顯微鏡圖像分析Dv50分別為10.41 μm 和19.54 μm(歐美克公司PIP8.1 圖像分析儀，重復(fù)性誤差<3%)。超聲測(cè)量誤差略大于圖像分析，除反演誤差外，還可能與硅顆粒的非球形度、超聲譜測(cè)量誤差有關(guān)。

圖6 樣品實(shí)驗(yàn)衰減譜Fig.6 Ultrasonic attenuation spectra of silicon particles

圖7 兩種方法結(jié)果對(duì)比Fig.7 Comparison of the results by Image analysis and Ultrasonic spectrum

5 結(jié)論

本文主要對(duì)和聲搜索算法在超聲衰減譜法顆粒粒徑測(cè)量的反演問(wèn)題進(jìn)行了研究和討論，得到了以下結(jié)論：

IHS迭代次數(shù)為200次時(shí)，BFGS開始介入的改進(jìn)算法試反演優(yōu)化結(jié)果精確度和穩(wěn)定性最好，但反演所需時(shí)間較長(zhǎng)；通過(guò)多次模擬，IHS 算法反演結(jié)果誤差均小于10%，優(yōu)于HS 算法(極端誤差可達(dá)170%)，其反演精度還可結(jié)合更寬頻帶的譜得以提高；通過(guò)實(shí)驗(yàn)及與圖像法對(duì)比進(jìn)一步驗(yàn)證表明IHS具有較強(qiáng)的實(shí)際應(yīng)用性。本文工作表明IHS 算法相較于HS 算法在反演準(zhǔn)確性和抗噪性方面均得到了明顯改善，在超聲衰減譜法粒徑表征具有較好的實(shí)際應(yīng)用價(jià)值。