基于最優化線性波數光譜儀的譜域光學相干層析成像系統?

吳彤 孫帥帥 王緒暉 王吉明 赫崇君 顧曉蓉 劉友文

(南京航空航天大學理學院應用物理系,南京 210016)

1 引 言

光學相干層析成像技術(optical coherence tomography,OCT)是一種高速、高分辨率、非接觸的生物醫學光學成像技術[1?8],能夠實現微米級分辨率和毫米級深度成像.該技術已經廣泛用于生物醫學研究、臨床診斷和材料特性的無損檢測等領域[9?15].在譜域OCT(spectral domain OCT,SDOCT)中,干涉光譜信號被光柵光譜儀探測[16].傳統的光柵光譜儀由衍射光柵、聚焦透鏡、線陣相機組成.根據光柵方程,光柵衍射角的正弦值與光波長呈線性關系,與光波數呈非線性關系,則光波數在線陣相機上呈非線性分布.傳統光柵光譜儀的非均勻波數分布特性對譜域OCT成像產生兩個影響:第一,直接對非線性波數光譜干涉信號進行傅里葉變換會導致軸向點擴散函數加寬;第二,線陣相機上每個像素所探測的帶寬不相同會導致靈敏度隨著成像深度增加而下降得很快[17].

為了解決非均勻波數分布問題,傳統的做法是在傅里葉轉換前對采集到的干涉光譜數據進行數值重采樣,把在波長空間均勻分布的數據轉換成波數空間均勻分布的數據[18].但是重采樣過程計算量大,靈敏度下降仍然不能得到改善.另一種做法是基于線性波數光譜儀在不需要大計算量插值處理的基礎上改善靈敏度下降性能.Hu和Rollins[19]第一次把線性波數光譜儀引進譜域OCT系統中,在衍射光柵和聚焦透鏡之間插入一個BK7玻璃材質的等腰色散棱鏡,以棱鏡出射角對波數的導數來衡量波數線性度.在它們的設計中只考慮了棱鏡頂角這一光學結構參數.Gelikonov等[20]設計的線性波數光譜儀同時考慮了棱鏡頂角、棱鏡入射面與光柵平面之間的夾角兩個結構參數,并提出以波數相對偏差為準則評估波數線性度.Watanabe和Itagaki[21]優化了棱鏡材料和光柵棱鏡之間的夾角,以位置-波數導數的標準差作為線性波數優化準則.Lee等[22]通過數值模擬優化光學參數以實現超高分辨率成像,這些參數包括光柵的光柵常數、色散棱鏡的材料、光柵和色散棱鏡間的旋轉角,其最優化角度和棱鏡材料通過最小化角度-波數導數獲得.最近,Lan和Li[23]設計的線性波數光譜儀用在譜域OCT系統中實現了超高分辨率成像,利用定量光線追跡技術計算了11個離散波數的均方根偏差,并以其為準則評估線性波數光譜儀中波數的線性度.

本文提出一個新的優化準則來優化線性波數光譜儀的結構參數,包括棱鏡頂角、棱鏡材料和光柵與棱鏡之間的旋轉角.把考慮的結構參數整合到波數和像素位置的理論關系式中以數值模擬光譜干涉信號及其點擴散函數,以點擴散函數半峰全寬值的倒數作為評價因子來優化這些結構參數.按照優化的結構參數搭建了線性波數光譜儀,實驗測量了軸向分辨率和系統靈敏度,并與理論值比較以評估搭建的線性波數光譜儀的性能.最后用基于搭建的線性波數光譜儀的譜域OCT系統對人的手指指甲進行了在體層析成像.

2 方 法

線性波數光譜儀的結構示意圖見圖1,由準直鏡、衍射光柵、色散棱鏡、聚焦透鏡和線陣相機組成.定義線性波數光譜儀的結構參數為:色散棱鏡頂角α、棱鏡入射面與光柵面之間的旋轉角β、棱鏡材料折射率n.為了使寬帶光源帶寬范圍內的波數均勻分布在線陣相機像素上,需要推導像素位置坐標和波數之間的函數關系.如圖1所示,令像素陣列的一端為坐標原點,中心像素位置坐標為x0,任一像素位置坐標為x.設計中首先應使寬帶光源發出的中心波數激光聚焦在中心像素位置上.當滿足閃耀條件時,中心波數激光的光柵衍射角θd0可以表示為

其中,λkc為中心波數對應的波長,d為光柵常數,θin為衍射光柵的入射角.如圖1所示,定義棱鏡出射角為棱鏡右端面的出射光與法線的夾角.則中心波數的棱鏡出射角θ0可以表示為

通過光線追跡,對應線陣相機上某個特定像素位置坐標x的棱鏡出射角θout可以表示為

圖1 線性波數光譜儀的結構示意圖Fig.1.Schematics of the linear-k spectrometer.

其中f是聚焦透鏡的焦距.給定一組線性波數光譜儀的結構參數(α,β,n),某波數對應的光柵衍射角θd可以由相應的棱鏡出射角θout推出:

把(1),(3),(4)式代入光柵方程,可以得到像素位置坐標x和波數k之間的函數關系:

根據(5)式,可以數值模擬對應某組線性波數光譜儀結構參數的波數分布函數k(x),進而模擬對應樣品某深度的干涉光譜信號2(k),表達式如下:

其中z代表樣品深度位置坐標,k代表對應某組結構參數的波數分布函數,S(k)代表光源的光譜密度函數,理想情況下OCT光源的光譜分布可取為高斯分布,即有S(k)=S0exp(?(k?kc)2/?k2),kc為光源的中心波數,?k為光源的波數帶寬.由干涉光譜信號的傅里葉變換可以得到對應于該樣品深度的點擴散函數.點擴散函數即OCT成像系統的軸向分辨率.定義點擴散函數的半峰全寬(FWHM)值的倒數1/FWHM作為優化設計線性波數光譜儀結構參數的評價因子c.通過精細改變結構參數得到最高的c值即對應于優化的線性波數光譜儀設計.

可以通過計算相對偏差函數ε(x)評估優化線性波數光譜儀的波數線性度.偏差函數描述特定的波數分布函數k(x)和理想的線性波數分布函數k1(x)的偏差[20].理想的線性波數分布函數定義如下:

其中kc代表中心波數,?x代表可探測波數范圍?k所對應的像素陣列空間范圍.

相對偏差函數定義如下:

因此波數分布函數k(x)的非線性度可以用度量δε來數值估計,其表達式如下:

其中i和j代表像素序數(0 6i,j

3 模擬與實驗

優化線性波數光譜儀需要考慮的結構參數有棱鏡頂角α、旋轉角β和棱鏡材料折射率n.在線性波數光譜儀的數值模擬優化中,首先根據商品化的光學元器件的典型參數指標設定模擬參數.其中透射式衍射光柵的光柵常數設定為1200線對每毫米,聚焦透鏡焦距設定為100 mm,線陣掃描相機的單個像素尺寸設定為10μm,可探測的光譜帶寬波數范圍設定為7.189—7.854μm?1.色散棱鏡的材料設定為F2,BK7和SF11,由于這三種光學玻璃材料在此可探測光譜范圍內的折射率改變很小,模擬中材料折射率n設定為中心波數所對應的折射率.商品化的標準色散棱鏡通常為等邊棱鏡,因此棱鏡頂角的取值范圍設定在60?附近.旋轉角的取值范圍設定為21?—31?度.根據(5)和(6)式分別模擬波數分布函數以及對應的干涉光譜信號.深度z設置為典型值1.6 mm.經過對數值模擬的干涉光譜信號進行傅里葉變換可以計算對應的點擴散函數半峰全寬值及其評價因子c值.不同的α,β和n值可以取得不同的c值.一組使c值最大的α,β和n值即為優化的線性波數光譜儀結構參數.為了在理論上評價優化線性波數分布函數的非線性度,根據(7)和(8)式可得優化波數分布函數和理想線性波數分布函數之間的相對偏差函數ε(x),根據(9)式可得優化波數分布函數的非線性度δε.

圖2是實驗室搭建的譜域OCT系統結構示意圖.光源為超發光二極管(SLD)寬帶光源(SLD-371,Superlum公司),中心波長為840 nm,帶寬為74 nm,最大輸出功率約為11 mW.從光源發出的寬帶激光經分光比為50/50的光纖耦合器分為兩束后分別進入樣品臂和參考臂.樣品臂由偏振控制器、光纖準直鏡、二維掃描振鏡(GVS002,Thorlabs Inc.)、聚焦物鏡和被測樣品組成.參考臂由偏振控制器、光纖準直鏡、聚焦透鏡和平面反射鏡(PF10-03-M03,Thorlabs Inc.)組成,其中偏振控制器用于調節參考臂的光波與樣品臂的光波偏振方向一致.從參考鏡和樣品返回的光匯合后形成的干涉光譜條紋信號由優化的線性波數光譜儀探測.線性波數光譜儀由光纖準直鏡、透射式衍射光柵、等邊色散棱鏡、聚焦透鏡和一個線陣掃描互補金屬氧化物半導體(CMOS)相機(sp2048-70 km,Basler Inc.)組成.為了校準和調節線性波數光譜儀中的光路,制作了光柵槽和棱鏡槽用來固定光柵和棱鏡,光柵槽與棱鏡槽之間的夾角設定為優化的旋轉角β.由線陣掃描CMOS相機探測到的干涉光譜條紋信號經圖像采集卡(IMAQ PCIe-1430,National Instruments)采集進入計算機內存,經數據處理后得到OCT圖像.

為了實現對干涉光譜信號的實時處理,搭建的譜域OCT成像系統使用支持通用并行計算性能的圖形處理器GPU對采集到的譜域OCT干涉光譜信號進行處理.實時信號處理和圖像顯示流程如圖3所示.首先,實時采集到的一幀干涉光譜信號數據從主機內存轉移到GPU板載內存;干涉光譜信號數據在數據類型轉換之后減除預先存儲在GPU中的參考光強度;繼而對一幀干涉光譜信號數據進行并行傅里葉變換(FFT)和對數變換,得到一幀OCT圖像數據;再經過一次類型轉換后OCT圖像數據從GPU板載內存轉移至計算機主機內存中,并在顯示器顯示OCT圖像.包括干涉光譜信號采集、GPU并行處理以及OCT圖像顯示的圖形用戶界面使用基于Visual C++2008程序實現.

圖2 基于最優化線性波數光譜儀的譜域OCT系統示意圖Fig.2.Schematic of the SD-OCT system with the optimized linear-k spectrometer.

圖3 譜域OCT實時干涉光譜信號處理和顯示流程圖Fig.3.Flowchart of the real-time signal processing and displaying for the SDOCT system.

4 結 果

4.1 模擬結果

圖4為基于Matlab程序優化線性波數光譜儀的模擬結果.其中圖4(a)為固定旋轉角的情況下,改變棱鏡頂角時數值模擬的點擴散函數曲線,從圖中可知點擴散函數的FWHM隨著棱鏡頂角的變化而變化.圖4(b)—(d)分別為對應三種典型光學材料的評價因子模擬結果.其中圖4(b)為采用F2玻璃色散棱鏡的評價因子c值的模擬結果,當棱鏡頂角為61.4?、光柵與棱鏡間的旋轉角為29.7?時,評價因子c值達到最大值0.1190μm?1,其最優化棱鏡頂角范圍是59.6?—61.8?.圖4(c)為采用BK7玻璃色散棱鏡的評價因子模擬結果,當棱鏡頂角為66.4?、光柵與棱鏡間的旋轉角為27.1?時,評價因子c值達到最大值0.1189μm?1,其最優化棱鏡頂角范圍是65.7?—66.8?.圖4(d)為采用SF11玻璃色散棱鏡的評價因子模擬結果,其材料特性和前兩者差異較大,當棱鏡頂角為51.4?、光柵與棱鏡間的旋轉角為14.9?時,評價因子c值達到最大值0.1181μm?1,其最優化棱鏡頂角范圍是49.5?—51.0?.綜合三種典型材料的數值模擬結果,在60?附近F2棱鏡的評價因子明顯大于BK7棱鏡相對應的評價因子,且F2棱鏡的最優化角度范圍比較大;而BK7和SF11棱鏡情況下的評價因子對角度變化比較敏感.考慮到商品化的色散棱鏡通常都為等邊棱鏡,根據數值模擬結果,實驗中選擇F2等邊色散棱鏡用于搭建線性波數光譜儀.由圖4(b)可知當F2色散棱鏡頂角角度為60?時,光柵與棱鏡間的最優化旋轉角為21.8?.

在最優化結構參數情況下,理論評估線性波數光譜儀波數線性度的結果如圖5所示.當深度坐標z取0.2 mm,線陣相機上探測的光譜信號分別為波長線性分布和波數線性分布情況下,波數與像素位置關系曲線以及根據(5)式模擬的干涉光譜信號如圖5(a)所示,紅色曲線表示波長呈線性分布時的波數-像素關系曲線及對應的干涉光譜信號,藍色曲線表示波數呈線性分布時的波數-像素關系曲線及對應的干涉光譜信號.圖5(b)所示為模擬的基于兩種光譜儀即線性波長光譜儀和線性波數光譜儀的譜域OCT系統在不同成像深度處的點擴散函數對比圖.如圖所示,在較淺成像深度處,非線性波數分布對軸向分辨率的影響不大,但隨著成像深度的增加,非線性波數情況下的點擴散函數顯著加寬.根據(8)式可以計算優化的波數分布函數和理想的波數分布函數之間的相對偏差函數ε(x).根據計算結果δε=0.045%低于臨界值1/N=0.1%,表明優化的線性波數光譜儀的波數線性度表現良好.

圖4 (a)不同角度組合下的點擴散函數模擬圖;F2棱鏡(b),BK7棱鏡(c)以及SF11棱鏡(d)情況下評價因子數值模擬三維輪廓圖Fig.4.(a)The simulated PSF under different combinations of angles;3D numerical simulation results of the cost factor for F2(b),BK7(c)and SF11(d)prism.

圖5 (a)深度坐標為0.2 mm處的干涉光譜信號的模擬結果和波數線性度模擬結果,紅色曲線對應于線性波長分布,藍色曲線對應于線性波數分布;(b)線性波長干涉信號與線性波數干涉信號經FFT后的點擴散函數圖像,彩色曲線對應線性波長情況,黑色曲線對應線性波數情況Fig.5.(a)Comparison of the simulated spectral interference signal at the depth coordinate of 0.2 mm(red,linearin-wavenumber case;blue,linear-in-wavelength case);(b)point spread functions directly Fourier transformed from the linear-in-wavelength signal(color curves)and linear-in-wavenumber signal(black curves).

4.2 實驗結果

為了通過實驗評估搭建的線性波數光譜儀的性能,在譜域OCT系統中以平面反射鏡作為樣品采集并處理了幾組干涉光譜信號.在進行直流項處理和FFT后,可以測得系統靈敏度、靈敏度下降曲線以及軸向分辨率.這些實驗結果顯示在圖6中.從圖6(a)可以看出,基于譜域OCT線性波數光譜儀的靈敏度在零光程附近處約為91 dB.在深度范圍1.2 mm內,靈敏度下降低于6 dB.

圖6 (a)基于譜域OCT系統的線性波數光譜儀靈敏度;(b)三種靈敏度下降曲線對比;(c)軸向分辨率的對比;(d)手指皮膚與指甲接縫處橫斷面OCT圖像Fig.6.(a)The sensitivity of linear-k spectrometer based SDOCT system;(b)sensitivity roll-of fcomparison;(c)the comparison of axial resolution;(d)cross-sectional OCT images of human finger skin obtained by the constructed linear-k spectrometer based SDOCT system.

為了把靈敏度下降的實驗結果和理論模型、模擬結果相比較,引用文獻[17]提出的譜域干涉光譜信號解析模型.此解析模型可以表達如下:

其中A是干涉光譜信號強度,p是像素尺寸,R是線色散率,即光譜聚焦平面上1μm長度所占的光譜帶寬,a是光譜聚焦光斑大小,z是成像深度.靈敏度下降對比曲線如圖6(b)所示,其中紅線代表實驗測得的系統靈敏度下降曲線,藍線代表在相機上光斑尺寸為6.8μm時用解析模型計算的理論靈敏度下降曲線,黑線代表數值模擬情況下最優化線性波數光譜儀靈敏度下降曲線,由代入最優化結構參數的波數函數求得.如圖所示,實驗測得的靈敏度下降性能和解析模型、數值模擬的結果基本一致.軸向分辨率的實驗結果及其與理論值的對比如圖6(c)所示.在圖中,理論值是由OCT軸向分辨率公式δz=0.44λ2/?λ決定的數值;優化值代表模擬點擴散函數的半峰全寬值;理想值代表在波數函數和像素位置成線性關系情況下模擬的點擴散函數的半峰全寬值;實驗值代表用實驗室搭建的基于線性波數光譜儀的譜域OCT系統測得的點擴散函數的半峰全寬值.在實驗中,光源的半峰全寬為38 nm,理論軸向分辨率為8.14μm.當平面鏡位于1.6 mm深度處,優化的軸向分辨率為8.45μm,等波數分布理想軸向分辨率為8.36μm,實驗測得軸向分辨率9.15μm.由圖中曲線可以看出,當深度較小時,實際軸向分辨率值接近于理論值;當深度較大時,兩者存在一定差距,推測主要是由于線性波數光譜儀搭建誤差引起.最后,對搭建的基于線性波數光譜儀的譜域OCT系統進行活體生物組織成像性能測試.圖6(d)展示的是對人手指指甲皮膚接縫處組織進行成像的OCT圖像結果,其橫向和軸向視場范圍分別為5.6 mm和2.4 mm.從圖中可見,手指指甲組織、皮膚組織、指甲和皮膚接縫處組織結構清晰可見.高分辨率手指皮膚組織的活體成像結果證明了搭建的基于線性波數光譜儀的譜域OCT系統的成像性能.

5 結 論

提出了一個優化準則來優化線性波數光譜儀的結構參數.需要考慮的結構參數有棱鏡材料、棱鏡頂角和光柵與棱鏡間的旋轉角.根據優化結果,選取F2玻璃等邊色散棱鏡、以光柵-棱鏡間旋轉角角度為21.8?搭建了最優化線性波數光譜儀,并引入譜域OCT成像系統.實驗測得成像系統的成像深度達到了2.42 mm,軸向分辨率達到8.52μm,靈敏度達到91 dB,6 dB成像深度達到1.2 mm.該線性波數光譜儀大大改善了光譜線性度,并在不需要重新采樣的情況下即能獲得高分辨率OCT圖像.利用GPU的并行運算能力實時顯示了橫斷面OCT圖像.

[1]Huang D,Swanson E A,Lin C P,Schuman J S,Stinson W G,Chang W,Hee M R,Flotte T,Gregory K,Puliafito C A,Fujimoto J G 1991Science254 1178

[2]Fercher F A,Hitzenberger C K,Kamp G,Elzaiat S Y 1995Opt.Commun.117 43

[3]Hausler G,Lindner M W 1998J.Biomed.Opt.3 21

[4]Cho H S,Jang S J,Kim K,Dan A V,Shishkov M,Bouma B E,Oh W Y 2013Biomed.Opt.Express5 223

[5]Wang R K,Zhang A Q,Choi W J,Zhang Q Q,Chen C L,Miller A,Gregori G,Rosenfeld P J 2016Opt.Lett.41 2330

[6]Liang Y M,Zhou D C,Meng F Y,Wang M W 2007Acta Phys.Sin.56 3246(in Chinese)[梁艷梅,周大川,孟凡勇,王明偉2007物理學報56 3246]

[7]Jia Y Q,Liang Y M,Zhu X N 2007Acta Phys.Sin.56 3861(in Chinese)[賈亞青,梁艷梅,朱曉農2007物理學報56 3861]

[8]Huang L M,Ding Z H,Hong W,Wang C 2011Acta Phys.Sin.60 023401(in Chinese)[黃良敏,丁志華,洪威,王川2011物理學報60 023401]

[9]Leitgeb R,Hitzenberger C K,Fercher A F 2003Opt.Expresss11 889

[10]Choma M A,Sarunic M V,Yang C,Izatt J A 2003Opt.Express11 2183

[11]Brauer B,Murdoch S G,Vanholsbeeck F 2016Opt.Lett.41 5732

[12]Zhang M,Hwang T S,Campbell J P,Bailey S T,Wilson J D,Huang D,Jia Y 2016Biomed.Opt.Express7 816

[13]Photiou C,Bousi E,Zouvani I,Pitris C 2017Biomed.Opt.Express8 2528

[14]Chen J B,Zeng Y G,Yuan Z L,Tang Z L 2018Acta Opt.Sin.38 0111001(in Chinese)[陳俊波,曾亞光,袁治靈,唐志列2018光學學報38 0111001]

[15]Gao W R,Chen Y D,Liu C,Zhang T Q,Zhu Y 2016Acta Opt.Sin.45 0611001(in Chinese)[高萬榮,陳一丹,劉暢,張秋庭,朱越2016光學學報45 0611001]

[16]Bao W,Ding Z H,Wang C,Mei S T 2013Acta Phys.Sin.62 114202(in Chinese)[鮑文,丁志華,王川,梅勝濤2013物理學報62 114202]

[17]Hu Z L,Pan Y S,Rollins A M 2007Appl.Opt.46 8499

[18]Dorrer C,Belabas N,Likforman J P,Jo ff re M 2000J.Opt.Soc.Am.B17 1795

[19]Hu Z L,Rollins A M 2007Phys.Opt.Lett.32 3525

[20]Gelikonov V M,Gelikonov G V,Shilyagin P A 2009Opt.Spectrosc.106 459

[21]Watanabe Y,Itagaki T 2009J.Biomed.Opt.14 48

[22]Lee S W,Kam H,Joo H P,Tae G L,Eun S L,Jae Y L 2015J.Opt.Soc.Korea19 55

[23]Lan G P,Li G Q 2017Sci.Rep.7 75