張俊友，齊宏，*，王一飛，任亞濤，阮立明

（1.哈爾濱工業(yè)大學(xué) 能源科學(xué)與工程學(xué)院，哈爾濱150001；2.哈爾濱工業(yè)大學(xué) 空天熱物理工業(yè)和信息化部重點實驗室，哈爾濱150001）

碳黑是一種由碳?xì)浠衔锊煌耆紵a(chǎn)生的含碳物質(zhì)，納米級的碳黑初級顆粒會因布朗運動碰撞而互相附著，通常會形成具有分形結(jié)構(gòu)的團(tuán)聚體，即碳黑團(tuán)聚體［1］。在航空發(fā)動機(jī)燃燒室中，碳?xì)淙剂系娜紵伯a(chǎn)生碳黑團(tuán)聚體。高溫高壓下，碳黑團(tuán)聚體不但會磨損機(jī)體，而且燃燒室內(nèi)的火焰?zhèn)鳠嵋暂椛鋼Q熱為主，碳黑團(tuán)聚體的輻射特性會影響燃燒室的輻射換熱［2-3］。因此，碳黑團(tuán)聚體的粒徑和輻射特性等性質(zhì)對于航空發(fā)動機(jī)燃燒室的壽命和性能有著重要意義。此外，航空發(fā)動機(jī)的高空排放也是大氣中碳黑團(tuán)聚體的重要來源之一。大氣中的碳黑顆粒吸收太陽輻射，加熱大氣并冷卻地表，因此具有碳黑團(tuán)聚體被認(rèn)為是工業(yè)時代氣候變化的第二重要的人類因素［4］。同時，含有有毒物質(zhì)的碳黑團(tuán)聚體對人類健康有害，可能導(dǎo)致慢性肺部疾病、肺癌、哮喘等疾?。?］。由此可見，碳黑團(tuán)聚體對于火焰中輻射換熱、氣候變化和空氣質(zhì)量都至關(guān)重要［6］。因此，碳黑團(tuán)聚體的性質(zhì)測量研究吸引了大量海內(nèi)外學(xué)者的關(guān)注。

本文提出了一種利用光學(xué)信號間接測量碳黑團(tuán)聚體結(jié)構(gòu)特征和光學(xué)特性的方法。將光學(xué)方法用于火焰中碳黑團(tuán)聚體的性質(zhì)測量研究具有以下特點：高靈敏度，原位測量，不對樣本產(chǎn)生干擾等。目前，已有學(xué)者開展了利用不同光學(xué)信號反演得到火焰中碳黑物理性質(zhì)的研究工作。其中有代表性成果包括：1992年Sorensen等［7］提出了一種利用光的散射-消光信號實現(xiàn)原位光學(xué)測定碳黑團(tuán)簇的基本粒子直徑、每簇的基本粒子數(shù)和分形維數(shù)的新方法。2007年Iyer等［8］同樣采用散射-消光信號來實現(xiàn)光學(xué)參數(shù)的重建。2011年Link等［9］使用3個角度的散射光信號確定多分散團(tuán)聚體的參數(shù)包括粒子尺寸分布和分形維數(shù)。2016年Amin和Roberts［10］使用瑞利-德拜-甘斯多分散分形團(tuán)聚體（RDG-PFA）散射理論計算2個角度的散射-消光信號，分別反演碳黑的多種物理性質(zhì)，包括碳黑體積分?jǐn)?shù)、基本粒子直徑、團(tuán)聚體平均回轉(zhuǎn)半徑和基本粒子數(shù)量密度。2017年Moghaddam等［11］利用多體T矩陣（MSTM）模型的精度和瑞利-德拜-甘斯分形團(tuán)聚體（RGD-FA）模型的計算速度，實現(xiàn)從散射光的角度分布反演氣溶膠中碳黑團(tuán)聚體的粒徑分布和結(jié)構(gòu)特征。但是上述研究存在一個共性問題，就是在精確已知很多物性的前提下進(jìn)行反演過程，這大大提高了反演測量準(zhǔn)備工作的難度，使反演方法偏離真實情況。實際上，很多物性是很難提前獲得或者精確測量的，這些已知參數(shù)的不確定度會嚴(yán)重影響反演方法的精度，使得提出的方法具有很大的局限性，甚至在實際測量中不具備可行性。本文為了克服這一問題，使用2種不同類型的光學(xué)信號，增加反問題輸入信息，實現(xiàn)了在關(guān)鍵物性參數(shù)都未知和測量誤差等諸多干擾下，碳黑團(tuán)聚體形狀特征參數(shù)和光學(xué)性質(zhì)參數(shù)精確且穩(wěn)定的同時反演。

1 正問題

1.1 分形團(tuán)聚體

本文研究對象為碳黑顆粒系統(tǒng)。如圖1所示，單個碳黑團(tuán)聚體由許多個近似球形的基本粒子組成，基本粒子間互相吸附，形成鏈狀的分形結(jié)構(gòu)。如果假設(shè)所有基本粒子為球形且具有相同的粒徑，就可以使用分形理論對單個碳黑團(tuán)聚體的幾何特征進(jìn)行參數(shù)化描述［12］：

圖1 碳黑團(tuán)聚體分形示意圖Fig.1 Schematic of fractalmorphology of soot aggregate

式中：Np為單個碳黑團(tuán)聚體中基本粒子總數(shù)；kf為分形前置因子；Rg為碳黑團(tuán)聚體的回轉(zhuǎn)半徑；dp為碳黑團(tuán)聚體中基本粒子的直徑；Df為分形維數(shù)。

測量的目標(biāo)區(qū)域內(nèi)有大量不同的分形碳黑團(tuán)聚體。如果以回轉(zhuǎn)半徑表征不同碳黑團(tuán)聚體的尺寸，回轉(zhuǎn)半徑的分布函數(shù)近似滿足對數(shù)正態(tài)分布［12］：

式中：μg和σg為分布函數(shù)的特征參數(shù)，分別為Rg的平均值和標(biāo)準(zhǔn)差。

1.2 RDG-PFA散射理論

RDG-PFA模型中，基本粒子的微分散射截面為［12］

單個碳黑團(tuán)聚體的微分散射截面為［12］

其中：q＝（4π／λ）sin（θ／2）為RDG-PFA散射理論中一個重要的物理量，與散射角θ和入射激光波長λ有關(guān)；M＝4；C1＝2M／（3Df）；C2＝C3＝0，C4＝1。

因此，整個測量體積內(nèi)所有碳黑團(tuán)聚體的角度散射光強(qiáng)度為［12］

式中：α為測量系統(tǒng)的效率，介于0與1之間；Iinc為入射光強(qiáng)度；nagg為激光探測體積內(nèi)的碳黑團(tuán)聚體數(shù)量密度，即單位體積內(nèi)碳黑團(tuán)聚體的個數(shù)。

將式（1）、式（3）和式（4）代入式（6），可以整理成如下形式：

式中：C為一個與復(fù)折射率m和碳黑團(tuán)聚體數(shù)量密度nagg有關(guān)的函數(shù)，C＝naggF（m）；c1為多個常數(shù)參數(shù)的函數(shù)，c1＝f（α，Iinc，dp，k，kf）。

此外，單個碳黑團(tuán)聚體的吸收截面為［12］

光譜準(zhǔn)直透射率為［12］

式中：L為激光在測量體積內(nèi)的行程長度。

將式（1）和式（8）代入式（9），可以整理成如下形式：

式中：c2＝f（r，L，dp，k，kf），r＝F（m）／E（m）。

1.3 實驗裝置

如圖2所示，本文采用了文獻(xiàn)［14］中提出的廣角光散射（WALS）測量裝置，用于在10°～170°的寬角度范圍內(nèi)收集散射光。通過透鏡，將散射光成像到ICCD相機(jī)檢測器的芯片上，并允許以約0.6°的角度分辨率同時采集全散射圖像［15］，因此通過散射圖像可以得到散射光在不同角度上的強(qiáng)度。獲取某一散射角強(qiáng)度的具體實驗操作請參考文獻(xiàn)［14-15］。反射鏡上有2條相對的狹縫，保證激光進(jìn)出，進(jìn)出的激光都會被束流拘束器收集，因此可以得到光譜準(zhǔn)直透射率。本文提出的反演方法中，在10°～165°范圍內(nèi)每隔5°取一個散射角，共計32個散射角，作為反演信號。

圖2 廣角光散射測量裝置［14］Fig.2 W ide angle light scattering measurement system［14］

2 反問題

反問題的目標(biāo)是同時反演4個目標(biāo)參數(shù)，分別為C、分形維數(shù)Df、粒徑分布特征參數(shù)μg和σg。如圖3所示，整個反演過程如下：

圖3 反演過程Fig.3 Inversion procedure

步驟2 目標(biāo)參數(shù)的預(yù)測值更新，代入RDGPFA模型計算得到不同散射角度的散射光強(qiáng)度的預(yù)測值。）和）代入適應(yīng)度函數(shù)計算適應(yīng)度值，根據(jù)適應(yīng)度值反演算法更新目標(biāo)參數(shù)的預(yù)測值，使測量值）和預(yù)測值）逐漸接近。

步驟3 當(dāng)適應(yīng)度值Fobj小于目標(biāo)值eps或者迭代次數(shù)T達(dá)到最大迭代次數(shù)maxgens時，反演過程停止并輸出目標(biāo)參數(shù)的最終預(yù)測結(jié)果。

解決反問題的過程就是使目標(biāo)函數(shù)值最小化的過程。按照步驟2中描述的過程，可以將目標(biāo)函數(shù)定義為

式中：下標(biāo)i代表第i個散射角的散射光強(qiáng)度，共計n個散射角。

我國高校肩負(fù)著為社會建設(shè)發(fā)展培養(yǎng)人才的歷史使命，近年來隨著我國社會及經(jīng)濟(jì)的騰飛，對于人才的需求在不斷增加，高校的地位得到了進(jìn)一步的提升，因此，必須要有效地做好高校的建設(shè)及發(fā)展工作。

在最初的反演中，僅僅采用了多角度的散射信號。但隨著測量誤差的增大，4個目標(biāo)參數(shù)的反演誤差無法同時滿足精度要求。因此，需要增加更多的光學(xué)信號來反映碳黑團(tuán)聚體的內(nèi)部信息，僅僅增加散射角數(shù)量效果微弱。加入光譜準(zhǔn)直透射信號是一個很好的選擇，在已有的設(shè)備條件下，光譜準(zhǔn)直透射率是方便實現(xiàn)的，并且可以實現(xiàn)與散射信號的同時測量。在后面的分析中也證明了采用散射與透射信號結(jié)合的方法要比僅使用散射信號病態(tài)性更弱，有利于在有誤差干擾下反演過程向全局最優(yōu)解收斂。此時，目標(biāo)函數(shù)變?yōu)?/p>

本文使用的CMA-ES算法的具體原理和源代碼在文獻(xiàn)［16］有詳細(xì)介紹，不再贅述。

圖4 沒有噪聲與10%高斯噪聲下的角度散射光強(qiáng)度對比Fig.4 Comparison of angular light scattering intensity under no noise and 10% Gaussian noise

3 結(jié)果和討論

3.1 殘余適應(yīng)度分析

為了說明本文采用的多角度散射-準(zhǔn)直透射率組合信號的優(yōu)勢所在，對比目標(biāo)參數(shù)分布范圍下的殘余適應(yīng)度分布情況。由于受限于圖片表達(dá)的維度，最多只能分析2個參數(shù)的殘余適應(yīng)度值的分布情況，為了展現(xiàn)問題的收斂特性，圖片中采用3D云圖與2D投影結(jié)合的展示方式。

圖5（a）、（b）分別為復(fù)折射率m采用2種信號組合時的殘余適應(yīng)度分布?？梢园l(fā)現(xiàn)，在僅使用散射信號的情況下，復(fù)折射率反演結(jié)果不唯一，體現(xiàn)在圖5（a）中成谷狀的收斂區(qū)間，這意味著反演存在嚴(yán)重的多峰情況，使反問題精確求解變得非常困難。加入透射信號后，有效改善了這一現(xiàn)狀，使適應(yīng)度函數(shù)收斂于一點，不僅使目標(biāo)參數(shù)的精確反演變得易于實現(xiàn)，也會使反演收斂速度明顯提高。

圖5 兩種信號方案下的殘余適應(yīng)度分布Fig.5 Residual fitness distribution contour under two signal schemes

3.2 碳黑團(tuán)聚體參數(shù)的同時反演結(jié)果及分析

本文中入射激光的波長是806.5 nm，從文獻(xiàn)［17］中查得在此波長下乙炔火焰碳黑的復(fù)折射率m＝1.57-0.46i。2個粒徑分布函數(shù)的特征參數(shù)的真值定為μg＝90 nm和σg＝1.6，相應(yīng)的回轉(zhuǎn)半徑范圍為0～500 nm。目標(biāo)參數(shù)C根據(jù)所含各個參數(shù)計算，C＝0.8。分形維數(shù)Df根據(jù)其數(shù)學(xué)意義應(yīng)該為1～3，本文定義Df＝1.65，可以產(chǎn)生合理的分形碳黑團(tuán)聚體。各個參數(shù)在反演時需要給定一個搜索范圍，為了證明本文方法在大搜索范圍的適用性和穩(wěn)定性，所選取的搜索范圍都盡可能得大，初始值就在此范圍內(nèi)線性隨機(jī)產(chǎn)生。如表1所示，每個參數(shù)的搜索范圍都是遠(yuǎn)遠(yuǎn)大于可能存在的區(qū)域，如μg不可能超出最大500 nm或小于最小0 nm的粒徑范圍。

實際上，CMA-ES算法除了初始化過程外并不使用此搜索范圍，而是向全域范圍擴(kuò)展式地搜索，因此實際的搜索范圍要比給定的搜索范圍更大。此外CMA-ES算法雖然有許多算法參數(shù)，但因為其引入了參數(shù)自適應(yīng)策略，算法參數(shù)會在算法進(jìn)程中自行調(diào)整，因此并不需要為尋找合適的算法參數(shù)費心。但算法的收斂精度和最大迭代次數(shù)需要聲明，如表2所示。

表1 目標(biāo)參數(shù)的真實值和搜索范圍Tab le 1 O riginal value and search range of target param eters

eps是目標(biāo)收斂精度，設(shè)eps＝10-10，只有在沒有噪聲的情況下，反演才會因目標(biāo)函數(shù)值滿足目標(biāo)精度而停止。而有噪聲存在時，最終的目標(biāo)函數(shù)值都會降到一個相對最低值。maxgens是最大迭代次數(shù)，設(shè)為maxgens＝1 000。CMA-ES算法是一種隨機(jī)算法，具有一定的隨機(jī)性。因此，每種情況都重復(fù)50次以減少隨機(jī)性的影響。限于篇幅，本文只展示了所有算例中的一個，其余基于不同目標(biāo)參數(shù)值的算例都能取得類似精度的反演結(jié)果。

評價反演結(jié)果的相對誤差絕對值定義如下：

式中：εrel為相對誤差；zori為每個參數(shù)的真實值；zest為預(yù)測值。相對誤差的最大可接受值定為1%。本文使用標(biāo)準(zhǔn)差來評價反演結(jié)果的集中程度，即數(shù)據(jù)集的標(biāo)準(zhǔn)差越大數(shù)據(jù)越分散，分布范圍越大，反之，數(shù)據(jù)集越集中，分布范圍越小。

表3列出了在不同高斯噪聲下使用不同信號組合的目標(biāo)參數(shù)反演結(jié)果?？梢园l(fā)現(xiàn)，在沒有噪聲的情況下，2種信號方案都可以非常準(zhǔn)確地反演4 個參數(shù)，每個參數(shù)的相對誤差都小于0.005%，而且標(biāo)準(zhǔn)差也都達(dá)到很小的數(shù)量級，說明50次反演結(jié)果都很好地收斂到全局最優(yōu)解（目標(biāo)參數(shù)真值）。隨著噪聲的增大，4個目標(biāo)參數(shù)的相對誤差都不可避免地增加。僅采用散射信號，4個參數(shù)在6%和10%的高斯噪聲下都超出可接受范圍（1%），其中μg更是達(dá)到了20%以上，而且標(biāo)準(zhǔn)差的數(shù)量級對比沒噪聲的情況驟增。不過對比6%，在10%下誤差及標(biāo)準(zhǔn)差變化趨勢不是很明顯。

表2 CM A-ES算法的參數(shù)設(shè)定值Tab le 2 Param eter value setting of CM A-ES algorithm

表3 不同高斯噪聲下使用不同信號方案的目標(biāo)參數(shù)反演結(jié)果Table 3 Inversion results of target param eters ob tained by different signal schem es under differen t Gaussian noise

采用散射＋透射的信號組合，除10%高斯噪聲情況下μg的相對誤差超出1%外，其余情況下，各個參數(shù)的相對誤差都小于1%，達(dá)到可接受范圍。隨著噪聲增大，反演結(jié)果的相對誤差和標(biāo)準(zhǔn)差基本上呈增加趨勢，不過此信號組合有效限制了增加的幅度。對比相同噪聲下僅采用散射信號的反演結(jié)果，可以發(fā)現(xiàn)同時采用散射和透射信號明顯優(yōu)化了反演結(jié)果，成功地將糟糕的反演結(jié)果提升到精確反演的程度，尤其是測量誤差不超過6%的情況下。各個參數(shù)的反演誤差和標(biāo)準(zhǔn)差數(shù)量級都顯著下降。這說明在噪聲干擾下，引入透射信號有效改善了問題的病態(tài)性，提高了反演的穩(wěn)定性與抗噪性。

此外還發(fā)現(xiàn)不論什么情況下，μg的參數(shù)標(biāo)準(zhǔn)差和相對誤差都是4個參數(shù)中最大的，這與3.1節(jié)的殘余適應(yīng)度云圖的分析結(jié)論相吻合，換言之就是4個參數(shù)中最難實現(xiàn)精確反演的，這也是散射＋透射信號在10%高斯噪聲下，相對誤差唯一超出范圍的參數(shù)。繪出并對比各個情況下的粒徑分布函數(shù)后（見圖6）發(fā)現(xiàn)，μg對分布函數(shù)的影響沒有σg大，即使其相對誤差達(dá)到3.75%，但只要σg反演結(jié)果足夠精確，反演得到的分布函數(shù)與真實的分布函數(shù)也幾乎相同。而如果σg反演結(jié)果不夠準(zhǔn)確（僅使用散射信號的情況），分布函數(shù)的曲線就很明顯偏離了真實的分布曲線。因此可以認(rèn)為，在10%高斯噪聲下，同時使用散射和透射信號得到的反演結(jié)果也達(dá)到了精確級別。

值得一提的是，CMA-ES算法為本文方法帶來收斂速度上的優(yōu)勢。圖7展示了不同高斯噪聲下，采用散射＋透射信號的反演過程中適應(yīng)度函數(shù)值的下降過程（分別是50次計算中具有代表性的一次，其余49次的下降情況也基本相同）。迭代上限是1 000，無高斯噪聲的情況下因達(dá)到目標(biāo)精度而提前結(jié)束反演過程。而在6%和10%高斯噪聲的干擾下，500代之后也基本平穩(wěn)不再明顯下降，因此只展示了前450代的收斂過程。從中可以發(fā)現(xiàn)，3種高斯噪聲下，250代以內(nèi)都很迅速地完成收斂過程，250代時平均耗時為3.5 s。250代之后，對于無干擾的情況，曲線會迅速降到目標(biāo)精度使反演結(jié)束。而高斯噪聲的存在使另外2種情況只能降到一個相對最小值，并在該值附近微弱地波動，實際上有效的反演過程已經(jīng)完成，余下迭代過程只是在等待迭代次數(shù)達(dá)到上限。由于算法的特點，計算時間是與迭代次數(shù)成線性正比，450代時平均計算耗時6.3 s。在大量數(shù)值模擬的結(jié)果基礎(chǔ)上，認(rèn)為可以將最大迭代次數(shù)縮小到450代，以減少無效的計算時間提高效率。雖然3.5 s的計算用時還不足以到達(dá)在線測量的要求，但是也為接下來研究工作提供了良好的基礎(chǔ)。

圖6 不同高斯噪聲下使用不同信號方案反演得到的粒徑分布曲線Fig.6 Particle size distribution profiles obtained by inversion results of different signal schemes under different Gaussian noise

圖7 不同高斯噪聲下散射-透射組合信號的收斂過程Fig.7 Convergence process of scattering-transm ittance combined signal under different Gaussion noise

4 結(jié) 論

數(shù)值計算結(jié)果證明了本文提出的在一定高斯噪聲下精確并且穩(wěn)定地用于同時重建碳黑團(tuán)聚體粒徑分布、分形維數(shù)和常數(shù)參數(shù)方法的可行性。這一結(jié)論在于：對比在相同高斯噪聲下只使用散射信號的4個參數(shù)的反演結(jié)果，采用散射-透射信號組合的反演結(jié)果更好，即4個參數(shù)各自的反演相對誤差明顯降低，各自的標(biāo)準(zhǔn)差縮小。尤其是隨著高斯噪聲的增加，反演的結(jié)果精度達(dá)到了可以忽略的范圍（小于1%），實現(xiàn)精確反演。

1）散射-透射信號改善了反問題的病態(tài)性，并在10%的高斯噪聲下實現(xiàn)碳黑團(tuán)聚體的形態(tài)和粒徑分布參數(shù)的精確同時反演。

2）使用CMA-ES算法使得反演過程迅速收斂，初步達(dá)到3.5 s內(nèi)基本完成反演過程的效果。

3）本文得到的數(shù)值結(jié)果是在很大的目標(biāo)參數(shù)搜索空間下得到的，這些范圍在其數(shù)學(xué)意義和參數(shù)自身性質(zhì)上都足夠大，因此該方法的大搜索范圍普適性可以得到驗證。