分?jǐn)?shù)評(píng)估法實(shí)現(xiàn)光通過(guò)散射介質(zhì)聚焦

方龍杰,程科,何修軍

(成都信息工程大學(xué)光電工程學(xué)院,四川 成都 610225)

1 引言

光散射成像技術(shù)是生物醫(yī)學(xué)領(lǐng)域一項(xiàng)非常重要的技術(shù)。隨著空間光調(diào)制器的迅速發(fā)展,波前調(diào)制在光通過(guò)散射介質(zhì)聚焦和成像領(lǐng)域得到了廣泛的研究[1-3]。這些研究成果在光聲成像[4]、顯微操作[5]、生物醫(yī)學(xué)檢測(cè)[6]、三維全息顯示[7]、光纖通信[8]、超衍射成像[9]等領(lǐng)域具有廣闊的應(yīng)用前景。為了實(shí)現(xiàn)光通過(guò)散射介質(zhì)聚焦,需要尋找與散射介質(zhì)相匹配的輸入光場(chǎng),其中用到的算法包括:數(shù)字相位共軛法[10]、傳輸矩陣算法[2,11]和迭代反饋算法[12]。數(shù)字相位共軛方法利用共軛相位來(lái)抵消散射介質(zhì)的隨機(jī)擾動(dòng)效應(yīng),但其光學(xué)系統(tǒng)較為復(fù)雜,且需精確校準(zhǔn)[10]。傳輸矩陣方法通過(guò)建立輸入光場(chǎng)和輸出光場(chǎng)的線性關(guān)系實(shí)現(xiàn)聚焦,但難以獲得高維度的傳輸矩陣,并且抗干擾能力較弱[13]。迭代反饋算法根據(jù)散射光場(chǎng)的實(shí)時(shí)反饋不斷調(diào)節(jié)入射光場(chǎng),使焦點(diǎn)的光強(qiáng)逐漸增大。這種算法操作簡(jiǎn)單,在實(shí)際應(yīng)用當(dāng)中易于實(shí)現(xiàn)[14]。常見(jiàn)的迭代反饋算法有蟻群算法[15]、麻雀搜索算法[16]等。然而,散射光場(chǎng)容易受到環(huán)境擾動(dòng)而發(fā)生畸變,在噪聲環(huán)境下難以實(shí)現(xiàn)輸出光場(chǎng)的精準(zhǔn)調(diào)控。

近年來(lái),深度學(xué)習(xí)算法在散射成像領(lǐng)域得到了廣泛的應(yīng)用。該算法基于數(shù)據(jù)集建立輸入光場(chǎng)和輸出光場(chǎng)之間的映射關(guān)系,具有強(qiáng)大的擬合能力,在抗散射過(guò)程中展示出其獨(dú)特的優(yōu)勢(shì)。例如,Horisaki等利用神經(jīng)網(wǎng)絡(luò)算法對(duì)散斑圖像進(jìn)行分類和重構(gòu)[17]。Li等使用基于密集塊的Unet網(wǎng)絡(luò)實(shí)現(xiàn)相干光通過(guò)毛玻璃成像[18]。Unet網(wǎng)絡(luò)具有架構(gòu)簡(jiǎn)單、精度高、易于實(shí)現(xiàn)等優(yōu)點(diǎn),在散射成像領(lǐng)域得到了廣泛的研究,但僅僅適用于薄散射介質(zhì)的圖像重建。對(duì)于厚散射介質(zhì)的成像,這種方法將會(huì)失效。為了解決這一問(wèn)題,Lyu等提出深度神經(jīng)網(wǎng)絡(luò)和混合神經(jīng)網(wǎng)絡(luò)結(jié)構(gòu),成功實(shí)現(xiàn)通過(guò)厚度為3 mm的聚苯乙烯散射介質(zhì)成像,成像深度為介質(zhì)平均自由程的13.4倍[19]。Guo等構(gòu)建的抗散射卷積神經(jīng)網(wǎng)絡(luò)將散射成像的記憶效應(yīng)范圍擴(kuò)大了40倍[20]。

深度學(xué)習(xí)算法雖然能夠?qū)崿F(xiàn)厚散射介質(zhì)的聚焦與成像,但是要求散射介質(zhì)為靜態(tài)散射介質(zhì)。而對(duì)于動(dòng)態(tài)散射介質(zhì),深度學(xué)習(xí)算法則失效。因此,本文提出分?jǐn)?shù)評(píng)估法實(shí)現(xiàn)光通過(guò)散射介質(zhì)聚焦。根據(jù)線性傳輸光學(xué)系統(tǒng)的輸入通道理論,散射聚焦的物理本質(zhì)是保留對(duì)聚焦有益的輸入通道,并且去掉對(duì)聚焦不利的輸入通道。而對(duì)于聚焦不利的通道,如果將其附加π的相位值,則不利通道會(huì)變?yōu)橛欣ǖ?。我們采取二元相位?yōu)化的方式提高迭代程序運(yùn)行的速度,縮短聚焦時(shí)間,便于將其運(yùn)用于生物成像。在實(shí)驗(yàn)中,將分?jǐn)?shù)評(píng)估法與蟻群算法[15]對(duì)比,結(jié)果發(fā)現(xiàn),在強(qiáng)噪聲環(huán)境下,分?jǐn)?shù)評(píng)估法仍能實(shí)現(xiàn)光通過(guò)散射介質(zhì)聚焦,而蟻群算法則失效,說(shuō)明分?jǐn)?shù)評(píng)估法具有更強(qiáng)的抗噪聲能力。

2 算法描述

分?jǐn)?shù)評(píng)估算法的流程圖如圖1所示。首先,在空間光調(diào)制器上加載大量的隨機(jī)相位分布(設(shè)總共加載N個(gè)相位分布)。為了實(shí)現(xiàn)快速聚焦,我們采用的相位板為二元相位板,即只含有0和 π 兩個(gè)相位的取值。對(duì)于每一個(gè)隨機(jī)相位分布,都用CCD相機(jī)記錄經(jīng)過(guò)散射介質(zhì)后的散斑光強(qiáng)分布(總共有N個(gè)散斑光強(qiáng)分布)。本實(shí)驗(yàn)的目標(biāo)是實(shí)現(xiàn)光通過(guò)散射介質(zhì)聚焦,因此目標(biāo)函數(shù)為

圖1 分?jǐn)?shù)評(píng)估算法示意圖Fig.1 The schematic diagram of the score assessment algorithm

ftarget(S)=Ifocus

(1)

其中,S為空間光調(diào)制器的相位分布矩陣,Ifocus為焦點(diǎn)處的光強(qiáng)值。

接下來(lái)構(gòu)造分?jǐn)?shù)評(píng)估矩陣Xi(i= 1, 2, …,N)。Xi的維數(shù)是S矩陣的二倍,設(shè)S矩陣總共含有n個(gè)元素,則Xi矩陣含有n×2個(gè)元素,其中第一列代表相位 π 的分值,第二列代表相位0的分值。給分的原則按照焦點(diǎn)光強(qiáng)的大小賦分。例如,對(duì)于第一個(gè)隨機(jī)相位分布,測(cè)得焦點(diǎn)光強(qiáng)為6,此時(shí)第一個(gè)空間光調(diào)制器單元取的是相位0,因此相位0給6分,相位 π 給0分,所以矩陣X1的第一行為0和6;第二個(gè)空間光調(diào)制器單元取的是相位 π ,因此相位 π 給6分,相位0給0分。以此類推,則可以找出每一個(gè)隨機(jī)相位分布對(duì)應(yīng)的分?jǐn)?shù)評(píng)估矩陣。

最后,計(jì)算總的分?jǐn)?shù)評(píng)估矩陣,方法是把所有的分?jǐn)?shù)評(píng)估矩陣相加,即總的分?jǐn)?shù)評(píng)估矩陣

(2)

矩陣X表示相位0和 π 獲得的總分。比較0和 π 兩個(gè)相位獲得的總分,則優(yōu)化位相分布Xopt應(yīng)當(dāng)選取總分較高的位相。例如,對(duì)于空間光調(diào)制器的第一個(gè)像素單元,相位 π 的得分是82分,而相位0的得分是48分,因此第一個(gè)像素單元的相位應(yīng)該是 π 。分?jǐn)?shù)評(píng)估法的物理實(shí)質(zhì)是按照對(duì)焦點(diǎn)光強(qiáng)貢獻(xiàn)的份額定量評(píng)估應(yīng)采用同相還是反相。當(dāng)焦點(diǎn)光強(qiáng)越高時(shí),說(shuō)明該位相為正確位相的概率越大,其得分應(yīng)該越高;反之,如果焦點(diǎn)光強(qiáng)較低,說(shuō)明該點(diǎn)位相正確的可能性更小,應(yīng)當(dāng)以較小的分?jǐn)?shù)評(píng)估該相位值。

3 數(shù)值模擬

為了驗(yàn)證分?jǐn)?shù)評(píng)估算法在噪聲擾動(dòng)下的抗散射能力,我們通過(guò)數(shù)值模擬使光通過(guò)散射介質(zhì)聚焦,并將分?jǐn)?shù)評(píng)估法與蟻群算法比較。數(shù)值模擬的光路示意圖如圖2所示。激光器發(fā)出的平面波通過(guò)空間光調(diào)制器進(jìn)行相位調(diào)制,經(jīng)過(guò)透鏡匯聚到散射介質(zhì)的前表面,然后經(jīng)介質(zhì)散射的光波用CCD接收。為了簡(jiǎn)化計(jì)算,我們假設(shè)CCD和空間光調(diào)制器的像素個(gè)數(shù)相同,散射介質(zhì)為服從均勻分布的隨機(jī)矩陣,聚焦的位置為CCD的中心處。根據(jù)光場(chǎng)傳輸?shù)膹?fù)振幅理論,散射輸出光場(chǎng)Eout(x,y)可以表示為

圖2 數(shù)值模擬光路示意圖Fig.2 Schematic of optical path in numeric simulation

Eout(x,y)=F{exp[iφSLM(x,y)]·exp[iφmedium(x,y)]}

(3)

其中,F表示透鏡的傅里葉變換作用,φSLM(x,y)表示空間光調(diào)制器的相位調(diào)制作用,φmedium(x,y)表示散射介質(zhì)的隨機(jī)相位調(diào)制作用。CCD接收到的圖像應(yīng)為散斑光強(qiáng)分布,可以用下式表示:

Iout(x,y)=|Eout(x,y)|2

(4)

其中|Eout(x,y)|表示輸出光場(chǎng)的振幅,Iout(x,y)表示CCD相機(jī)接收到的散斑光強(qiáng)分布。在每次數(shù)值模擬中加入隨機(jī)噪聲來(lái)比較算法的抗噪聲能力。隨機(jī)噪聲大小與焦點(diǎn)光強(qiáng)有關(guān),可以表示為

Inoise(x,y)=Iout(x,y)+μ·Iout(x,y)·δ

(5)

其中,μ是一個(gè)介于[-0.5, 0.5]區(qū)間的隨機(jī)數(shù);δ為比例常數(shù),在模擬中取四個(gè)值:0、0.01、0.05和0.1,分別表示無(wú)噪聲、加入方差為1%、5%和10%的高斯白噪聲。

數(shù)值模擬中,所使用的空間光調(diào)制器像素單元個(gè)數(shù)是64×64=4096,仿真程序的迭代次數(shù)為10000次迭代。蟻群算法中采用的參數(shù)為:螞蟻個(gè)體數(shù)量為100,初始轉(zhuǎn)移概率p0為0.2,信息素衰減系數(shù)ρ為0.8,信息素濃度閾值(算法終止條件)為0.0001。這些參數(shù)的確定都是通過(guò)大量的模擬來(lái)優(yōu)化參數(shù),使模擬程序在當(dāng)前條件下達(dá)到最佳狀態(tài),從而使比較更加客觀。

圖3展示了分?jǐn)?shù)評(píng)估算法和蟻群算法在不同噪聲情況下的增強(qiáng)因子隨迭代次數(shù)的變化曲線。其中SAA表示分?jǐn)?shù)評(píng)估算法,ACO表示蟻群算法。從圖3(a)可以看出,在不引入噪聲的情況下,分?jǐn)?shù)評(píng)估算法在迭代10000次后的增強(qiáng)因子與蟻群算法的增強(qiáng)因子基本相同。需要指出的是,分?jǐn)?shù)評(píng)估算法的初始增長(zhǎng)斜率明顯高于蟻群算法,但在之后的迭代中,增長(zhǎng)斜率逐漸放緩,這主要是因?yàn)橄辔粌?yōu)化已經(jīng)接近于飽和狀態(tài)。當(dāng)引入隨機(jī)噪聲時(shí),從圖3(b)-(d)可以看出,在相同的噪聲環(huán)境下,分?jǐn)?shù)評(píng)估算法得到的焦點(diǎn)增強(qiáng)因子高于蟻群算法,從而獲得最佳信噪比。同時(shí)也證實(shí)了分?jǐn)?shù)評(píng)估算法在散射聚焦和成像方面具有較好的抗噪能力。值得注意的是,分?jǐn)?shù)評(píng)估算法不僅在增強(qiáng)速率上比蟻群算法更快,而且在相同條件下,與蟻群算法相比,獲得最優(yōu)解所需的迭代次數(shù)更少。造成這種結(jié)果的主要原因是兩種算法不同的工作機(jī)制。蟻群算法在優(yōu)化過(guò)程中具有隨機(jī)性,可能會(huì)使得當(dāng)前解偏離最優(yōu)位置附近。而分?jǐn)?shù)評(píng)估法是通過(guò)在現(xiàn)有最佳相位分布的基礎(chǔ)上探索更優(yōu)的解決方案。

圖3 不同算法在不同噪聲下的增強(qiáng)因子。其中SAA表示分?jǐn)?shù)評(píng)估算法,ACO表示蟻群算法。(a) 無(wú)噪聲情形;(b) 加入1%的噪聲;(c) 加入5%的噪聲;(d) 加入10%的噪聲Fig.3 Enhancement factors of different algorithms under different noise level. SAA means score assessment algorithm, and ACO means ant colony algorithm. (a) no noise; (b) 1% noise; (c) 5% noise; (d) 10% noise

在數(shù)值模擬中,分?jǐn)?shù)評(píng)估算法采用的迭代次數(shù)是10000次。原則上講,實(shí)驗(yàn)次數(shù)越多,焦點(diǎn)光強(qiáng)越大。而從模擬結(jié)果可以看出,迭代10000次時(shí)焦點(diǎn)光強(qiáng)已達(dá)到飽和,因此選取10000次迭代是可行的。在實(shí)際應(yīng)用中,迭代次數(shù)由介質(zhì)的特性確定。如果散射介質(zhì)是動(dòng)態(tài)的,則需要在介質(zhì)的去相關(guān)時(shí)間內(nèi)完成迭代優(yōu)化,否則會(huì)因?yàn)榻橘|(zhì)的變化導(dǎo)致聚焦失敗。在實(shí)際應(yīng)用過(guò)程中,也會(huì)存在大量的隨機(jī)誤差,實(shí)驗(yàn)次數(shù)越多,累積的隨機(jī)誤差也越多。因此,需要在焦點(diǎn)光強(qiáng)達(dá)到預(yù)期值時(shí)立刻終止實(shí)驗(yàn)程序,防止誤差累積太大。

4 實(shí)驗(yàn)

聚焦實(shí)驗(yàn)采用的實(shí)驗(yàn)裝置圖如圖4所示。波長(zhǎng)為532.8 nm的激光器(MSL-FN-532.8 nm-200 mw,CNI)發(fā)出圓偏振光,經(jīng)過(guò)偏振片P1產(chǎn)生線偏振光,然后通過(guò)由透鏡L1(f=150 mm)和透鏡L2(f=250 mm)組成的4f系統(tǒng)擴(kuò)束。擴(kuò)束后的光通過(guò)偏振片P2改變偏振方向,并且偏振片P2的偏振方向與空間光調(diào)制器(Holoeye PLUTO-TEL-CO,1920×1080像素)的純相位調(diào)制工作模式的偏振方向一致。旋轉(zhuǎn)偏振片P1和P2可以調(diào)節(jié)入射光的強(qiáng)度。經(jīng)過(guò)空間光調(diào)制器反射的光波通過(guò)透鏡L3(f=100 mm)和L4(f=50 mm)縮束,并使用10×顯微物鏡(OL1,NA=0.25)將光聚焦到散射介質(zhì)的前表面。本實(shí)驗(yàn)使用的散射介質(zhì)為毛玻璃。經(jīng)過(guò)介質(zhì)散射的光通過(guò)20 ×顯微物鏡(OL2,NA=0.4)成像到CCD(Point Grey Research,Blackfly BFLY-U3-23S6M)上,并將接收到的圖像輸入計(jì)算機(jī),可獲取焦點(diǎn)光強(qiáng)。

圖4 散射聚焦的實(shí)驗(yàn)裝置圖。其中P1、P2為偏振片;L1、L2、L3、L4為透鏡;SLM為空間光調(diào)制器;OL1、OL2為顯微物鏡;S為散射介質(zhì);CCD為電荷耦合器件Fig.4 Experimental setup for focusing light through scattering media. P1, P2: polarizer. L1, L2, L3, L4: lenses. SLM: spatial light modulator. OL1, OL2: objective lenses. S: scattering medium. CCD: charge coupled device

實(shí)驗(yàn)得到的聚焦效果如圖5所示。為了比較在強(qiáng)噪聲和弱噪聲情況下兩種算法的抗噪聲能力,我們通過(guò)調(diào)節(jié)激光器的輸出功率產(chǎn)生兩種不同的噪聲環(huán)境。激光器的最大輸出功率是200 mw,將激光器輸出功率分別調(diào)為200 mw和100 mw,對(duì)應(yīng)的噪聲大小為10%的噪聲和30%的噪聲。此百分比為測(cè)量30分鐘內(nèi)輸出光強(qiáng)的均方差占輸出功率均值的比例,噪聲來(lái)源包括光學(xué)平臺(tái)的振動(dòng)、激光器輸出功率的不穩(wěn)定、實(shí)驗(yàn)室的雜散光等等。圖5a為分?jǐn)?shù)評(píng)估算法和蟻群算法的聚焦結(jié)果,其中SAA為分?jǐn)?shù)評(píng)估算法,ACO為蟻群算法?？梢钥闯?當(dāng)噪聲為10%的情況下,兩種算法均能夠得到焦點(diǎn),分?jǐn)?shù)評(píng)估算法的焦點(diǎn)光強(qiáng)為優(yōu)化前散斑平均光強(qiáng)的71.4倍。而在強(qiáng)噪聲條件下(30%噪聲),蟻群算法失效,而分?jǐn)?shù)評(píng)估算法仍然能夠得到焦點(diǎn)。圖5b和圖5c分別為兩種算法在10%和30%噪聲下焦點(diǎn)增強(qiáng)因子隨迭代次數(shù)的演化曲線。可以看出,當(dāng)噪聲為10%時(shí),蟻群算法的最終增強(qiáng)因子為66.7,與分?jǐn)?shù)評(píng)估算法的增強(qiáng)因子基本相同。但是,當(dāng)噪聲為30%時(shí),蟻群算法的最終增強(qiáng)因子下降為4.1,而蟻群算法的增強(qiáng)因子為51.6。說(shuō)明分?jǐn)?shù)評(píng)估算法的抗噪聲能力高于蟻群算法。

圖5 實(shí)驗(yàn)得到的聚焦結(jié)果。(a)最終聚焦結(jié)果;(b)10%噪聲下的增強(qiáng)曲線;(c)30%噪聲下的增強(qiáng)曲線Fig.5 Experiment results of focusing. (a) Final focusing results. (b)Enhancement curves of 10% noise. (c) Enhancement curves of 30% noise

5 討論

在分?jǐn)?shù)評(píng)估算法的數(shù)值模擬部分,我們采用二元位相調(diào)制方法,獲取聚焦的優(yōu)化位相分布。假設(shè)采取連續(xù)位相調(diào)制方法,以空間光調(diào)制器的調(diào)制單元數(shù)為64×64=4096為例,迭代次數(shù)為10000次,在[0, 2π]區(qū)間等間距取16個(gè)相位,獲得的焦點(diǎn)增強(qiáng)倍數(shù)如圖6所示。為了討論散射的隨機(jī)性對(duì)最終焦點(diǎn)增強(qiáng)倍數(shù)的影響,我們總共進(jìn)行30次模擬。二元位相調(diào)制方法同樣也進(jìn)行30次模擬,獲得的最終焦點(diǎn)增強(qiáng)倍數(shù)也展示在圖6中。對(duì)于每一種方法,所有的增強(qiáng)因子按照升序排列。從圖6可以看出,連續(xù)調(diào)制方法的增強(qiáng)因子高于二元位相調(diào)制方法的增強(qiáng)因子。對(duì)于30次數(shù)值模擬,連續(xù)調(diào)制算法的平均增強(qiáng)因子是1101.4,二元調(diào)制方法的平均增強(qiáng)因子是1004.7。連續(xù)調(diào)制算法的增強(qiáng)倍數(shù)僅僅比二元調(diào)制算法高出9.7%。而從測(cè)量時(shí)間的角度來(lái)看,用二元調(diào)制方法需要的測(cè)量次數(shù)為2×10000=20000次,而用連續(xù)調(diào)制方法,需要的測(cè)量次數(shù)是16×10000=160000次,是二元調(diào)制方法耗時(shí)的8倍。因此,采用二元調(diào)制方法完全是可行的。

圖6 連續(xù)調(diào)制方法和二元調(diào)制方法的焦點(diǎn)增強(qiáng)因子。其中CA表示連續(xù)調(diào)制方法,BA表示二元優(yōu)化算法Fig.6 Enhancement factors of continuous algorithm and binary algorithm. CA means continuous algorithm, and BA means binary algorithm

分?jǐn)?shù)評(píng)估算法的前提假設(shè)是,相位調(diào)制單元對(duì)聚焦的貢獻(xiàn)是確定的。但是,由于尋找優(yōu)化位相分布這一問(wèn)題是多解問(wèn)題,每一個(gè)調(diào)制單元的相位其實(shí)與其他單元的相位有關(guān),它們是相互干涉后的結(jié)果。因此,通過(guò)分?jǐn)?shù)評(píng)估算法得到的優(yōu)化位相分布其實(shí)是眾多優(yōu)化解當(dāng)中的一個(gè)。在相位優(yōu)化過(guò)程中,由于隨機(jī)噪聲的引入,會(huì)使得局域相位值發(fā)生改變,因此對(duì)于同一散射介質(zhì),不同的實(shí)驗(yàn)次數(shù)得到的優(yōu)化結(jié)果會(huì)不同。我們對(duì)比了同一散射介質(zhì),調(diào)制單元數(shù)64×64=4096,加入10%的噪聲情況下,5次模擬得到的優(yōu)化位相分布,結(jié)果如表1所示。從表1可以看出,由于噪聲的影響,每一次獲得的相位分布均不相同,但是都能夠獲得聚焦效果,說(shuō)明這是一個(gè)多解問(wèn)題,并且隨機(jī)噪聲的擾動(dòng)會(huì)對(duì)實(shí)驗(yàn)結(jié)果產(chǎn)生影響。

表1 五次數(shù)值模擬得到的優(yōu)化位相數(shù)量統(tǒng)計(jì)(調(diào)制單元數(shù)為64×64=4096)Table 1 Statistics of the number of optimized phases obtained by five numerical simulations (number of modulation units 64×64=4096)

6 結(jié)論

本文提出了一種基于分?jǐn)?shù)評(píng)估的方式控制入射光場(chǎng)來(lái)實(shí)現(xiàn)光通過(guò)散射介質(zhì)聚焦。與蟻群算法相比,分?jǐn)?shù)評(píng)估算法具有較強(qiáng)的抗噪聲能力,能夠?qū)崿F(xiàn)強(qiáng)噪聲環(huán)境下的散射聚焦。分?jǐn)?shù)評(píng)估算法利用二元相位優(yōu)化,提升了散射聚焦優(yōu)化過(guò)程的速度,縮短了聚焦時(shí)間,可將該方法用于動(dòng)態(tài)散射介質(zhì)的聚焦和成像。通過(guò)實(shí)驗(yàn)獲得分?jǐn)?shù)評(píng)估法的焦點(diǎn)光強(qiáng)是優(yōu)化前散斑平均光強(qiáng)的71.4倍?？傮w而言,分?jǐn)?shù)評(píng)估法有效地實(shí)現(xiàn)了散射光場(chǎng)的控制,在生物成像、醫(yī)學(xué)診斷、生物組織實(shí)時(shí)監(jiān)測(cè)等領(lǐng)域具有廣闊的應(yīng)用前景。