基于改進差分進化算法的超聲衰減譜反演計算

蔣 瑜， 賈 楠， 蘇明旭

（上海理工大學 能源與動力工程學院，上海 200093）

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

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

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

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

當超聲波通過含顆粒物介質(zhì)時，由于顆粒與聲波的作用而產(chǎn)生衰減，聲波在兩相之間的分界面處，會在固體顆粒的內(nèi)部和外部產(chǎn)生壓縮波，同時在內(nèi)部產(chǎn)生熱波、剪切波。在球坐標下，運用邊界條件得到在顆粒內(nèi)部和周圍介質(zhì)的3 種波動的勢函數(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 模型計算理論聲衰減譜，同時將其同實驗聲衰減譜比較，定義誤差函數(shù)

式中：αs為理論聲衰減預測譜；αm為實驗聲衰減譜。

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

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

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

2 差分進化算法

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

2.1 標準差分進化算法

基本的差分進化算法主要由4 個步驟組成。

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

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

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

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

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

d. 選擇操作。差分進化算法采用貪婪算法來選擇進入下一代種群的個體。

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

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

2.2 算法改進

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

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

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

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

圖 2 適應值隨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ù)，隨機產(chǎn)生初始種群并計算個體的適應度（fitness）。

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

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

3 結(jié)果與分析

3.1 數(shù)值模擬

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

3.1.1 單峰分布

按照表1 設定參數(shù)，對服從R-R、正態(tài)和對數(shù)正態(tài)分布的顆粒系進行反演計算。R，K（或σ）分別為特征尺寸以及分布特征參數(shù)?？梢钥闯觯瑯藴蔇E 算法和改進算法反演得到的結(jié)果均與原始分布相同，改進算法IDE 計算一次需要時間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 時，反演結(jié)果和設定值吻合。增加頻寬的范圍對于反演更有利，但頻寬增加對于實驗技術提出了更高要求。

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

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

表 2 R-R 分布函數(shù)下的5 個參數(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 抗噪性分析

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

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

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

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

4 結(jié) 論

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