结合条件预优的流动气溶胶动态光散射正则化反演

韩锦壮，李鑫强，申晋，王保珺，刘伟，王雅静

（山东理工大学电气与电子工程学院，山东 淄博 255049）

0 引言

动态光散射（Dynamic Light Scattering，DLS）技术是测量亚微米及纳米颗粒粒度及其分布的有效方法［1］，该方法通过测量悬浮颗粒散射光强信号的自相关函数（Autocorrelation Function，ACF），来获取颗粒的粒度分布（Particle Size Distribution，PSD）信息［2-4］。常规的DLS测量是针对非流动样品在固定容器中进行，样品中除布朗运动外的其它宏观运动被认为是不存在的。与稳态悬浮颗粒的DLS测量不同，对于流动颗粒，由于附加平移运动导致的对传统DLS理论的修正，在理论与实践中都极大地增加了动态光散射测量的难度。

对流动颗粒进行DLS测量的探索可追溯到20世纪80年代，CHOWDHURY D P等［5］推导了单分散流动颗粒的ACF模型，将原表征布朗运动的扩散项改为其与表征定向运动的平移项的乘积。1986年，TAYLOR T W和SORENSEN C M［6］通过研究入射光束的高斯形状特征，发现平移项所对应的平移特征时间是由流速和入射光束交点处的束腰半径所决定，与散射体相较在入射光束焦点的位置无关，他们因此将ACF表达式中的高斯光束半径改为高斯光束焦点处的束腰半径。1998年，WEBER R等建立了多分散流动颗粒的ACF理论模型［7］，进而提出可通过拟合ACF获取气溶胶的颗粒粒度、浓度和流速［8］。此后，用于自混频相干测量的流动布朗运动颗粒的ACF和功率谱模型也建立起来［9］。2021年，牟彤彤等基于WEBER R的ACF模型进行了流动气溶胶PSD的实测［10］，进而分析了流动对气溶胶PSD测量的制约机制［11］：流动气溶胶粒度反演困难的一个重要原因是，流速的增加导致ACF反演方程的病态性加剧，表现为反演方程核矩阵条件数增大。

正则化方法是改善病态方程核矩阵条件数的常用方法，其基本原理是通过对核矩阵的奇异值进行加法修正来实现核矩阵条件数的降低［12-13］，已在常规DLS粒度反演中得到广泛应用。对于流动颗粒，反演方程的病态程度随着流速的增加而显著增大，反演过程中正则化方法对病态程度的降低受到流速增加的制约。为在较高流速下得到更为准确的PSD，本文采用条件预优［14-15］结合Tikhonov正则化方法进行流动气溶胶的PSD反演，通过条件预优利用先验流速信息降低核矩阵条件数，显著改善了颗粒流动引起反演方程病态性增高导致的PSD准确性降低的问题。

1 流动气溶胶DLS测量原理

散射场中，光强ACF与电场ACF满足Siegert关系［16］

式中，g(2)(τ)为归一化的光强ACF，τ为延迟时间，β为相干因子，g(1)(τ)为归一化电场ACF。

对于单分散颗粒系，电场ACF为指数衰减函数

式中，Γ为衰减线宽，其表达式为

式中，d为颗粒粒度大小，kB、T、η、λ0、n和θ分别为玻尔兹曼常数、绝对温度、介质的粘滞系数、入射光波长、介质的折射率和散射角度。

对于流动颗粒，光强ACF衰减由扩散运动与平移运动两部分共同作用，若散射体内的颗粒数足够多，则其电场ACF为［10］

式中，v为气溶胶颗粒流速为聚焦束腰，ω0为入射光的束腰，f1为透镜的焦距。对于多分散气溶胶，以衰减线宽为权重来表示g(1)(τ)，其离散化形式为

式中，N为离散的颗粒点数，G(Γi)为归一化的衰减线宽分布，满足=1。将式（3）带入式（5），并直接以颗粒粒度为权重表示电场ACF，可以得到

式中，f(di)为离散PSD，满足=1。式（6）的向量形式为

式中，g是由g(1)(τj)组成的列向量，j=1，2，…，M，M为延迟通道数，f是 由f(di)组成的列向量，i=1，2，…，N，核矩阵A可表示为

式（7）为一典型的病态方程，g的微小改变便会引起解的巨大变化。Tikhonov正则化通过正则参数和正则矩阵作用于原病态核矩阵来改善方程的病态性，其表达式为

式中，||·||2为欧式范数，α为正则参数，采用L-curve准则选取，起到平衡病态核矩阵与正则化矩阵的作用，反映了正则化矩阵的权重大小。L为正则矩阵，选用二阶差分矩阵，起到修正病态核矩阵的作用。f为颗粒粒度分布，需要添加非负约束。式（9）的等价形式可转化为

不同于非流动悬浮颗粒粒度分布的常规DLS反演，由于流动气溶胶存在定向运动，其光强ACF为表征布朗运动的扩散项与表征定向运动的平移项间的乘积，导致反演方程的病态性增加。设核矩阵A中相应参数为分别为：λ0=532 nm，T=298.15 K，n=1.0003，θ=90°，kB=1.3807×10-23J·K-1，η=18.37×10-3mPa·s，入射光束腰ω0=0.409 mm，透镜焦距f1=175 mm。光强ACF的延迟时间为1 μs到4000 μs，按比例间隔分25组，每组内延迟时间按线性间隔。反演的PSD取值按8.32 nm间隔，在2.01 nm到1000.01 nm之间取120个点。可计算出流速为0 m/s、0.4 m/s、0.8 m/s、1.2 m/s、1.6 m/s和2.0 m/s时的条件数分别为3.63×1019、3.93×1022、4.08×1031、3.38×1048、1.15×1070和8.07×1097，可见，流速的增加导致病态性的显著增大。采用Tikhonov正则化后，方程的条件数分别降低为3.04×107、3.94×107、4.04×107、4.91×107、5.87×107和7.86×107，可以看出，相对于原方程，正则化方法可显著降低条件数。但是，流动条件下的条件数明显高于非流动情况，随着颗粒流速的增加正则化方法的条件数也随之增加，直接的后果是，流速的增加导致反演结果失真。

为降低流速增加对反演结果的影响，本文对反演方程进行条件预优化处理，即，构造与式（7）等价的方程组

式中，P＞0，为条件预优矩阵，满足

可利用先验流速信息和延迟时间构造为M×M维方阵

将式（13）带入式（11）

得条件预优结合Tikhonov正则化方法求取f(di)的表达式

不难看出，经过条件预优处理后的核矩阵条件数不再受流速的影响。

2 模拟数据反演

模拟的多分散气溶胶PSD采用正态分布［17］

式中，a和b为PSD峰值权重的比例参数，且满足a+b=1，μ1和μ2为PSD峰值对应的粒度，σ1和σ2为颗粒粒度分布的标准差。表1为模拟PSD的分布参数，模拟时通过分段积分对连续分布离散化。模拟实验条件为λ0=532 nm，ω0=0.409 mm，f1=175 mm，T=298.15 K，n=1.0003，θ=90°，kB=1.3807×10-23J·K-1，η=18.37×10-3mPa·s，β=0.7。模拟中对应光强ACF的延迟时间按1 μs到4000 μs分25组，取150个值，分组按组间比例间隔、组内线性间隔。单峰粒度范围从208 nm到1000 nm，双峰粒度范围从8 nm到1000 nm，颗粒粒度间隔均为4 nm。

表1 模拟PSD的参数Table 1 PSD simulation parameters

为与实测情况吻合，模拟中数据加入高斯随机噪声

式中，n（τ）为高斯随机噪声。

流速取值范围通过雷诺数限制，以保证流体处于层流状态。雷诺数为

式中，ρ0和D0分别为流体介质的密度和流体所经圆管的当量直径。取ρ0=1.293 kg/m3，D0=12 mm，由层流条件下雷诺数范围（Re≤2000）得气溶胶流速范围为0～2.3 m/s。以0.4 m/s为间隔，从0.4 m/s到2.0 m/s选取5个速度值作为已知条件参数。

粒度反演时，根据式（1）由光强ACF求得电场ACF，再通过电场ACF反演颗粒粒度。反演效果通过峰值位置相对误差EP和分布误差EF评估

式中，P表示峰值位置处的粒径，下标true和measure分别表示真实值和反演值。

图1和图2分别为600 nm单峰分布和200 nm/700 nm双峰分布气溶胶颗粒在不同流速下的反演结果，相应的性能参数如表2和表3所示。图中的“true PSD”表示真实的粒度分布，“Tik”代表Tikhonov正则化，“Pre-Tik”代表结合条件预优的Tikhonov正则化，表中的PP表示峰值位置，Cond表示反演方程核矩阵的条件数。

表2 600 nm单峰气溶胶反演结果的性能指标Table 2 The performance indices of the recovery of 600 nm unimodal PSDs

图1 600 nm单峰气溶胶不同速度下的反演结果Fig.1 The recovery of 600 nm unimodal aerosol PSDs at different velocities

图2 200 nm/700 nm双峰气溶胶不同速度下的反演结果Fig.2 The recovery of 200 nm/700 nm bimodal aerosol PSDs at different velocities

从图1和表2可以看出，对于600 nm单峰气溶胶，随着流速的增加，反演结果中PSD出现展宽，峰值位置左移，Tik方法的EP和EF分别由0和0.0020增大到0.0467和0.0069。在流速相同的条件下，Pre-Tik方法的反演结果优于Tik方法，并且随着流速的增大，Pre-Tik方法的优化效果愈趋明显。颗粒流速为0时，两种方法的EP和EF相同，当流速达2.0 m/s时，Pre-Tik方法得到的EP和EF分别降低到0.0267和0.0057。

从图2和表3可以看出，在200 nm/700 nm双峰气溶胶反演结果中，两个峰均向粒度减小的方向偏移，且偏移随流速增加逐渐增大，峰值位置误差和分布误差随之增大，但两种方法反演结果的性能指标差异显著，当流速为2 m/s时，Tik方法反演结果的EP和EF分别为0.6600/0.2743和0.0121，Pre-Tik方法的EP和EF则分别为0.4600/0.0914和0.0091，反演结果的准确性显著提高。

表3 200 nm/700 nm双峰气溶胶反演结果的性能指标Table 3 The performance indices of the recovery of 200 nm/700 nm bimodal PSDs

3 实验

实验数据取自本课题组研制的流动气溶胶DLS在线测量实验平台（图3（a）），实验装置主要由气溶胶发生器、导管、比色皿、激光、光子计数探测器、光子相关器等组成。单、双峰分布气溶胶分别由3321PLUS和SX-Q5气溶胶发生器产生，激光器采用MGL-III-532nm-10mW，光子计数探测器选用CH326、光子相关器选用Brookhaven BI-9000AT（图3（b）），流速通过热线风速仪SAR886A测量，其他实验参数与模拟参数相同。相关器取150个延迟通道，通道从1μs到4000μs分25组，组间比例间隔、组内线性间隔。颗粒粒度以10 nm间隔从10 nm到1000 nm取100个点。粒度反演时，通过来自光子相关器的光强ACF，计算出电场ACF，再由电场ACF反演PSD。反演结果如图4、图5、图6和表4、表5所示，反演结果的峰值位置误差和分布误差分别用EP0和EF0表示，0 m/s时的反演结果作为对比参照的气溶胶分布。

图3 实验装置和相关器Fig.3 Experimental devices and Correlator

图4 不同流速下气溶胶颗粒的光强ACFFig.4 Measured light intensity ACFs at different velocities

从图5和表4可以看出，对于单峰气溶胶，在流速为1.09 m/s和1.77 m/s时，Tik方法反演结果的EP0分别为0.0526和0.0877，EF0分别为0.0063和0.0190，Pre-Tik方法得到的EP0分别为0.0351和0.0526，EF0分别为0.0045和0.0162，后一方法反演效果好于前者。

表4 单峰气溶胶反演结果的性能指标Table 4 The performance indices of the recovery of unimodal PSDs

图5 不同流速时单峰气溶胶的反演结果Fig.5 The recovery of unimodal aerosol PSDs at different velocities

图6和表5表明，双峰气溶胶的正则化反演对流速表现得更为敏感，与Tik方法相比，Pre-Tik方法反演的峰值位置误差和分布误差均显著下降。当速度为1.28 m/s时，Tik方法反演出的小峰峰高仅为0速对应PSD小 峰 峰高的0.0926倍，而Pre-Tik方法反 演出的小 峰 峰高为0速对应PSD峰 高的0.2472倍，Pre-Tik方法反演结果的大峰峰值位置误差约为Tik方法的1/3，Pre-Tik方法的小峰峰值位置误差仅为Tik方法的1/5，改善效果显著。

表5 双峰气溶胶反演结果的性能指标Table 5 The performance indices of the recovery of bimodal PSDs

图6 不同流速时双峰气溶胶的反演结果Fig.6 The recovery of bimodal aerosol PSDs at different velocities

4 结论

在流动气溶胶DLS测量中，被测颗粒的定向流动导致基于光强ACF构建的反演方程病态性增加，流速越快，方程的病态程度越高，表现为反演方程核矩阵的条件数增加。作为降低病态性常用方法的正则化反演方法，其在流动导致方程病态性增加时，不能完全消除流动导致的条件数增加问题。正则化方法在流动颗粒反演中的这种局限性，可通过条件预优改善。通过条件预优处理，以先验流速信息和延迟时间构建对角阵形式的条件预优矩阵，对病态方程进行乘法修正，可降低流速对反演方程病态性的加剧作用，从而降低流动气溶胶DLS测量时正则化反演方法对流速增加的敏感性。模拟与实测数据的反演结果表明，对于流动颗粒的DLS测量，采用结合条件预优的正则化反演，可以克服正则化方法在降低条件数过程中存在的局限性，改善PSD反演的性能指标，提高流动气溶胶DLS测量的准确性。