基于特征譜線與約束擬合相位的絕對波數標定方法*

王迪 韓濤 錢黃河 劉智毅 丁志華

(浙江大學光電科學與工程學院,現代光學儀器國家重點實驗室,杭州 310027)

譜域光學相干層析成像(spectral-domain optical coherence tomography,SD-OCT)系統中普遍存在波數域的非線性采樣問題.為實現常規快速傅里葉變換算法下離散界面的精確定位與OCT 圖像的高質量重建,需要解決光譜儀中離散采樣點絕對波數的精確標定問題.本文提出了一種基于精確光程差下特征譜線與約束擬合相位的絕對波數標定方法,在譜域OCT 系統的樣品臂中,使用具有精確厚度差異的金屬量規,獲得特征譜線對應的絕對相位值,進一步精確求解出特征譜線對應的相位包裹次數,克服了常規干涉光譜相位方法中普遍存在的2π 歧義,結合窗口約束條件下高信噪比區域的擬合相位,實現光譜儀采樣點絕對波數的精確標定.通過全面比較本文方法與傳統插值重采樣方法在離散界面定位、軸向分辨率以及圖像重建質量等方面的差異,驗證了本方法的顯著優勢.

1 引言

光學相干層析成像(OCT)[1,2]是一種基于低相干干涉原理的三維斷層成像技術,可以實現生物組織內部微米量級的高分辨率實時成像,也可用于材料的無損檢測與精細結構的三維重建[3-7].譜域OCT 技術是當前應用最為廣泛的OCT 成像技術,通過光柵光譜儀探測并采集干涉光譜信息,相對于時域OCT 技術具有顯著的信噪比、靈敏度和成像速度優勢[8-10].

譜域OCT 技術中,波數域與深度域是傅里葉變換關系.對干涉光譜信號進行傅里葉變換,可以得到沿深度方向分布的反射系數,反映了樣品深度方向的結構信息.由于光柵光譜儀色散分光原理,譜域系統中普遍存在非線性采樣問題,具體表現為在線陣探測器等像素間隔位置上,采集的干涉光譜在波數域上的分布不均勻.為滿足傳統快速傅里葉變換(fast Fourier transform,FFT)算法對波數域均勻采樣的要求,需要對干涉光譜信號進行插值重采樣,以獲得波數域上均勻分布的干涉光譜[11],從而避免軸向點擴散函數的展寬以及靈敏度快速下降[12].而且,獲得均勻分布光譜的絕對波數值,可以實現深度域位置信息的精確定位.因此,光譜儀絕對波數的標定是實現譜域OCT 系統中離散界面精確定位及高信噪比、高軸向分辨率OCT 圖像重建的前提.Hu 和Rollins[13]首次在譜域系統中引入線性波數光譜儀,Wu 等[14]在此基礎上進行了光譜儀結構參數的優化,從而實現了均勻波數采樣,但不可避免地提升了硬件結構復雜度.Hyle Park 等[15]提出了基于光譜儀物理模型的標定方法,結合具體幾何參數計算各像素點的對應波長,但理論推導過程極為復雜.特征譜線多項式擬合方法是經典的波長標定方法,這類方法基于標準定標光源或物質的特征吸收光譜,可以建立波長-采樣像素對應關系,如Perret 等[16]利用法布里-珀羅干涉濾光片的等間距譜線,Eom 等[17]利用布拉格光柵的多條特征譜線,最終通過多項式擬合獲得采樣空間的波數標定結果,但這類方法仍受制于特征譜線數量、線寬以及分布情況.基于干涉光譜相位的波數標定方法無需額外硬件,通過計算特定光程差下干涉光譜相位與波數間的線性關系,可以實現光譜儀的標定,如Wang 和Ding[18]基于兩次干涉光譜相位相減方法消除色散影響,實現了光譜儀均勻分布波數的重采樣,但尚未考慮相位噪聲對波數標定的負面影響,而且,該方法無法直接確定采樣點的絕對波數.Wu 等[19]在共路譜域OCT 系統樣品臂中,使用不同厚度的金屬量規,求解第一個采樣像素的初始相位值,結合相對相位分布,恢復采樣空間的絕對相位,但其初始相位值的求解存在較大誤差,同時絕對相位中混入大量的相位噪聲,無法實現絕對波數的求解.

本文提出了一種基于精確光程差下特征譜線與約束擬合相位的絕對波數標定方法,在譜域OCT系統的樣品臂中使用具有精確厚度差異的金屬量規,準確求解特征譜線對應的絕對相位值的包裹次數,克服了常規干涉光譜相位方法中普遍存在的2π 歧義,并結合窗口約束條件下高信噪比區域的相位擬合,實現了采樣點空間絕對波數的精確標定.通過全面比較所提方法與傳統重采樣方法在離散界面定位、軸向分辨率以及圖像重建質量等方面的差異,驗證了本方法的顯著優勢.

2 理論方法

不失一般性,譜域OCT 系統探測得到單層離散反射面的干涉光譜信號,去除直流項和自相關項后可表示為

式中,ki為光譜儀采集的離散波數,i=1,2,···,N表示采樣序列;rS和rR分別表示參考反射面和樣品反射面的反射系數;S(ki)為光源功率譜分布函數;z為反射面與參考面光程差的一半;φ(ki,z)為對應的光譜域絕對相位離散分布,可以用非線性采樣函數g(k)=ki與兩干涉臂間不平衡色散項h(ki)來表示

其中,k為波數.對(1)式實施希爾伯特變換,可以得到限制在 [-π,+π] 主值區間內的包裹相位分布φw(ki,z)[20,21]:

式中,floor 表示向負無窮方向取整運算,不同的離散采樣波數ki有不同次數的相位包裹,存在 2π 歧義,即相位卷繞問題[22].

以起始采樣點為基準對包裹相位φw(ki,z)做連續化解包裹處理,可以得到相對相位分布,但一般的干涉光譜中心區域信噪比高(相位噪聲小)、邊緣區域信噪比低(相位噪聲大),邊緣區域的相位噪聲將在連續化解包裹過程中累積放大,從而影響恢復相位的準確性.因此,可采用以臨近中心位置kc的相位主值作為起點的雙向連續化解包裹處理,得到相對相位分布φr(ki,z):

式中,N=為固定值,表示臨近中心位置(對應波數kc)的相位包裹次數.為消除(2)式中兩干涉臂間的不平衡色散,采集不同光程下兩組干涉信號并分別提取相位,進而獲得精確光程差(2Δz=2(z2-z1))下對應的相位差分布:

實際情況下,當特征譜線處于光譜儀臨近中心區域時,可進行光譜儀精確標定.以He-Ne 激光特征譜線為例,對應的絕對相位為 2kHeNeΔz,可以計算出對應相對相位的包裹次數N2-N1,表示為

式中,round 運算表示取最接近的整數值,當特征譜線kHeNe與光程差 2Δz引入的誤差相位不超過π時,可以求得準確的絕對相位包裹次數.由此恢復得到光程差 2Δz下所有采樣光譜點的絕對相位分布為

考慮到相位噪聲分布 δφ(i)對光譜采樣的影響,實際系統得到的相對相位分布為=φr(ki,Δz)±δφ(i),與理論值具有一定偏差,最終影響絕對相位恢復結果以及波數標定的準確性.常見的基于全采樣序列光譜相位φ-波數k的最小二乘多項式擬合方法,可以在一定程度上去除相位噪聲的干擾,但不能避免中心區域相位的整體偏移.Han 等[23]在最小二乘擬合過程中,加入過中心采樣點擬合的約束條件,改善了中心區域的相位偏移情況,但仍無法剔除邊緣噪聲的影響.因此,在最小二乘擬合過程中,加入擬合窗口約束條件,可以充分結合采樣中心區域高信噪比、低相位噪聲的特點,得到更接近高信噪比相位理論分布的擬合相位.在采樣中心區域加入窗口寬度為 2L的約束條件,進行絕對相位φa(ki,Δz)的部分擬合,可以求得最小二乘多項式的擬合系數{c0,c1,c2,c3,...},并以此恢復得到對應光程差 2Δz下的全采樣序列擬合相位:

最小二乘擬合算法獲得的均方偏差會隨著窗口約束條件 [N/2-L,N/2+L] 的改變而改變,通過比較不同約束條件下,高信噪比區域的相位均方偏差值,最終可以確定最佳窗口約束擬合條件.此時獲得的擬合相位在高信噪比區域貼合于采樣恢復相位,在采樣邊緣位置則符合光譜域相位平滑變化趨勢,相比于恢復絕對相位φa(ki,Δz),更接近于理想相位分布.

在獲得特征絕對相位的基礎上,結合特征譜線波長數值,可以精確求解實際的精確光程差=φa(kHeNe,Δz)/kHeNe,繼而用于精確絕對波數分布的標定.在最佳窗口約束相位擬合條件下,得到絕對波數分布,均勻化處理后獲得重采樣向量分布k＇[i],依次表達為

3 實驗系統與標定過程

如圖1 所示,在已有超高分辨譜域OCT 系統[24]基礎上,本實驗系統進行部分改進.超連續光源發出光譜帶寬約為 200 nm 的可見低相干光,經2 倍擴束器擴束后,平行入射至分光比為1∶1 的分束器上(BS,Thorlabs),透射光部分經 30 mm 透鏡聚焦到參考臂平面鏡上,反射光部分則通過掃描振鏡(GVS012,Thorlabs)掃描,經消色差透鏡(f=30 mm,Edmund optics)聚焦到樣品上;最終,兩束返回光重新匯聚,經擴束、50 μm 小孔濾波后入射至衍射光柵上(600 線/mm,Wasatch Photonics),形成寬光譜色散,最終被像素數2048、行頻250 kHz的線陣互補金屬氧化物半導體相機(CMOS,OCT OPLUS,Teledyne E2v)采集,實現了最高250 fps(frames/s)成像速度(對應于橫向1000 個采樣點數)的快速成像.此外,引入基于He-Ne 激光器的特征譜線(632.8 nm)的光譜標定模塊,并在樣品臂成像位置處設置兩組金屬量規(量規厚度分別為d1=1350 μm ,d2=1250 μm,厚度偏差|δd|＜0.1 μm)作為高反射樣品來實施光譜儀的波數標定.兩組金屬量規對應的雙臂光程差分別為 2z1和2z2,且兩量規的上表面位置對稱分布于譜域OCT 系統的零光程附近,滿足關系z1＜0＜z2.分別采集兩組量規上表面與參考臂反射鏡返回光的干涉信號,提取相位并通過相減獲取光譜位相差,以消除色散失配.該光譜相位差分布可認為是由虛擬“空氣隙”間的干涉所致,且虛擬“空氣隙”的光程為量規的厚度差(2Δz=2(d1-d2))所導致.

圖1 用于光譜儀標定實驗的譜域OCT 系統Fig.1.Schematic diagram of the spectral domain optical coherence tomography system for spectrometer calibration.

根據(7)式求解光譜絕對相位及(6)式求解相位包裹次數的方法,在高信噪比相對相位值φr(kHeNe,Δz)確定的條件下,“空氣隙”對應的絕對位相差分布的準確求解,依賴于設定光程差 2Δz和特征譜線值kHeNe的精確性.本實驗選用國標量規標稱值與實際值間存在不大于 0.1 μm 的厚度偏差,對應“空氣隙”光程差小于 0.2 μm,由此造成的相位誤差kHeNe·2δd≈0.63π＜π ;同時,由于本系統光柵光譜儀具有較高的光譜分辨能力(約 0.2 nm),由特征譜線偏差造成的相位偏差 δkHeNe·2Δz可以忽略不計.因此,結合特征譜線的絕對相位差標定方法,可以準確求解絕對相位差的包裹次數,并準確恢復“空氣隙”對應的絕對相位差分布.同時,根據特征絕對相位恢復結果φa(kHeNe,Δz),結合特征譜線波長對“空氣隙”實際光程進行反解,得到與標稱值的光程偏差,滿足標定要求.

圖2 展示了利用He-Ne 特征譜線恢復“空氣隙”絕對相位差的過程.圖2(a)中He-Ne 激光特征譜線出現在CMOS 相機像素序列i=1025 處,其半高全寬小于一個像素間隔;圖2(b)展示了兩組金屬量規樣品對應的互相關干涉信號,紅色與藍色曲線分別對應了負光程差z1與正光程差z2下的干涉信號;依照前述方法處理獲得正負光程下的連續化相對相位φr(ki,zi)如圖2(c),兩者呈現相反的變化趨勢;最終,根據構造的“空氣隙”相對相位φr(ki,Δz),結合(6)式計算得到“空氣隙”對應的相位包裹次數為N2-N1=315,由此恢復的絕對相位曲線如圖2(d)所示.

圖2 He-Ne 特征譜線恢復“空氣隙”絕對相位的過程(a)He-Ne 特征譜線標定光譜;(b)兩組金屬量規的互相關干涉信號;(c)連續化處理后的相對相位分布;(d)恢復的“空氣隙”絕對相位分布Fig.2.Process of recovering the absolute phase of “air gap” from He-Ne characteristic spectral line:(a)Spectrum of He-Ne characteristic line calibration;(b)cross correlation interference signals of two groups of metal gauges;(c)relative phase distribution after continuity;(d)recovered absolute phase of “air gap”.

圖3 展示了上述絕對相位φa(ki,Δz)的相位差波動,其中左邊緣區域、中心區域以及右邊緣區域的局域放大圖分別用不同顏色的線框圈出;采樣邊緣區域相位波動大(可達 0.1 rad),具有較高相位噪聲,而采樣中心區域相位波動小(僅約 0.01 rad),對應了更高的信噪比,印證了本方法的可行性.

圖3 “空氣隙”采樣絕對相位的相位差波動及不同區域相位差波動情況Fig.3.Phase difference fluctuations of absolute phase of “air gap” and phase difference fluctuations in different sampling regions.

在恢復采樣絕對相位的基礎上,控制窗口約束條件進行絕對相位的多項式擬合,獲得平滑的擬合絕對相位分布.圖4 展示了最佳窗約束條件下的擬合絕對相位,以及不同窗約束條件下的相位均方偏差大小.以中心采樣點為擬合窗口中心位置,依次改變窗口寬度進行絕對相位最小二乘多項式擬合,得到了不同窗約束條件下擬合相位和恢復相位之間的相位均方偏差,圖4(b)描述了相位均方偏差隨約束窗口變化的規律,其最小值對應了最佳約束條件和最佳擬合窗口位置(擬合區域[216,1832]),由此得到的窗口約束下最佳擬合相位如圖4(a)中所示,在邊緣位置相位擬合曲線偏離采樣相位值,反映算法對低信噪比區域的相位進行了擬合校正.最終,根據(9)式的計算過程,獲得光譜儀采樣空間的絕對波數分布.

圖4 最佳窗口約束條件下的擬合絕對相位及不同窗約束條件下相位標準偏差分布曲線(a)最佳擬合絕對相位分布;(b)不同窗約束下相位標準偏差變化曲線Fig.4.Absolute phase fitting curve under optimum window constraints and phase standard deviation curve under different window constraints:(a)The optimal fitting of the absolute phase distribution;(b)phase standard deviation curve under different window constraints.

4 結果與討論

4.1 離散反射界面的成像實驗結果

基于最佳窗口約束條件獲得的擬合絕對相位,可以用來計算采樣軸上的絕對波數分布,繼而得到線性化后的波數序列,以此為基礎進行干涉光譜信號的重采樣、消除色散及快速傅里葉變換,可以獲得相比傳統重采樣方法更佳的成像結果.

圖5 展示了譜域OCT 系統使用不同處理算法獲得的軸向點擴散函數(PSF),包括直接快速傅里葉變換方法(direct FFT)、傳統波數域插值重采樣方法(k-domain spline interpolation without fit,KDSI),以及基于不同窗約束條件下相位擬合的插值重采樣方法(window-constrainedk-domain spline interpolation,WC-KDSI),所有的點擴散函數均以本文提出方法(最佳窗約束條件下的WC-KDSI算法,WC-KDSI(Best fitting))的點擴散函數峰值為基準進行了歸一化處理.其中,直接FFT 方法的軸向點擴散函數極大展寬,證明重采樣以及消除色散過程具有必要性;傳統的KDSI 算法進行了波數域的直接插值以及色散去除,其峰值信號強度及點擴散函數的半高全寬較直接FFT 方法有很大提高,但與WC-KDSI 方法相比則遜色不少;在最佳窗約束條件下,WC-KDSI 方法的軸向點擴散函數質量相較于其他窗約束條件也略有所提高.表1詳細表述了上述各方法獲得的離散界面點擴散函數的半高全寬及歸一化強度差異,本文提出方法獲得點擴散函數的半高全寬即軸向分辨率可達 1.6 μm,高于傳統的KDSI 方法的 2.1 μm 及直接FFT方法的 5.3 μm,同時其歸一化強度最大,深度域信號的信噪比更高.

圖5 不同方法獲得的軸向點擴散函數比較Fig.5.Axial PSFs obtained by different methods.

表1 不同方法獲得的軸向分辨率和峰值歸一化強度Table 1. Axial resolution and peak intensity resulted from different methods.

研究發現,傳統KDSI 方法不僅在軸向點擴散函數的半高全寬、峰值強度上與本文提出方法具有一定差距,針對單反射面峰值位置的精確定位,也存在離散采樣位置差異.為了驗證所提出方法相較傳統KDSI 方法在單界面精確光程定位上具有一定優勢,對標定過程中金屬量規“空氣隙”的峰值離散位置和峰間離散采樣間隔進行了定位.圖6 中分別使用KDSI 方法和WC-KDSI 方法處理兩金屬量規離散界面的干涉信號,獲得KDSI 方法下兩量規峰間像素采樣間隔為n3-n1=163 ,根據dz=π/(Ndk)=0.62 μm,可求實際光程差為OPDKDSI=2(n3-n1)dz=201.1 μm;同時,WC-KDSI 方法下兩量規峰間采樣間隔為n4-n2=162,對應實際光程差為 OPDWC-KDSI=2(n4-n2)dz=199.9 μm,實驗中使用的兩組量規具有精確的光程差為 200 μm,其預設誤差小于 2δd=0.2 μm .顯然,本文提出方法相較于傳統的KDSI 方法在離散界面的精確定位方面,具有突出優勢;進行光譜域信號多倍補零處理后,可以獲得更小的深度域采樣間隔,兩種方法間的定位精度偏差將更加凸顯.

圖6 經KDSI 方法和WC-KDSI 方法處理得到的量規反射面定位結果Fig.6.Positioning results of gauges’ interfaces processed by KDSI method and WC-KDSI method.

4.2 鏡面腐蝕部位的成像實驗結果

為了實際說明提出方法在實際三維形貌重建方面具有一定優勢,本文首先針對鏡面腐蝕位置的微形貌特征進行了實際成像,分別使用直接FFT方法、KDSI 方法及WC-KDSI 方法處理了三維圖像數據.處理獲得的OCT 重建圖像如圖7 所示,圖中選用了同一掃描位置的B-scan 進行展示,為了提高圖像的信噪比,使用相鄰10 幅B-scan 層析圖像進行了實際處理[25].圖7(a)中,鏡面點擴散函數展寬,所得重建圖像模糊,且三維重建結果不能分辨瑕疵細節;圖7(b)及圖7(c)分別使用KDSI方法及WC-KDSI 方法對表面腐蝕的微形貌進行了重建,腐蝕部分以及鏡面部分均具有較高的信噪比,所得C-scan 均可以清晰分辨鏡面腐蝕情況,但在腐蝕邊緣WC-KDSI 方法擁有更好的細節分辨能力,更能展現邊緣的凹凸以及腐蝕中心的起伏情況;選擇兩處橫向掃描位置(X1=250,X2=125)并采集一維軸向信號,得到如圖7(d)中的A-scan信號分布,可見在不同橫向掃描位置,WC-KDSI方法獲得的重建圖像相比傳統KDSI 方法獲得的重建圖像,均具有更高的分辨率和信噪比;同時通過讀取其PSF 的半高全寬,可以知道在探測深度100 μm 左右,使用WC-KDSI 方法可以分辨深度小于 4 μm 的鏡面瑕疵.

圖7 使用(a)直接FFT,(b)KDSI,(c)WC-KDSI 算法獲得的鏡面腐蝕微形貌OCT 重建圖像;(d)不同掃描位置軸向PSF 比較Fig.7.Reconstructed OCT images of mirror corrosion by(a)direct FFT,(b)KDSI,(c)WC-KDSI method;(d)comparison of axial PSFs at different scanning points.

4.3 橘子果肉樣本的成像實驗結果

針對植物或生物等高散射組織,本文提出方法在三維圖像高質量重建方面同樣具有一定優勢.選用橘子果肉樣本,依次使用KDSI 方法及WC-KDSI方法進行相關處理,獲得成像結果見圖8.其中,圖8(a)和圖8(b)展示了KDSI 方法和WC-KDSI方法處理的果肉二維層析圖像;由于系統光源譜段處于可見光波段,針對橘子果肉樣本等高散射樣本的穿透能力不足,因此獲得的層析圖像成像深度有限;但是,在果肉淺層深度,重建圖像仍能較為清晰地分辨橫向及縱向的果肉細胞壁組織,整體視野下WC-KDSI 方法擁有更高的對比度、分辨和細節分辨能力.圖8(c)選取不同方法處理下、同一掃描位置的A-line 信號進行對比,在圖中框選部分信號,相比于KDSI 算法處理結果(藍色曲線),WC-KDSI 算法處理結果(紅色信號)在多個果肉細胞壁位置均擁有更高的分辨率及信號強度,信號峰更突出,同時KDSI 處理信號的噪聲普遍大于WC-KDSI 方法.綜上,以上實驗說明本文提出方法在高質量圖像重建方面具有一定優勢.

圖8 使用(a)KDSI 算法及(b)WC-KDSI 算法處理獲得的橘子果肉OCT 重建圖像;(c)特定掃描位置的軸向PSF比較Fig.8.Reconstructed OCT images of orange pulp by(a)KD SI and(b)WC-KDSI algorithm;(c)comparison of axial PSFs at certain scanning position.

5 結論

在譜域OCT 系統中,光柵光譜儀的精確標定是快速傅里葉變換算法處理干涉光譜數據的前提,基于光譜域相位的光譜儀標定方法可以獲得采樣軸上的波數相對分布,因此被廣泛采用.針對譜域干涉信號在采樣邊緣區域的信噪比不高、相位噪聲較大的特點以及采樣過程中的非線性特征,本文提出了基于特征譜線及約束擬合相位的光譜絕對波數標定方法,利用高信噪比區域的部分相位來擬合整個采樣軸上的相位分布,結合特征譜線的特征相位值精確計算相位包裹次數,從而實現了采樣軸上絕對相位的求解及絕對波數的精準標定.本文通過反射界面點擴散函數實驗比較說明本文提出方法相比傳統方法具有更高分辨率和更高信噪比,通過對比不同方法下兩組量規離散界面的深度域采樣位置和信號峰之間的采樣間隔,說明了本文提出方法在離散界面精確定位方面具有優勢;通過對反射鏡鏡面腐蝕微形貌、橘子果肉樣本進行三維圖像重建,詳細分析軸向點擴散函數及層析圖像信息,說明了本文提出方法可以實現更高質量的圖像重建.