估计点扩展函数对内源性光学成像信号的提取

尹 愚 宋雪梅 尧德中（电子科技大学生命科学与技术学院，成都 60054）

2（中国科学院上海生命科学研究院，上海 200031）

引言

脑内源性信号的光学成像（optical imaging based on intrinsic signals）是从20世纪70年代发展起来的一种脑成像技术，是研究脑的工作机制的重要手段之一。同其它研究手段相比，内源性光学成像在时间和空间分辨率上都有较大优势。在哺乳动物的大脑内部，执行特定功能或有相同特性的细胞聚集在一起［1-2］。这些细胞发生兴奋活动时，皮层神经活动的改变会导致相应区域的血流速度、去氧血红蛋白（Hbr）浓度等血液参数的变化，进而引起其所在皮层对特定波长光吸收率的变化。这些变化是依赖于脑活动的内源性信号，显示了大脑皮层的功能构筑［3］。在活体脑组织中基于这种血氧动力学而产生的光学信号，其强度变化非常微弱，相对于全部反射光大约占0.1％～0.2％（605nm光照条件）［4］。

因为进行研究的脑功能区域一般在皮层以下（深度为200～800μm），光学信号在透射出皮层的过程中，大脑组织会对光产生散射作用。这使得记录到的功能区域会比真实活动的皮层范围有所扩散［4］。另外在图像生成和传输过程中的各种干扰也会对图像产生模糊作用，而这种模糊过程一般可以看作是系统点扩展函数与图像卷积的结果。因此寻找一个对点扩展函数准确的估计，并利用反卷积方法就能对这种退化图像进行恢复。对于内源性光记录信号的去噪处理，主要还是采用空间平滑滤波。本研究根据内源性光学信号图像二阶统计量的分析，利用其基于生理特性的分布形式对点扩展函数进行合理的估计，并恢复出功能柱图像。相对于恢复后的图像，图像噪声得到了有效抑制，功能柱结构清晰，而且功能区域更为集中。

1 实验方法和数据处理

1.1 动物准备和光学成像

实验采用体重2.0～3.0kg的成年猫3只，性别不限。肌肉注射30mg/kg盐酸氯胺酮诱导麻醉，施行气管插管和股静脉插管手术。将动物俯卧固定于立体定位仪（SN-3，NARISHIGE）；持续静脉注射3mg/（kg·h）戊巴比妥钠和10mg/（kg·h）三碘季胺酚，以维持实验过程中动物的麻醉和麻痹状态；150mg/（kg·h）葡萄糖溶液为动物提供必要的营养物质。人工呼吸机（CIV-100，Cloumbus Instruments）设置为20次/min，潮气量约为12.3mL/kg。动物体温维持仪（SS20-1，蚌埠市大江电子公司）将被试动物体温控制在38℃～38.5℃。实验过程中监视动物心电（ECG）和脑电（EEG），以了解动物的麻醉状态和生理情况。根据图谱在Horsley-Clarke坐标系的P0-P10、L0-L6处，开颅暴露皮层17区，面积约为5mm×5mm，剥离硬脑膜。在暴露皮层的颅骨窗外围用牙科水泥封固一个金属小室，使用甲基硅油填充满小室后加玻璃盖片封闭以维持适当压力，从而减少心跳和呼吸导致的皮层表面机械波动［5］。动物眼睛加滴1％阿托品扩瞳，5％新福林收缩瞬膜。然后佩戴接触镜，用以矫正视力和防止角膜干燥。用检眼镜将动物眼底血管投影于视觉刺激显示屏上，以确定视盘位置。在接触镜外放置直径为3mm的人工瞳孔。

高灵敏度数字CCD摄像机［6］（iXon897，ANDOR）在波长550nm（半波宽为30nm）的绿光照射条件下，将镜头聚焦于皮层表面，拍摄出皮层表面的血管图像。然后改为使用波长为605nm（半波宽为10nm）的红光进行照明，并通过对CCD电动支架调节，使得镜头聚焦于皮层表面以下约500μm处，以减弱血管的干扰，进行脑功能光学图像的采集。

1.2 视觉刺激和信号采集

视觉刺激通过一台刺激服务器连接一台显示器给出。显示器置于猫眼前57cm处，刺激范围为30°×40°视角［7］。视觉刺激的平均亮度为20cd/m2，对比度为100％。运动的正弦光栅作为视觉刺激，其时间频率为2Hz，空间频率为0.8c/deg，刺激朝向为0°～157.5°，间隔22.5°共八种。高灵敏度CCD摄像机正对暴露皮层区域，拍摄速度为20帧/s，空间分辨率为16μm，景深约为0.1mm。每个刺激周期内视觉刺激呈现前持续拍摄2s（共计拍摄40帧），刺激呈现后持续拍摄8s（共计拍摄160帧）。不同朝向的光栅在不同的记录周期内随机分别呈现，实验重复记录10次。视觉刺激为单眼刺激，即仅暴露一只眼，遮盖另一只眼进行刺激。

1.3 数据预处理

由于被探测皮层不可能处于理想的均匀光照下。为获得皮层活动的图像，首先将皮层受到刺激时获取的图像与基准图像进行对比，以修正光照的不均匀性。将实验中全部刺激呈现前所拍摄的图像进行叠加平均作为皮层的基线图像（baseline image），这被称为空白本底（blank）。将多次刺激的全部图像进行叠加平均的目的是减弱生物体生理活动和光照系统的低频波动所引入的噪声，多次试验的叠加平均可以将图像中每个像素的标准差减小为噪声标准差的1/（K是进行叠加平均的图像数目）。对内源性光学信号图像进行标准分析，将某种刺激条件下，获得的全部拍摄图像进行叠加平均，然后与本底图相除。而后使用“第一帧分析”技术（first frame analysis），去除低于0.3Hz的低频噪声（主要是呼吸等生物体活动引入的干扰）。即将每帧图像都减去未刺激时所采集的第一帧图像。图1（b）显示了实验动物1受朝向为0°刺激时，通过上述两项预处理后所得到的结果。对照于图1（a）的直接拍摄所得图像，从中可见存在以簇状成团方式呈现的功能活动区。在图像中暗色区域是有着较强的光吸收区域（相对于本底图），即为有着较强活动的功能柱区域。通过多次重复实验的叠加平均，能将各种低频周期波动所引入的噪声进行一定程度的压制。但是因为大脑组织对光的散射，内源信号会发生一定程度的扩散。所以图中所显示的暗色区域边缘部分并非完全是皮层的相关活动功能区，有部分是因为光学弥散而产生的伪迹［4］。

图1 内源光学信号。（a）CCD采集图像；（b）预处理结果Fig.1 Opticalimaging based on intrinsic signals.（a）CCD image；（b）function image

1.4 二阶统计分析

二阶统计量既方差所表示的是偏离平均值的误差情况，因此可以通过计算像素点（m，n）邻域S内的二阶统计量来评价像素点间的差异情况。

邻域S范围选择192μm×192μm的空间大小，实验动物1的处理结果如图2（a）所示。其中暗色区域为f*（m，n）值较小的区域，在这些位置上像素点（m，n）与其相邻区域S内的像素间的差异较小。可以观察到图1（b）中最暗和最亮的区域都对应着图2（a）中的最暗区域。这说明在功能图1（b）中，因视觉刺激而激活的皮层区域内所产生的光学信号的相对变化程度的一致性较高。这是由于激活的皮层区域内大量功能相同的神经元发生放电活动，放电活动代谢能量而引起血氧浓度变化，使得对照射光的吸收程度增强。在有着较强活动的功能柱区域内，功能相似的神经元的活动所产生的光学信号也相似，因此这些位置上图像二阶统计处理结果的值也较小。这样的结果是由功能柱的结构和生理活动所决定的。图2（a）结果中显示在未激活的皮层区域也同样保持着较小的二阶统计量。这是因为在非活动的皮层区域神经元未被激活，其生理活动也一致的保持为较低水平，所以这个区域的皮层对照射光吸收程度的变化差异不大。因此在功能区域和非功能区域内的二阶统计量较为一致，其差异主要体现为均值的不同。白色区域表示这些区域的f*（m，n）值较大，说明像素点之间的差异较大。对应于图1（b），这些位置是处于功能柱和未激活皮层之间。这些区域并非与刺激相关的神经元集中的位置，因此并没有较强的放电活动代谢能量而引起血氧浓度变化。在这些位置光学信号的变化更可能是由于附近位置血氧浓度扩散和不完全透明的脑组织造成的弥散结果。在这些部位，靠近功能柱区域的颜色偏暗，而靠近未激活皮层区域的颜色偏亮。图2（b）为图像二阶统计量归一化值的分布直方图。

图2 功能图二阶统计结果。（a）二阶统计结果；（b）二阶统计量分布直方图Fig.2 Second-order statistics results.（a）the second-order statistics image；（b）the secondorder statistics histogram.

1.5 点扩展函数与估计

内源光学信号十分微弱，而且预处理后的图像对比度不高。为尽可能利用所采集的图像信息，对功能柱区域的信号部分与弥散部分进行分离，需要对所采集的图像进行恢复处理。弥散信号部分可以认为是信号在皮层中的扩散结果，这被视为是图像的一种模糊过程，使得图像产生了退化。通常使用以下卷积模型来模拟模糊过程

式中，g为模糊图像；h为点扩展函数（PSF）；f为原图像；n为随机噪声（通常被假设为具有均值为0，方差为σ2的高斯分布）；*表示卷积运算。如果可以得知点扩展函数的形式，就能通过反卷积恢复出真实的图像。

因为无法得知点扩展函数的真实情况，为恢复出原始图像，需要在合理的假设条件下对点扩展函数进行估计。在许多光学测量系统和成像系统模型中，一般都假定点扩展函数形式为高斯函数［8］。一种典型的高斯点扩展函数模型为

式中，σ为高斯函数的方差。要直接对二维PSF进行估计比较困难和复杂［9］，根据1.3中光记录预处理的描述可知经过预处理的图像在照明上可以认为是大体上均匀的，因此该PSF可以假设为具有圆对称特性［10］。设x方向上的线扩展函数（LSF）为h*（x），并且采用式（3）的形式，于是有

这样将二维PSF的估计转化为一维LSF的估计。

对于PSF的估计需要选取合理的考察区域，因为功能柱图中没有明确的边缘特征，传统的阶跃边缘法［9，11］和边缘检测法［12］等利用边界信息进行PSF估计的方法并不适用。根据上文对内源性光学图像的二阶统计分析，在图像中共有三种活动模式：功能柱区、弥散区和非功能区。假设活动模式是具有高斯分布的形式而进行的，并且认为在预处理中其它噪声已被得到有效的处理，可以暂不考虑其影响。所以将每个区域活动形式与LSF相卷积可以得到

式中，fa（x）为功能柱区的活动函数，fu（x）为非功能区的活动函数，fd（x）为弥散区的活动函数，其形式均为高斯函数。根据高斯函数的卷积规则，由式（5）直接可得

式中，σa、σu分别是功能柱区、和非功能区的真实活动方差。由于弥散区中既有功能柱区的活动也有非功能区的活动，反映的是两个活动之和，所以可以用两者的高斯函数相加来描述。根据高斯函数相加的原则，可知相加后的函数方差应为（σa+σu）2。σh为需要估计的LSF的方差。σga、σgd、σgu为模糊图像中这三个区域的方差。从功能图像中这三个区域取若干样本，即可以由式（6）解出σh。实际上此选取不同区域的样本有所不同，得到的σh也不会完全相同。可以根据如下准则：在线扩展函数h*（x）定义域的任何子集Ω中，整合后h*（x）的能量应该与所有样本得到的（x）的能量的平均值相同［13］，即

基于上述准则，可以从整个图像的全部样本点中估计出与图像唯一对应的h*（x）。

2 结果

2.1 点扩展函数估计结果

利用对内源信号功能图二阶统计分析的结果，结合预处理功能图的信号与非信号区域，针对实验动物1，共选取200组特征差异明显的像素点作为样本。通过式（6）和式（7）得到对该组数据的PSF的估计，结果如图3所示。

图3 点扩展函数估计结果Fig.3 Estimation of point spread function

2.2 反卷积恢复

根据对点扩展函数的估计结果，对实验动物1的预处理图像结果图1（b）进行反卷积图像恢复，图4（a）为功能图像恢复结果。相比预处理图像，图中的噪声得到了有效的抑制，同时保留的图像细节。更为重要的是二阶统计量分布直方图发生了明显变化。预处理图像的二阶统计量的分布（见图2（b））没有明显的峰值，很难从中选取一个阈值以划分二阶统计量较小区域和较大区域。在反卷积恢复后图像的二阶统计量分布图中（见图4（c）），其分布形式发生明显变化，左端的分布变得更为集中。可见在二阶统计值的较小端形成明显的峰值（图中用“*”标出），并且在峰的左右分布形状对称。而在其附近同时存在另外一个较为微弱的分布峰。相对于恢复前图像的直方图，图4（c）的主峰明显，之前未见的次峰也能可见。根据图像阈值划分的原理，可以在这两个分布峰之间的谷值处选取一个值以作为二阶统计量的大小的划分阈值，如图中“↓”标出。因为从预处理功能图和功能柱的生理机制可知，功能柱区域是一个功能相似的神经元集中区域，根据1.4节可知，共同相似的神经功能活动所引起的光学信号在二阶统计量上应表现为较小和分布集中。因此有效的功能区域在二阶统计量的分布上对应于其二阶统计值较小端的峰值。

图4 反卷积处理。（a）功能图恢复；（b）二阶统计量划分模板；（c）恢复图像的二阶统计量分布直方图Fig.4 Deconvolution results.（a）function image restoration；（b）edge template based on second-order statistics；（c）second-order statistics histogram of deconvolution result.

图4（b）中显示了利用该选取的阈值从反卷积恢复图像中得到分割模板结果。因为功能柱区域和非功能柱区域均属于二阶统计量值小的区域，所以该模板图中既包含功能柱区域，又包含非功能柱区域。结合反卷积恢复图可以准确划分激活功能柱区域。

2.3 功能柱区域分割

激活区域应满足两点：位于由二阶统计分析得到的分割模板的白色区域（二阶统计量小的区域）；位于反卷积恢复图中灰度低于全图平均值区域（激活区域）。据此标准对激活功能柱区域进行了划分，3只实验动物通过该方法所得结果分别表示于图5中。

图5 功能柱区域划分。（a）～（c）分别为3只实验动物结果Fig.5 Edge segmentation.（a）～（c）the results of the three cats，respectively

3 讨论

3.1 反卷积图像增强

光学信号在传递、采集过程中都会受到各种因素的干扰，使得图像产生模糊。这种模糊可以看作是系统点扩展函数与图像卷积的结果。因而如果可以预先知道系统的点扩展函数，就可以通过传统的图像复原算法对模糊图像进行逆运算，得到较清晰的复原图像。脑内源性光记录信号不仅非常微弱而且噪声成分复杂。传统处理手段主要以空间滤波为主，虽然处理结果使得图像更加平滑，图像的颗粒噪声也能得到抑制。内源信号光学成像去噪使用的平滑滤波是利用滤波核对图像进行卷积，实质是个模糊图像的过程。脑内源性光记录信号又具有弥散性的特点，生理结构的组织上没有明显边界。用空间滤波去增强图像会使得原本弥散的结果更加模糊，对功能柱区域划分产生干扰。利用脑内源性光记录信号的生理机制特点，以其二阶统计量的分布关系为基础所估计的点扩展函数，反卷积图像增强的结果良好。功能柱区域不再是不可分辨的一团模糊的暗色区域，而是形成具有类似等高线结果的阶梯区域。这样的结果更准确的反应了功能柱区域发生活动最强烈的位置。由于反卷积图像增强是去除图像的模糊过程，在恢复结果中使得功能柱区域更为集中，使弥散区域能产生一定程度的复原。

3.2 内源光学信号功能柱定位

脑内源性光记录研究工作需要对脑活动的功能柱区域进行定位和提取。特别是针对某些需要利用功能柱边界或者明确的功能柱区域信息所进行的研究。这类研究希望能通过内源性光学记录结果的引导，迅速准确的获取相关区域的信息。传统上多采取图像灰度阈值划分，但由于图像对比度低，不同情况下采集数据结果均有不同，而阈值选择主观性比较强。另外更是因为功能柱区域的弥散，生理结构无明显边界，这些都使得直接采用图像灰度阈值划分的结果无法准确定量，很难准确分割出功能柱区域。本研究中，利用功能柱活动的生理结构基础，根据相关活动区域的二阶方差应当较小的这一结果，通过反卷积恢复后图像弥散区域被有效去除的条件下，给出阈值选取的合理标准。最后结果显示，划分结果清楚，对血管等干扰的屏蔽作用强，基于这个方法下的功能柱划分是具有实用意义。

4 结论

本研究基于估计点扩展函数对内源性光学成像信号的提取，在对本实验室采集的内源性光记录信号的处理应用中得到了较好结果，证明了此方法适用于内源性光记录信号的图像增强和功能柱区域提取。并在一定程度上比直接采用图像灰度阈值划分更为合理。此方法对于内源性光记录信号在脑功能成像的研究，不仅提升了图像质量，更为明确的确定了功能区域。对于利用内源性光记录信号在神经活动的探索中具有进行挖掘信号的实用价值。在此后进一步研究中，计划结合电生理记录对本方法进行的功能柱区域划分结果进行单细胞记录的验证。

（致谢：感谢李朝义院士的大力支持与指导。感谢徐杏珍老师在动物实验中的认真指导与热情帮助。）

［1］Mountcastle VB.Modality and topographic properties of single neurons of cat＇s somatic sensory cortex［J］.Journalof Neurophysiology，1957，20（4）：408-434.

［2］HubelDH，Wiesel TN.Receptive fields and functional architecture in two nonstriate visual areas（18 and 19）of the cat［J］.J Neurophysiol，1965，28：229-289.

［3］Grinvald A.Functional architecture of cortex revealed by optical image of intrinsic signals［J］.Nature，1986，324：361-364.

［4］Grinvald A.In-vivo optical imaging of cortical architecture and dynamics［M］//Windhorst U，Johansson H，eds.Modern Techniques in Neuroscience Research.Berlin：Springer-Verlag，1999：893-969.

［5］俞洪波，邢大军，寿天德.脑内源信号光学成像术：猫视皮质方位功能柱的活体显示［J］.中国神经科学杂志，2000，16（4）：355-359.

［6］White BR，Bauer AQ，Snyder AZ，et al.Imaging of functional connectivity in the mouse brain［J］.PloS I，2011，6（1）：e16322.

［7］俞洪波，寿天德.用脑光学成像术研究不同空间拓扑位置猫初级视皮层的空间频率反应特性［J］.生理学报，2000，52：411-4151.

［8］吕成淮，何小海，陶青川，等.图像复原中高斯点扩展函数参数估计算法研究［J］.计算机工程与应用，2007，43（10）：31-34.

［9］李东兴，赵剡，许东.基于参数估计的降晰函数辨识及图像复原算法［J］.红外与激光工程，2010，1：166-172.

［10］Claxton CD，Staunton RC.Measurement of the point spread function of a noisy imaging system［J］.Journal of the Optical Society of America A：Optics，Image Science，and Vision，2008，25（1）：1592170.

［11］Vishwakumara K，Martens J B.Estimation of perceived image blur using edge features［J］.International Journal of Imaging Systems and Technology，1996，7：102-109.

［12］Elder JH，Zucker SW.Local scale control for edge detection and blur estimation［J］.IEEE Trans on Pattern Analysis and Machine Intelligence，1998，20（7）：699-716.

［13］吴俊芳，刘桂雄，韦宁.散焦含噪图像的点扩散函数估计与边缘检测［J］.华南理工大学学报（自然科学版），2010，38（6）：118-121.