基于散射光強度的碳黑團聚體分形結構和粒徑分布同時反演

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

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

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

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

1 正問題

1.1 分形團聚體

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

圖1 碳黑團聚體分形示意圖Fig.1 Schematic of fractalmorphology of soot aggregate

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

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

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

1.2 RDG-PFA散射理論

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

單個碳黑團聚體的微分散射截面為［12］

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

因此，整個測量體積內所有碳黑團聚體的角度散射光強度為［12］

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

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

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

此外，單個碳黑團聚體的吸收截面為［12］

光譜準直透射率為［12］

式中：L為激光在測量體積內的行程長度。

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

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

1.3 實驗裝置

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

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

2 反問題

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

圖3 反演過程Fig.3 Inversion procedure

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

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

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

式中：下標i代表第i個散射角的散射光強度，共計n個散射角。

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

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

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

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

3 結果和討論

3.1 殘余適應度分析

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

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

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

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

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

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

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

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

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

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

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

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

表3 不同高斯噪聲下使用不同信號方案的目標參數(shù)反演結果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%，達到可接受范圍。隨著噪聲增大，反演結果的相對誤差和標準差基本上呈增加趨勢，不過此信號組合有效限制了增加的幅度。對比相同噪聲下僅采用散射信號的反演結果，可以發(fā)現(xiàn)同時采用散射和透射信號明顯優(yōu)化了反演結果，成功地將糟糕的反演結果提升到精確反演的程度，尤其是測量誤差不超過6%的情況下。各個參數(shù)的反演誤差和標準差數(shù)量級都顯著下降。這說明在噪聲干擾下，引入透射信號有效改善了問題的病態(tài)性，提高了反演的穩(wěn)定性與抗噪性。

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

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

圖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 結 論

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

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

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

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