馬光耀，曹茂永,，馬鳳英，紀 鵬

(1.山東科技大學 電氣與自動化工程學院，山東 青島 266590；2.齊魯工業(yè)大學 電氣工程與自動化學院，山東 濟南 250353)

顆粒粒徑及其分布情況是衡量顆粒產品優(yōu)良性最為重要的參數(shù)，粒徑的大小及分布情況直接決定了由顆粒生產所得產品的性能以及環(huán)境污染程度的判別[1-2]。動態(tài)光散射技術(dynamic light scattering,DLS)具有非接觸、測量速度快等諸多特性，在顆粒粒度分布(particle size distribution,PSD)檢測方面得到了廣泛應用[3-4]。然而，此方法在檢測過程中存在比較大的缺陷：當測量數(shù)據含有微小擾動或噪聲時，會導致最終檢測結果與真實結果有較大差距。所以在運用該方法時需要盡可能地去除噪聲的影響，將病態(tài)不適定方程問題轉換成一個與之近似的適定方程問題求解[5-7]。

目前，已有很多關于動態(tài)光散射反問題的優(yōu)秀研究方法被提出，最為常用且擁有較好反演精度和抗噪能力的是Tikhonov正則化[8]和截斷奇異值(truncated singular value decomposition,TSVD)正則化[9]方法。其中，Tikhonov正則化法是通過引入正則化參數(shù)和穩(wěn)定泛函來設定濾子函數(shù)，從而改善矩陣的病態(tài)性；截斷奇異值正則化法是通過設定濾子函數(shù)把容易造成觀測數(shù)據噪聲放大的小奇異值進行截斷處理[10-13]，從而將不適定的病態(tài)問題轉換成適定問題并得到結果。為了讓Tikhonov和TSVD算法反演得到更好的效果，文獻[10]將Morozov偏差確定正則參數(shù)的迭代原理運用到DLS含噪數(shù)據反演中，在一定程度上優(yōu)化了低噪聲單峰顆粒粒度的反演精度，但是該迭代算法的終止參數(shù)卻難以確定；林東方等[11]、郭家橋等[12]通過對奇異值矩陣的小奇異值部分進行修正減小了一定的正則化誤差引入量，雖然該方法引入的是小奇異修正值，但引入量中仍存在部分未能及時修正的奇異值，導致結果出現(xiàn)偏差；王雅靜等[14]提出以小波正則化與對反演問題進行分解，用Tikhonov正則化、TSVD正則化分別對粗尺度、細尺度進行尋優(yōu)并處理，得到了具有較好抗噪性以及準確性的結果，但該算法計算量大且尋優(yōu)過程中需保證非負性，難度較大。在模擬實驗過程中，發(fā)現(xiàn)Tikhonov相對于TSVD而言，反演結果的平滑性更好；而TSVD的抗噪性、精確性更好一些[15]。

本研究在Tikhonov正則化法和TSVD正則化法的基礎上，通過對兩者濾子函數(shù)的研究[16-17]，提出一種將TSVD和Tikohonv相結合的新正則化方法TTSVD，并對動態(tài)光散射含噪數(shù)據進行模擬反演實驗。通過將反演結果與真實值、Tikhonov算法結果、TSVD算法結果進行比對，發(fā)現(xiàn)新算法有較好的精度和抗噪性。

1 動態(tài)光散射及Tikhonov、TSVD反演原理

1.1 動態(tài)光散射原理

在動態(tài)光散射測量技術中，粒度分布函數(shù)是根據歸一化光強自相關函數(shù)[18]獲得的：

(1)

其中，τ是延遲時間，Γ為Rayleigh線寬，又稱衰減率。若想求PSD結果，需要在求得衰減線寬分布函數(shù)G(Γ)之后，根據平移擴散系數(shù)公式和Stokes-Einstein關系進行求解。為求解衰減線寬函數(shù)G(Γ)，將式(1)進行離散化得：

(2)

式中：j表示相關器的采樣通道序列，i是樣品顆粒粒徑分級數(shù)。式(2)可以用矩陣形式進行簡寫：

Ax=b，

(3)

其中：xi=G(Γi)，bj=g(τj)；矩陣A的元素ai,j=exp(-Γiτj)。對方程(3)進行求解即可得到PSD結果。

1.2 Tikhonov正則化和TSVD反演理論

DLS反演問題即是對病態(tài)不適定方程(3)的求解問題，理論上可由最小二乘的求解問題來替換：

(4)

假設矩陣A∈Rm×n(m≥n)，根據奇異值分解理論，對矩陣A進行奇異值分解得：

(5)

式中，U=(u1,u2,…,um)∈Rm×m，V=(v1,v2,…,vn)∈Rn×n，Σ=diag(σ1,σ2,…,σn)∈Rn×n，σ1≥σ2≥…≥σn為奇異值，ui、vi分別是矩陣A的左、右奇異向量。

所以，理論上的最小二乘解為：

(6)

由于奇異值矩陣中含有非常小的奇異值，當測量數(shù)據b包含微小噪聲或受到輕微擾動時，xLS就會產生非常大的隨機波動，求取的結果會非常不穩(wěn)定[15]。

Tikhonov正則化是通過引入展平泛函將不適定問題近似轉換程適定問題，并進行近似求解：

(7)

式中，L是正則算子，x0為解的初始估計，λ為正則化參數(shù)。當x0=0，L取單位矩陣時，Tikhonov正則化可表示成：

(8)

在選取適當?shù)恼齽t參數(shù)后，即可求得正則解：xTikhonov=min{Mλ[x,b]}。

TSVD正則化是通過將奇異值矩陣中使得數(shù)據中噪聲進行放大的小奇異值部分進行截斷，從而消除觀測數(shù)據中噪聲對結果的影響。為了避免小奇異值對正則解的影響，當任意k

(9)

2 新正則化算法TTSVD

根據Lanczos奇異值分解理論可知，在方程病態(tài)時，奇異值σn是一個非常小的近似為0的值，且奇異值σ1遠大于σn。矩陣A中較大的奇異值及其相應的左右奇異向量構成了其數(shù)值模型中的可靠部分，而小奇異值及其相應的左右奇異向量構成了數(shù)值模型中的不可靠部分。矩陣中的不可靠部分會將式(6)中數(shù)據b所包含的噪聲進行放大，導致結果出現(xiàn)較大偏差。Tikhonov正則化通過引入濾子函數(shù)對矩陣的奇異值進行修正，在修正不可靠部分的小奇異值過程中也不可避免地將可靠部分的大奇異值進行變換，導致反演結果對噪聲數(shù)據較為敏感，穩(wěn)定性較差；TSVD正則化通過將矩陣不可靠部分也就是奇異值矩陣中的小奇異值部分進行完全截斷，在消除矩陣不可靠部分中對噪聲放大部分的同時造成不可靠部分中所含信息的缺失，降低了解的準確性[14-15]。

在對Tikhonov正則化和TSVD正則化兩種算法進行綜合考慮后，結合這兩種算法的優(yōu)點，通過設定不同的閾值將矩陣按照大奇異值、小奇異值以及剩余奇異值分成非?？煽坎糠?、不確定部分以及非常不可靠部分，讓濾子函數(shù)僅對中間不確定部分進行修正調整，對非?？煽坎糠植蛔鲎兓?，并將非常不可靠部分進行完全截斷，以期實現(xiàn)在保留不可靠部分細節(jié)性信息的同時降低其中非常不可靠成分對噪聲的放大作用。

Tikhonov正則化算法的濾子函數(shù)是:

(10)

根據濾子函數(shù)式(10)及選取的正則參數(shù)，可求得解為：

(11)

TSVD正則化算法的濾子函數(shù)是：

(12)

根據濾子函數(shù)式(11)及其選取的截斷正則參數(shù)，可求得解為：

(13)

由式(11)與(13)可以看出，Tikhonov正則化會將所有矩陣的所有奇異值均進行變化，在保留小奇異修正值的同時也將大奇異值進行了變化；而TSVD正則化是直接將其所認定的小奇異值進行了完全截斷，導致數(shù)值模型的部分信息缺失。

改進的正則化算法TTSVD的濾子函數(shù)是：

(14)

式中，k代表非?？煽坎糠值脑O定閾值，其值為截斷奇異值的序列號；λ為選取的正則化參數(shù)，t代表的是不確定部分的設定閾值，其值是根據正則化參數(shù)λ和截斷奇異值序列號k來確定的。當正則化參數(shù)為λ時，奇異值序列中必定存在一個與之相近的奇異值σλ1滿足以下條件：

(15)

t是以奇異值σλ1的序列號λ1為中心，從λ1向右擴充截斷參數(shù)k后的奇異值總數(shù)，即：

t=λ1+k。

(16)

根據濾子函數(shù)式(14)及其選取的參數(shù)k和t，可求得解為：

(17)

在以上三種求解算法中，均需要確定正則化參數(shù)，本研究采用應用廣泛的L曲線法進行選取。

3 模擬實驗分析

3.1 模擬實驗及參數(shù)設定

為了更全面地展現(xiàn)并研究TTSVD算法對于動態(tài)光散射噪聲數(shù)據反演結果的影響，本研究通過建立模擬動態(tài)光散射含噪數(shù)據來代替實際觀測數(shù)據進行全面分析。在實際顆粒粒度檢測過程中，樣品的雜質與濃度情況、實驗溫度的變化、光電倍增管的背景雜散光等諸多因素的存在，會使得實際觀測數(shù)據中包含大量由不同因素所導致的隨機噪聲。依據中心極限定理思想，通過以正態(tài)分布白噪聲對模擬的無噪數(shù)據進行激勵，讓模擬產生的觀測數(shù)據能夠達到與真實結果較為接近的效果。具體模擬實現(xiàn)過程：首先通過模型分別建立單峰窄分布顆粒系、單峰寬分布顆粒系、雙峰窄分布顆粒系、雙峰寬分布顆粒系四種不同的粒度分布函數(shù)作為模擬的真實顆粒粒度分布結果，然后通過公式分別得到其自相關函數(shù)數(shù)據，也就是在顆粒粒度檢測過程中所得到的不含噪聲污染的理論觀測數(shù)據，最后在相關函數(shù)數(shù)據中分別加入不同強度水平的正態(tài)分布白噪聲，得到噪聲干擾下的實際觀測數(shù)據來進行研究。

模擬實驗中，使用Johnson’s SB[19]分布作為模擬粒徑分布模型(即真實粒度分布)，其表達式是：

(18)

式中，t=(d-dmin)/(dmax-dmin)是粒徑大小歸一化的結果，dmax和dmin分別是模擬粒度分布中粒子的最大粒徑和最小粒徑；u和δ表示模型的分布參數(shù)，這兩個參數(shù)值將決定模擬的粒度分布情況。模擬實驗參數(shù)分別為：溫度25 ℃，分散介質(水)的折射率1.331，入射光的波長632.8 nm，散射角90°，水的黏度系數(shù)0.89×10-3N·s/m2，波爾茲曼常數(shù)1.380 7×10-23J/K。通過對自相關函數(shù)g(τ)添加不同水平的正態(tài)分布白噪聲，來模擬在實際測量過程中可能出現(xiàn)的不同程度的干擾情況。

為了能夠更直觀地展示算法對模擬反演結果的影響，將新算法所得的結果與傳統(tǒng)的Tikhonov正則化、TSVD正則化結果進行對比，同時引入峰值粒徑(peak particle size,PPS)、峰值相對誤差(peak particle size relative error,PPSRE)、粒度分布誤差(particle size distribution error,PSDE)三個重要指標參數(shù)并以數(shù)據形式展現(xiàn)出來。其中峰值粒徑指的是粒度分布峰值點所對應的粒徑值，峰值相對誤差指的是算法求取峰值與模擬峰值的相對誤差大小

RPPSRE=(|dpps_inv-dpps_true|/dpps_true)×100%。

(19)

其中，dpps_inv和dpps_true分別代表算法峰值粒徑和模擬峰值粒徑。粒度分布誤差指的是算法求取的粒度分布與模擬的粒度分布之間的誤差：

(20)

其中，xpps_inv和xpps_true分別代表算法求取的粒度分布和模擬的粒度分布。

3.2 結果對比與數(shù)據分析

3.2.1 單峰顆粒系的反演結果與數(shù)據對比

單峰窄分布顆粒系對應的模擬參數(shù)為：u=3.8，sigma=2.0，dmin=1，dmax=800；單峰寬分布對應的模擬參數(shù)為：u=0.8，sigma=1.8，dmin=1，dmax=800。兩種單峰分布顆粒系的的反演結果和數(shù)據如圖1～2以及表1～2所示。

圖1 單峰窄分布顆粒系在不同噪聲強度下的反演結果Fig.1 Inversion results of unimodal narrow distribution particles under different noise intensities

由圖1(a)與表1數(shù)據所示，對于在噪聲強度為0.000 1時的單峰窄分布顆粒系，TTSVD算法和TSVD算法均能反演出與真實結果相同的PPS，比Tikhonov算法的結果更精確；同時，TTSVD結果的PSDE比Tikhonov的低0.001 2，比TSVD低0.000 1。這也就是說，在低噪聲的情況下TTSVD所選取的可靠部分加上不確定部分已經近似替代了原矩陣，使得TTSVD反演結果與真實結果之間的擬合程度更優(yōu)于其他兩種算法；而非常不可靠部分的截斷使得修正后的矩陣獲得比Tikhonov更精確的結果。同樣的，對于單峰寬分布顆粒系也印證了這一解釋。如圖2(a)所示，在噪聲為0.000 1時的單寬峰顆粒系，TTSVD的PPSRE比Tikhonov低3.35%、比TSVD低0.67%；在PSDE方面，TTSVD結果也表現(xiàn)出與單窄峰顆粒系相類似的情況。這也說明對選取的不確定部分進行修正以及對非?？煽坎糠直３植蛔兊奶幚硎欠浅Ｓ斜匾模瑑?yōu)化了反演的精度以及平滑性。

圖2 單峰寬分布顆粒系在不同噪聲強度下的反演結果Fig.2 Inversion results of unimodal wide distribution particles under different noise intensities

表1 單峰窄分布顆粒系在不同噪聲強度下的反演數(shù)據Tab.1 Inversion data of unimodal narrow distribution particles under different noise intensities

表2 單峰寬分布顆粒系在不同噪聲強度下的反演數(shù)據Tab.2 Inversion data of unimodal wide distribution particles under different noise intensities

如圖1和圖2(b)、2(c)所示，隨著噪聲強度不斷增強，三種算法的結果均出現(xiàn)了不同程度的變化，其中以Tikhonov算法受噪聲的干擾最為嚴重，而TTSVD和TSVD算法受噪聲干擾的影響相對較小。相較于Tikhonov算法，隨著噪聲的不斷增強，TTSVD的PPSRE相較真實值變化不大，而Tikhonov算法的PPSRE也隨著噪聲強度增強而變大，在噪聲強度為0.010 0時，兩種算法的單寬峰顆粒系峰值粒徑相對誤差達到4.36%，單窄峰顆粒系峰值粒徑相對誤差達到5.62%。，Tikhonov算法已經不能得到準確的PPS。這說明，相較于將奇異值完全修正的Tikhonov算法，TTSVD的部分截斷與部分修正是完全有必要的。相較于TSVD算法，TTSVD和TSVD的PPSRE相差不大，且在單寬峰分布顆粒系TTSVD結果的PPSRE要優(yōu)于TSVD結果，在噪聲強度為0.010 0時，單寬峰顆粒系TTSVD算法的PPSRE要比TSVD的低5%；與此同時，不管噪聲強度如何變化，TTSVD算法的PSDE始終保持著比TSVD算法更低的值。這是由于TTSVD算法在截斷非常不可靠部分的小奇異值后，相較于Tikhonov算法受噪聲的擾動影響變??；確定部分的完全保留與不確定部分的修正，使得TTSVD算法能夠具有比Tikhonov和TSVD算法更低的PSDE，提高了Tikhonov、TSVD算法與真實結果的擬合度。

綜上所述，改進的TTSVD算法能夠在單峰分布顆粒系的反演中得到比較精確且穩(wěn)定的結果。相較于Tikhonov算法，TTSVD具有更好的抗噪性與精確性，特別是噪聲強度為0.010 0時，峰值誤差至少降低了3%，粒度分布結果也更接近真實粒度分布情況；相較于TSVD算法，TTSVD算法反演結果具有更高的擬合度，在低噪聲的情況下具有更好的精確性。

3.2.2 雙峰顆粒系的反演結果與數(shù)據對比

本實驗中雙峰顆粒系是由兩種單峰顆粒系按照1∶1的光強比進行混合而成，雙峰窄分布對應的模擬系數(shù)為：u1=4.9，sigma1=3.1，u2=-3.9，sigma2=4.0，dmin=1，dmax=800；雙峰寬分布顆粒系對應的模擬系數(shù)為：u1=3.0，sigma1=2.1，u2=-2.1，sigma2=2.2，dmin=1，dmax=800。兩種由不同顆粒系混合而成的雙峰分布顆粒系反演結果和數(shù)據如圖3～4以及表3～4所示。

圖3 雙峰窄分布顆粒系在不同噪聲強度下的反演結果Fig.3 Inversion results of bimodal narrow distribution particles under different noise intensities

表3 雙峰窄分布顆粒系在不同噪聲強度下的反演數(shù)據Tab.3 Inversion data of bimodal narrow distribution particles under different noise intensities

如圖3和圖4(a)所示，在噪聲強度為0.000 1時，三種算法均能得到比較精準的雙峰顆粒系的峰值粒徑，其中又以TTSVD算法結果的PPSRE最小，擬合效果最好。由表3數(shù)據可知，在雙窄峰顆粒系中三種算法得到的第一個峰值粒徑完全一樣，在第二個峰值上TTSVD和TSVD算法比Tikhonov算法的PPSRE低0.85%；根據表4數(shù)據，TTSVD算法在雙寬峰顆粒系的PPSRE比TSVD算法低0.5%左右，比Tikhonov算法低1.5%～5.8%。這說明，在外部噪聲較小的情況下，TTSVD算法和TSVD算法對非常可靠部分奇異值的保留是很有必要的，提高了反演的精度。從粒度分布反演誤差PSDE的角度以及圖像雙峰的分辨率來看，寬峰顆粒系和窄峰顆粒系的規(guī)律基本一致，TTSVD算法結果的粒度分布誤差均低于TSVD和Tikhonov算法，與真實結果的擬合度較高，能夠較好地反映真實結果的粒度分布情況。

圖4 雙峰寬分布顆粒系在不同噪聲強度下的反演結果Fig.4 Inversion results of bimodal wide distribution particles under different noise intensities

表4 雙峰寬分布顆粒系在不同噪聲強度下的反演數(shù)據Tab.4 Inversion data of bimodal wide distribution particles under different noise intensities

如圖3和圖4(b)、4(c)所示，隨著噪聲強度的不斷增強，三種算法對顆粒粒度反演的結果不斷變差，特別是在噪聲強度為0.010 0時，三種算法已經基本失去了第二種顆粒系的分布情況。相比較而言，TTSVD算法具有比Tikhonov算法和TSVD算法更好的結果。如圖3(b)所示，在噪聲強度為0.001 0時，TSVD算法反演結果在第一個顆粒系上的峰值準確性要優(yōu)于TTSVD和Tikhonov算法，但其同時也丟失掉第二顆粒系信息；而TTSVD算法能夠同Tikhonov算法一樣獲取第二個峰的顆粒信息，且得到了比Tikhonov算法在第一峰值PPSRE上低1.54%、第二峰值上低5.63%的反演結果。與之相類似的是，在圖4(b)中，Tikhonov算法在0.001 0噪聲時取得了較噪聲強度0.000 1時更好的結果。出現(xiàn)這種情況均是因為施加噪聲隨機性較大以及L曲線選取的參數(shù)未能達到最優(yōu)從而導致Tikhonov和TTSVD算法反演結果受到影響，但與Tikhonov算法不同的是，TTSVD對非常不可靠部分奇異值的截斷處理使得該算法能夠取得比Tikhonov算法更好的結果，進一步說明TTSVD算法對非常不確定部分奇異值的截斷處理是很有必要的。如圖3和圖4(c)所示，在噪聲強度為0.010 0時，TTSVD算法得到的第一個峰值顆粒系結果比Tikhonov和TSVD算法更準確，同時在雙窄峰顆粒系上，能夠分辨出第二個顆粒系的存在情況。造成這種結果的原因是TSVD將不確定部分與不可靠部分的奇異值進行了截斷，使得在反演中缺失部分有效數(shù)據，隨著噪聲強度增加截斷的奇異值就越多，所保留的有效信息也就越少，最終導致反演結果在高噪聲情況下失去第二種顆粒系的粒度分布情況；Tikhonov正則化算法雖然保留了大量有效信息，但是其濾子函數(shù)對矩陣數(shù)學模型中大小奇異值的處理不充分，使得小奇異值對噪聲變化非常敏感，大奇異值的變化使得結果對可用信息的處理效率降低。而TTSVD算法將小奇異值部分進行劃分并分別處理，在盡可能保留有效信息的同時，降低了小奇異值對噪聲的放大作用，在一定程度上對Tikhonov和TSVD算法進行了優(yōu)化。

綜合上述分析得出：改進的TTSVD算法具有比TSVD算法、Tikhonov算法與原顆粒粒度分布更高的擬合度，對雙峰顆粒系的反演能夠得到比TSVD算法與Tikhonov算法更好的結果。TTSVD算法在低噪聲的情況下反演得到的結果精確度要高于TSVD算法與Tikhonov算法；在高噪聲的情況下，具有比Tikhonov算法更強的抗噪性、比TSVD算法更好的雙峰分辨率。

4 實驗數(shù)據的反演分析

為驗證該算法在顆粒粒度檢測應用中的有效性，使用DLS實驗裝置獲取nm單峰分布的標準聚苯乙烯顆粒相關函數(shù)數(shù)據并進行反演驗證。DLS裝置主要由波長為532 nm的綠色激光器、端窗式光電倍增管以及128通道的高速數(shù)字相關器組成，測試實驗溫度為27℃，散射角為90°，分散介質選擇純凈水。根據動態(tài)光散射技術理論可知，該技術通過采用相關去噪法對獲取的光強信號進行去噪，去噪效果與測量時間成正相關，測量時間越長去噪效果越好。而測量時間是采樣時間與數(shù)據長度的乘積，這就要求實驗檢測時的數(shù)據量越大越好。為了得到不同程度噪聲情況下的觀測數(shù)據，在采樣數(shù)約為2×106和2×105時分別獲取相關函數(shù)數(shù)據作為低噪聲與高噪聲下的觀測數(shù)據，并以TSVD、Tikhonov和TTSVD算法分別對實驗數(shù)據進行反演，結果和數(shù)據如圖5和表5所示。

由圖5和表5可以看出，在采樣數(shù)量為2×106時，三種算法均能比較準確地反演出顆粒的峰值粒徑情況，其中Tikhonov算法結果的峰值誤差為0.41%，粒度分布較寬；TTSVD和TSVD算法結果的峰值誤差為0.21%，粒度分布較窄。相較于采樣數(shù)量為2×106時的反演結果，在采樣數(shù)量為2×105時，TTSVD和TSVD算法仍能夠比較準確地反演出顆粒粒徑峰值結果，其中TTSVD算法結果的峰值誤差保持不變，TSVD算法結果的峰值誤差增大到0.47%，而Tikhonov算法結果的粒徑峰值出現(xiàn)較大變化，峰值誤差增大5.65%。因此，TTSVD算法能夠更好反演出實際的顆粒粒度分布結果，驗證了模擬數(shù)據的結論。

圖5 350 nm顆粒在不同采樣數(shù)量下的反演結果Fig.5 Inversion results of 350 nm particles under different sampling numbers

表5 350 nm顆粒在不同采樣數(shù)量下的反演數(shù)據Tab.5 Inversion data of 350 nm particles under different sampling numbers

5 結論

在動態(tài)光散射數(shù)據的反演實驗的基礎上，結合Tikhonov正則化算法和TSVD正則化算法，提出一種將奇異值序列進行區(qū)域劃分并分別處理的改進正則化方法TTSVD。將奇異值序列劃分為非?？煽坎糠?、不確定部分以及非常不確定三部分，對三部分以不同的濾子函數(shù)進行處理，解決了Tikhonov算法對所有奇異值進行修正以及TSVD算法將所有小奇異值進行完全截斷的缺陷。通過將原始顆粒粒度分布作為基準，以改進的正則化方法TTSVD和Tikhonov、TSVD算法分別在強度為0.000 1、0.001 0、0.010 0下的DLS數(shù)據反演進行結果比對，確定新算法結果能夠有較強的抗噪性、較好的精確性以及與真實結果之間更好的擬合性。當然由于計算存在隨機性，可能偶爾會出現(xiàn)其他方法某個參數(shù)更好的情況，但總體而言，TTSVD方法是一種可行的顆粒粒度反演方法。