本发明涉及一种利用短脉冲激光测量球形颗粒分形聚集特征参数的方法，属于颗粒聚集特征测量技术领域。
背景技术：
分形聚集形态广泛存在于颗粒态物质中，如大气气溶胶、航空发动机燃烧室中产生的碳烟颗粒等。颗粒之间的分形聚集影响了颗粒系自身光辐射特性(如吸收特性、散射特性等)，改变了光辐射在颗粒系中的传输过程，最终导致光辐射传输与颗粒系之间的相互作用机制发生改变。因此，开展颗粒分形聚集特征参数研究对分析航空发动机中碳黑颗粒生产、大气光辐射传输等均有重要的指导意义。
非接触式测量方法和接触式测量方法是两种常见的颗粒分形聚集特征参数测量方法。其中，非接触式测量通常借助获得颗粒系外部辐射信号和辐射逆问题求解技术间接获得颗粒分形聚集特征参数。相比于接触式测量，非接触式测量方法能实现在线监测，能获得具有时间和空间分辨能力的测量结果，且不会干扰测量对象，因此更受青睐。常见的非接触式测量方法有光谱消光法、红外发射ct法、双色法等。
与以往非接触式方法不同，本项发明利用短脉冲激光辐照球形颗粒系的一侧，然后用探测器获得另一侧的多角度时域散射信号，并用粒度分析仪测量颗粒系的粒径分布情况，用傅里叶光谱分析仪、k-k关系式、mie散射理论及人工蚁群优化算法等反演获得颗粒光学常数，最后基于测量获得的多角度时域散射信号结合逆问题求解技术间接得到球形颗粒分形聚集特征参数。本项发明充分借助短脉冲激光提供的丰富的时间分辨辐射测量信息，能有效解决以往研究中测量精度不高、抗干扰性差等问题。
技术实现要素：
本发明的主要目的在于：提供一种实验结合反演算法求解球形颗粒分形聚集特征参数的方法，基本思路是通过实验测得颗粒系的多角度时域散射信号，然后结合多体t矩阵理论模型和人工蜂群算法反演得到颗粒分形聚集特征参数。
本发明的方法是：首先将要测量的球形颗粒聚集体做成悬浮颗粒系样本，然后用短脉冲激光照射颗粒系一侧，用探测器测量颗粒系另一侧的多角度时域散射信号，再用粒径分析仪测量得到颗粒粒径分布情况，用傅里叶光谱分析仪结合k-k关系式、mie散射理论及人工蜂群算法获得颗粒光学常数。将探测器所测得多角度时域散射信号作为反问题的初始条件，结合多体t矩阵理论模型和人工蜂群优化算法反演得到颗粒分形聚集特征参数。
本发明的具体步骤如下：
一种利用单频调制激光辐照技术的球形颗粒分形聚集特征参数的测量方法，包括如下步骤：
步骤1：将待测颗粒装在有机玻璃样本容器中，使样本容器中的样本颗粒系处于悬浮流动状态；
步骤2：利用脉冲宽度为tp、波长为λ的短脉冲激光沿着与样本容器表面法线成θc角的方向入射到样本容器左侧表面，其中0＜θc＜π/2；用探测器在样本容器右侧表面测量不同角度方向θj上的时域散射信号，获得样本容器右侧表面角度时域散射信号强度
步骤3：测量样本颗粒系中颗粒聚集体总数n、球形颗粒总数n0及分形聚集体中颗粒数ns，同时测量样本颗粒系颗粒粒径分布情况，得到粒径分布函数f(d)，d为颗粒聚集体中单个颗粒直径；
步骤4：制作球形颗粒与溴化钾混合压片，测量压片的光谱透射率τλ,meas，结合kramers-kranigs关系式、mie散射理论及逆问题求解方法反演得到球形颗粒光学常数mλ＝nλ+ikλ，i为虚数单位，nλ为折射指数，kλ为吸收指数；
步骤5：计算颗粒分形聚集体的光谱吸收截面cabs,λ,pred和光谱散射截面csca,λ,pred；具体步骤如下：
步骤5.1：采用逆问题求解方法在颗粒聚集体特征参数可能的取值范围内随机假设一组颗粒聚集体特征参数值，即分形维数df、回转半径rg、前向因子kf；
步骤5.2：根据假设的颗粒聚集体特征参数值重构颗粒聚集体的微观几何结构；
步骤5.3：结合光学常数mλ＝nλ+ikλ和多体t矩阵模型，求得颗粒分形聚集体的光谱吸收截面cabs,λ,pred和光谱散射截面csca,λ,pred；
步骤6：根据光谱吸收截面cabs,λ,pred和光谱散射截面csca,λ,pred，计算得出样本颗粒系的光谱吸收系数κa,λ和光谱散射系数κs,λ；
步骤7：根据光谱吸收系数κa,λ和光谱散射系数κs,λ求解辐射传输方程，获得计算域内的任意位置x在θj方向上的时域散射信号强度iλ(t,x,θj)；
步骤8：利用步骤7获得的计算域内的时域辐射强度场，结合公式：
sλ(t,l,θj)pred＝id,λ(t,l,θj)+ic,λ(t,l,θj)(1)
(0≤θj<π/2,j＝1,2,...,n)
获得样本容器右侧表面角度θj上估计的时域散射信号强度sλ(t,l,θj)pred；式中，id,λ(t,l,θj)为估计的样本容器右侧表面角度θj上扩散光时域信号强度；ic,λ(t,l,θj)为估计的样本容器右侧表面角度θj上平行光时域信号强度；
步骤9：利用步骤2中探测器获得的测量样本容器右侧表面角度θj方向上的时域散射信号强度sλ(t,l,θj)exp(j＝1,2,…,n)，与步骤8中估计的时域散射信号强度sλ(t,l,θj)pred(j＝1,2,…,n)，再结合公式：
获得适应度函数fit；n为测量角度个数；
步骤10：判断步骤9中适应度函数值fit是否小于设定阈值ξ，若是，则将步骤5中假设的颗粒聚集体特征参数值作为结果，完成利用单频调制激光辐照技术的球形颗粒分形聚集特征参数的测量，否则重复步骤5至步骤10。
进一步的，步骤4包括如下步骤：
步骤4.1：破坏样本颗粒系的颗粒聚集特征，使其呈现单个分散状态，制作球形颗粒与溴化钾混合压片，测量压片的光谱透射率τλ,meas；
步骤4.2：采用逆问题求解方法在颗粒光学常数可能的取值范围内随机产生一组球形颗粒的光学常数：折射指数nλ和吸收指数kλ；
步骤4.3：根据光学常数方程mλ＝nλ+ikλ，计算出球形颗粒的光学常数mλ；
步骤4.4：根据mλ，结合mie散射理论计算压片的模拟光谱透射率τλ,pred；
步骤4.5：根据τλ,meas和τλ,pred构造并计算适应度函数fit′，判断fit′是否小于设定阈值η，若fit′<η，则步骤4.2产生的球形颗粒的光学常数即为真实的光学常数，否则继续采用逆问题求解方法随机产生一组颗粒光学常数，重复步骤4.2至步骤4.5。
进一步的，其特征在于，所述逆问题求解方法为人工蜂群算法。
更进一步的，模拟光谱透射率τλ,pred通过公式(3)计算得出，适应度函数fit′通过公式(4)计算得出，折射指数nλ和吸收指数kλ的关系式如(5)、(6)所示；
式中，l′为压片厚度；p表示柯西主值积分；π表示圆周率；qext(d,m,λ)为单个颗粒的衰减因子，由mie散射理论获得；λ0为真空中的波长。
进一步的，步骤5，基于分形理论和团簇-团簇分形聚集动力学仿真模型重构颗粒聚集体的微观几何结构，其中分形聚集特征参数之间满足关系式：
式中，rj为颗粒聚集体中第j个粒子中心到颗粒聚集体重心的几何距离。
进一步的，步骤5中颗粒分形聚集体的光谱吸收截面cabs,λ,pred和光谱散射截面csca,λ,pred的计算公式为：
cabs,λ,pred＝cext,λ,pred-csca,λ,pred(11)
式中，nmax为最大截断级数，k0为光在真空中的波矢，m、n、m'、n'表示球形谐波函数多项式展开的项数，实际计算过程中根据所需要的精度来选择这些参数的值。
进一步的，步骤5中样本颗粒系的吸收系数κa,λ和散射系数κs,λ通过公式：
实现。
进一步的，步骤7中辐射传输方程为：
式中，λ为短脉冲激光波长；φ(θ',θ)为θ′方向入射并从θ方向散射出去的散射相函数；id,λ(t,x,θ)为t时刻入射波长为λ在θ方向x处的扩散光时域辐射强度；计算域内平行光时域辐射场强度，利用lambert-beer定律：
ic,λ(ω,l,θj)＝δ(θj-θc)i0,λ(ω,0,θc)exp[-(κa,λ+κs,λ)l](15)
实现；式中，i0,λ(t,0,θc)为t时刻入射波长为λ入射角为θc的入射边界处平行光时域辐射强度；ic,λ(t,l,θj)为t时刻入射波长为λ的平行光在l处θj方向上平行光时域辐射强度。
进一步的，步骤2中脉冲宽度tp、波长λ的取值范围为10-12≤tp≤10-9s，0.6μm<λ<2.0μm。
本发明通过建立球形颗粒系多角度时域散射信号测量的正问题和反问题求解模型，解决颗粒系颗粒分形聚集特征参数不能直接测量和测量结果不准确的问题，提出了利用短脉冲激光测量球形颗粒分形聚集特征参数的方法。优点在于：采用短脉冲激光，该激光能提供丰富的时变测量信号，能提高测量结果的鲁棒性和精度；采用多体t矩阵理论模型，该模型理论核心是将入射场、散射体内部场及散射体外部场均严格采用矢量球谐函数展开，获得待求散射场和已知入射场之间的线性矩阵关系，能很精确的反应出颗粒之间的电磁相互作用；采用人工蜂群优化算法，该算法求解优化问题时有高稳定性、高灵敏度、高鲁棒性等优点，可以使得测量精确度提高10％。该项发明为研究颗粒分形聚集特征参数提供一种快速准确的方法，对航天、国防和民用工业具有十分重要的意义。
附图说明
图1为样本颗粒系左侧受到入射方向为θc波长为λ的短脉冲激光照射时的多角度时域散射信号测量示意图。
具体实施方式
本实施方式所述一种利用短脉冲激光测量球形颗粒分形聚集特征参数的方法，该方法的具体步骤为：
步骤一，根据颗粒系物性测量要求，将待测颗粒装在有机玻璃样本容器中。在实验测量过程中通过搅拌方式使得颗粒系处于悬浮流动状态，以确保样本容器内颗粒系的均匀分布。
步骤二，如图1所示，利用脉冲宽度为tp、波长为λ的短脉冲激光沿着与样本容器表面法线成θc角的方向入射到样本颗粒系左侧表面，其中，10-12≤tp≤10-9s，0<θc<π/2，0.6μm<λ<2.0μm；并用探测器在样本颗粒系右侧l处测量不同散射角度方向θj上的时域散射信号，获得样本右侧l处角度时域散射信号强度sλ(t,l,θj)exp(j＝1,2,…,n)，n为测量角度个数，l为样本容器厚度。
步骤三，利用扫描电镜测量样本中颗粒聚集体总数n、球形颗粒总数n0及分形聚集体中颗粒数ns，同时采用粒径分析仪测量得到样本颗粒系颗粒粒径分布情况f(d)。
步骤四，通过搅拌破坏颗粒聚集特征，使其呈现单个分散状态，制作球形颗粒与溴化钾(kbr)混合压片，由傅里叶光谱分析仪测量压片的光谱透射率τλ，同时结合kramers-kranigs(k-k)关系式、mie散射理论及人工蚁群优化算法反演得到颗粒光学常数mλ＝nλ+ikλ，i为虚数单位。
步骤五，利用逆问题算法思路假设出待测颗粒系的颗粒分形聚集特征参数(即分形维数df、回转半径rg、前向因子kf)，并基于分形理论和团簇-团簇分形聚集动力学仿真模型重构颗粒聚集体的微观几何结构，然后结合光学常数mλ＝nλ+ikλ和多体t矩阵模型计算得到颗粒分形聚集体的光谱吸收截面cabs,λ,pred和光谱散射截面csca,λ,pred，最后根据已测得样本中颗粒聚集体总数n，计算得出颗粒系光谱吸收系数κa,λ和光谱散射系数κs,λ，并通过对时域辐射传输方程求解，获得计算域内任意位置x在θj方向上的时域散射信号强度iλ(t,x,θj)。
步骤六，利用步骤五获得计算域内时域辐射强度场，结合公式：
sλ(t,l,θj)pred＝id,λ(t,l,θj)+ic,λ(t,l,θj)(1)
(0≤θj<π/2,j＝1,2,...,n)
获得测量样本右侧l处角度θj方向上估计的时域散射信号强度sλ(t,l,θj)pred。式中，id,λ(t,l,θj)为估计获得的样本右侧l处角度θj上扩散光时域信号强度；ic,λ(t,l,θj)为估计获得的样本右侧l处角度θj上平行光时域信号强度；n为测量角度个数。
步骤七，利用步骤二中探测器获得测量样本右侧l处角度θj方向上的时域散射信号强度sλ(t,l,θj)exp(j＝1,2,…,n)，与步骤六中预测的时域散射信号强度sλ(t,l,θj)pred(j＝1,2,…,n)，再结合公式：
获得逆问题算法中的适应度函数fit。式中，n表示测量角度个数。
步骤八，判断步骤七中适应度函数值fit是否小于设定阈值ξ，若是，则将步骤五中获得的待测颗粒系颗粒分形聚集特征参数(即分形维数df、回转半径rg、前向因子kf)作为结果，完成利用短脉冲激光测量球形颗粒分形聚集特征参数的方法，否则返回步骤五，采用人工蜂群优化算法重新修正预测的颗粒分形聚集特征参数。
步骤四采用傅里叶光谱分析仪测量压片的光谱透射率τλ，同时结合kramers-kranigs(k-k)关系式、mie散射理论及人工蚁群优化算法反演得到颗粒光学常数mλ＝nλ+ikλ，其中光谱透射率τλ和k-k关系式为：
式中λ表示激光波长；l′为kbr混合压片厚度；p表示柯西主值积分；π表示圆周率；qext(d,m,λ)为单个颗粒的衰减因子，可由mie散射理论获得；d为颗粒粒径。
步骤五基于分形理论和团簇-团簇分形聚集动力学仿真模型重构颗粒聚集体微观几何结构，其中分形聚集特征参数之间满足关系式：
式中，ns为颗粒聚集体中颗粒总数；d为颗粒聚集体中单个颗粒直径；kf为颗粒团聚体的前向因子；rg为颗粒聚集体的回转半径，rj为颗粒聚集体中第j个粒子中心到聚集体重心的几何距离。
步骤五获得颗粒系的吸收系数κa,λ和散射系数κs,λ通过公式：
实现。式中，n表示所测颗粒聚集体总个数。
步骤五获得计算域内的扩散光时域信号强度方法为，利用辐射传输方程：
实现。式中，λ为短脉冲激光波长；φ(θ',θ)为θ′方向入射并从θ方向散射出去的散射相函数；id,λ(t,x,θ)为t时刻入射波长为λ在θ方向x处的扩散光时域辐射强度；κa,λ为颗粒系吸收系数；κs,λ为颗粒系散射系数。此外，计算域内平行光时域辐射场强度，利用lambert-beer定律：
ic,λ(ω,l,θj)＝δ(θj-θc)i0,λ(ω,0,θc)exp[-(κa,λ+κs,λ)l](11)
实现。式中，i0,λ(t,0,θc)为t时刻入射波长为λ入射角为θc的入射平行光时域辐射强度；ic,λ(t,l,θj)为t时刻入射波长为λ的平行光在l处θj方向上平行光时域辐射强度；l为样本厚度。
以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。