熒光光譜無用散射的自動去除算法

許佩軍, 周 末, 張明振, 郭秋含, 牟曉紅, 李 焱

(1.遼寧師范大學 物理與電子技術學院,遼寧 大連 116029;2.中國科學院 大連化學物理研究所 分子反應動力學國家重點實驗室,遼寧 大連 116023)

三維熒光光譜是近些年廣泛應用而發展起來的一種熒光分析技術, 它在食品科學[1-2]、分析化學[3-6]、生物化學[7-8]以及環境科學[9-10]等領域都有廣泛的應用. 三維熒光光譜是由激發波長(Emission, Em), 發射波長(Excitation, Ex), 熒光強度(熒光強度數值用Z表示)組成的三維矩陣光譜(Excitation-Emission Matrix, EEM), 它描述了熒光強度隨激發波長和發射波長的變化關系. 三維熒光光譜不僅能夠獲得激發波長與發射波長, 同時能夠獲取變化時的熒光強度信息. 在熒光光譜的生成過程中, 由散射機制引起的拉曼(Raman)散射或瑞利(Rayleigh)散射, 經常會單獨存在或同時存在于熒光光譜數據中. 這兩種散射的能量通常要高于目標物質對應的能量. 而且, 這兩種散射的能量范圍有可能干擾甚至部分覆蓋目標物質的能量, 于是, 去除熒光光譜中的無用散射就具有重要的理論研究意義.

無用散射的去除工作, 一直是熒光光譜分析工作中的重要步驟[11-13]. 2014年, Eilers和Kroonenberg通過高斯峰刻畫了散射的位置信息、高度以及寬度. 在平滑處理后, 實現拉曼和瑞利散射的去除[14]. 2015年, Rinnan和Anderson對熒光光譜數據中的拉曼和瑞利散射進行去除[15]. 2015年,Elcoroaristizabla和Bro等人對熒光光譜數據處理的相關工作進行了總結[16]. 2019年,Chiappini、Alcaraz和Goicorchea等人提出改進的信息保守方法用于處理熒光數據中的拉曼和瑞利散射[17]. 2019年, Chiappini等人基于已有的算法, 給出處理熒光光譜數據的圖形交互界面, 這使得即使沒有程序語言基礎, 也能夠輕松實現光譜數據的處理[18].

上述工作中, 對無用散射的判斷, 大部分是在Ex方向上基于高斯分布判斷散射的寬度, 高度等信息. 當散射能量不滿足對稱或高斯分布時, 散射的信息便有待進一步確定. 本文旨在提出一種基于數據自動判斷無用散射的機制, 在完成去除后, 通過剩余數據對刪除數據實現補全. 這一算法對熒光光譜數據并無高斯分布的要求.

1 去無用散射算法

在光譜數據中, 瑞利散射滿足

Em=n×Ex,

(1)

其中,n代表n階散射.

在光譜數據中, 拉曼散射滿足

(2)

其中,wnr是光譜數據中瑞利散射和拉曼散射的波數之差.

在無用散射去除過程中, 熒光光譜需要滿足如下兩條件要求：①瑞利和拉曼散射的能量大于目標物質的能量；②瑞利和拉曼散射至少存在一種. 在Em方向中，散射能量不要求為高斯分布. 最高值可以為定值. 無用散射的去除過程有3個主要步驟：①數據預處理及散射判斷；②散射信息處理以及散射去除；③去除補齊.

在對光譜數據實現預處理和散射判斷的過程中, 由于生成機制所致, 熒光光譜中會存在大量的散點噪聲, 這使得光譜的能量并不平滑, 會影響之后的無用散射信息的判斷準確性. 于是, 在本文提出的算法中, 通過2D的均值濾波實現平滑處理, 在基本保留熒光光譜的主要信息的前提下, 去除此類噪聲. 具體的, 如下式所示：

(3)

其中,ni、nj分別是在2D數據中沿Ex、Em方向上的平滑寬度.

在熒光光譜的數據中, 通過能量值判斷瑞利或拉曼散射是否存在. 根據兩種散射的理論基礎, 兩者在光譜中的位置基本固定, 需要確定的是必須去除的無用散射寬度.

本文的思路是通過迭代的方式, 在一定范圍不斷選擇能量高的數據, 通過最小二乘法擬合出直線方程.當直線方程在散射的范圍之內, 則認為存在散射;若直線方程不在散射的范圍之內, 有可能是因為選取能量高的數據點過少, 則需放寬能量點的選擇范圍, 篩選更多的能量點, 通過上述的方法直至擬合出在散射范圍內的直線方程. 具體的判斷流程如圖1所示.

圖1 無用散射的確定過程流程圖Fig. 1 The flowchart of detecting the useless scatter

對于光譜數據, 通過針對每個Em值上不同Ex對應的強度Z, 實現散射寬度的確定. 首先, 通過一維數據說明無用散射的信息確定過程. 如圖2(A)顯示的是特定Em值, 強度Z隨Ex的變化的數據. 由于已經經過2D的平滑處理, 強度曲線中已經不存在散點噪聲, 但曲線依舊不平滑. 本文采用1D的平滑處理, 使曲線保持其原有的變化趨勢, 但能平滑掉細小的凹凸不平, 具體的如圖2(B)所示. 因為Em已知, 根據之前確定的散射方程, 可以得到無用散射的峰值的粗略范圍,Sp在本論文中同時考慮曲線的一階導數和二階導數判斷散射峰的范圍.

針對平滑后的數據求解能量強度Z, 對于Ex的一階導數和二階導數, 分別在圖2(C)和圖2(D)所示. 其中，距離Sp最近的一階導數為0的點為無用散射的寬度范圍,二階導數為0的點用于進一步判斷寬度范圍是否準確.

圖2 無用散射信息確定過程的示意圖Fig.2 The illustration of deciding the range of useless scatterA.隨Ex的變化的2D平滑后的強度值(特定Em);B.1D平滑后的強度值;C.強度的一階導數變化;D.強度的二階導數變化

在確定無用散射的寬度后, 將在確定的無用散射范圍內的強度值設為空(NaN).

在去除無用散射的能量后, 通過插值的方式, 將NaN區域進行補齊, 以確保熒光光譜的完整性. 本文中采用的是三次多項式的插值實現.

對于一個特定Em值對應的能量曲線, 分別針對拉曼和瑞利散射實現上述過程, 完成單條能量曲線上的兩種無用散射位置的補齊. 對于每個Em值均完成上述處理, 即完成整個光譜的無用散射的處理工作.算法的流程圖如圖3所示.對數據的處理示意圖如圖4所示.

圖3 三維熒光光譜數據去除無用散射的實現步驟Fig. 3 The steps of deleting the useless scatter in 3D EEM fluorescence

圖4 三維熒光光譜數據去除無用散射的實現流程圖Fig. 4 The illustrations of deleting the useless scatter in 3D EEM fluorescence

2 測試光譜生成

針對新的去散射算法, 完全基于實驗所采集的光譜數據進行驗證需要消耗大量的資源. 并且由于真實的實驗環境所限, 絕大部分的實驗過程都無法保證實驗所使用數據的完備性, 于是本文嘗試使用依據理論生成的熒光光譜數據來驗證所提算法針對不同特點數據的魯棒性.

首先, 對于要生成的光譜數據, 確定其中目標能量的位置. 如前所述, 在光譜數據中拉曼瑞利兩種散射生成機制相對固定, 其中,不同譜之間無用散射改變的是其對應的長度、寬度、強度和分散程度等性質. 為了與現實中的光譜數據更加接近, 在生成譜的過程中還在數據中加入散點噪聲和半徑小于Rnoise的圓形噪聲. 為了驗證所提算法對無用散射和目標能量的區分能力, 針對上述信息隨機生成具有不同性質的無用散射的光譜數據.

將目標物對應的能量生成在光譜數據的有效范圍中, 其中, 拉曼散射和瑞利散射包圍范圍之外的部分不生成目標物的能量, 目標物的能量采用ntag個二維高斯分布代替. 其中，σ1,σ2在一定的范圍內選取在本數值實驗中選擇ntag=1,2,3.

目標散射和無用散射基本存在三種位置關系, 分別是相距較遠, 部分相交和部分覆蓋.

在圖5(A)中顯示的是目標物的能量與兩種散射均不相交的情況.在圖5(B)中顯示的是目標物的能量只部分與散射相交且散射的能量覆蓋目標物能量. 在圖5(C)中顯示的是目標物的能量散射能量覆蓋目標物.圖5(D)～圖5(F)分別對應著子圖5(A)～圖5(C)中固定Em值能量強度對應Ex的變化過程. 基于生成的熒光光譜數據,采用本文提出的去噪算法,會得到去除無用散射后的光譜數據.針對三種不同的情況,分別生成若干個測試用光譜數據,圖6～圖8分別針對每種情況選擇三組數據用于顯示和說明.每個圖的子圖中,左側子圖是依據理論生成的光譜數據,側子圖是去散射之后,目標能量的分布.由于目標能量相比于無用散射低所以一般的右圖中的能量范圍要小于左圖中的能量范圍.三種情況中,本論文算法均能很準確地判斷無用散射的位置和能量,實現無用散射的去除和補全.

圖5 目標能量和無用散射的位置關系Fig. 5 The position relationship of target energy and useless scatterA.目標物的能量與無用散射均不相交的三維光譜圖；B.目標物的能量只部分與無用散射相交的三維光譜圖；C.無用散射能量覆蓋大部分目標物能量的三維光譜圖；D.基于A圖數據下固定Em值能量強度對應Ex的變化過程；E.基于B圖數據下固定Em值能量強度對應Ex的變化過程；F.基于C圖數據下固定Em值能量強度對應Ex的變化過程

具體如圖6所示, 目標能量在散射之間, 且分別距離兩種散射較遠, 這使得對無用散射的判斷不受目標能量影響, 容易確定其覆蓋范圍. 于是, 具有較高能量的無用散射被剔除后, 基本沒有影響目標物的能量. 在三組生成的光譜數據中, 所提算法均很準確地實現了無用散射的剔除.

圖6 去無用散射前后光譜數據的對比圖. 目標能量和無用散射相距較遠Fig. 6 The comparison between the fluorescence data before and after deleting the useless scatter. Large distance between target energy and useless scatter

在下面的各個光譜數據圖像中, 橫軸為Em, 縱軸為Ex, 光譜強度以等高線的形式顯示, 激發波長和發射波長單位為nm,光譜強度單位為a.u.,其中，等高線顏色對應的光譜強度由每個子圖右側的colorbar顯示. 之后生成的光譜數據、實驗數據的圖像均以此種方式表達.

如圖7所示, 當目標能量存在一部分與無用散射相交的情況時, 目標能量對無用散射的判斷也產生影響, 但這種影響相對于無用散射的高能量, 基本可以忽略. 目標能量中, 能量較低的部分由插值得到. 于是, 本文所提算法依舊可以準確地實現無用散射的剔除.

圖7 去無用散射前后光譜數據的對比圖.目標能量和無用散射部分相交Fig.7 The comparison between the fluorescence data before and after deleting the useless scatter. Partly intersection of target energy and useless scatter

如圖8所示, 部分的目標能量與無用散射重疊, 這使得對無用散射的范圍的判斷難度提升. 在目標能量中, 存在一部分通過插值得到的數據, 但根據光譜數據進行插值, 可以確保保留目標物的能量變化趨勢, 于是, 所提算法還是可以準確實現無用散射的剔除.

圖8 去無用散射前后光譜數據的對比圖.目標能量和無用散射部分重疊Fig.8 The comparison between the fluorescence data before and after deleting the useless scatter. Partly overlapping of target energy and useless scatter

3 實驗光譜的去無用散射

通過生成的熒光光譜可以驗證所提算法的可行性. 下面將該算法用于實驗過程中得到的真實的光譜數據中.

對比于生成的光譜數據, 真實的熒光光譜數據具有更強的隨機性, 其中目標函數的范圍也更加的不規則, 這給去除無用散射增加了難度.

如圖9所示, 在光譜數據中, 目標能量覆蓋的范圍更大, 且只存在一組無用散射的影響. 對于每個Em值, 判斷無用散射的難度加大. 但由于散射與目標能量相互獨立, 所以所提算法很準確地去除了無用散射, 并保留下大范圍的目標能量.

圖9 去無用散射前后光譜數據的對比圖. 目標能量和無用散射相距較遠(實驗數據)Fig.9 The comparison between the fluorescence data before and after deleting the useless scatter. Large distance between target energy and useless scatter(experiment)

如圖10所示, 目標能量有一部分與無用散射相交, 而且無用散射以很高的能量影響目標能量. 即便如此, 目標能量也被準確的保留, 得到較好的無用散射去除結果.

圖10 去無用散射前后光譜數據的對比圖.目標能量和無用散射部分相交(實驗數據)Fig.10 The partly intersection of target energy and useless scatter (experiment). Partly intersection of target energy and useless scatter

4 結 論

本文提出一種基于光譜數據判斷無用散射的去除算法. 針對生成的光譜數據和實驗得到的光譜數據, 均得到了較好的無用散射去除效果, 本算法為進一步的高通量自動去除無用散射提供可能. 本文所提出的算法, 暫時只針對一階瑞利和拉曼散射. 在未來的工作中, 如何基于現有機制實現高階無用散射的去除, 是我們要進一步研究和探討的主要內容.