岳樸杰, 張世瑋, 孟 磊, 寧 翔, 俞天陽

(1. 大唐環(huán)境產(chǎn)業(yè)集團(tuán)股份有限公司, 北京 100097; 2. 上海理工大學(xué) 上海市動(dòng)力工程多相流動(dòng)與傳熱重點(diǎn)實(shí)驗(yàn)室, 上海 200093)

0 引 言

近年來, 兩相流中顆粒粒徑測量問題受到國內(nèi)外廣泛關(guān)注和研究, 如能源、 化工、 制藥和環(huán)保等行業(yè)中各種懸濁液和乳劑的顆粒在線測量問題[1-3]。 其中顆粒的形狀、 粒徑分布及存在狀態(tài)對(duì)工業(yè)產(chǎn)品的性能質(zhì)量、 生產(chǎn)效率和安全均有重要影響, 因此, 對(duì)兩相介質(zhì)中顆粒物特征參數(shù)的測量顯得尤為重要[4-5]。 在眾多顆粒尺寸測量技術(shù)中, 超聲波具有較強(qiáng)的穿透性、 測量范圍寬、 適合高濃度體系的非接觸式測量等特點(diǎn), 同時(shí)有測量裝置簡單、 易實(shí)現(xiàn)在線檢測等優(yōu)勢[6-8]。 超聲波在兩相體系中傳播與顆粒發(fā)生相互作用并產(chǎn)生衰減, 得到的超聲衰減譜包含被測顆粒系的粒徑信息, 結(jié)合ECAH(Epstein-Carhart-Allegra-Hawley), McC(McClements), BLBL(Bouguer-Lambert-Beer-Law)等模型預(yù)測, 構(gòu)造目標(biāo)函數(shù)進(jìn)一步經(jīng)反演計(jì)算可得到顆粒系的粒徑分布[9-11]。

實(shí)驗(yàn)衰減譜的測量、 理論模型的建立和反演計(jì)算是超聲衰減譜法顆粒測量的3個(gè)主要過程, 其中, 反演算法直接關(guān)系到顆粒粒徑分布求解的精度。 反演方法可以分為非獨(dú)立模式算法和獨(dú)立模式算法。 非獨(dú)立模式算法通常需事先假定顆粒系的粒徑分布遵循某種函數(shù)形式, 然后以理論預(yù)測和實(shí)際測量的衰減譜誤差作為目標(biāo)函數(shù)來構(gòu)造最優(yōu)化問題[12-14]。 但大多數(shù)情況下并不知道顆粒系粒徑分布情況, 通常采用無函數(shù)約束的獨(dú)立模式求解方法, 利用正則化技術(shù)克服方程的病態(tài)問題, 引入正則化因子和光順矩陣, 采用非負(fù)最小二乘算法進(jìn)行超聲譜的反演得到顆粒系粒徑分布[15-16]。

在使用超聲系統(tǒng)測量懸濁液或乳劑的實(shí)驗(yàn)衰減譜過程中, 不可避免地會(huì)帶來測量誤差, 但非負(fù)最小二乘算法并未考慮該誤差對(duì)求解粒徑分布的影響。 文章提出一種基于最優(yōu)正則化的加權(quán)非負(fù)最小二乘OrtWtls(Optimum regularization technique Weighted total least squares)算法, 反演顆粒系的粒徑分布, 以不同粒徑和信噪比的理論衰減譜對(duì)該算法進(jìn)行驗(yàn)證。 此外, 將此算法應(yīng)用于不同濃度石灰石漿液粒度分布的反演計(jì)算, 并采用圖像法對(duì)其結(jié)果進(jìn)行驗(yàn)證。

1 加權(quán)最小二乘粒徑反演

1.1 獨(dú)立模式算法求解粒徑分布

顆粒粒徑分布測量原理是基于不同粒徑產(chǎn)生的超聲譜特性差異, 通過尋找理論和試驗(yàn)數(shù)據(jù)之間的最佳擬合值確定。 借助波的散射理論聯(lián)系單個(gè)顆粒的散射效應(yīng)和整個(gè)離散顆粒群的宏觀效果, 得到懸濁液中的復(fù)波數(shù)[17], 進(jìn)一步求解得到的超聲衰減系數(shù)

(1)

式中:αs表示超聲衰減系數(shù);An為散射系數(shù);φ為體積濃度;kc為連續(xù)介質(zhì)波數(shù);R為顆粒半徑;ω為角頻率;qi為顆粒半徑在第i個(gè)尺寸區(qū)間[Ri,Ri+1]的體積分?jǐn)?shù);n為顆粒粒徑分散等級(jí); Re為取實(shí)部運(yùn)算。

為獲得準(zhǔn)確的粒徑分布, 將模型改寫為標(biāo)準(zhǔn)數(shù)學(xué)求解問題, 改變系數(shù)矩陣和超聲譜向量, 將式(1)離散并在多頻條件下轉(zhuǎn)化為線性方程形式[18]

AF=G,

(2)

式中:A是系數(shù)矩陣;F是待求的離散化顆粒尺寸體積分布;G是由不同頻率下的實(shí)驗(yàn)超聲譜組成的向量。

該矩陣方程非適定且高度病態(tài), 常用矩陣求逆無法得到準(zhǔn)確結(jié)果, 通常由最優(yōu)正則化(ORT, Optimum Regularization Technique)算法結(jié)合非負(fù)最小二乘求解[19-20], 通過引入正則化因子γ和光順矩陣H, 對(duì)病態(tài)方程進(jìn)行改善, 提高解的穩(wěn)定性。 表達(dá)式為

F=(ATA+γH)-1ATG,

(3)

式中:AT表示A的轉(zhuǎn)置矩陣;γ稱正則化因子, 表示矩陣的權(quán)重。

1.2 加權(quán)最小二乘求解粒徑分布

為減小兩相介質(zhì)中超聲衰減譜的測量誤差對(duì)粒徑分布求解精度的影響, 提出一種最優(yōu)正則化算法結(jié)合加權(quán)最小二乘, 求解顆粒系粒徑分布。 引入權(quán)重矩陣W, 粒徑分布求解公式表達(dá)式為

F=(ATWA+γH)-1ATWG。

(4)

在OrtWtls算法中, 正則化參數(shù)的選擇較為關(guān)鍵, 該參數(shù)的確定取決于系數(shù)矩陣A, 測量向量G的誤差程度以及顆粒粒徑分布的光滑程度。 鑒于在顆粒測量問題中很難估計(jì)測量向量G的誤差水平, 因此常使用廣義交叉檢驗(yàn)(Generalized Cross Validation, GCV)準(zhǔn)則和L曲線準(zhǔn)則(L-Curve)對(duì)正則化因子進(jìn)行優(yōu)化[21]。 其中GCV準(zhǔn)則考慮測量向量中各元素都應(yīng)被相應(yīng)的正則解預(yù)測, 且正則化參數(shù)與G的正交變換無關(guān)。 經(jīng)推導(dǎo),γ可歸結(jié)為對(duì)式(5)求最小值來確定[22]。

(5)

式中:I為單位矩陣;m為所用頻率數(shù)目, 也是單位矩陣的階數(shù), 上標(biāo)T為矩陣轉(zhuǎn)置運(yùn)算, 符號(hào)Trace是求矩陣的跡。

式(6)為超聲衰減譜權(quán)重矩陣W的計(jì)算公式

(6)

式中:ej為超聲衰減譜與其擬合譜的殘差絕對(duì)值;σ為ej的標(biāo)準(zhǔn)差。 根據(jù)式(6)將超聲衰減譜的殘差分為3類:① 殘差不是粗大誤差, 殘差ej小于k0σ, 權(quán)重因子為1;② 殘差可能是粗大誤差, 此時(shí)殘差ej在k0σ和k1σ之間, 權(quán)重因子為k0/(|ej/σ|);③ 殘差是粗大誤差, 此時(shí)殘差ej大于k1σ, 權(quán)重因子為0。 通常取k0為1.5,k1為2.5[23-24]。

圖1 為OrtWtls算法流程圖, 基于ORT算法結(jié)合加權(quán)最小二乘反演顆粒系粒徑分布。

圖1 OrtWtls算法流程圖Fig.1 Flow chart of OrtWtls algorithm

由權(quán)重系數(shù)dw構(gòu)建權(quán)重矩陣W, 并根據(jù)權(quán)重系數(shù)對(duì)實(shí)驗(yàn)超聲衰減譜進(jìn)行修正, 直至權(quán)重系數(shù)dw中無粗大誤差或殘差的標(biāo)準(zhǔn)差小于0.01結(jié)束迭代。 與最小二乘算法相比, 加權(quán)最小二乘使超聲衰減譜中具有粗大誤差的權(quán)重系數(shù)置零, 能夠有效降低實(shí)驗(yàn)衰減譜測量誤差對(duì)反演結(jié)果的影響。

2 數(shù)值模擬

2.1 OrtWtls算法驗(yàn)證

為驗(yàn)證OrtWtls算法的性能, 以體積濃度1%的水-玻璃微珠懸濁液為研究對(duì)象, 采用McC+BLBL理論模型分別計(jì)算顆粒直徑10 μm～80 μm懸濁液中的理論超聲衰減譜。 使用GCV準(zhǔn)則對(duì)γ參數(shù)進(jìn)行優(yōu)化,V(γ)函數(shù)取最小值時(shí)對(duì)應(yīng)的γ即為最優(yōu)正則化因子, 不同粒徑下的計(jì)算結(jié)果如表1 所示。 經(jīng)分析計(jì)算可知, 對(duì)上述直徑10 μm～80 μm的理論衰減譜進(jìn)行三階多項(xiàng)式擬合處理, 其相關(guān)系數(shù)均值為0.999 8, 因此, 不同粒徑下的超聲衰減譜可用三階多項(xiàng)式近似表達(dá)。

表1 不同粒徑的正則化因子Tab.1 Regularization factors of different particle sizes

基于上述分析, 使用OrtWtls反演顆粒系的粒徑分布, 計(jì)算不同算例粒徑分布的體積中位徑DV50(小于該直徑的顆粒體積各占顆粒總體積的50%), 并與設(shè)定直徑進(jìn)行對(duì)比, 如圖2 所示。 反演結(jié)果相對(duì)誤差均小于5%, 驗(yàn)證了OrtWtls算法的計(jì)算準(zhǔn)確性。

圖2 設(shè)定直徑與OrtWtls反演結(jié)果對(duì)比Fig.2 Comparison between set diameter and inversion result of OrtWtls

2.2 抗噪性分析

在實(shí)驗(yàn)超聲衰減譜獲取過程中, 不可避免會(huì)引入測量誤差, 為模擬實(shí)驗(yàn)過程中帶來的隨機(jī)干擾和測量誤差, 在理論衰減譜中加入隨機(jī)噪聲, 使其信噪比分別為10 dB, 20 dB, 使用OrtWtls算法反演顆粒系的粒徑分布。 為分析提出算法的魯棒性, 統(tǒng)計(jì)連續(xù)20組含噪衰減譜的反演結(jié)果, 并與預(yù)設(shè)直徑和Twomey算法[25]計(jì)算結(jié)果進(jìn)行對(duì)比。 圖3 為信噪比10 dB時(shí)不同粒徑的理論和含噪衰減譜, 圖4 為含噪衰減譜分別由Twomey和OrtWtls算法計(jì)算得到的顆粒系粒徑分布。

圖3 信噪比10 dB超聲衰減譜Fig.3 Ultrasonic attenuation spectrum of signal noise ratio 10 dB

圖4 不同算法反演結(jié)果對(duì)比Fig.4 Comparison of calculation results of different algorithms

由圖3 中超聲衰減譜可知, 含噪衰減系數(shù)隨超聲頻率的變化趨勢與理論衰減譜基本一致, 但個(gè)別頻率下衰減系數(shù)的殘差為粗大誤差, 致使圖4 中Twomey算法得到的DV50與預(yù)設(shè)直徑差異較大, 顆粒直徑10 μm, 20 μm, 30 μm, 40 μm時(shí)的相對(duì)誤差分別為55.3%, 16.4%, 23.8%, 14.6%。 而對(duì)超聲衰減譜做加權(quán)處理的OrtWtls算法得到的特征參數(shù)DV50與設(shè)定直徑基本一致, 不同粒徑下的相對(duì)誤差均小于5%。 表2 為信噪比10 dB超聲衰減譜連續(xù)20次計(jì)算粒徑分布DV50的統(tǒng)計(jì)結(jié)果。

表2 10 dB衰減譜粒徑分布統(tǒng)計(jì)結(jié)果Tab.2 Statistical results of particle size distribution of 10 dB attenuation spectrum

圖5(a) 為信噪比10 dB衰減譜的粒徑分布統(tǒng)計(jì)結(jié)果。 相比Twomey算法, 提出的OrtWtls算法得到的特征參數(shù)DV50與設(shè)定直徑更接近, 不同粒徑下相對(duì)誤差分別為4.0%, 2.1%, 1.4%, 1.2%。 為分析OrtWtls算法在不同信噪比下的抗噪性, 表3 給出了信噪比20 dB衰減譜連續(xù)20次計(jì)算粒徑分布DV50的均值、 方差和標(biāo)準(zhǔn)差, 圖5(b) 為該信噪比下的粒徑分布統(tǒng)計(jì)結(jié)果。

表3 20 dB衰減譜粒徑分布統(tǒng)計(jì)結(jié)果Tab.3 Statistical results of particle size distribution of 20 dB attenuation spectrum

(a) 信噪比10 dB衰減譜結(jié)果統(tǒng)計(jì)

由圖5 可知, 對(duì)于信噪比10 dB和20 dB的衰減譜, 不同粒徑時(shí)OrtWtls算法得到的DV50與設(shè)定直徑誤差均小于5%, 且多次計(jì)算的標(biāo)準(zhǔn)差均小于1 μm。 設(shè)定顆粒直徑10 μm時(shí), OrtWtls算法得到的特征參數(shù)DV50與設(shè)定值基本一致, 且標(biāo)準(zhǔn)差為0.43 μm, 而Twomey算法的計(jì)算結(jié)果偏差較大。 當(dāng)顆粒直徑大于20 μm時(shí), 兩種算法得到的DV50與設(shè)定值基本一致, 但OrtWtls算法標(biāo)準(zhǔn)差明顯小于Twomey算法。 提出的OrtWtls算法對(duì)于含噪超聲衰減譜的反演具有較好的魯棒性, 能夠準(zhǔn)確計(jì)算顆粒系的粒徑分布。

3 實(shí)測顆粒反演計(jì)算

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

圖6 為基于超聲法原理測量漿液粒度的實(shí)驗(yàn)系統(tǒng)。 樣品池寬度12.5 mm, 在厚度10 mm的有機(jī)玻璃緩沖塊兩端分別安裝超聲發(fā)射、 接收換能器(Olympus V310, 探頭頻率5 MHz), 使用CTS-8077PR信號(hào)發(fā)射接收器激勵(lì)超聲換能器產(chǎn)生超聲波, 另一換能器接收經(jīng)被測漿液后的超聲透射信號(hào), 經(jīng)信號(hào)放大后選用NI-PCI-5114高速A/D卡(采樣率250 MHz)完成對(duì)超聲信號(hào)的采集。 對(duì)超聲信號(hào)進(jìn)行頻譜分析, 得到被測漿液的衰減譜, 最終經(jīng)數(shù)值反演得到被測漿液的粒徑分布。

圖6 超聲測量系統(tǒng)示意圖Fig.6 Schematic diagram of ultrasonic measurement system

以石灰石粉末分別配置質(zhì)量濃度為10%, 20%的漿液, 使用BT-1600光學(xué)顯微鏡觀察漿液顆粒的大小和形態(tài), 圖7 為經(jīng)稀釋后的石灰石漿液光學(xué)顯微鏡圖片。 使用ImageJ圖像分析軟件對(duì)相同試樣的10張顯微鏡圖片進(jìn)行分析, 統(tǒng)計(jì)得到石灰石漿液的粒徑分布。

圖7 石灰石漿液顯微鏡圖片F(xiàn)ig.7 Microscope image of limestone slurry

3.2 實(shí)驗(yàn)數(shù)據(jù)與分析

超聲衰減系數(shù)可以通過超聲波穿過漿液和背景介質(zhì)的信號(hào)幅值差異來計(jì)算, 采用透射法測量漿液衰減系數(shù)的表達(dá)式為[26]

(7)

式中:Aw為背景介質(zhì)的透射信號(hào)幅值;As為被測漿液的透射信號(hào)幅值,l2為測量區(qū)的寬度。 多個(gè)超聲頻率的衰減系數(shù)即可組成實(shí)驗(yàn)超聲衰減譜。

使用圖6 測量系統(tǒng), 分別采集質(zhì)量濃度10%, 20%漿液超聲透射信號(hào)并計(jì)算超聲衰減譜, 結(jié)果如圖8 所示。 超聲衰減系數(shù)均隨超聲頻率的增加而增大, 近似呈線性變化趨勢。 另外, 由于漿液中單位體積內(nèi)的顆粒數(shù)目隨漿液濃度的增加而增大, 對(duì)超聲波能量的耗散增強(qiáng), 故質(zhì)量濃度20%的衰減系數(shù)較10%時(shí)的大。 根據(jù)不同濃度漿液的實(shí)驗(yàn)超聲衰減譜, 分別使用OrtWtls和Twomey算法求解漿液的粒徑分布, 并與圖像法的統(tǒng)計(jì)結(jié)果做對(duì)比, 結(jié)果如圖9 所示。

圖8 石灰石漿液超聲衰減譜Fig.8 Ultrasonic attenuation spectrum of limestone grout

圖9 超聲法與圖像法粒徑分布對(duì)比Fig.9 Comparison of particle size distribution between ultrasonic and image method

由圖9 中的粒徑分布可知, 漿液質(zhì)量濃度為10%和20%時(shí), OrtWtls算法得到的粒徑分布DV50分別為9.50 μm, 9.67 μm, 與圖像法體積中位徑8.82 μm的測量結(jié)果基本吻合, 且不同濃度下的相對(duì)誤差分別為7.71%和9.64%。 但Twomey算法并未得到確切的粒徑分布, 這是因?yàn)樵撍惴]有對(duì)超聲衰減譜的測量誤差進(jìn)行處理。 對(duì)于實(shí)驗(yàn)得到的含噪超聲衰減譜, 文章提出的OrtWtls算法通過加權(quán)處理能有效減小反演誤差, 準(zhǔn)確計(jì)算石灰石漿液的粒徑分布。

4 結(jié) 論

鑒于超聲衰減譜的測量誤差對(duì)粒徑反演結(jié)果的影響, 文章提出了OrtWtls算法求解顆粒系的粒徑分布, 主要結(jié)論如下:

1) 根據(jù)設(shè)定單一顆粒系直徑為10 μm～80 μm的理論超聲衰減譜, 使用GCV準(zhǔn)則優(yōu)化正則化因子, OrtWtls算法計(jì)算結(jié)果與設(shè)定直徑基本相同, 相對(duì)誤差小于5%。

2) 以10 dB和20 dB含噪衰減譜為例, 分析OrtWtls算法的抗噪性, 該算法得到的粒徑與設(shè)定直徑相對(duì)誤差小于5%, 且多次測量的標(biāo)準(zhǔn)差小于1 μm。

3) 使用OrtWtls算法修正質(zhì)量濃度10%和20%的石灰石漿液實(shí)驗(yàn)衰減譜, 并計(jì)算漿液的粒徑分布, 測量結(jié)果與圖像法吻合, 其相對(duì)誤差小于10%。