针灸疗效评估中激光散斑衬比图像稳定性研究

仇福康，陈晓荣，刘宏业，卫晓翔，张凯凯，孔 平

（1．上海理工大学 光电信息与计算机工程学院，上海 200093；2．上海健康医学院 文理教学部, 上海 201318）

引 言

中医理疗如今被越来越多的人关注与接受。中医理疗的关键是通过外在物理刺激，如针灸、拔罐、推拿等，作用于人体气血淤积之处，进而使得人体的气血得以通畅，疼痛得以缓解。中医不同于以生物、物理、化学等为基础的西医，中医的治疗效果不能通过一些生理指标来直接评估，一般是通过患者的反馈加之主治医师的医疗经验得出综合评判。近期的研究表明，中医理疗的效果可以从皮肤表层的微循环的变化得以观测[1]，而皮肤微循环的具体情况需要现代医学仪器来探究。激光多普勒[2]是常用的一种应用于监测微循环系统的技术手段，现在市面上已经存在开发成熟的激光多普勒血流仪。激光多普勒的作用原理是：在激光的单点照射下，生物的被测区域中流动的红细胞会改变散射回去的激光的频率，这些回到接收端的激光的频率存在频率偏移[3]，通过数学分析这些偏移量可以推导出被测血流的速度分布情况。激光多普勒具有侵入深度深、成像清晰的优势，但不足之处是单点扫描采集相应的频移信息，其成像速度过于缓慢[4]。

随着研究的进一步深入，检测效果优于激光多普勒单点扫描成像的激光散斑衬比成像[5](laser speckle contrast imaging, LSCI)的技术标准得以制定。相对于激光多普勒，LSCI具有成像速度快、成像维度为二维以及成像空间分辨率高等优点[6]，这些成像优势使得LSCI可用来实时监测针灸过程中被测对象皮肤表层的微循环。LSCI的原理是：在全场照明下，被测区域微循环中流动的红细胞会产生散射激光的相移[7]，在一定曝光时间内不同相移的激光在接收端的像平面产生不同程度的干涉，分析不同程度的干涉所形成的散射图像即可得出被测微循环区域的血流速度分布。LSCI的分析是基于散斑图像空间信息的数理统计，定义了反映血液中红细胞流动速度信息的衬比值K[8]。在实验研究中发现，激光散斑衬比成像的衬比值K的稳定性会随着光源调整而变化，为了提高激光散斑衬比成像中流速信息反馈的准确性，对影响散斑成像的衬比值及光源进行测试与分析，并进一步改进散斑成像的衬比值的处理算法，以提高激光散斑成像在针灸效果评估中的稳定性。

1 激光散斑衬比成像理论与处理算法改进

1.1 成像理论基础

LSCI的理论基础是随机波动的激光光强在CCD相机设定的曝光时间内在探测器上形成的不同的光强积分，即成像面上对应位置产生随机分布的灰度值。在设定的曝光时间内，激光光源保持稳定输出，成像面的灰度值分布反映了被测区域的运动情况，静止的部分显示为黑白分明的散斑图案，运动部分则根据速度的快慢显示为模糊程度不一的灰度图案，该灰度图案称之为散斑图[9]。定义散斑衬比K的数学表达式为

式中：β是一个常系数，其值由光强探测器探测单元的尺寸与散斑尺寸的比值决定；σI和表示空间滑动窗口内的光强分布的标准差与均值；x是CCD相机曝光时间T与光强波动的自相关时间τ的比值，其中光强波动的自相关时间τ与被测区域内散射介质运动速度v成反比。

1.2 处理算法改进

激光散斑成像理论是基于散斑图像灰度值的空间数理统计，散射激光的干涉分布的散斑衬比K可表示为

式中Ii为空间窗口内某点的光强。散斑衬比K的定义是，在N×N的空间滑动窗口内，激光光强的标准差与激光光强的均值的比值。图1显示了激光散斑空间算法的基本操作流程：滑动窗口（一般选取5×5大小）即在散斑图像上自左至右扫描，分别计算扫描过的区域内的图像灰度值的标准差与均值的比值，分别记为K11、K12、K13、…；以此类推，扫描完散斑图像的第一行，图1中的滑动窗口下移一行，完成对散斑图像第二行的扫描，所得衬比分别记为K21、K22、K23、…。

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

由式(2)以及图1所示的激光散斑空间算法，处理激光在CCD像面上形成的包含被测区域内散射介质的运动信息的灰度值分布图，即在一定曝光时间内，运动的散射介质形成散射激光相移，进而干涉形成光强的积分分布图。基于空间散斑算法的基本公式，计算衬比分布的耗时随着CCD采集的图片像素的增大而增加。针灸疗效的检测系统需具备快速数据处理功能，Tom等发表的快速处理算法能有效缩短散斑运算的耗时，记为Roll算法[10]，算法的基本表达式为

采用Roll算法能显著减少运算时间，例如采用 Intel® CoreTMi7-4790K CPU@4.00 GHz、内存16 GB、Microsoft windows 7、操作系统为64位的计算机，在散斑图像大小为1040×1392、滑动窗口设置为5×5的情况下，激光散斑基本算法的耗时为76.512 s，而Roll算法耗时为6.027 s，运算效率明显提升。因此在采用Roll算法的基础上，可以实现针灸过程中皮肤表层血流微循环的实时监测。研究发现，虽然运算速度得到提升，但是由于激光光源的非理想性，以及激光散斑本身具有随机波动的特性[11]，因此待测区域内的静止部分仍会随着时间进行随机的、细微的波动。将激光散斑应用到针灸治疗效果的评估时，被测区域的相对静止部分能得到一个比较稳定的散斑衬比，使得针灸刺激区域的细微反应可以被观察、记录与分析。

为使针灸疗效检测系统达到实时显示、输出稳定的需要，在Roll算法的基础上，添加了基于时间的数理统计分析，以期得到稳定的衬比输出。改进的流程如图2所示，记此算法为Roll with time-sampling，简称RT。采用与Roll算法相同的计算机配置进行实验，在散斑图像大小为1040×1392、滑动窗口设置为5×5的情况下，RT算法耗时为 6.192 s。相对于 Roll算法，RT算法的运算效率存在细微的下降，这是因为载入多组数据进行对比运算会产生延时。

图2 RT 提升散斑图像衬比稳定的流程Fig.2 RT enhances the stability of the speckle image contrast

RT的处理思路是，通过基于快速数据处理的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，这样既可以得到稳定化的衬比值，又可以满足测量系统对实时性的要求。

2 针灸测量系统与系统参数优化

2.1 针灸测量系统

基于LSCI成像系统的针灸测量系统主要由激光光源输出、图像采集、数据分析3个模块构成，如图3所示，其中，区域1为选自背景部分的静态对照区域，区域2为以毫针刺入的位置为中心的刺激区域。激光经光束扩展装置照射到针灸作用的区域，图中椭圆形区域为CCD相机的有效视场。激光经被测区域后，以不同角度散射并被成像器件捕获，成像器件包括一定带宽的滤波片以及面阵CCD相机。数据分析用LabVIEW与MATLAB混合编制的分析程序来完成。

图3 基于激光散斑衬比成像的针灸测量系统Fig.3 Acupuncture measurement system based on laser speckle contrast imaging

2.2 针灸测量系统参数优化

针灸测量系统所得衬比稳定性的影响因素主要为测量系统参数以及散斑成像系统参数，其中：散斑成像系统参数包括散斑大小[12]、散斑空间处理的滑动窗口大小[13]；成像系统参数包括CCD的靶面大小[14]、激光光源的功率[15]。随着实验的进一步深入，激光的波长作为一个新的影响因素被列入考察范围。

实验是基于国内外认可的空间散斑衬比分析方法，测量被测试者手部微循环在针灸刺激下的变化情况。测量系统激光发生器选用了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)，通过对比实验筛选出能产生最佳衬比的激光光源。激光发射的光束经由两片具有一定间隔的透镜所组成的光束扩散系统，发散形成激光光柱照射到被测区域，实验对象为平整A4纸。

实验分为4组进行，通过变量控制及添加激光光强衰减器，将激光功率统一调制为30 mW。通过分析衬比的数学统计特性来评判衬比的稳定性，其中，选取衬比的方差为稳定性的判断指标，实验结果见表1。

表1 激光光源测试结果归纳Tab.1 Summary of laser source test results

通过分析对比各光源下衬比值的方差大小，可以得出：采用半导体激光器1可使衬比稳定性最佳。因此在针灸效果评估的实验中，采用波长为785.0 nm的半导体激光器作为光源。

3 实验结果与讨论

图4是针灸实验激光散斑图，实验环境温度为26 ℃，测试者静坐于实验台前，右手放在被测区域内，保持静止状态。实验时长共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个抽样时刻的原始散斑图所得衬比伪彩图。伪彩图区间设置为0(蓝色)～70(红色)。

图5 Roll和 RT 算法处理所得衬比伪彩图对比Fig.5 Comparison pseudo color pictures of Roll and RT algorithm

图5中：标记“2”的区域为针灸刺入人体后人体皮肤表层血流微循环的观察区域；标记“1”的区域为静态参考点，用来评价算法处理原始散斑图所得衬比图分布的稳定程度；经过Roll算法处理所得的(a)、(b)、(c)3幅伪彩图之间出现了比较明显的整体的数据波动；相对而言，经过RT算法处理所得的(d)、(e)、(f)3幅伪彩图之间波动比较平稳。区域2的变化情况能直观表现针灸的治疗效果，在RT算法的处理下，将图像的整体抖动进行了优化，进而可以比较准确地反映出区域2的衬比值波动情况。

图6给出了时长为22 min的手背针灸测试结果稳定性对比图，其中：“+”表示RT算法处理原始散斑图像区域1数据后所得衬比均值随时间的分布曲线；“□”表示Roll算法处理原始散斑图像区域1数据后所得衬比均值随时间的分布曲线；“◇”表示Roll算法处理原始散斑图像区域2数据后所得衬比均值随时间的分布曲线；“○”表示RT算法处理原始散斑图像区域2数据后所得衬比均值随时间的分布曲线。

理想状态下，激光光源保持持续的稳定输出，测试系统所用的背景板是处于绝对静止状态，此时，静止背景部分的衬比值为1(实际上由于被测试人员的手部在测试过程中不能保持长时间的相对静止，即存在肌肉的无意识舒张运动，会使底端的背景板产生细微的抖动)。衬比的分布区间为0～1，数值越小表示对应区域内运动的散射介质运动越快；同理，数值越大表示对应区域的散射介质运动越慢。经过RT算法优化处理之后，由图6可见：2 min之后(0～2 min为针灸前)针灸区域2开始呈现衬比减小的趋势；3～12 min为针灸中，针灸区域2保持衬比减小趋势并达到极小值；13～22 min为针灸后，针灸区域2出现衬比增加的趋势，并最终回归到静止状态时的衬比分布值。相对而言，Roll算法处理区域2的原始数据后，衬比分布中存在明显的波动，在针灸前第2个时间点，针灸中第4、6、10个时间点，针灸后的第4个时间点均产生异常衬比，这表明了Roll算法处理原始散斑数据所得衬比的稳定性欠佳。参考静态对照区域1的衬比分布，相对于Roll算法，RT算法处理所得衬比的波动程度明显小很多，相对平稳，衬比波动值为1.73%(衬比随时间分布的标准差)，小于3.47%(激光散斑近期研究中定义的衬比波动稳定临界值[16])，因此相对于Roll算法，RT算法能在一定程度上改善衬比分布的稳定性。

图6 Roll和 RT 处理所得衬比随时间分布Fig.6 Distribution of speckle contrast over time obtained by Roll and RT algorithm

4 结 论

通过光源对照组的选择以及相应的散斑衬比值K的方差分析比较可知，波长785.0 nm、功率50 mW的半导体激光器产生的散斑衬比稳定性最优。采用对散斑衬比稳定性能最优的近红外激光光源进行实验，对时长为22 min的实验数据分别进行Roll算法和RT算法处理。通过对比分析Roll算法、RT算法处理静态背景区域1和针灸区域2的原始散斑图像数据所得的衬比随时间分布状况，可以得出：RT算法处理所得的衬比分布具有更佳的去散斑图像抖动的效果。