蔣 瑜， 賈 楠， 蘇明旭

（上海理工大學(xué) 能源與動(dòng)力工程學(xué)院，上海 200093）

近年來隨著與顆粒緊密關(guān)聯(lián)的行業(yè)如能源、化工、食品、環(huán)境、醫(yī)藥等的飛速發(fā)展，顆粒兩相流測量技術(shù)獲得了廣泛的關(guān)注[1?5]。超聲顆粒測量技術(shù)具有穿透力強(qiáng)、非浸入式、易實(shí)現(xiàn)在線檢測、可測粒徑范圍寬等特點(diǎn)，其中，超聲衰減譜法[6?7]是一種具有廣闊發(fā)展前景的在線顆粒測量方法，超聲波在顆粒兩相體系中傳播并與顆粒之間發(fā)生物理作用，得到超聲衰減譜中包含的大量顆粒粒徑信息，結(jié)合BLBL（Bouguer-Lambert-Beer-Law）、ECAH（Epstein-Carhart-Allegra-Hawley）模型[8?9]，進(jìn)一步反演計(jì)算，可以快速分析顆粒粒徑分布。

顆粒粒徑測量的反演問題得到了廣泛關(guān)注。欒志超[10]利用擬牛頓算法模擬反演粒徑分布。蘇明旭等[11]曾通過非線性最小二乘優(yōu)化LM 算法對(duì)目標(biāo)函數(shù)進(jìn)行反演計(jì)算?？酌鞯萚12]用自動(dòng)基線擬合的累計(jì)量算法對(duì)光子相關(guān)光譜函數(shù)進(jìn)行反演。汪雪等[13]采用改進(jìn)遺傳算法（genetic algorithm）反演顆粒分布參數(shù)。與上述方法相比，差分進(jìn)化算法（differential evolution, DE）[14]具有很好的全局搜索能力和魯棒性。與標(biāo)準(zhǔn)遺傳算法相比，差分進(jìn)化算法所需控制參數(shù)較少，原理簡單易行，已在許多優(yōu)化問題中獲得了良好效果[15?17]。

針對(duì)超聲法顆粒測量，本文研究了一種非獨(dú)立模式條件下的超聲衰減譜反演差分進(jìn)化改進(jìn)算法，針對(duì)標(biāo)準(zhǔn)差分進(jìn)化算法容易陷入局部最優(yōu)以及早熟收斂問題，提出了自適應(yīng)控制參數(shù)因子并進(jìn)行優(yōu)化，通過對(duì)高斯分布、R-R（Rosin-Rammler）分布以及對(duì)數(shù)正態(tài)分布函數(shù)的反演過程模擬，將改進(jìn)差分進(jìn)化算法和其他算法比較，分析了改進(jìn)差分進(jìn)化算法的特點(diǎn)，尤其是不同設(shè)定峰值數(shù)情況下噪聲水平影響及抑制。

1 超聲衰減譜法及反演優(yōu)化問題

當(dāng)超聲波通過含顆粒物介質(zhì)時(shí)，由于顆粒與聲波的作用而產(chǎn)生衰減，聲波在兩相之間的分界面處，會(huì)在固體顆粒的內(nèi)部和外部產(chǎn)生壓縮波，同時(shí)在內(nèi)部產(chǎn)生熱波、剪切波。在球坐標(biāo)下，運(yùn)用邊界條件得到在顆粒內(nèi)部和周圍介質(zhì)的3 種波動(dòng)的勢函數(shù)，其線性方程組可簡寫為[18]

式中：M為6 階方陣；An，Bn，Cn分別表示壓縮波、熱波以及剪切波；e為散射系數(shù)向量。

求解方程組（1）可得到系數(shù)An，聲波總衰減由系數(shù)An決定，最終的衰減系數(shù)

式中：kc為連續(xù)介質(zhì)中的波數(shù)；φ為顆粒相體積濃度；R為顆粒半徑。

式中：α為聲衰減分布列向量；S為衰減系數(shù)矩陣；W為被測顆粒的粒徑分布。

為了便于采用最優(yōu)化理論求解，往往根據(jù)超聲頻率、粒徑分布以及模型參數(shù)結(jié)合ECAH 模型計(jì)算理論聲衰減譜，同時(shí)將其同實(shí)驗(yàn)聲衰減譜比較，定義誤差函數(shù)

式中：αs為理論聲衰減預(yù)測譜；αm為實(shí)驗(yàn)聲衰減譜。

通過反演算法最小化誤差函數(shù)，即可求解得最優(yōu)的顆粒粒徑分布。對(duì)于最優(yōu)化求解，通常顆粒系被描述為粒徑頻度分布服從一定的函數(shù)形式，例如，高斯分布、對(duì)數(shù)正態(tài)分布或R-R（Rosin-Rammler）分布。其中，R-R 分布函數(shù)為

對(duì)不同的優(yōu)化算法適用性進(jìn)行研究，設(shè)計(jì)算例，對(duì)半徑為5 μm 的單峰R-R 分布玻璃微珠顆粒和設(shè)定參數(shù)雙峰分布系K2=25，W1=0.5）進(jìn)行數(shù)值反演。W1為第一個(gè)峰所占的比例。按式（6）形式構(gòu)造目標(biāo)函數(shù)，分別采用DFP（Davidon-Fletcher-Powell），LM（Levenberg-Marquard）、BFGS（Broyden-Fletcher-Goldfarb-Shanno）和基本差分進(jìn)化算法（DE）在0.1～50 μm粒徑區(qū)間尋優(yōu)。從圖1 中可以看出，無論是給定單峰還是雙峰分布，3 種局部最優(yōu)化算法DFP，LM，BFGS的結(jié)果均明顯偏離真值，而DE 算法的結(jié)果總體上符合對(duì)真值范圍的預(yù)期，雙峰反演結(jié)果存有一定偏差，但明顯優(yōu)于DFP，LM，BFGS 算法，由此可見，基于DE 原理的反演算法完全可行。V（D）為體積頻度分布。

圖 1 4 種算法反演計(jì)算的顆粒粒徑分布Fig.1 Inversion results by four algorithms for the given particle size distribution

2 差分進(jìn)化算法

差分進(jìn)化算法是一種基于群體智能理論的優(yōu)化算法，通過群體內(nèi)個(gè)體間的合作與競爭來對(duì)目標(biāo)進(jìn)行優(yōu)化。近年來差分進(jìn)化算法在最優(yōu)化問題求解方面得到了廣泛的應(yīng)用。

2.1 標(biāo)準(zhǔn)差分進(jìn)化算法

基本的差分進(jìn)化算法主要由4 個(gè)步驟組成。

式中：N為種群大??；分別為xj取值的下限和上限，rand（0,1）為在（0,1）區(qū)間均勻分布的隨機(jī)數(shù)。

b. 變異操作。通過如下差分策略，產(chǎn)生變異中間體

式中：F為縮放因子，本文取F=0.35。

c. 交叉操作。對(duì)第g代種群{xi（g）}及其變異的中間體{vi（g+1）}進(jìn)行個(gè)體間的交叉操作

式中：CR為交叉因子，本文取CR=0.85；jrand為隨機(jī)取出解向量的位置。

d. 選擇操作。差分進(jìn)化算法采用貪婪算法來選擇進(jìn)入下一代種群的個(gè)體。

式中，F(xiàn)為目標(biāo)函數(shù)值。

反復(fù)執(zhí)行式（9）～（12），直到滿足終止條件為止。

2.2 算法改進(jìn)

DE 算法主要有3 個(gè)控制參數(shù)，即種群大小N、縮放因子F和交叉因子CR。縮放因子F控制偏差變量的放大作用，改變搜索的方向，若種群過早收斂，則F應(yīng)增加，反之亦然。交叉因子CR改變種群多樣性，較大的CR會(huì)加速收斂。標(biāo)準(zhǔn)DE 算法中的F，CR為固定值，在算法迭代后期會(huì)因其種群個(gè)體聚集，造成種群間差異性變小，使算法易陷入局部最優(yōu)解、出現(xiàn)早熟收斂等問題。因此，在算法中引入自適應(yīng)控制參數(shù)因子加以改進(jìn)，改進(jìn)后的縮放因子和交叉因子為

式中：τ1，τ2為調(diào)節(jié)因子；FL，F(xiàn)U分別為縮放因子的下限和上限。

圖2 給出了不同的FL，F(xiàn)U對(duì)應(yīng)的目標(biāo)函數(shù)值，為了提高優(yōu)化效率，給定FL（0.1,0.4），F(xiàn)U（0.5,0.9）的范圍，在多種粒徑分布參數(shù)下進(jìn)行反演尋優(yōu)，選取目標(biāo)函數(shù)值最小處對(duì)應(yīng)的FL，F(xiàn)U。圖2（b）為FL，F(xiàn)U對(duì)應(yīng)的目標(biāo)函數(shù)值等高線圖，由圖2（b）可確定FL，F(xiàn)U的可選范圍，對(duì)應(yīng)目標(biāo)函數(shù)值最小，將FL，F(xiàn)U分別選取為0.3，0.8。

將改進(jìn)后的算法稱為IDE（improved differential evolution）算法，算法流程如圖3 所示。

圖 2 適應(yīng)值隨FU，F(xiàn)L 的變化Fig.2 Variation of the best fitness value along with FU and FL

圖 3 IDE 算法流程Fig. 3 Flowchart of the improved differential evolutionalgorithm

a. 確定DE 算法參數(shù)，隨機(jī)產(chǎn)生初始種群并計(jì)算個(gè)體的適應(yīng)度（fitness）。

b. 判斷是否最佳適應(yīng)度基本不變。若是，則輸出最佳個(gè)體為最優(yōu)解；若否，結(jié)合自適應(yīng)控制參數(shù)產(chǎn)生縮放因子F、交叉因子CR，進(jìn)行變異和交叉操作，得到中間種群，在原種群和中間種群選擇個(gè)體，得到新種群。

c. 若滿足停止條件，算法結(jié)束；反之，進(jìn)化代數(shù)t=t+1，轉(zhuǎn)步驟b。

3 結(jié)果與分析

3.1 數(shù)值模擬

為了驗(yàn)證改進(jìn)算法的性能，分別采用標(biāo)準(zhǔn)DE 算 法 和 改 進(jìn) 算 法IDE（ improved differential evolution）對(duì)滿足單峰、雙峰分布的顆粒群進(jìn)行數(shù)值模擬。

3.1.1 單峰分布

按照表1 設(shè)定參數(shù)，對(duì)服從R-R、正態(tài)和對(duì)數(shù)正態(tài)分布的顆粒系進(jìn)行反演計(jì)算。R，K（或σ）分別為特征尺寸以及分布特征參數(shù)?？梢钥闯觯瑯?biāo)準(zhǔn)DE 算法和改進(jìn)算法反演得到的結(jié)果均與原始分布相同，改進(jìn)算法IDE 計(jì)算一次需要時(shí)間10 s 左右，比DE 算法多2 s 左右。

表 1 3 種分布函數(shù)下的單峰參數(shù)反演結(jié)果Tab.1 Inversion results of three kinds of distribution functions by DE and IDE

3.1.2 雙峰分布

從表2 數(shù)據(jù)還可以看出，在上述算例將頻譜范圍從之前1～10 MHz 增至1～15 MHz 時(shí)，反演結(jié)果和設(shè)定值吻合。增加頻寬的范圍對(duì)于反演更有利，但頻寬增加對(duì)于實(shí)驗(yàn)技術(shù)提出了更高要求。

根據(jù)反演結(jié)果求得的顆粒特征尺寸如表3 所示，Dv10，Dv50，Dv90分別表示顆粒體積累積概率達(dá)到10%，50%，90%時(shí)對(duì)應(yīng)的顆粒粒徑，可以看出，改進(jìn)差分進(jìn)化算法反演得到的特征粒徑更接近設(shè)定分布值，其相比于設(shè)定分布的誤差均小于5%，而原始反演算法得到的誤差甚至可以達(dá)到30%。

對(duì)雙峰分布的反演結(jié)果如圖4 所示，可以看出，IDE 反演結(jié)果均接近原始分布，而DE 算法反演分布曲線相比原始分布均出現(xiàn)一定偏離。如圖4（a）所示，其第一個(gè)峰的位置向右產(chǎn)生了偏移，同時(shí)兩個(gè)峰的峰值均小于原始設(shè)定值。同步獲得目標(biāo)函數(shù)反演殘差，原始DE 算法譜誤差ESSD=0.011，IDE 算法反演得到的ESSD=1.09×10?14。圖4（b）中第二個(gè)峰的結(jié)果與設(shè)定基本吻合，但第一個(gè)峰的反演粒徑偏小且峰寬偏窄，對(duì)應(yīng)原始DE 算法反演得到的ESSD=0.0027，IDE 算法反演得到的ESSD=1.04×10?12，可見IDE 求解精度得到明顯提高。

表 2 R-R 分布函數(shù)下的5 個(gè)參數(shù)反演值Tab.2 Inversion results of R-R distribution functions by DE and IDE

表 3 R-R 分布函數(shù)下反演得到的顆粒特征尺寸Tab.3 Inversion results of R-R distribution functions

圖 4 雙峰顆粒粒徑分布反演結(jié)果Fig.4 Inversion results of bimodal particle size distribution

3.2 抗噪性分析

在超聲衰減譜的實(shí)際測量過程中會(huì)不可避免地引入測量誤差（如實(shí)驗(yàn)過程噪聲或譜分析誤差），為了模擬實(shí)際測量中的隨機(jī)噪聲和其他干擾，在超聲衰減譜中人為加入3%，5%，10%的隨機(jī)噪聲，進(jìn)行反演計(jì)算。為了討論算法的隨機(jī)性影響和抑制，分析了25 組連續(xù)反演結(jié)果，并將其均值作為最終結(jié)果與設(shè)定顆粒粒徑分布和無噪聲情形進(jìn)行對(duì)比。

不同噪聲幅度下單峰R-R 分布反演結(jié)果如圖6（a）所示，可以看出，采用改進(jìn)差分進(jìn)化算法的反演曲線光滑平順，未加入噪聲時(shí)反演分布和設(shè)定曲線重合；超聲譜加入3%，5%，10%的噪聲后，反演曲線稍有偏移，但其值偏差始終在5%以內(nèi)，且基本保持了原有分布寬度。R-R 雙峰分布顆粒系的超聲譜加入1%，2%噪聲后的反演結(jié)果如圖6（b），其反演出兩個(gè)峰的分布寬度相對(duì)偏窄，但偏差均在5%及以內(nèi)；當(dāng)加入3%噪聲后，偏差變大，第一個(gè)峰分布過寬，第二個(gè)峰也向左偏移且分布寬度變窄，其偏差在15%及以內(nèi)。考慮到實(shí)驗(yàn)室條件下超聲譜的測量隨機(jī)誤差為1%～2%，上述結(jié)果表明，IDE 對(duì)于單峰和雙峰分布顆粒系的反演均是有效的，對(duì)反演結(jié)果求平均值（本文采用了25 次），可以較好地抑制此類噪聲的影響，相對(duì)誤差小于5%。

圖 5 在不同噪聲幅度下K，值的變化Fig.5 Changes of the values of and K with different noise magnitudes

圖 6 R-R 分布函數(shù)加入噪聲的反演計(jì)算結(jié)果Fig.6 Inversion results of R-R distribution functions with noise

4 結(jié) 論

超聲衰減譜法粒徑反演時(shí)，通過對(duì)差分進(jìn)化算法優(yōu)化，自適應(yīng)控制縮放因子F、交叉因子CR，由本文數(shù)值模擬結(jié)果可知，改進(jìn)差分進(jìn)化算法反演得出的分布參數(shù)值，K誤差小于5%，而標(biāo)準(zhǔn)差分進(jìn)化算法反演得到的，K值誤差達(dá)到20%及以上。改進(jìn)差分進(jìn)化算法得到的Dv10，Dv50，Dv90與設(shè)定分布的相對(duì)誤差均小于5%，而原始差分進(jìn)化算法得到的誤差甚至可以達(dá)到30%。反演結(jié)果顯示，改進(jìn)差分進(jìn)化算法具有很強(qiáng)的抗干擾能力，反演峰值位置偏差很小，粒度相對(duì)誤差不超過5%。