光子学报  2018, Vol. 47 Issue (7): 0729001  DOI: 10.3788/gzxb20184707.0729001
0

引用本文  

黄钰, 申晋, 徐敏, 孙成, 刘伟, 孙贤明, 王雅静. 基于核矩阵扩展的动态光散射截断奇异值分解反演[J]. 光子学报, 2018, 47(7): 0729001. DOI: 10.3788/gzxb20184707.0729001.
HUANG Yu, SHEN Jin, XU Min, SUN Cheng, LIU Wei, SUN Xian-ming, WANG Ya-jing. Truncated Singular Value Decomposition Inversion Based on Extended Nuclear Matrix in Dynamic Light Scattering[J]. Acta Photonica Sinica, 2018, 47(7): 0729001. DOI: 10.3788/gzxb20184707.0729001.

基金项目

山东省自然科学基金(Nos.ZR2018MF032, ZR2016EL16, ZR2017LF026, ZR2017MF009)资助

第一作者

黄钰(1993-), 女, 硕士研究生, 主要研究方向为动态光散射测量. Email: yuhuang_93@126.com

通讯作者

申晋(1962-), 男, 教授, 博士, 主要研究方向为光电精密测试技术. Email: shenjin@sdut.edu.cn

文章历史

收稿日期:2018-02-04
录用日期:2018-03-25
基于核矩阵扩展的动态光散射截断奇异值分解反演
黄钰 , 申晋 , 徐敏 , 孙成 , 刘伟 , 孙贤明 , 王雅静     
(山东理工大学 电气与电子工程学院, 山东 淄博 255049)
摘要:针对截断奇异值分解方法进行动态光散射反演存在的颗粒粒度信息丢失问题,本文在分析自相关函数不同衰减时段粒度信息分布差异的基础上,提出利用每一角度核矩阵与对应角度下粒度信息在自相关函数不同延迟时刻的分布构建扩展矩阵的扩展截断奇异值分解方法.该方法通过用自相关函数中每一延迟时刻的粒度信息,调节同一时刻原核矩阵数据对信噪比的贡献,进而保留了更多的有效奇异值,减少了由于奇异值截断引起的信息丢失,在保证抗噪性的基础上,提高了自相关函数的信息利用率.在1×10-3噪声水平下,对一组单峰宽分布(260 nm)和三组双峰颗粒分布(250/750 nm)、(270/800 nm)以及(306/974 nm)的模拟动态光散射数据,进行了单角度、3角度和6角度反演.结果表明,与截断奇异值分解方法相比,采用扩展截断奇异值分解方法反演得到的峰值粒度误差和分布误差均明显减小.对306/974 nm颗粒体系的6角度实测数据的反演表明,采用扩展截断奇异值分解法得到的颗粒峰值粒度误差由非扩展方法的0.032/0.016降至0.029/0.006,且得到的峰值比更接近真实值.
关键词动态光散射    粒度分布    反演    颗粒测量    奇异值分解    
中图分类号:TN911.74      文献标识码:A      文章编号:1004-4213(2018)07-0729001-10
Truncated Singular Value Decomposition Inversion Based on Extended Nuclear Matrix in Dynamic Light Scattering
HUANG Yu , SHEN Jin , XU Min , SUN Cheng , LIU Wei , SUN Xian-ming , WANG Ya-jing     
(School of Electrical and Electronic Engineering, Shandong University of Technology, Zibo, Shandong 255049, China)
Foundation item: The Natural Science Foundation of Shandong Province (Nos. ZR2018MF032, ZR2016EL16, ZR2017LF026, ZR2017MF009)
Abstract: In order to solve the problem of particle size information loss in dynamic light scattering inversion by truncated singular value decomposition method, on the base of analyzing the difference of particle size information distribution in autocorrelation function at different delay time, extended truncated singular value decomposition method was proposed which uses nuclear matrix and particle size information distribution in autocorrelation function at different delay time in the same angle to build the extended nuclear matrix. The method adjusts the contribution of the original nucleus matrix datas to the SNR at the same time using the particle size information at each delay time in the autocorrelation function, thereby preserves more effective singular values and reduces information loss due to singular value truncation. On the basis of ensuring noise reduction, the information utilization for the autocorrelation function is improved. At noise level 1×10-3, the datas a set of single peak width distribution (260 nm) and three sets of bimodal particle distributions (250/750 nm), (270/800 nm) and (306/974 nm) simulated in dynamic light scattering in single angle, 3 angles and 6 angles was inverted. The results show that, compared with truncated singular value decomposition method, the peak size errors and distribution errors obtained by extended truncated singular value decomposition method are obviously reduced. The inversion of the measured datas for 306/974 nm particle system in 6-angles shows that the peak particle size error obtained by the extended truncated singular value decomposition method is reduced from 0.032/0.016 by the truncated singular value decomposition method to 0.029/0.006, and peak ratio is closer to actual value.
Key words: Dynamic light scattering    Particle size distribution    Inverse problems    Particle size    Singular value decomposition    
OCIS Codes: 290.0290;290.3200;290.5820;290.5850
0 引言

动态光散射(Dynamic Light Scattering, DLS)技术是测量亚微米及纳米颗粒的有效方法,已广泛应用于材料、化工、医药等领域,具有检测方便快捷、不破坏样品原状态等优点.在DLS颗粒测量技术中,粒度分布(Particle Size Distribution, PSD)的获取需要求解第一类Fredholm积分方程,该方程是典型的病态方程,即在测量过程的微小干扰或噪声都会导致粒度分布反演结果的严重失真,甚至造成反演失败.为得到准确的反演结果,很多反演算法被提出,包括累积分析法[1]、双指数法[2]、NNLS法[3]、CONTIN算法[4]、贝叶斯算法[5]、截断奇异值分解[6-7](Truncated Singular Value Decomposition,TSVD)和正则化方法[8]等,这些算法各有优点与不足,迄今为止,尚无一种算法具有绝对的优势,TSVD方法因其良好的抗噪性和重复性成为得到较多应用的算法之一.

1987年,HANSEN P C[6]采用TSVD方法求取粒度分布,对核矩阵进行奇异值分解,并将对扰动起放大作用的奇异值进行截断,进而求取粒度分布.2006年,ARIAS M L和FRONTINI G L[9]利用TSVD求取粒度分布,在此基础上添加了一个约束项保证了粒度分布的非负性.2010年,ZHU X[10]等采用非负最小二乘TSVD,对核矩阵进行奇异值分解,构造粒度分布的约束优化问题,同时根据粒度分布相对误差和解范数存在弱对称关系求取粒度分布.2013年,于向飞[11]等利用光强自相关函数(Autocorrelation Function, ACF)构造Hankel矩阵,根据分解后奇异值的分布确定噪音级别,重建矩阵的秩,通过对新的矩阵进行特征值分解得到颗粒的粒度分布.2017年,胡华[12]等利用奇异值分解对光能矩阵进行分解,发现灵敏度参数超过某个临界值时,会随着上边界的增大而急剧增大,推导出测量上限与仪器物理参数之间关系的解析表达式.

对于含噪声ACF,TSVD是以牺牲解的准确性来保证其抗噪性的,反演结果不可避免地会出现与真实粒度分布间的偏差[9],因而,难以获得准确的宽分布和双峰分布测量结果.本文对大、小两种颗粒组成的双峰及构成双峰的单峰颗粒体系的ACF衰减特性的量化分析表明,在双峰颗粒ACF衰减速度显著增加段,小颗粒对衰减起主要作用,在衰减速度显著减小段,大颗粒对衰减起主要作用.据此提出利用核矩阵与ACF不同延迟时刻的粒度信息构造扩展核矩阵的扩展截断奇异值分解(Extended TSVD,ETSVD)方法,通过用ACF每一延迟时刻的粒度信息,调节同一时刻原核矩阵数据对信噪比的贡献,从而在奇异值截断后保留更多的有效奇异值,减少了奇异值截断引起的信息丢失,进而减小反演偏差,获得更为准确的粒度分布结果.

1 动态光散射测量与TSVD反演

光强ACF与电场ACF满足Siegert关系

$ G_{{\theta _{\rm{r}}}}^{\left( 2 \right)}\left( {{\tau _j}} \right) = G_{\infty ,{\theta _{\rm{r}}}}^{\left( 2 \right)}\left( {{\tau _j}} \right)\left( {1 + \beta {{\left| {g_{{\theta _{\rm{r}}}}^{\left( 1 \right)}\left( {{\tau _j}} \right)} \right|}^2}} \right)\;\;\;\;\left( {{\theta _{\rm{r}}} = {\theta _{\rm{1}}}, \cdots ,{\theta _m},j = 1,2, \cdots ,M} \right) $ (1)

式中,θr为散射角,Gθr(2)(τj)为散射角θr处的光强ACF,τj为延迟时间,M为相关器的通道数,G∞, θr(2)(τj)为实验基线,β(β≤1)为相干因子,gθr(1)(τj)为电场ACF,m为散射角个数.

离散电场ACF为

$ g_{{\theta _{\rm{r}}}}^{\left( 1 \right)}\left( {{\tau _j}} \right) = {k_{{\theta _{\rm{r}}}}}\sum\nolimits_{i = 1}^{{N_{\rm{s}}}} {\exp \left( { - {\mathit{\Gamma }_i}{\tau _j}/{D_i}} \right){C_{{\theta _{\rm{r}}}}}\left( {{D_i}} \right)f\left( {{D_i}} \right)} $ (2)

式中,kθr为角度权重系数,可通过光强均值法[13]求取,衰减线宽Γi=DΓiq2,颗粒平移扩散系数DΓi=KBT/3πηDi,散射矢量的模值q=4πnsin (θr/2)/λKB为玻尔兹曼常数,T为绝对温度,η为分散介质粘度系数,n为介质的折射率,λ为真空中激光的波长,Ns表示粒度分布的离散点数,Cθr(Di)代表粒度为Di的颗粒在散射角θr的散射光强分数,可通过Mie理论计算,f(Di)表示待求粒度分布.

多角度电场ACF的矩阵形式为

$ \mathit{\boldsymbol{g}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{g}}_{{\theta _1}}}}\\ {{\mathit{\boldsymbol{g}}_{{\theta _2}}}}\\ \cdots \\ {{\mathit{\boldsymbol{g}}_{{\theta _{\rm{m}}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{A}}_{{\theta _1}}}}\\ {{\mathit{\boldsymbol{A}}_{{\theta _2}}}}\\ \cdots \\ {{\mathit{\boldsymbol{A}}_{{\theta _{\rm{m}}}}}} \end{array}} \right]\mathit{\boldsymbol{f}} = \mathit{\boldsymbol{Af}} $ (3)

式中,g=[gθ1; gθ2; …; gθm]为多角度电场ACF,A=[Aθ1; Aθ2; …; Aθm]为相应的核矩阵,$ {\mathit{\boldsymbol{A}}_{{\theta _{\rm{m}}}}} = {k_{{\theta _{\rm{m}}}}}\sum\nolimits_{i = 1}^{{N_{\rm{s}}}} {\exp \left( { - {\mathit{\Gamma }_i}{\tau _j}/{D_i}} \right)} {C_{{\theta _{\rm{m}}}}}\left( {{D_i}} \right)f\left( {{D_i}} \right) $f为待求粒度分布.

对于式(3),其最小二乘解可表示为

$ \mathit{\boldsymbol{f}} = \min {\left\| {\mathit{\boldsymbol{Af}} - \mathit{\boldsymbol{g}}} \right\|_2}\;\;\;\;\left( {0 \le \mathit{\boldsymbol{f}} \le 1} \right) $ (4)

A进行奇异值分解,得

$ \mathit{\boldsymbol{A}} = \mathit{\boldsymbol{U}}\sum {{\mathit{\boldsymbol{V}}^{\rm{T}}}} = \sum\nolimits_{i = 1}^n {{\sigma _i}{\mathit{\boldsymbol{\mu }}_i}\mathit{\boldsymbol{v}}_i^{\rm{T}}} $ (5)

式中,uivi分别为矩阵UV的列向量,σiA的奇异值,且σ1σ2…>σn>0,n=min(m×MNs).

对于含噪声ACF,即gθr-noise=g+egθr-noise为含噪声电场ACF,e为噪声.其粒度分布为

$ \mathit{\boldsymbol{f}} = \sum\nolimits_{i = 1}^n {\frac{{\mathit{\boldsymbol{\mu }}_i^{\rm{T}}\mathit{\boldsymbol{g}}}}{{{\sigma _i}}}{\mathit{\boldsymbol{v}}_i}} + \sum\nolimits_{i = 1}^n {\frac{{\mathit{\boldsymbol{\mu }}_i^{\rm{T}}\mathit{\boldsymbol{e}}}}{{{\sigma _i}}}{\mathit{\boldsymbol{v}}_i}} $ (6)

从式(6)第二项可以看出,小的奇异值会放大噪声e的影响,为减小这种奇异值的负面作用,通常对其进行截断,则截断正则矩阵和最小二乘的解为

$ \left\{ \begin{array}{l} {\mathit{\boldsymbol{A}}_p} = \sum\nolimits_{i = 1}^n {{\sigma _i}{\mathit{\boldsymbol{\mu }}_i}\mathit{\boldsymbol{v}}_i^{\rm{T}}} \\ {\mathit{\boldsymbol{f}}_p} = \sum\nolimits_{i = 1}^p {\frac{{\mathit{\boldsymbol{\mu }}_i^{\rm{T}}\mathit{\boldsymbol{g}}}}{{{\sigma _i}}}{\mathit{\boldsymbol{v}}_i}} + \sum\nolimits_{i = 1}^p {\frac{{\mathit{\boldsymbol{\mu }}_i^{\rm{T}}\mathit{\boldsymbol{e}}}}{{{\sigma _i}}}{\mathit{\boldsymbol{v}}_i}} \end{array} \right. $ (7)

式(7)中的有效奇异值(σ1, σ2, …, σp)由正则参数决定,正则参数由解范数‖fp2和残差范数‖gθr-noiseAfp2构造的L曲线[14]确定,大于正则参数的奇异值为有效奇异值,p为有效奇异值个数.

2 ACF的衰减特性分析与核矩阵扩展机理

在反演过程中,抑制噪声影响的同时提高粒度信息的利用率是准确获取粒度分布的关键[15],常规的TSVD方法均等地处理核矩阵数据,未对各延迟时刻ACF数据对粒度分布反演的贡献加以区分.为了分析ACF的衰减特性和不同延迟时刻的数据对粒度分布反演的贡献,以200 nm、900 nm和200/900 nm双峰颗粒体系为例,计算相应的ACF(图 1(a))、ACF的导数(图 1(b))及双峰颗粒ACF与等效平均粒度重构ACF的差值(图 1(c)).由图 1(b)可以看出,同一ACF在不同延迟时刻的衰减速度是不同的,呈现先增加后减小的趋势,不同ACF达到最大衰减速度的延迟时刻也不一样.为区分ACF在不同延迟时刻的衰减特性,按双峰颗粒衰减特性将延迟时间分为4段(图 1(b)),衰减速度无显著变化段(时段Ⅰ),衰减速度显著增加段(时段Ⅱ),衰减速度显著减小段(时段Ⅲ)和最小衰减结束段(时段Ⅳ).

图 1 200/900 nm颗粒体系衰减特性和粒度信息分布 Fig.1 Attenuation characteristic andparticle size information of 200/900 nm particles

图 1(b)不难看出,在时段Ⅰ,双峰及构成双峰的两个单峰颗粒对应的ACF衰减速度变化均很小,衰减速度依次为200 nm>200/900 nm>900 nm.在时段Ⅱ,三种颗粒ACF衰减速度显著增加,但仍保持200 nm>200/900 nm>900 nm的相对值,衰减的增加由小颗粒主导.在时段Ⅲ,ACF衰减速度变为900 nm>200/900 nm>200 nm,即衰减的增加由大颗粒主导.在时段Ⅳ,衰减速度均趋近于0.可见,小颗粒对双峰ACF的影响主要在时段Ⅱ,在该时段其ACF衰减速度较快,大颗粒对双峰ACF的影响则主要在时段Ⅲ,其ACF在此时段衰减速度较快.

为定量描述被测颗粒在不同延迟时刻对粒度分布反演的贡献,计算了粒度信息分布,即双峰颗粒ACF与其等效平均粒度(通过累积分析法[1]得到)重构的ACF的差值,见图 1(c).可以看出,反映双峰粒度特征的信息集中分布在ACF的时段Ⅱ和时段Ⅲ,而在时段Ⅰ和时段Ⅳ,粒度信息分布相对较少.文献[15]为了提高ACF数据的利用率,构建了一个对角矩阵作为正则化残差和的权重,其中对角元素采用ACF为底数、粒度信息分布为指数.从而改善了对ACF的信息利用,但未能解决TSVD方法中奇异值截断引起的信息丢失问题.

为了充分利用不同延迟时刻的数据对粒度分布的贡献差异,改善由于奇异值减少导致的常规TSVD方法反演准确性降低,本文在原待求方程gθr=Aθrf的基础上,利用每一角度核矩阵与对应角度下粒度信息在不同延迟时刻分布的结合,构建扩展矩阵Ar=[AθrIθr],从而,利用每一延迟时刻的粒度信息调节同一时刻原核矩阵数据对信噪比的贡献.其中,IθrM×1的向量, 代表粒度信息分布,ArM×(Ns+1)矩阵,表示由核矩阵和粒度信息分布构成的扩展核矩阵.则多角度电场ACF变为

$ \mathit{\boldsymbol{g}} = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{g}}_{{\theta _1}}}}\\ {{\mathit{\boldsymbol{g}}_{{\theta _2}}}}\\ \cdots \\ {{\mathit{\boldsymbol{g}}_{{\theta _{\rm{m}}}}}} \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{\mathit{\boldsymbol{A}}_{{\theta _1}}}{I_{{\theta _1}}}}\\ {{\mathit{\boldsymbol{A}}_{{\theta _2}}}{I_{{\theta _2}}}}\\ \cdots \\ {{\mathit{\boldsymbol{A}}_{{\theta _{\rm{m}}}}}{I_{{\theta _{\rm{m}}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} f\\ \delta \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} {{A_{I{\theta _1}}}}\\ {{A_{I{\theta _2}}}}\\ \cdots \\ {{A_{I{\theta _{\rm{m}}}}}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} f\\ \delta \end{array}} \right] = {A_I}\left[ {\begin{array}{*{20}{c}} f\\ \delta \end{array}} \right],\delta \to 0 $ (8)

式中,g为多角度电场ACF,f粒度分布,δ为一个无限趋于0的实数,AI为多角度扩展核矩阵.

对扩展核矩阵进行奇异值分解,再按照式(6)和式(7)即可实现核矩阵扩展后的TSVD反演——ETSVD反演.

为比较奇异值截断对TSVD法和ETSVD法保留奇异值的影响,本文在5种噪声水平下,对200/900 nm双峰颗粒体系的3角度(40°、60°和90°)含噪声ACF数据进行了基于正则参数的奇异值截断,正则参数值和截断后剩余的有效奇异值数量见表 1,其中,噪声水平为1×10-4和1×10-3下有效奇异值的分布及局部放大图见图 2.λ1λ2分别为TSVD方法和ETSVD方法的正则参数,正则参数由L曲线选取,n1n2分别为TSVD方法和ETSVD方法的有效奇异值个数.

表 1 200/900 nm颗粒在不同噪声水平下的正则参数及有效奇异值个数 Tab.1 Regularization parameters and numbers of effective singular values for 200/900 nm particles at different noise levels
图 2 噪声水平在1×10-4和1×10-3下的有效奇异值 Fig.2 Effective singular values at noise levels 1×10-4 and 1×10-3

表 1图 2可以看出,在不改变L曲线截断规则的前提下(即保证抗噪性的基础上),采用两种方法所得的较大有效奇异值没有显著差异,在接近正则参数值的小奇异值部分,ETSVD方法得到的奇异值略大于TSVD方法的奇异值,且保留了多于TSVD方法的有效奇异值.

3 数值模拟

本文采用半对数函数和两个半对数函数的组合形式分别模拟一组宽单峰S1(260 nm)和三组双峰颗粒体系S2(250/750 nm)、S3(270/800 nm)以及S4(306/974 nm)在6个角度(40°、50°、60°、70°、90°和110°)的光强ACF,同时采用TSVD方法及ETSVD方法对电场ACF进行反演,分析反演结果.半对数函数表达式[16]

$ f\left( {{D_i}} \right) = \frac{a}{{{D_i}{\sigma _1}\sqrt {2{\rm{ \mathsf{ π} }}} }}\exp \left[ { - \frac{{{{\left( {\ln \left( {{D_i}/{D_1}} \right)} \right)}^2}}}{{2\sigma _1^2}}} \right] $ (9)

式中,f(Di)是颗粒粒度分布,a为分布系数,D1为峰值处颗粒粒度,σ1为标准偏差,通过调整半对数函数参数可以获得相应的粒度分布.四组颗粒体系的参数如表 2所示.模拟实验条件满足:KB=1.380 7×10-23 J/K,分散介质折射率n=1.33,入射光在真空中的波长λ=632.8 nm,绝对温度为T=298.15 K,介质粘度系数η=0.89×10-9 g/nm·s.

表 2 4组颗粒的模拟参数 Tab.2 Simulation parameters of 4 groups particles

为与实测数据得到的ACF相吻合,通过式(10)添加随机噪声,即

$ G_{\theta {\rm{r}} - {\rm{noise}}}^{\left( 2 \right)}\left( {{\tau _j}} \right) = G_{\theta {\rm{r}}}^{\left( 2 \right)}\left( {{\tau _j}} \right) + \delta n\left( {{\tau _j}} \right) $ (10)

式中,Gθr()(τj)为无噪声的光强ACF,Gθr-noise()(τj)为含噪声的光强ACF,δ为噪声水平取1×10-3n(τj)为高斯随机噪声.

为量化比较两种方法的反演结果,引入峰值粒度误差E、峰值比R和分布误差V三个性能指标

$ \left\{ \begin{array}{l} E = \frac{{\left| {{D_{{\rm{sim}}}} - {D_{{\rm{inv}}}}} \right|}}{{{D_{{\rm{sim}}}}}}\\ R = {P_1}:{P_2}: \cdots {P_n}\\ V = {\left( {\sum\nolimits_{i = 1}^{{N_{\rm{s}}}} {{{\left[ {{f_{{\rm{sim}}}}\left( {{D_i}} \right) - {f_{{\rm{inv}}}}\left( {{D_i}} \right)} \right]}^2}} /\sum\nolimits_{i = 1}^{{N_{\rm{s}}}} {{{\left[ {{f_{{\rm{sim}}}}\left( {{D_i}} \right)} \right]}^2}} } \right)^{1/2}} \end{array} \right. $ (11)

式中,DsimDinv分别表示模拟分布和反演分布中峰值粒度,P表示峰值,fsimfinv分别表示模拟分布和反演分布.峰值粒度误差越小、模拟分布的峰值比与反演分布的峰值比越接近、分布误差越小,表示反演粒度分布与模拟粒度分布越接近.

图 3~7分别表示对单峰宽分布S1(260 nm)和三组双峰颗粒分布S2(250/750 nm)、S3(270/800 nm)以及S4(306/974 nm),采用1个角度(90°)、3个角度(40°、60°、90°)和6个角度(40°、50°、60°、70°、90°和110°)的TSVD方法和ETSVD方法的反演结果,表 3~6分别表示四组颗粒分布体系反演的性能参数.其中,'Sim PSD'表示模拟粒度分布,TSVD和ETSVD分别表示采用TSVD方法和ETSVD方法.

图 3 260 nm单峰颗粒采用TSVD方法和ETSVD方法的反演效果 Fig.3 PSDs of 260 nm unimodal particles using TSVD and ETSVD methods
图 4 250/750 nm双峰颗粒采用TSVD方法和ETSVD方法的反演效果 Fig.4 PSDs of 250/750 nm bimodal particles using TSVD and ETSVD methods
图 5 270/800 nm双峰颗粒采用TSVD方法和ETSVD方法的反演效果 Fig.5 PSDs of 270/800 nm bimodal particles using TSVD and ETSVD methods
图 6 306/974 nm双峰颗粒采用TSVD方法和ETSVD方法的反演效果 Fig.6 PSDs of 306/974 nm bimodal particles using TSVD and ETSVD methods
图 7 306/974 nm双峰颗粒采用TSVD方法和ETSVD方法的反演效果 Fig.7 PSDs of306/974 nm bimodal particles using TSVD and ETSVD methods
表 3 260 nm颗粒采用TSVD方法和ETSVD方法的性能参数 Tab.3 Performance parameters of 260 nm using TSVD and ETSVD methods
表 4 250/750 nm颗粒采用TSVD方法和ETSVD方法的性能参数 Tab.4 Performance parameters of 250/750 nm using TSVD and ETSVD methods
表 5 270/800 nm颗粒采用TSVD方法和ETSVD方法的性能参数 Tab.5 Performance parameters of 270/800 nm using TSVD and ETSVD methods
表 6 306/974 nm颗粒采用TSVD方法和ETSVD方法的性能参数 Tab.6 Performance parameters of 306/974 nm using TSVD and ETSVD methods

图 3表 3可以看出,对于260 nm宽分布颗粒体系,在单角度测量中,采用TSVD方法得到峰值粒度误差为0.076、分布误差为0.211的分布,采用ETSVD方法反演结果中峰值粒度误差减小到0.038、分布误差减小到0.198.无论是单角度测量还是多角度测量,与TSVD方法相比,ETSVD方法反演结果中的分布误差都有明显减小,特别是6个角度测量时,分布误差由0.409减小到0.175,峰值粒度误差由0.069减小到0.019.

对于250/750 nm双峰颗粒体系,从图 4表 4可以看出,在单角度测量中,采用TSVD方法虽能得到双峰分布,但峰值粒度误差达到0.496/0.052,分布误差为0.801,采用ETSVD方法所得分布中峰值粒度误差显著减小到0.056/0.013,分布误差减小到0.601.随着散射角度的增多,ETSVD方法的优势更为明显,当采用6个散射角时,其峰值粒度误差仅为0/0.006.

对于270/800 nm双峰颗粒体系,从图 5表 5可以看出,在单角度测量中,ETSVD方法避免了TSVD方法中“翘尾”现象,且无论是单角度还是多角度测量,与TSVD方法相比,采用ETSVD方法均能得到更为准确的粒度分布反演结果.

图 6表 6可以看出,对于306/974 nm双峰颗粒体系,随散射角个数的增加,采用TSVD方法和ETSVD方法反演得到的分布误差均有所减小,且ETSVD方法的分布误差均小于TSVD方法,同时,峰值粒度误差也逐渐减小,在6个角度测量时,采用ETSVD方法的峰值粒度误差仅为0.013/0.004.

4 实验验证

为了进一步验证ETSVD方法的反演性能,采用TSVD方法和ETSVD方法对峰值粒度为306/974 nm、峰值比为1:1的标准聚苯乙烯乳胶颗粒在1个角度(90°)、3个角度(40°、60°、90°)和6个角度(40°、50°、60°、70°、90°和110°)的ACF数据进行了反演.实验测量装置采用BI-2030AT相关器、He-Ne激光器和步进电机控制的测角仪.反演结果和性能参数见图 7表 7.

表 7 306/974 nm颗粒采用TSVD方法和ETSVD方法的性能参数 Tab.7 Performance parameters of 306/974 nm using TSVD and ETSVD methods

图 7表 7可以看出,在单角度测量中,采用TSVD方法得到的峰值粒度误差为0.091/0.315、峰值比为1: 3.162,采用ETSVD方法反演的峰值粒度误差为0.094/0.030,峰值比为1: 1.757.随着测量角度的增加,峰值粒度误差逐渐减小,与TSVD方法相比,ETSVD方法得到的峰值粒度和峰值比更趋接近实际情况,在6个散射角测量情况下,采用ETSVD得到了接近峰值粒度306/974 nm和峰值比1: 1的分布.

5 结论

为解决奇异值截断造成的粒度分布反演偏差,通过量化分析粒度信息在ACF不同延迟时段分布特性发现:双峰及构成双峰的两个单峰颗粒对应的ACF衰减速度变化在衰减速度无显著变化段均很小;在衰减速度显著增加段,三种颗粒ACF衰减速度显著增加,衰减的增加由小颗粒主导;在衰减速度显著减小段,衰减的增加由大颗粒主导;在最小衰减结束段,衰减速度均趋近于0.本文据此提出利用核矩阵与ACF不同延迟时刻上的粒度信息构建扩展矩阵的ETSVD方法,该方法以扩展矩阵作为电场ACF的新核矩阵,通过用ACF每一延迟时刻的粒度信息,调节同一时刻原核矩阵数据对信噪比的贡献,保留了更多的有效奇异值,减少了由于奇异值截断引起的信息丢失,在保证抗噪性的基础上,提高了ACF的信息利用率.在相同噪声水平下,对单峰宽分布和双峰分布颗粒体系的单角度和多角度动态光散射模拟数据反演表明,与TSVD方法相比,采用ETSVD方法反演得到的峰值粒度误差和分布误差均明显小,对双峰颗粒实测数据的反演,得到相同的结论.

宽分布和双峰分布颗粒体系的准确测量受制于测量噪声和反演算法的信息利用率,由于噪声是DLS测量中无法完全去除的因素,在不降低抗噪性的前提下提高测量准确性,是提高和改进粒度反演方法的主要目标之一,进行核矩阵扩展进而提出相应的改进反演方法,提高了宽分布和双峰分布颗粒体系的测量准确性,也将为多峰和复杂分布颗粒体系的测量提供新途径.

致谢 本文所用双峰聚苯乙烯乳胶颗粒实验数据由阿根廷Jorge R. Vega教授课题组提供,在此致谢.

参考文献
[1] LIU Wei, WANG Ya-jing, SHEN Jin. Optimal fitting cumulants method for dynamic light scattering[J]. Acta Optica Sinica, 2013, 33(12): 1229001.
刘伟, 王雅静, 申晋. 动态光散射最优拟合累积分析法[J]. 光学学报, 2013, 33(12): 1229001.
[2] INANOV D A, WINKELMANNJ. Multiexponential decay autocorrelation function in dynamic light scattering in near-critical ternary liquid mixture[J]. Journal of Chemical Physics, 2006, 125(10): 104507. DOI:10.1063/1.2338312
[3] ANSARI R R, NYEO S L. Submicron particle size distributions by dynamic light scattering with non-negative least-squares algorithm[J]. Chinese Journal of Physics, 2012, 50(3): 459-477.
[4] ZHU X, SHEN J, THOMAS J C, et al. Analysis of noisy dynamic light scattering data using constrained regularization techniques[J]. Applied Optics, 2012, 51(31): 7537-7548. DOI:10.1364/AO.51.007537
[5] CLEMENTI L A, VEGA J R, GUGLIOTTA L M, et al. A bayesian inversion method for estimating the particle size distribution of latexes from multiangle dynamic light scattering measurements[J]. Chemometrics & Intelligent Laboratory Systems, 2011, 107(1): 165-173.
[6] HANSEN P C. The truncated SVD as a method for regularization[J]. Bit Numerical Mathematics, 1987, 27(4): 534-553. DOI:10.1007/BF01937276
[7] DOU Zhen-hai, WANG Ya-jing, SHEN Jin, et al. A hybrid non-negative inversion of dynamic light scattering based on truncated singular value decomposition[J]. Chinese Journal of Lasers, 2013, 40(6): 0608001.
窦震海, 王雅静, 申晋, 等. 动态光散射混合非负截断奇异值反演[J]. 中国激光, 2013, 40(6): 0608001.
[8] XU Min, SHEN Jin, ZHU Xin-jun, et al. Recovery of bimodal particle size distributions with multiangle dynamic light scattering[J]. Acta Photonica Sinica, 2017, 46(2): 0229001.
徐敏, 申晋, 朱新军, 等. 双峰分布颗粒体系的多角度动态光散射数据反演[J]. 光子学报, 2017, 46(2): 0229001.
[9] ARIAS M L, FRONTINI G L. Particle size distribution retrieval from elastic light scattering measurements by a modified regularization method[J]. Particle & Particle Systems Characterization, 2006, 23(5): 374-380.
[10] ZHU X, SHEN J, LIU W, et al. Nonnegative least-squares truncated singular value decomposition to particle size distribution inversion from dynamic light scattering data[J]. Applied Optics, 2010, 49(34): 6591-6596. DOI:10.1364/AO.49.006591
[11] YU Xiang-fei, YANG Hui, YANG Hai-ma, et al. An inversion algorithm of dynamic light scattering based on singular value decomposition[J]. Acta Photonica Sinica, 2013, 42(11): 1324-1328.
于向飞, 杨晖, 杨海马, 等. 基于奇异值分解的动态光散射反演算法[J]. 光子学报, 2013, 42(11): 1324-1328.
[12] HU Hua, ZHANG Fu-gen, LV Qie-ni, et al. Study on the upper limit of laser particle size analyzer[J]. Acta Optica Sinica, 2018, 38(04): 0429001.
胡华, 张福根, 吕且妮, 等. 激光粒度仪测量上限研究[J]. 光学学报, 2018, 38(04): 0429001.
[13] GAO Ming-liang, WANG Xue-min, SHEN Jin, et al. Influence path of angular error on multiangle dynamic light scattering measurement[J]. Acta Photonica Sinica, 2017, 46(10): 1029002.
高明亮, 王雪敏, 申晋, 等. 多角度动态光散射角度误差影响测量的途径分析[J]. 光子学报, 2017, 46(10): 1029002.
[14] LIU X, SHEN J, THOMAS J C, et al. Multiangle dynamic light scattering analysis using angular intensity weighting determined by iterative recursion[J]. Applied Optics, 2012, 51(7): 846-854. DOI:10.1364/AO.51.000846
[15] XU M, SHEN J, THOMAS J C, et al. Information-weighted constrained regularization for particle size distribution recovery in multiangle dynamic light scattering[J]. Optics Express, 2018, 26(1): 15-31. DOI:10.1364/OE.26.000015
[16] BERMEO L A, CAICEDO E, CLEMENTI L, et al. Estimation of the particle size distribution of colloids from multiangle dynamic light scattering measurements with particle swarm optimization[J]. Ingeniería E Investigación, 2015, 35(1): 49-54. DOI:10.15446/ing.investig.v35n1.45213