針灸療效評估中激光散斑襯比圖像穩(wěn)定性研究

仇福康，陳曉榮，劉宏業(yè)，衛(wèi)曉翔，張凱凱，孔 平

（1．上海理工大學 光電信息與計算機工程學院，上海 200093；2．上海健康醫(yī)學院 文理教學部, 上海 201318）

引 言

中醫(yī)理療如今被越來越多的人關注與接受。中醫(yī)理療的關鍵是通過外在物理刺激，如針灸、拔罐、推拿等，作用于人體氣血淤積之處，進而使得人體的氣血得以通暢，疼痛得以緩解。中醫(yī)不同于以生物、物理、化學等為基礎的西醫(yī)，中醫(yī)的治療效果不能通過一些生理指標來直接評估，一般是通過患者的反饋加之主治醫(yī)師的醫(yī)療經(jīng)驗得出綜合評判。近期的研究表明，中醫(yī)理療的效果可以從皮膚表層的微循環(huán)的變化得以觀測[1]，而皮膚微循環(huán)的具體情況需要現(xiàn)代醫(yī)學儀器來探究。激光多普勒[2]是常用的一種應用于監(jiān)測微循環(huán)系統(tǒng)的技術手段，現(xiàn)在市面上已經(jīng)存在開發(fā)成熟的激光多普勒血流儀。激光多普勒的作用原理是：在激光的單點照射下，生物的被測區(qū)域中流動的紅細胞會改變散射回去的激光的頻率，這些回到接收端的激光的頻率存在頻率偏移[3]，通過數(shù)學分析這些偏移量可以推導出被測血流的速度分布情況。激光多普勒具有侵入深度深、成像清晰的優(yōu)勢，但不足之處是單點掃描采集相應的頻移信息，其成像速度過于緩慢[4]。

隨著研究的進一步深入，檢測效果優(yōu)于激光多普勒單點掃描成像的激光散斑襯比成像[5](laser speckle contrast imaging, LSCI)的技術標準得以制定。相對于激光多普勒，LSCI具有成像速度快、成像維度為二維以及成像空間分辨率高等優(yōu)點[6]，這些成像優(yōu)勢使得LSCI可用來實時監(jiān)測針灸過程中被測對象皮膚表層的微循環(huán)。LSCI的原理是：在全場照明下，被測區(qū)域微循環(huán)中流動的紅細胞會產(chǎn)生散射激光的相移[7]，在一定曝光時間內(nèi)不同相移的激光在接收端的像平面產(chǎn)生不同程度的干涉，分析不同程度的干涉所形成的散射圖像即可得出被測微循環(huán)區(qū)域的血流速度分布。LSCI的分析是基于散斑圖像空間信息的數(shù)理統(tǒng)計，定義了反映血液中紅細胞流動速度信息的襯比值K[8]。在實驗研究中發(fā)現(xiàn)，激光散斑襯比成像的襯比值K的穩(wěn)定性會隨著光源調(diào)整而變化，為了提高激光散斑襯比成像中流速信息反饋的準確性，對影響散斑成像的襯比值及光源進行測試與分析，并進一步改進散斑成像的襯比值的處理算法，以提高激光散斑成像在針灸效果評估中的穩(wěn)定性。

1 激光散斑襯比成像理論與處理算法改進

1.1 成像理論基礎

LSCI的理論基礎是隨機波動的激光光強在CCD相機設定的曝光時間內(nèi)在探測器上形成的不同的光強積分，即成像面上對應位置產(chǎn)生隨機分布的灰度值。在設定的曝光時間內(nèi)，激光光源保持穩(wěn)定輸出，成像面的灰度值分布反映了被測區(qū)域的運動情況，靜止的部分顯示為黑白分明的散斑圖案，運動部分則根據(jù)速度的快慢顯示為模糊程度不一的灰度圖案，該灰度圖案稱之為散斑圖[9]。定義散斑襯比K的數(shù)學表達式為

式中：β是一個常系數(shù)，其值由光強探測器探測單元的尺寸與散斑尺寸的比值決定；σI和表示空間滑動窗口內(nèi)的光強分布的標準差與均值；x是CCD相機曝光時間T與光強波動的自相關時間τ的比值，其中光強波動的自相關時間τ與被測區(qū)域內(nèi)散射介質(zhì)運動速度v成反比。

1.2 處理算法改進

激光散斑成像理論是基于散斑圖像灰度值的空間數(shù)理統(tǒng)計，散射激光的干涉分布的散斑襯比K可表示為

式中Ii為空間窗口內(nèi)某點的光強。散斑襯比K的定義是，在N×N的空間滑動窗口內(nèi)，激光光強的標準差與激光光強的均值的比值。圖1顯示了激光散斑空間算法的基本操作流程：滑動窗口（一般選取5×5大小）即在散斑圖像上自左至右掃描，分別計算掃描過的區(qū)域內(nèi)的圖像灰度值的標準差與均值的比值，分別記為K11、K12、K13、…；以此類推，掃描完散斑圖像的第一行，圖1中的滑動窗口下移一行，完成對散斑圖像第二行的掃描，所得襯比分別記為K21、K22、K23、…。

圖1 激光散斑空間算法的基本原理Fig.1 Basic principle of laser speckle space algorithm

由式(2)以及圖1所示的激光散斑空間算法，處理激光在CCD像面上形成的包含被測區(qū)域內(nèi)散射介質(zhì)的運動信息的灰度值分布圖，即在一定曝光時間內(nèi)，運動的散射介質(zhì)形成散射激光相移，進而干涉形成光強的積分分布圖?；诳臻g散斑算法的基本公式，計算襯比分布的耗時隨著CCD采集的圖片像素的增大而增加。針灸療效的檢測系統(tǒng)需具備快速數(shù)據(jù)處理功能，Tom等發(fā)表的快速處理算法能有效縮短散斑運算的耗時，記為Roll算法[10]，算法的基本表達式為

采用Roll算法能顯著減少運算時間，例如采用 Intel? CoreTMi7-4790K CPU@4.00 GHz、內(nèi)存16 GB、Microsoft windows 7、操作系統(tǒng)為64位的計算機，在散斑圖像大小為1040×1392、滑動窗口設置為5×5的情況下，激光散斑基本算法的耗時為76.512 s，而Roll算法耗時為6.027 s，運算效率明顯提升。因此在采用Roll算法的基礎上，可以實現(xiàn)針灸過程中皮膚表層血流微循環(huán)的實時監(jiān)測。研究發(fā)現(xiàn)，雖然運算速度得到提升，但是由于激光光源的非理想性，以及激光散斑本身具有隨機波動的特性[11]，因此待測區(qū)域內(nèi)的靜止部分仍會隨著時間進行隨機的、細微的波動。將激光散斑應用到針灸治療效果的評估時，被測區(qū)域的相對靜止部分能得到一個比較穩(wěn)定的散斑襯比，使得針灸刺激區(qū)域的細微反應可以被觀察、記錄與分析。

為使針灸療效檢測系統(tǒng)達到實時顯示、輸出穩(wěn)定的需要，在Roll算法的基礎上，添加了基于時間的數(shù)理統(tǒng)計分析，以期得到穩(wěn)定的襯比輸出。改進的流程如圖2所示，記此算法為Roll with time-sampling，簡稱RT。采用與Roll算法相同的計算機配置進行實驗，在散斑圖像大小為1040×1392、滑動窗口設置為5×5的情況下，RT算法耗時為 6.192 s。相對于 Roll算法，RT算法的運算效率存在細微的下降，這是因為載入多組數(shù)據(jù)進行對比運算會產(chǎn)生延時。

圖2 RT 提升散斑圖像襯比穩(wěn)定的流程Fig.2 RT enhances the stability of the speckle image contrast

RT的處理思路是，通過基于快速數(shù)據(jù)處理的Roll算法，加入5個抽樣時刻t1、t2、t3、t4、t5，對應的襯比值為k1、k2、k3、k4、k5，抽樣間隔分別為 Δt1、Δt2、Δt3、Δt4，總間隔為ΔT=Δt1+Δt2+Δt3+Δt4。

針灸效果評估要求實時顯示，故設定Δt1=5 s、Δt2=10 s、Δt3=15 s、Δt4=25 s，這樣既可以得到穩(wěn)定化的襯比值，又可以滿足測量系統(tǒng)對實時性的要求。

2 針灸測量系統(tǒng)與系統(tǒng)參數(shù)優(yōu)化

2.1 針灸測量系統(tǒng)

基于LSCI成像系統(tǒng)的針灸測量系統(tǒng)主要由激光光源輸出、圖像采集、數(shù)據(jù)分析3個模塊構成，如圖3所示，其中，區(qū)域1為選自背景部分的靜態(tài)對照區(qū)域，區(qū)域2為以毫針刺入的位置為中心的刺激區(qū)域。激光經(jīng)光束擴展裝置照射到針灸作用的區(qū)域，圖中橢圓形區(qū)域為CCD相機的有效視場。激光經(jīng)被測區(qū)域后，以不同角度散射并被成像器件捕獲，成像器件包括一定帶寬的濾波片以及面陣CCD相機。數(shù)據(jù)分析用LabVIEW與MATLAB混合編制的分析程序來完成。

圖3 基于激光散斑襯比成像的針灸測量系統(tǒng)Fig.3 Acupuncture measurement system based on laser speckle contrast imaging

2.2 針灸測量系統(tǒng)參數(shù)優(yōu)化

針灸測量系統(tǒng)所得襯比穩(wěn)定性的影響因素主要為測量系統(tǒng)參數(shù)以及散斑成像系統(tǒng)參數(shù)，其中：散斑成像系統(tǒng)參數(shù)包括散斑大小[12]、散斑空間處理的滑動窗口大小[13]；成像系統(tǒng)參數(shù)包括CCD的靶面大小[14]、激光光源的功率[15]。隨著實驗的進一步深入，激光的波長作為一個新的影響因素被列入考察范圍。

實驗是基于國內(nèi)外認可的空間散斑襯比分析方法，測量被測試者手部微循環(huán)在針灸刺激下的變化情況。測量系統(tǒng)激光發(fā)生器選用了4種，分別為He-Ne激光器(波長632.8 nm，功率30 mW)、半導體激光器1(波長785.0 nm，功率50 mW)、半導體激光器2(波長532.0 nm，功率300 mW)、半導體激光器3(波長671.0 nm，功率300 mW)，通過對比實驗篩選出能產(chǎn)生最佳襯比的激光光源。激光發(fā)射的光束經(jīng)由兩片具有一定間隔的透鏡所組成的光束擴散系統(tǒng)，發(fā)散形成激光光柱照射到被測區(qū)域，實驗對象為平整A4紙。

實驗分為4組進行，通過變量控制及添加激光光強衰減器，將激光功率統(tǒng)一調(diào)制為30 mW。通過分析襯比的數(shù)學統(tǒng)計特性來評判襯比的穩(wěn)定性，其中，選取襯比的方差為穩(wěn)定性的判斷指標，實驗結果見表1。

表1 激光光源測試結果歸納Tab.1 Summary of laser source test results

通過分析對比各光源下襯比值的方差大小，可以得出：采用半導體激光器1可使襯比穩(wěn)定性最佳。因此在針灸效果評估的實驗中，采用波長為785.0 nm的半導體激光器作為光源。

3 實驗結果與討論

圖4是針灸實驗激光散斑圖，實驗環(huán)境溫度為26 ℃，測試者靜坐于實驗臺前，右手放在被測區(qū)域內(nèi)，保持靜止狀態(tài)。實驗時長共22 min，分為針灸前2 min、針灸進行中10 min、針灸后10 min。實驗采取的抽樣間隔為1 min，每次采集5張圖片，圖片采集間隔依次分為5 s、10 s、15 s、25 s。

圖4 針灸實驗激光散斑圖Fig.4 Original laser picture of acupuncture experiment

分別通過Roll算法及RT算法處理測得的激光散斑圖片，所得的襯比偽彩圖如圖5所示。圖5中，(a)表示Roll算法處理針灸前第1個抽樣時刻的原始散斑圖所得襯比偽彩圖，(d)表示RT算法處理針灸前第1個抽樣時刻的原始散斑圖所得襯比偽彩圖；(b)表示Roll算法處理針灸過程中第4個抽樣時刻的原始散斑圖所得襯比偽彩圖，(e)表示RT算法處理針灸過程中第4個抽樣時刻的原始散斑圖所得襯比偽彩圖；(c)表示Roll算法處理針灸后第5個抽樣時刻的原始散斑圖所得襯比偽彩圖，(f)表示RT算法處理針灸后第5個抽樣時刻的原始散斑圖所得襯比偽彩圖。偽彩圖區(qū)間設置為0(藍色)～70(紅色)。

圖5 Roll和 RT 算法處理所得襯比偽彩圖對比Fig.5 Comparison pseudo color pictures of Roll and RT algorithm

圖5中：標記“2”的區(qū)域為針灸刺入人體后人體皮膚表層血流微循環(huán)的觀察區(qū)域；標記“1”的區(qū)域為靜態(tài)參考點，用來評價算法處理原始散斑圖所得襯比圖分布的穩(wěn)定程度；經(jīng)過Roll算法處理所得的(a)、(b)、(c)3幅偽彩圖之間出現(xiàn)了比較明顯的整體的數(shù)據(jù)波動；相對而言，經(jīng)過RT算法處理所得的(d)、(e)、(f)3幅偽彩圖之間波動比較平穩(wěn)。區(qū)域2的變化情況能直觀表現(xiàn)針灸的治療效果，在RT算法的處理下，將圖像的整體抖動進行了優(yōu)化，進而可以比較準確地反映出區(qū)域2的襯比值波動情況。

圖6給出了時長為22 min的手背針灸測試結果穩(wěn)定性對比圖，其中：“+”表示RT算法處理原始散斑圖像區(qū)域1數(shù)據(jù)后所得襯比均值隨時間的分布曲線；“□”表示Roll算法處理原始散斑圖像區(qū)域1數(shù)據(jù)后所得襯比均值隨時間的分布曲線；“◇”表示Roll算法處理原始散斑圖像區(qū)域2數(shù)據(jù)后所得襯比均值隨時間的分布曲線；“○”表示RT算法處理原始散斑圖像區(qū)域2數(shù)據(jù)后所得襯比均值隨時間的分布曲線。

理想狀態(tài)下，激光光源保持持續(xù)的穩(wěn)定輸出，測試系統(tǒng)所用的背景板是處于絕對靜止狀態(tài)，此時，靜止背景部分的襯比值為1(實際上由于被測試人員的手部在測試過程中不能保持長時間的相對靜止，即存在肌肉的無意識舒張運動，會使底端的背景板產(chǎn)生細微的抖動)。襯比的分布區(qū)間為0～1，數(shù)值越小表示對應區(qū)域內(nèi)運動的散射介質(zhì)運動越快；同理，數(shù)值越大表示對應區(qū)域的散射介質(zhì)運動越慢。經(jīng)過RT算法優(yōu)化處理之后，由圖6可見：2 min之后(0～2 min為針灸前)針灸區(qū)域2開始呈現(xiàn)襯比減小的趨勢；3～12 min為針灸中，針灸區(qū)域2保持襯比減小趨勢并達到極小值；13～22 min為針灸后，針灸區(qū)域2出現(xiàn)襯比增加的趨勢，并最終回歸到靜止狀態(tài)時的襯比分布值。相對而言，Roll算法處理區(qū)域2的原始數(shù)據(jù)后，襯比分布中存在明顯的波動，在針灸前第2個時間點，針灸中第4、6、10個時間點，針灸后的第4個時間點均產(chǎn)生異常襯比，這表明了Roll算法處理原始散斑數(shù)據(jù)所得襯比的穩(wěn)定性欠佳。參考靜態(tài)對照區(qū)域1的襯比分布，相對于Roll算法，RT算法處理所得襯比的波動程度明顯小很多，相對平穩(wěn)，襯比波動值為1.73%(襯比隨時間分布的標準差)，小于3.47%(激光散斑近期研究中定義的襯比波動穩(wěn)定臨界值[16])，因此相對于Roll算法，RT算法能在一定程度上改善襯比分布的穩(wěn)定性。

圖6 Roll和 RT 處理所得襯比隨時間分布Fig.6 Distribution of speckle contrast over time obtained by Roll and RT algorithm

4 結 論

通過光源對照組的選擇以及相應的散斑襯比值K的方差分析比較可知，波長785.0 nm、功率50 mW的半導體激光器產(chǎn)生的散斑襯比穩(wěn)定性最優(yōu)。采用對散斑襯比穩(wěn)定性能最優(yōu)的近紅外激光光源進行實驗，對時長為22 min的實驗數(shù)據(jù)分別進行Roll算法和RT算法處理。通過對比分析Roll算法、RT算法處理靜態(tài)背景區(qū)域1和針灸區(qū)域2的原始散斑圖像數(shù)據(jù)所得的襯比隨時間分布狀況，可以得出：RT算法處理所得的襯比分布具有更佳的去散斑圖像抖動的效果。