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

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

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

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

1 正問(wèn)題

1.1 分形團(tuán)聚體

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

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

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

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

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

1.2 RDG-PFA散射理論

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

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

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

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

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

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

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

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

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

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

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

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

1.3 實(shí)驗(yàn)裝置

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

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

2 反問(wèn)題

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

圖3 反演過(guò)程Fig.3 Inversion procedure

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

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

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

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

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

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

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

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

3 結(jié)果和討論

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

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

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

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

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

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

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

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

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

評(píng)價(jià)反演結(jié)果的相對(duì)誤差絕對(duì)值定義如下：

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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

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