基于贗熱光照明的單發光學散斑成像*

肖曉 杜舒曼 趙富 王晶 劉軍? 李儒新

1) (中國科學院上海光學精密機械研究所, 強場激光物理國家重點實驗室, 上海 201800)

2) (中國科學院大學, 北京 100049)

(2018 年 9 月 17 日收到; 2018 年 10 月 30 日收到修改稿)

散射介質對光的散射是當前限制光學成像深度或距離的一個嚴重的問題. 本文首先數值模擬比較了光透過隨機散射介質成像研究中常用的基于光學記憶效應(memory effect, ME)和自相關(autocorrelation, AC)方法的HIO&ER算法和乒乓(Ping-Pang, PP)算法的優缺點. 通過對HIO&ER算法和PP算法的恢復效果和迭代次數進行比較, 發現PP算法在保持較高恢復效果的前提下擁有更快的運行速度. 實驗中, 利用連續He-Ne激光器和旋轉毛玻璃產生贗熱光源, 通過物鏡對隨機散射介質后數毫米距離內的不同形狀物體進行了單幀成像, 并采用PP算法成功地恢復出微米量級物體的實際圖像. 這一研究結果將進一步促進ME和AC方法在深層生物組織醫學成像研究上的應用. 最后, 實驗研究了不同的物鏡和散射介質的間距對成像恢復的放大率、分辨率和圖像強度的影響特性, 并進行了詳細研究.

1 引 言

隨著納米技術、光子學技術和生物醫學技術的發展, 活體光學成像在生物光子及現代醫療科學領域的應用越來越廣泛[1?3]. 然而, 由于生物組織對光具有吸收和散射作用, 激發光和信號光都難以穿透足夠深度的生物組織, 并且探測光在生物組織內的散射會引起其相位信息的破壞, 因而難以直接通過相機來對深層生物組織進行高分辨成像[4?6]. 散射介質除了引入以上問題, 研究也發現光經過隨機散射介質可以增加實際成像系統的數值孔徑、增大成像系統的視場角、接收來自物體表面的倏逝波并將其散射為行波在遠場傳播, 因而散射介質也被一些研究小組用來進行光學超衍射成像[7?9]. 多年來,國內外的一些研究小組利用波前調制等反饋控制調節法[10,11]及相位共軛時間反演法[12?14], 來實現光經過散射介質的深度成像. 此外, 在早期鬼成像的研究中, 就已經通過光學互相關的原理來減小散射介質對成像的影響[15?17]. 但這些方法的光路系統相對復雜, 容易受環境影響, 難以適應生物組織深度成像的實際應用.

最近的研究驚喜地發現, 利用贗熱光源基于“光學記憶效應”(memory effect, ME)[18?22]和自相關 (autocorrelation, AC)方法[23?25], 可以從探測器探測到的雜亂無章的散斑圖中直接恢復出物體的二維圖像信息. 其中, 非常重要的一點是用到了Gerchberg-Saxton (GS)算法[26?29]來進行物體的相位恢復. 本文對GS算法中的HIO&ER算法[30?32]和 PP (Ping-Pang, PP)算法[33]的恢復效果和迭代次數進行了比較, 結果表明PP算法在保持較高恢復效果的前提下擁有更快的運行速度. 實驗上, 采用PP算法, 并通過前置物鏡將贗熱光源聚焦到特定的成像目標上, 通過后置物鏡將經過散射介質的散射光收集到相機上, 其物體和散射介質間距可達毫米量級. 同時, 通過平移臺移動成像物體即可對其不同部位進行掃描成像, 本文成功恢復出了不同形狀的微米量級物體的實際圖像.

2 理論與算法

2.1 隨機散射介質特性

光波經過隨機散射介質會發生散射效應, 即一部分光在散射介質中隨機傳播, 導致其傳播的相位面發生畸變, 不能清晰成像, 最終形成散斑圖, 如圖1所示. 散斑指光束經過隨機散射介質, 并在介質內發生多重散射后, 產生的隨機無序的顆粒狀圖樣, 本質上是一種干涉現象. 因為觀察點的光場是散射介質顆粒上各點發出的相干子波的疊加, 且光波波長小于散射介質顆粒尺寸, 所以到達觀察點的各個相干子波的相位是隨機分布的, 隨機相位的相干疊加就產生了散斑的隨機強度圖樣.

圖1 散斑產生示意圖Fig.1. Schematic of speckle generation.

但實質上, 攜帶物體信息的光波經過隨機散射介質作用后, 物體信息并沒有完全丟失, 它們只是在散射過程中進行了隨機組合. 因此, 本文假設散射介質是一個光學透鏡, 利用它的光學傳播特性,從散射光中重組并恢復出物體的原有信息. 例如波前整形方法便可有效地補償扭曲的波陣面, 但該方法需要空間光調制器對扭曲波形面的相位進行逐一調制, 其操作過程耗時復雜, 不適于實時成像.

2.2 記憶效應

當入射光照射到散射介質上時, 其空間相位信息將會被打亂因而具有隨機性, 但其中也會包含相應的相關性. 對于入射光束, 其相對散射介質的入射角度可控, 物體經過散射介質后所成的 像與入射角度有關,

式中I(θ) 表示相機探測到的強度圖,O(r) 表示物體的實際像,S(r) 表示散斑強度,r表示空間坐標矢量,θ表示入射光束的入射角,d表示物體與散射介質的間距,?表示卷積操作, 即物體透過散射介質的像為物體實際像和散斑強度的卷積.

為了將物體實際像信息從隨機散斑中提取出來 , 對探測到的光強信息做AC操作,

式中?表示互相關操作,〈·〉表示取平均操作.

假定入射光束寬度為w, 則

式中, J1是一階貝塞爾函數;L是散射介質厚度,遠大于平均自由程;k=2π/λ是波數. (3)式中第一部分可以用來表示圖案相關性程度C(θ)=對于 ME, 其小角度條件需滿足第二部分中表示平均散斑大小, 通過增大w可無限接近于衍射極限.

假定散斑在空間分布是隨機的, 則散斑的AC 結果是一個 δ 函數, 即由此可將散斑信息從成像AC中消除, 從而得到物體實際像AC的近似值根據卷積定理 , 可得

2.3 相位恢復算法

根據相機采集到的光強圖像信息, 首先對其進行濾波等預處理操作, 得到相機光強圖像根據卷積定理, 對其做AC處理可得

根據維納-辛欽定理, 物體的能量譜等于物體AC的傅里葉變換振幅大小, 因此可以通過矩形窗口截取AC的中心部分, 再對其進行二維傅里葉變換, 即可得到物體的能量譜

由于前面所述的圖像信息均是強度信息, 因此在計算推導傅里葉變換過程中丟失了圖像的相位信息. 由此, 將采用GS算法對丟失的相位信息進行恢復, 在迭代過程中需要將得到的物體能量譜開根作為替代模量不斷替換傅里葉變換的模量. 基本的GS算法分為如下五步:

1)對隨機相位值gk(x,y) 做傅里葉變換

2)求得傅里葉變換域的角度

3)用測量得到的物體能量譜開根替換傅里葉變換的模量

4)對(9)式做傅里葉逆變換

5)根據物理約束條件, 生成新的迭代相位值

根據(11)式的不同迭代公式, 可將其分為Error Reduction (ER)和 Hybrid Input-Output(HIO)兩種算法:

1) ER 算法

2) HIO 算法

這里Γ表示不滿足非負的實數的集合.

2.4 Ping-Pang算法

由上可見, ER和HIO方法各有優缺點, ER方法是目前唯一在數學上被證明的解決相位問題的方法, 但該方法收斂速度慢、易受噪聲影響;HIO方法則是目前應用最廣泛的一種方法, 其算法簡單、運行效率高. 鑒于上述兩種算法的優缺點,將HIO和ER兩種算法結合, 衍生出更具潛力的PP算法. 如圖2所示, 首先運用HIO算法對隨機預測的初始相位進行迭代, 將HIO算法迭代的結果作為ER算法的輸入, 利用ER算法進行迭代, 最后根據物理約束條件來確定是否需要繼續迭代. 在PP算法中, 每次迭代均先后進行HIO算法和ER算法, 利用兩種算法的優點, 使得總迭代次數較小, 可以極大地縮短恢復算法運行時間.

3 數值模擬結果

圖3是按照圖2所示的相位恢復算法進行數值模擬的過程, 該方法截取 2 0×20 像素大小的數字圖像和 3 600×3600 像素大小的散射介質點擴散函數進行卷積, 模擬物像經過散射介質后的散斑圖. 圖 3(a)為標準數字“5”圖像; 圖 3(b)為模擬的散射介質的點擴散函數; 圖3(c)為用1000×1000像素大小窗口截取點擴散函數后的卷積合成的帶有數字“5”的散斑圖; 圖3(d)為散射介質的點擴散函數的三維AC結果, 可以看出其是一個 δ 函數,說明在AC操作下散射介質對物體圖像信息沒有影響; 圖 3(e)為標準數字“5”的 AC 結果; 圖 3(f)為散斑圖的AC結果, 與標準數字的AC相同;圖3(g)為優化的能量譜開根; 圖3(h)為最終重建出來的數字“5”.

圖3的成像恢復數值模擬使用的是基于GS算法的HIO算法和ER算法的順序疊加, 其最終的迭代參數設定為物理約束次數N=30 , HIO約束回歸系數β=1:-0.04:0 , 即先對得到的能量譜開根圖像進行 3 0×{(1/0.04)+1}=780 次的HIO迭代, 再對 HIO的迭代結果進行 30次的ER迭代, 最終結果如圖3(h)所示.

為了進一步探究迭代次數以及不同迭代算法對恢復結果的影響, 圖4還利用matlab對其不同情形進行仿真. 圖 4(a)—(c)是 HIO & ER 算法恢復結果, 圖4(d)—(f)是PP算法恢復結果, 參數設定如表1所列.

圖2 相位恢復算法框圖Fig.2. Schematic of phase retrieval algorithm.

圖3 成像過程的數值模擬 (a) 物體; (b) 點擴散函數; (c) 散斑圖; (d) 點擴散函數 AC; (e) 物體 AC; (f) 散斑 AC; (g) 能量譜開根; (h) 重建結果Fig.3. Simulations of imaging process: (a) Object; (b) point diffusion function; (c) speckle pattern; (d) AC of point diffusion function; (e) AC of object; (f) AC of speckle pattern; (g) square root of power spectrum; (h) result of reconstruction.

通過對比圖4(a)—(c)和表1可以看出, 隨著約束回歸系數β步長的增加, 總迭代次數減少, 表示其運行速度加快, 但圖像恢復效果也逐漸下降.圖4(d)是采用PP算法進行恢復的結果, 可以看出其恢復效果接近圖4(a)的效果, 但總迭代次數僅為202次, 相較于圖4(a)有近8倍的速度提升,進一步證明PP算法的優勢. 同時, 隨著約束回歸系數步長的增加, 迭代次數進一步減少, 但恢復效果更差, 圖像出現模糊, 如圖4(e)和圖4(f)所示, 因而PP算法也難以無限制地縮短運行時間.

4 實驗結果

本文通過散射介質分別對標準分辨率板上的不同數字進行成像, 并通過PP算法對其進行圖像恢復, 驗證其對透過散射介質的成像能力.

圖4 不同迭代次數下的恢復效果 (a)—(c) HIO&ER 算法的恢復結果, 其中, (a) β =1:-0.02:0 , (b) β =1:-0.04:0 ,(c) β =1:-0.05:0 ; (d)—(f) PP 算法的恢復結果, 其中, (d) β =3:-0.02:1 , (e) β =3:-0.05:1 , (f) β=3:-0.1:1Fig.4. Retrieval results in different interation times: (?a)?(c) Retrieval results of HIO&ER algorithm when (a) β =1:-0.02:0 ,(b) β =1:-0.04:0 , (c) β =1:-0.05:0 ; (d)(f) Retrieval results of PP algorithm when (d) β =3:-0.02:1 ,(e) β =3:-0.05:1 , (f) β =3:-0.1:1 .

表1 不同情形下算法迭代次數Table 1. Interation times of algorithm in different conditions.

4.1 不同物體的恢復效果

實驗光路圖如圖5(a)所示, 實驗中采用He-Ne 激光器 (632.8 nm,?0.48 mm, ThorLabs)作為光源, 經空間擴束系統(透鏡1焦距為25 mm,透鏡2焦距為200 mm)進行8倍的擴束后入射到旋轉的散射片 (600 砂,?2 口徑, ThorLabs)上, 該散射片利用電機以20 Hz的轉速驅動會將擴束之后連續光的空間相干性打亂, 產生實驗所需的贗熱光, 效果如圖5(b)所示.

產生的贗熱光經由一個前置顯微物鏡(× 20,0.4NA, 工作距 1.2 mm, Olympus)聚焦到標準分辨率板 (1951 USAF 負測試靶,?1 口徑, ThorLabs)的特定數字上, 其中, 分辨率板可以利用一手動平移臺移動位置, 實現前置顯微鏡對分辨率板上不同數字圖像的聚焦提取. 透過分辨率板上不同數字的光, 又經過一片距離分辨率板8 mm遠的散射片(220 砂,?1 口徑, ThorLabs), 使得物體像信息模糊紊亂形成散斑圖, 經過散射片后的實際效果如圖5(c)所示. 這些散斑圖最后通過后置顯微物鏡 2 (× 10, 0.25NA, 工作距 10.6 mm, Olympus)收集到相機sCMOS(Quantalux?黑白相機, 1920×1080 像素, USB 3.0 接口, ThorLabs)上采集成像,獲得的散斑數據輸送到PC機上進行算法恢復.

圖5 通過散射介質成像的光學裝置 (a)實驗光路圖; (b)贗熱光的產生; (c)散斑的產生Fig.5. Optical setups used for imaging through the scattering media: (a) Optical path in experiment; (b) generation of pseudothermal light source; (c) generation of speckle pattern.

圖6是根據圖5所述光學裝置進行實驗的結果, 其中圖 6(a)—(e)分別為 (a)原始物像, (b)相機接收到的散斑圖, (c)散斑圖的AC結果, (d)光譜能量的根值, (e)通過PP相位恢復算法恢復出來的物體幅度信息. 在恢復過程中, 首先截取相機采集到的散斑圖中心 1 050×1050 像素大小的區域,對該區域圖像做預處理(濾波、歸一化處理等, 第二列圖), 然后對該區域進行AC運算操作(第三列圖)并取得近似的物體能量譜信息(第四列圖), 最后通過PP算法進行相位恢復(第五列圖), 為了得到較好的恢復效果, 在此設定的物理約束條件為N=30,β=3:-0.01:1 , 總 迭 代 次 數 為 402 次.圖 6(f)—(t)是對數字“3”,“5”,“6”的恢復過程, 這些數據也驗證了該方法對不同形狀物體成像具有普適性.

4.2 物鏡和散射介質間距的影響

本實驗系統采用前置物鏡將贗熱光源聚焦到特定的成像目標上, 并通過后置物鏡將經過散射介質的散射光收集到相機上. 為了探究物鏡和散射介質間距對物體成像質量的影響, 實驗調節后置顯微鏡焦平面與散射介質的間距, 利用PP算法對sCMOS采集到的散斑圖進行恢復.

如圖7所示, 其中圖7(a)—(f)表示物鏡焦平面與散射介質的間距分別為 7 00 , 9 00 , 1 100 , 1 300 ,1500和 1 700 μm 時對應的sCMOS采集到的散斑圖的AC結果, 圖7(g)—(l)表示相對應的最終恢復效果. 通過對比可以發現, 隨著物鏡焦平面和散射介質間距的增加, 所得到的AC圖和恢復結果圖都有相應的放大, 這表明物鏡對成像具有放大作用,且放大程度與物鏡焦平面和散射介質間距有關. 此外, 隨著間距的增加, sCOMS上收集到的物體采樣點信息增加, 分辨率得到進一步提高. 但由于物鏡的遠離, 其所收集到的經過散射介質后的散射光強變小, 最后恢復出來的圖像強度有所減弱. 因此,實際情況中, 應根據需要來選取合適的散射片與物鏡的距離, 以平衡成像分辨率和圖像強度, 實現最好的效果.

圖6 不同數字的實驗結果 (a)—(e)數字“1”的恢復過程, 其中, (a)物體, (b) sCOMS 成像, (c)散斑 AC, (d) 能量譜開根,(e)重建結果; (f)—(t)數字“3”, “5”, “6”的恢復過程Fig.6. Experimental results for different numbers: (a)?(e) Retrieval process of number “1”, namely, (a) object, (b) sCOMS image,(c) autocorrelaction of speckle pattern, (d) square root of power spectrum, (e) result of reconstruction; (f)?(t) retrieval processes of number “3”, “5” and “6”.

圖7 不同物鏡和散射介質間距對成像效果的影響 (a)—(f)不同間距下的散斑AC結果; (g)—(l)不同間距下的恢復結果Fig.7. Effects of different distance between objective and diffuser: (a)?(f) AC results of speckle pattern in different distance;(g)?(l) retrieval results in different distance.

5 結 論

本文研究了一種基于贗熱光照明的散斑成像方法, 該方法利用ME和AC原理消除了散射介質對物體成像的影響, 利用相位恢復算法通過單幅散斑圖即可實現物體成像的恢復. 同時, 針對傳統基于GS算法的HIO和ER算法迭代次數多的缺陷,設計出一種快速高效的PP算法, 該算法在保持較高分辨率的同時可以有效縮短運行時間, 有利于實現生物組織的實時成像.

采用 2 0×20 像素大小的數字“5”圖像和3600×3600像素大小的散射介質點擴散函數進行卷積模擬散斑圖, 通過數值模擬驗證了基于相位恢復算法成像的可行性. 對HIO&ER算法和PP算法的成像效果進行了模擬對比, 結果表明PP算法相對于HIO&ER算法有近8倍的速度提升, 證明PP算法具有更高的效率. 通過實驗對不同形狀的物體進行了散斑成像, 利用PP算法獲得微米量級物體的圖像恢復結果. 同時, 通過改變物鏡焦平面與散射介質的間距, 使其從 700 μ m 增加到 1700 μ m ,發現最終物體散斑成像的放大率和分辨率有了相應的提高, 而圖像強度有所下降, 這對今后實際生物組織醫學成像的發展具有重要指導意義.