基于局部Radon变换和卡尔曼滤波的超声图像肌束方向自动跟踪方法

温慧莹 杨晓娟 郭燕荣 张 帅 张志国 汪天富 陈思平 陈 昕*

1(深圳大学医学院生物医学工程系 医学超声关键技术国家地方联合工程实验室，广东省生物医学信息检测与超声成像重点实验室，深圳 518060)2(香港大学电气与电子工程系，香港 999077)

基于局部Radon变换和卡尔曼滤波的超声图像肌束方向自动跟踪方法

温慧莹1杨晓娟1郭燕荣1张 帅2张志国2汪天富1陈思平1陈 昕1*

1(深圳大学医学院生物医学工程系 医学超声关键技术国家地方联合工程实验室，广东省生物医学信息检测与超声成像重点实验室，深圳 518060)2(香港大学电气与电子工程系，香港 999077)

超声成像是一种无创、实时和便捷的成像方法，已被广泛用于人体肌肉运动分析。临床上目前主要采用手工处理的方法来提取超声图像中的形态结构，不仅依赖于操作者的主观经验，且耗时严重、重复性低。提出一种自动计算与跟踪超声图像序列中肌束方向的新方法，利用基于约束互信息的自由变换算法自动提取感兴趣区域，然后基于局部Radon变换对超声图像进行投影，并利用相邻帧间肌束的连贯性，采用卡尔曼滤波的图像处理方法。对6名试验对象进行步行实验，共采集小腿内侧腓肠肌的动态超声影像1 080帧，并用所提出的方法计算肌束方向。实验结果表明，所提出的自动检测方法有较好的鲁棒性。与人工检测结果相比，两者的平均相关系数为0.92±0.02，平均误差为0.30°±0.62°。该方法可以一定程度地替代人工测量，在处理大量超声数据时具有较大的优势。

超声成像；肌束方向；Radon变换；卡尔曼滤波

引言

肌肉是构成人体的重要组织，其主要功能为产生收缩和控制身体的运动。肌电图学(electromyography, EMG)是最早用于分析人体肌肉运动的方法，其原理是记录神经肌肉收缩和静止状态时的生物电信号。其中，表面肌电图因其实时、无创等优点，被广泛用于肌肉运动研究中[1-2]。但是，表面肌电图估计也存在一些不足，如对体内深层肌肉难以区别、对干扰比较敏感等，这制约了肌电信号在肌肉运动分析中的应用。另外，从上世纪90年代起，有学者开始利用超声评估肌肉的功能状态，并将分析结果应用在生物力学的研究领域[2-8]。超声成像是一种无创、实时和便捷的成像方法。二维超声图像可提供肌肉收缩时肌束的实时动态图像，已被广泛用于研究肌肉运动中形态结构的改变，为临床诊断和康复评估提供了途径。

利用超声可测量肌肉厚度(muscle thickness)[3, 6, 8-9]、羽状角(pennation angle)(或称肌束方向)[3, 5-8]、肌束长度(fascicle length)[3, 5-8]和肌肉横截面积(cross-section area)[4]等结构性参数来描述肌肉的状态变化，这些参数和肌肉组织的力学特性直接相关。其中，肌束方向描述肌纤维的排列，是决定肌肉力学性能的主要特征之一。在肌肉收缩过程中，羽状角越大，更多的平行排列且可收缩的肌纤维束将能量传递到肌腱上，使肌肉产生更大的力量。肌肉处于最大收缩状态时，羽状角的变化能达到静息时的120%～170%[10-11]。最近的研究还表明，肌肉肌束特性的改变与神经肌肉疾病相关[12-14]。目前，临床上主要是通过手工处理的方法来提取肌束方向信息。这不仅依赖于操作者的主观经验，在处理实时超声视频数据时耗时严重，而且需通过人工多次测量来提高重复性。

多种自动提取肌束方向的算法已被提出，根据方法大致可分为两类：Hough变换和Radon变换。Zhou等[15]首先提出了用迭代Hough变换(revoting hough transform, RVHT)，计算图像中的肌束方向。该方法通过将超声的边缘图像映射到霍夫空间，并将局部极值对应于该图像的直线位置，在检测线设定线宽范围内的所有特征点将被移除，重复上述过程，多条直线被依次提取，直至低于自定义阈值。然而，该方法极大地依赖边缘检测子的性能，容易受斑点噪声的干扰。此外，RVHT的检测结果高度依赖于所涉及迭代过程中的线宽参数，稍微不同的参数大小可能会导致完全不同的结果。临床中的肌肉超声图像，通常很难找到一个合适的线宽代表所有的肌肉纤维。Zhao[16]采用局部Radon变换(localized radon transform, LRT)，自动计算和跟踪肌束方向，不需要边缘检测子，受斑点噪声干扰小。然而，相同的迭代策略用于依次提取直线，仍存在相同的问题，如果初始检测不正确，后续难以恢复正确的肌束线结构[17]。最近文献[18]提出，利用局部梯度图来估计主要的肌束方向，而不是单独检测几根肌束方向，因此迭代策略是不必要的，并且更具稳健性，但是对局部方向上的每个像素的计算复杂度比[15-17]高很多。此外，上述方法存在两个缺点：一是需要在超声图像序列的第一帧图像上设置感兴趣区域(reading of interest, ROI)，对所有图像的肌束方向跟踪都是在该ROI范围中进行。然而，在真实的肌肉收缩过程中，深浅层筋膜区域是动态变化的，可能超过了固定ROI的范围，这将大幅度降低自动跟踪肌束方向的精度；二是在追踪过程中没考虑运动模型，因此容易受干扰，得到的结果曲线常出现毛刺。以往的研究都是针对静态或者准静态运动，干扰相对比较小。对步行运动，超声探头相对滑动等干扰因素比较严重，必须设法消除。

根据对现有方法的分析，本研究提出一种新的自动化计算与跟踪肌肉肌束方向的方法。与传统检测肌束直线的方法相比，该方法自动检测与跟踪ROI，利用Gabor滤波来自适应地增强超声图像，并应用局部Radon变换和卡尔曼滤波(Kalman filter, KF)来自动计算和跟踪肌束方向，提高了算法的鲁棒性和准确性。

1 方法

1.1 方法流程

研究方法的流程如图1所示。对连续的超声影像，首先基于约束互信息的自由变换算法(constrained mutual information-based free-form deforamtion, C-MI-FFD)，实现ROI区域的自动化选取，采用Gabor滤波抑制图像斑点噪声，并纵向增强图像肌束区域。此外，使用迟滞阈值法(hysteresis thresholding, HT)和图像细化(image thinning, IT)进一步增强图像，从而解决迭代策略中的线宽问题。利用文献[16]所述的方法来捕获肌束方向，并为KF提供校正测量。由于两个超声图像相邻帧之间的肌束具有连贯性，可取前一帧KF预测的肌束方向变化作为先验知识来追踪当前图像内的肌束。

图1 方法流程Fig.1 The diagram of the proposed method

1.2 基于约束互信息的自由变换算法

文献[4]首次提出基于约束互信息的自由变换算法(C-MI-FFD)，即在连续的超声图像序列中自动提取肌肉的横截面积。本研究扩展了该算法，自动选取连续超声图像序列深浅层筋膜间的感兴趣区域。通过最小化以互信息为基础的目标函数，确定描述连续两帧图像的转换函数。为了精确地匹配感兴趣区域，C-MI-FFD算法分为两个步骤。

步骤1:全局变换。通过全局匹配两帧图像确定全局变换参数，对它们相对比例、平移和旋转建模。

步骤2：局部变换。局部变换参数由二维样条函数定义，并被细化为规则网络下进行插值样条函数的位移值。图2(a)、(b)分别为典型的原始和自动获取ROI后的超声图像。

1.3 Gabor滤波器

由于肌肉超声图像具有连贯性，采用Gabor滤波可以抑制斑点噪声，并在纵向上增强肌束图像[19]。Gabor滤波器(Gabor filter, GF)为带通滤波器，目前广泛应用于图像特征提取、边缘检测、图像纹理分析等。二维Gabor函数的实数部分定义如下：

(1)

式中,λ为波长，θ为Gabor函数并行条纹的方向，φ为相位偏移，σ为高斯因子标准差，γ为空间纵横比，(x,y)为图像坐标。

图2(b)、(c)分别表示Gabor滤波前后的超声图像。可以看出，使用Gabor滤波后的图像斑点噪声被抑制，肌束区域图像被增强。值得注意的是，GF中的θ是自适应更新的、由前一帧KF预测的方向作为处理当前帧的并行条纹的方向。

为解决LRT迭代策略中的线宽问题和提升其检测精度，采用HT、IT两个技术进一步增强GF滤波后的图像信息。

步骤1：首先在HT中设有两个阈值，像素响应高于阈值上界的输出为1，而低于阈值下界的输出为0，像素响应位于阈值上界与下界之间的，且如果它们通过一条边缘链连接到其他像素响应均大于下界阈值的输出为1。

步骤2：图像细化。通过抑制非极大值部分的输出，可减少输出线结构的单像素边缘宽度，从而得到非常小的线宽，且LRT中重新表决的线宽问题得到有效解决。图2(d)、(e)分别表示HT和IT后的效果。

图2 方法过程追踪。(a)原始超声图像；(b)提取ROI(绿色线表示)；(c)Gabor滤波效果；(d)HT效果；(e)IT效果LRT得到的主要肌束方向用蓝色线表示Fig.2 Illustration of the proposed muscle fiber orientation tracking approach. (a) A raw ultrasound video frame; (b) Extracted ROI (green line surrounded area in (a)); (c) GF result; (d) HT result; (e) IT result. The dominant orientation obtained by using LRT is labeled by blue line

1.4 局部Radon变换和线提取

Radon变换已广泛用于直线特征的提取，标准Radon变换在2D变换空间内的定义为

(4)

式中,f(x,y)为在(x,y)位置的图像灰度，δ为狄拉克δ函数，ρ为直线到图像中心的距离，θ为x轴和线的夹角，D为图像网格。

标准Radon变换表示整个图像平面内的灰度积分。在具体应用中，常采用矩形来定义ROI，无法真实反映浅层和深层筋膜之间的区域的真实轮廓。如果只对预定义范围内部的位置和方向积分，标准Radon变换则可修改为LRT。运用先验知识如追踪ROI和肌束方向角度，LRT可以显著增强捕获肌束的准确度，降低边缘效应。在本研究中，采用C-MI-FFD算法可排除有强回声强度的浅层与深层筋膜，且同时最大程度地跟踪其余区域。当检测到指定的直线数量，或是低于峰值振幅与平均振幅的比值，表明这帧图像中所有直线已检测完毕。

1.5 卡尔曼滤波

卡尔曼滤波是一种高效率的递归滤波器，它能够从一系列的不完全及包含噪声的测量中估计动态系统的状态。目前广泛用于工程领域，如雷达、计算机视觉，近来更被应用于计算机图像处理，如人脸识别、图像分割、图像边缘检测等[20]。

正如上文描述，方向参数是Gabor滤波器和后续LRT直线特征提取的关键参数。笔者提出用KF自适应地更新方向的变化，即用LRT获得当前帧的肌束方向作为KF的测量值，KF输出的最优估计值用于Gabor滤波器处理下一帧图像中的方向参数的反馈输入，使达到更好的图像处理效果。KF的离散时间线性状态空间模型如下：

(5)

zk=Ckxk+Dkyk+vk

(6)

相对传统的LRT，使用KF后，跟踪肌束方向的精度得到提高，运行时间显著降低。

2 实验

选取试验对象6名，均为男性。年龄(24.6±1.5)岁，体重(69.4±6.3)kg，身高(170±5.6)cm。所有试验对象身体健康，均没有肌肉、神经类疾病病史。受试者均签署知情同意书。

使用迈瑞DC-6B超实时采集肌肉的动态超声图像，配备线阵超声探头，中心频率为10MHz。用一个定制的可调节的支架固定超声探头，使探头的长轴与试验对象小腿内侧腓肠肌的表面走向平行，并通过大量超声耦合剂使探头紧贴在目标肌肉表面。实验开始前，试验对象被要求在跑步机上做适应性步行，一般为2～3min。实验中，当试验对象准备好后，开启跑步机并且调节到合适的速度1.0km/h，试验对象在跑步机上做正常的步行运动。超声扫描仪的视频输出通过视频采集卡(大恒公司，CG400)来采集，视频帧率为25Hz。实验中持续采集30s的视频数据，对起始和结尾时间段的视频进行剪切，去除运动以及不稳定因素的影响，最终对每名试验对象提取180帧连续的超声图像序列。

3 结果

利用所提出的方法，对6名试验对象(共6×180=1 080帧)进行检测。以人工测量的结果作为参考标准，在对所有方法的结果进行比较和分析后，进一步证明本方法的正确性。所提方法在Matlab2013b版本下实现。

在本实验中，手动过程由两位临床专家完成。人工测量时，每一位临床专家在每帧图像的深浅层筋膜的感兴趣区域内，选择标记三段显著的肌束，取平均后得到最终的人工测量结果，如图3所示。

图3 手工测量示意图Fig.3 The diagram of the manual method

在步行运动中，肌肉收缩引起肌肉快速运动，并产生了较大的运动模糊。此时，人工和自动检测都很难准确地检测肌束方向，因此数据分析中将先排除因运动模糊而无法检测的超声图。跟踪的典型结果如图4所示，其中取两名临床专家的测量结果的平均值作为人工测量的最终结果，用黑色虚线表示，本方法用黑色实线表示。总体来看，本方法能捕捉正常人步行中肌束方向的变化趋势，与人工测量结果相差较小。进一步分析同一组人工检测与自动方法检测结果的相关性，Pearson相关系数r=0.93(见图5)，即本方法与人工测量结果有高度相关性。

图4 肌束方向跟踪结果Fig.4 Tracking result of the muscle fascicle orientation

图5 人工与自动检测结果的相关性分析(实线是线性趋势线)Fig.5 The correlation analysis of the manual and automated test result, the sold line suggests the linear trend line

表1 自动测量方法和人工测量方法比较Tab.1 The comparison between automated and manual measurement

表1总结了6名试验对象的肌束方向测量结果，包括人工与笔者提出的自动测量结果的相关系数与误差。从整体上看，相关系数的范围是0.88～0.94，在统计学意义上的平均相关系数为0.92±0.02，即人工与自动检测结果有很好的相关性,证明本方法的测量精度较高。此外，自动和人工检测方法的平均误差为0.30°±0.62°，可以推断该方法的准确性能够满足临床需求。可见，本方法在连续检测正常步行中肌肉肌束的方向具有很强的鲁棒性。

4 讨论

实验表明，所提出的自动测量方法与人工测量的平均误差达到0.30°±0.62°，Pearson相关系数为0.92±0.02，性能优于或者接近其他的一些自动测量方法。例如，Zhou等提出的迭代Hough变换算法，平均误差是0.18°±2.41°，Pearson相关系数为0.98[15]。Zhao等提出的局部Radon变换算法，误差约为1°[16]。而且，上述两个研究均是在静态或准静态的情况下测得的，深浅层筋膜区域被设为固定。本研究是针对正常步行运动进行分析，运动模型和噪声等问题更加复杂，自动化计算与跟踪超声图像中肌束方向的难度增大，因此本方法的准确性和鲁棒性更高。此外，本研究对提出方法的验证更充分。在实验设计中，搜集了6名试验对象，对每名试验对象处理至少180帧连续超声图像，总计处理帧数为1 080帧，充分验证了测量肌束方向新方法的稳定性和准确性。其他的研究一般只采用了1～2名试验对象, 图像处理帧数有限[15-16]。综上所述，本方法在自动化计算与跟踪肌束方向上有较高精度，且可实现正常步行运动中对感兴趣区域的自动化选取，实现了肌束方向的全自动测量。

所提方法的计算复杂度主要集中在局部Radon变换的方向估计。若采用M表示ROI尺寸的规模， 表示Radon方向的总数规模，则LRT方法的计算复杂度为O(θ·M2)，而应用中ROI的尺寸常小于图像整体尺寸，因此该方法的实际运算效率较高。下一步工作将采用执行效率更高的语言(如C语言)来实现本方法，有望实现肌束角度的实时在线测量。

在实验中，误差产生的原因有多个方面。首先，较静态或准静态下超声在肌肉运动的评估，在正常步行运动分析中，运动模型和噪声等问题更加复杂，难以准确、动态地计算与跟踪肌束方向。本研究引入了一阶运动模型，今后可能考虑在步行中采用更高阶的运动模型，进一步提高测量精度。其次，步行运动中肌肉收缩引起肌束快速运动，产生较大的运动模糊，对自动和人工检测方法的准确性检测造成一定的难度。在以后的研究中，将采用更高的视频采样帧率和更准确的读取图像方法，获取更快步行速度和更清晰的超声图像序列。最后，实验中如何在步行运动中固定超声探头，防止探头不规则滑动影响超声数据的采集，也是急需解决的问题。要从根本上解决该问题，需要开发出更加轻便灵活的超声探头。

5 结论

综上所述，本研究实现了对感兴趣区域的自动提取，基于LRT和卡尔曼滤波，自动化计算与跟踪超声图像序列中的肌束方向。实验证明，在正常步行运动中，本研究提出的方法对连续检测肌束方向有较高的鲁棒性。选用卡尔曼滤波器作为辅助跟踪的运动模型，预测存在观测噪声情况下的运动估计，然后再对估计的方向参数进行修正，并利用Gabor滤波自适应地增强图像肌束区域，提高了测量结果的准确性。本方法与临床人工检测结果具有高度的一致性，可以一定程度地替代人工测量，在处理大量超声数据时具有较大的优势。随着超声在肌肉检测上的应用，利用重复性好的自动方法定量分析肌肉肌束方向参数，在临床诊断和康复评估上具有巨大潜力。

[1] Disselhorst-Klug C, Schmitz-Rode T, Rau G, et al. Surface electromyography and muscle force: Limits in semg-force relationship and new approaches for applications [J]. Clin Biomech (Bristol, Avon), 2009, 24(3): 225-35.

[2] 李乔亮，易万贯，陈昕，等. 融合实时超声影像的多模态运动特性研究 [J]. 中国生物医学工程学报, 2012, 31(4): 518-525.

[3] Blazevich AJ, Gill ND, Zhou S, et al. Intra- and intermuscular variation in human quadriceps femoris architecture assessed in vivo [J]. J Anat, 2006, 209(3): 289-310.

[4] Chen Xin, Zheng Yongping, Guo JingYi, et al. Sonomyographic responses during voluntary isometric ramp contraction of the human rectus femoris muscle [J]. Eur J Appl Physiol, 2012, 112(7): 2603-2614.

[5] Fukunaga T, Kubo K, Kawakami Y, et al. In vivo behaviour of human muscle tendon during walking [J]. Proc Biol Sci, 2001, 268(1464): 229-233.

[6] Hodges PW, Pengel LH, Herbert RD, et al. Measurement of muscle contraction with ultrasound imaging [J]. Muscle Nerve, 2003, 27(6): 682-692.

[7] Kawakami Y, Ichinose Y, Fukunaga T, et al. Architectural and functional features of human triceps surae muscles during contraction [J]. J Appl Physiol (1985), 1998, 85(2): 398-404.

[8] Maganaris CN, Baltzopoulos V, Sargeant AJ, et al. In vivo measurements of the triceps surae complex architecture in man: Implications for muscle function [J]. J Physiol, 1998, 512 (Pt 2): 603-614.

[9] 李乔亮，任盼盼，张会生，等. 基于光流的超声图像肌肉厚度自动测量方法 [J]. 中国生物医学工程学报, 2013, 32(2): 149-153.

[10] Herbert RD, Gandevia SC. Changes in pennation with joint angle and muscle torque-in-vivo measurements in human brachialis muscle [J]. Journal of Physiology-London, 1995, 484(2): 523-532.

[11] Narici MV, Binzoni T, Hiltbrand E, et al. In vivo human gastrocnemius architecture with changing joint angle at rest and during graded isometric contraction [J]. Journal of Physiology-London, 1996, 496(1): 287-297.

[12] Gao F, Zhang LQ. Altered contractile properties of the gastrocnemius muscle poststroke [J]. Journal of Applied Physiology, 2008, 105(6): 1802-1808.

[13] Mohagheghi AA, Khan T, Meadows TH, et al. In vivo gastrocnemius muscle fascicle length in children with and without diplegic cerebral palsy [J]. Developmental Medicine and Child Neurology, 2008, 50(1): 44-50.

[14] 刘芳，蒲传强，时宵冰. 肌炎特异性自身抗体在多发性肌炎/皮肌炎及其他神经肌肉疾病的表达 [J]. 临床神经病学杂志, 2009 (3): 219-220.

[15] Zhou Yongjin, Zheng Yongping. Estimation of muscle fiber orientation in ultrasound images using revoting hough transform (RVHT) [J]. Ultrasound In Medicine And Biology, 2008, 34(9): 1474-1481.

[16] Zhao Heng, Zhang Liqun. Automatic tracking of muscle fascicles in ultrasound images using localized radon transform [J]. IEEE Transactions on Biomedical Engineering, 2011, 58(7): 2094-2101.

[17] Rana M, Hamarneh G, Wakeling JM. Automated tracking of muscle fascicle orientation in b-mode ultrasound images [J]. Journal of Biomechanics, 2009, 42(13): 2068-2073.

[18] Zhou Yongjin, Li Jizhou, Zhou Guangquan, et al. Dynamic measurement of pennation angle of gastrocnemius muscles during contractions based on ultrasound imaging [J]. Biomedical Engineering Online, 2012, 11(17):63.

[19] Zhou Yongjin, Zheng Yongpin. Enhancement of muscle fibers in ultrasound images using Gabor filters[C]// 2009 IEEE International Ultrasonics Symposium (IUS). Roma:IEEE, 2009: 2296-2299.

[20] 牟锴钰，韦明，郭建平. 基于小波和卡尔曼平滑的事件相关电位单次提取 [J]. 中国生物医学工程学报, 2012, 31(2): 167-174.

Automatic Tracking of Muscle Fascicle Orientation in Ultrasound Images: A Novel Approach Based on Local Radon Transform and Kalman Filter

Wen Huiying1Yang Xiaojuan1Guo Yanrong1Zhang Shuai2Zhang Zhiguo2Wang Tianfu1Chen Siping1Chen Xin1*

1(DepartmentofBiomedicalEngineering,SchoolofMedicine,ShenzhenUniversity,National-RegionalKeyTechnologyEngineeringLaboratoryforMedicalUltrasound,GuangdongKeyLaboratoryforBiomedicalMeasurementandUltrasoundImaging,Shenzhen518060,Guangdong,China)2(DepartmentofElectricalandElectronicEngineering,TheUniversityofHongKong,HongKong999077,China)

Ultrasound imaging is a noninvasive, real-time and convenient imaging method, which has been widely applied in analysis of human muscle contraction. So far manual operation is commonly used in clinical applications to extract morphological parameters from ultrasound image, however, may bring outcomes such as subjective, time-consuming and low repetitive. In this paper, an automatic method is proposed to calculate and track the fascicle orientation in sequences of ultrasound images. We constrained mutual information-based free-form deformation algorithm to automatically detect reading of interest. Then the detected results were projected based on local Radon transform. Then the Kalman filter image processing technique was applied to predict and track the next frame change of fascicle orientation. Six subjects participated in the walking experiment. We collected the dynamic ultrasound images from the medial gastrocnemius muscle during the process and calculated the muscle fascicle orientation using the proposed method. Experimental results showed that the proposed automatic detection method achieved better robustness. Compared with the manual measurement, the average correlation coefficient was 0.92±0.02, and the average error was 0.30°±0.62°. To some extent, the proposed method can replace the manual measurement, presenting a significant advantage especially when dealing with a large number of ultrasound data.

ultrasonography; fascicle orientation; radon transform; Kalman filter

10.3969/j.issn.0258-8021. 2016. 02.003

2015-12-26， 录用日期:2015-11-26

国家自然科学基金(81000637)

R318

A

0258-8021(2016) 02-0141-07

*通信作者(Corresponding author)， E-mail: chenxin@szu.edu.cn