颗粒物在日常生活、工业生产和医药行业中普遍存在。颗粒粒度就是描述颗粒占据空间大小的尺度。颗粒粒度对产品性能有着重要影响，如在食品生产和药物制备中，其中间体大多是粉体状颗粒，合适的颗粒大小及分布是保证食品品质或药效特性的重要一步^[1]；煤炭燃烧效率与煤粉颗粒粒度大小密切相关^[2-3]；在医学中，人体血细胞中含有的红细胞、白细胞及血小板的大小与数量直接反映人体健康状况等^[4]。因此，进一步发展先进的、低成本的颗粒粒度测量技术及仪器具有现实意义与需求^[5]。

颗粒粒度测量的实现方法多种多样，如筛分法、沉降法、消光法、光衍（散）射法等^[6]。其中光衍（散）射法由于适用性广、可测范围大、速度快、精度高、非侵入性、可在线测量且结构简单等显著优点得到快速发展，现已成为颗粒粒度分析领域的主流技术之一^[7-8]。以此法为基础所研制的各种测量系统也得到了广泛的应用，如Polakowski等^[9]设计了光衍射测量系统，用于对土壤团聚体粒径分布进行测量。

传统激光粒度仪使用环形硅光电池或CCD线阵作为探测器，虽然技术成熟且工作稳定，但通常需要人工调整或专用软硬件进行对中，这使得操作不仅繁琐，而且容易受到外界环境的影响，同时分辨率也较低^[10]。随着CCD图像传感器技术的迅猛发展，其在替代传统光电探测器方面展现出优势，能够自动采集并记录高分辨率的散射光信号。面阵CCD探测器能记录颗粒的衍（散）射光全貌，不仅能反映颗粒的尺寸分布，还能反映颗粒的形貌特征^[11]，具备高自动化水平、高测量精度和较强适应性等优点^[12-13]。曲佩玙等^[14]基于面阵CCD，采用差分进化反演算法进行颗粒粒度测量，取得了一定效果。但这些反演算法存在收敛速度慢、计算资源消耗大，且针对双峰及多峰分布颗粒反演时，很可能陷入局部极值解等问题。因此，本文提出了利用光衍（散）射图样数据库，并结合衍射图样特征值的颗粒粒度反演方法。

本文搭建了激光颗粒衍（散）射测量系统。采用面阵CCD相机接收衍射光信号，利用计算机数据库的存储搜寻技术^[15]，预先编写程序计算，并建立了在被测颗粒测量范围内各种不同尺寸颗粒的归一化衍射图样数据库。采用数据库最小二乘搜寻法，并结合衍射图样本身特征对多种粒度颗粒进行了反演仿真和实验研究。此方法无需对衍射图样手动定中，可直接根据衍射图样特征，逐个比较数据库中所存数据与实测数据，二者最佳吻合时所对应库中的颗粒尺寸，即认为是被测颗粒尺寸。本文提出的数据库反演算法，区别于人工蜂群算法、差分进化算法等基于方程求解的反演算法，单次反演计算时间均在十几秒以内，计算速度较传统方程求解类算法有提高^[15-16]。原因在于：基于方程求解的反演算法计算量大，通常需要进行复杂的数值计算，如求解非线性方程、迭代计算等，这些过程往往需要大量的运算和较长的计算时间，且针对双/多峰分布的颗粒群的反演结果精度较低；而基于数据库的反演算法是利用已预先计算出的存储的光衍射曲线（作为数据库）和输入曲线进行比对，结合最小二乘法等快速匹配方法确定颗粒信息。该过程仅涉及曲线比对与误差计算，无需复杂的数学模型推导或数值求解。同时，算法通过对不同尺寸颗粒的逐个比较，可有效避免其他反演方法将局部极值误判为全局最优值的缺陷^[14,16]。除仿真计算外，本文还对5种不同尺寸的国家标准乳胶球微粒及组合进行了实测验证。

1 颗粒光衍射理论及数据库 1.1 颗粒光衍射理论

图1为激光颗粒衍射测量原理图。假定颗粒为球形，直径为D。根据光衍射（或散射）理论，在光强为$ {I}_{0} $、波长为$ \lambda $的平行单色光照射后，置于傅里叶透镜焦平面处的CCD探测器各环形单元上的衍射光强可表示为^[6]

$ I(\theta )=\frac{{\text{π}}^{2}{D}^{4}{I}_{0}}{16{f}^{2}{\lambda }^{2}}{\left(\frac{2{\mathrm{J}}_{1}(X)}{X}\right)}^{2} $

(1)

式中：$ \theta $为衍射角；$ f $为傅里叶接收透镜的焦距；$ {\mathrm{J}}_{1} $为1阶Bessel函数，$ X $为

$ X=\frac{\text{π}D\sin \theta }{\lambda } $

(2)

本文采用FP平面上的CCD相机面阵探测器。如图1（b）所示，将其等距分成若干个同心圆环。设第n环的面积为$ {H}_{n} $，则第n环（设环的半径从$ {S}_{n} $到$ {S}_{n+1} $，对应的衍射角从$ {\theta }_{n} $到$ \theta{{{_{{n+1}}}}} $）单位面积上的衍射光强^[6]为

$ {e}_{n}=\int \limits_{{S}_{n}}^{{S}_{n+1}}I(\theta )\frac{2\text{π}S}{{H}_{n}}\mathrm{d}S $

(3)

计算式（3）各环上的衍射光强$ {e}_{n} $，可得到单个颗粒在CCD探测器上的衍射光强分布。图2为计算得到的不同颗粒粒径的光衍射图样在某个截面上的衍射光分布曲线。其中横坐标为面阵CCD器件划分的等距同心圆环数，纵坐标为归一化后各环单位面积上的衍射光强。从图2可看出，不同粒径大小对应不同的衍射光曲线，颗粒粒径越大，衍射光相对变化越快，光强分布也越集中。

当被测颗粒为含有多个颗粒的颗粒群时，假设颗粒群的光散射符合不相关的单散射^[6]，则多个颗粒的衍射光强是单个颗粒衍射光强的N倍（N为颗粒数），故有

$ {E}_{n}=N\int \limits_{{S}_{n}}^{{S}_{n+1}}I(\theta )\frac{2\text{π}S}{{H}_{n}}\mathrm{d}S $

(4)

实际测量时，一般将被测颗粒尺寸区间划分为若干个子区间（若干个分档）$ \left[\left(D_{\min},D_1\right),\cdots, \left(D_i,D_{i+1}\right),\cdots,\left(D_{\mathrm{\mathit{m}}},D_{\max}\right)\right] $，其中$ {D}_{\min } $，$ {D}_{\max } $分别为最小和最大被测颗粒尺寸。设i代表颗粒尺寸的第i个小分档，则每档尺寸中颗粒数为$ {N}_{i} $（或重量为$ {W}_{i} $），这一组$ {N}_{i} $（或$ {W}_{i} $）就代表了颗粒的尺寸分布。$ {N}_{i} $和$ {W}_{i} $的关系为

$ {{N}}_{i}=\frac{6{W}_{i}}{\text{π}\rho {D}_{i}{}^{3}} $

(5)

式中，$ \rho $为颗粒的密度。通常，颗粒群的尺寸分布符合某个函数，当颗粒群的尺寸分布函数已知时，便可以求得各尺寸档中的颗粒数$ {N}_{i} $（或重量为$ {W}_{i} $）。鉴于实际情况下大多数颗粒的尺寸分布符合Rosin-Rammler（R-R）分布^[6,9]，因此本文用R-R分布函数作为颗粒群尺寸的分布函数。R-R函数的形式为

$ R(D)=\exp \left[-{\left(\frac{\;D\;}{\overline{D}}\right)}^{K}\right] $

(6)

式中：$ R(D) $为颗粒粒径大于D的颗粒累积重量分数；$ \overline{D} $为尺寸参量，其值大小反映颗粒群整体尺寸大小；$ K $为分布参量，反映颗粒群尺寸的分散程度。若颗粒群的尺寸分布为双峰或多峰分布，则也可由单个R-R函数的组合来模拟^[6,17]。如双峰分布形式可写为^[6]

$ {R}_{2}(D)=A{R}_{1}\left(D\right)+\left(1-A\right){R}_{1'}\left(D\right) $

(7)

式中：$ A $为比例系数，其值介于$ \left[0,1\right] $内，$ {R}_{1}(D) $和$ {R}_{1'}(D) $为两个不同的单峰分布，可由式（6）求得。当分布函数确定后，则各尺寸小档中的重量分数$ {W}_{i} $为

$ {W}_{i}=\left(\exp \left[-{\left(\frac{{D}_{i}}{\overline{D}}\right)}^{K}\right]-\exp \left[-{\left(\frac{{D}_{i+1}}{\overline{D}}\right)}^{K}\right]\right)\times 100\mathrm{\% } $

(8)

将式（1）代入式（4）积分并化简后，可得

$ \begin{split} {E}_{n}=&C\sum\limits_{i}\frac{{W}_{i}}{{D}_{i}{H}_{n}}\left[\mathrm{J}_{0}^{2}({X}_{i,n})+\mathrm{J}_{1}^{2}({X}_{i,n})-\right.\\& \left.\mathrm{J}_{0}^{2}({X}_{i,n+1})-\mathrm{J}_{1}^{2}({X}_{i,n+1})\right] \end{split} $

(9)

式中：C为常数；n为CCD探测器等距划分的环数。为了方便且不失一般性，本文设n的最大值$ {n}_{\max } $为30，也就是将CCD划分为等距的30个环。$ {\mathrm{J}}_{0} $、$ {\mathrm{J}}_{1} $分别是0阶、1阶Bessel函数，根据式（2）有

$ {X}_{i,n}=\frac{\text{π}{D}_{i}{\theta }_{n}}{\lambda },\quad{X}_{i,n+1}=\frac{\text{π}{D}_{i}{\theta }_{n+1}}{\lambda } $

(10)

式（9）为一方程组，也可写成矩阵形式

$ \left[\begin{array}{c}E_1 \\ E_2 \\ \vdots \\ E_{30}\end{array}\right]=\left[\begin{matrix}T_{1,1} & T_{2,1} & \cdots & T_{30,1} \\ T_{1,2} & T_{2,2} & \cdots & T_{30,2} \\ \vdots & \vdots & \vdots & \vdots \\ T_{1,30} & T_{2,30} & \cdots & T_{30,30}\end{matrix}\right]\left[\begin{array}{c}W_1 \\ W_2 \\ \vdots \\ W_{30}\end{array}\right] $

(11)

等式左边为光强分布列向量，等式右边的$ {T}_{i,j}=C\displaystyle\sum\limits_{i,j}\dfrac{[\mathrm{J}_{0}^{2}({X}_{i,n})+\mathrm{J}_{1}^{2}({X}_{i,n})-\mathrm{J}_{0}^{2}({X}_{i,n+1})-\mathrm{J}_{1}^{2}({X}_{i,n+1})]}{{D}_{i}{H}_{n}} $为方程系数，为一确定数。根据式（11），当测得光强分布$ {E}_{n} $后，便可基于该方程，反演计算得到尺寸分布$ {W}_{i} $。反之，当尺寸分布$ {W}_{i} $确定后，可计算出光强分布$ {E}_{n} $。反演计算过程简述如下：先设定R-R分布函数的2个参数$ \left(\overline{D},K\right) $，由式（8）计算各个$ {W}_{i} $，再由式（11）计算得到光强分布$ {E}_{n} $；将其与实测的由面阵CCD探测器获得的颗粒衍射图样，即光强分布$ {E}_{n}{}^{\prime} $作比较，若二者达到最佳拟合，则设定的$ \left(\overline{D},K\right) $就是所求的颗粒尺寸分布；若未达到拟合，则重新设定$ \left(\overline{D},K\right) $，如此循环，理论上总能找到一组与实测光强分布最佳拟合的$ \left(\overline{D},K\right) $，即为所求的颗粒尺寸分布。为提高反演算效率和准确度，本文采用PostgreSQL数据库技术。

1.2 数据库的建立

建立的数据库中存储了ID、颗粒尺寸分布参量、各环衍射光强$ {E}_{n} $数据及衍射图样（曲线）特征值4个部分，其中ID具有唯一性，用于标识每条数据。首先，根据式（11），通过MATLAB编程计算出对应的衍射光曲线，这些曲线由30个归一化后的单位面积上的衍射光强$ {E}_{n} $数据值构成。然后，将这些衍射光数据值和衍射图样特征值存入PostgreSQL数据库。PostgreSQL是功能强大的开源数据库管理软件，支持不同编程语言和自定义计算程序，如欧氏距离计算。

颗粒的衍射光曲线具有多种特征，如形状和数据特性，其中数据特性包括离散度、集中度和变化趋势等。本文采用颗粒衍射光的艾里斑直径（又叫圆孔的零级衍射斑，用于评价衍射角$ \theta $的弥散分布程度）和衍射曲线本身斜率作为特征值。利用合适的特征值，可缩小数据库的搜寻范围，从而极大减少计算时间。否则，需要将实际衍射光强$ {E}_{n}{}^{\prime} $数据值与库中大量数据逐一比较，增加计算时间。艾里斑直径$ {D}_{\text{airy}} $计算式^[5]为

$ {D}_{\text{airy}}=\frac{1.22\lambda }{\theta } $

(12)

斜率特征$ {S}_{\text{lope}} $计算式为

$ {S}_{\mathrm{lope}}=\sum\limits_{1}^{{n}_{\max }}\frac{{E}_{n}}{n} $

(13)

计算各个颗粒尺寸的衍射光曲线及特征值后，存入到PostgreSQL数据库中。

2 颗粒粒度的反演方法 2.1 反演算法

本文采用均方根误差（root mean square error, RMSE）为指标来评价数据库搜寻结果，表达式为

$ {R}_{\mathrm{mse}}=\sqrt{\frac{1}{{n}_{\max }}\sum\limits_{i=1}^{{n}_{\max }}{[{{L}_{\text{in}}}({{D}_{i}})-{{L}_{\text{out}}}({{D}_{i}})]}^{2}} $

(14)

式中：$ {L}_{\text{in}}(D) $表示数据库中理论上的颗粒衍射光强和特征值数据；$ {L}_{\text{out}}(D) $表示实测的颗粒衍射光强和特征值数据；$ {n}_{\max } $为30。数据库反演仿真计算的具体运行流程步骤如图3所示。首先，模拟实际颗粒的尺寸分布，设置颗粒参数$ \left(\overline{D},K\right) $，即可依据光衍射（或散射）理论，得到该颗粒群的光衍射曲线；随后，对艾里斑直径$ {D}_{\mathrm{airy}} $采用最小二乘法搜索^[15]，获得在设定的排名指标（即误差范围）内的数据集合$ {a}_{x} $，$ x\in \{1,2,\cdots \} $。接着在$ {a}_{x} $范围内，对斜率特征$ {S}_{\mathrm{lope}} $再次搜索，最终得到数据集合b。对集合$ b $中的数据，按$ b $中每条衍射光数据与实际衍射光数据的排名指标值的升序排序，进行重新排列，因此集合$ b $中第一行数据的颗粒分布参量$ \left(\overline{D},K\right) $，即可作为被测的颗粒尺寸分布。

由于排序时存在大量数据无法同时一次性存入内存，因此采用多路归并排序（multipath merge sort, MMS）。该方法将数据集分割成多个较小内存容量的子集合，对每个子集合执行快速排序（quick sort，QS），并将各个子集极小值存入内存；随后，对各个极小值再次采用快速排序获得最小值，并将该最小值写入存储单元。接着，从最小值所在的子集合中取出下一个数替换到内存，重复上述过程，将各个内存中最小的数存入存储单元，循环操作直至所有数据集排序完成。

2.2 反演算法的数值仿真验证

表1、表2、表3分别为单峰、双峰和三峰颗粒尺寸分布的仿真计算结果，表中的颗粒直径单位均为$ \text{μm} $。表中第一列为设定的颗粒尺寸参数$ \left(\overline{D},K\right) $，第二列为数据库方法反演得到的颗粒尺寸参数$ \left(\overline{D},K\right) $，第三列为二者之间的均方根误差值RMSE。

表4为对各环单位面积上归一化后的衍射光强$ {E}_{n} $，加入不同随机噪声后，取10次反演计算的平均值作为仿真结果。可以看到单峰、双峰、三峰的RMSE值均小于0.1，验证了数据库反演算法的可行性；加入随机噪声后其反演结果同样稳定，RMSE也小于0.1，验证了其抗干扰性和精确性。

3 颗粒粒度的实验测量 3.1 颗粒衍射测量系统

搭建的激光颗粒衍射测量系统的光路图，如图4所示。用市售的普通半导体激光器作为光源，半导体激光器发出功率为5 mW、波长为650 nm的激光束，经过滤光片后照射到样品池上。光通过装有待测颗粒的样品池后发生衍射（散射），衍射光透过焦距为63 mm的傅里叶透镜汇聚到焦平面上。在焦平面处，用面阵CCD相机记录衍射光信息并连接计算机分析。CCD相机感光面范围为36 mm × 24 mm，单个像素元大小约为$ 6.4\;\mu \mathrm{m}\times 6.4\;\mu \mathrm{m} $。这些参数能满足实验测量的要求。

3.2 实验结果与分析

为验证数据库技术反演方法的有效性，使用编号为GBW120021、GBW120134、GBW120164、GBW120024、GBW120041的乳胶微粒国家标准物质作为实验被测颗粒。其D50粒径标称值分别为2，3.5，7，10和15.3$ \text{μm} $。在实验过程中，将标准颗粒溶液与水混合均匀，制备成待测颗粒溶液。CCD相机记录的部分衍射光图样如图5。

对CCD面阵像素点上颗粒衍射光图样各环积分^[14]，并作归一化处理，得到的各环积分值作为实际衍射光强$ {E}_{n}{}^{\prime} $，输入数据库进行反演计算。本方法对衍射图样无需对中，也不需要进行专门的图像处理。图6表示了直径为10$ \text{μm} $的标准微粒的衍射光实测分布曲线和相应的理论分布曲线。对比两者，其均方根误差为0.0174，二者之间十分吻合。

分别对5种标准颗粒进行实测，表5为实测结果。从表5可以看出，反演实测值与D50标称值的相对误差均在5%以内。

进行双峰分布实验时，将D50粒径分别为7$ \text{μm} $、15.3$ \text{μm} $及10$ \text{μm} $、15.3$ \text{μm} $颗粒两两混合于蒸馏水中，实验过程与单峰颗粒测量一致，表6为实验结果。从表6可知，标准颗粒D50标称值与实测值的相对误差均在5%内。

4 结　论

提出了一种基于激光颗粒衍射图样预先计算并建立数据库的颗粒尺寸反演方法。预先建立的数据库中，存储了在傅里叶接收透镜焦距对应的被测颗粒尺寸范围内，与小步长间隔的不同颗粒尺寸及分布对应的衍射图样数据。通过数据库系统的高效管理与查询能力，并充分利用颗粒衍射图样曲线的几个特征参数来缩小搜寻范围，提高运算效率，减少计算时间，实现了精准的全局搜索。文中对算法进行了数值仿真验证，并对标准颗粒进行了实测验证，得到的中位粒径D50误差小于5%。实验表明，本方法测量精度较高，能满足工业生产和医药行业对颗粒粒度的实际测量要求。值得指出的是，本方法对基于Mie 散射理论的激光颗粒粒度测量，也同样具有参考借鉴价值。