动态骨骼肌超声图像处理研究进展

【作 者】汪灵梦 ，赵万明 ，邢运成 ，孙 通 ，陈 昕

1 深圳大学生物医学工程学院，深圳市，518060

2 医学超声关键技术国家地方联合工程实验室，深圳市，518060

3 广东省医学信息检测与超声成像重点实验室，深圳市，518060

0 引言

肌肉是构成人体的重要组织，与各项生命活动息息相关。人体总共有600多块肌肉，按结构和功能肌肉可以分为骨骼肌、心肌和平滑肌三种。肌肉通过收缩和舒张可以帮助我们完成各种复杂运动，例如在消化道运输食物，心脏搏动以及牵引骨骼运动。肌肉形态结构与功能密切关联，定量分析和评估肌肉的长度、厚度、羽状角、曲率以及横截面积等各项形态学参数是当前研究的热点，在临床医学、康复、运动学等领域得到了广泛应用[1-2]。

测量肌肉运动最主要方法是表面肌电信号(Surface Electromyography, SEMG)。肌电信号时频特性直接反映肌肉力的大小、肌肉的功能和状态、肌群的相互协作协调等特性。然而肌电信号容易受到各种潜在因素影响，如电极位置、肌肉类型、邻近肌肉干扰等，这些都制约了肌电信号在肌肉评估中的应用。

超声（Ultrasound, US）成像相对于X线成像、计算机断层成像（Computer Tomography, CT）以及磁共振成像（Magnetic Resonance Imaging, MRI）等其它成像方式，具有实时、快速、无辐射、价格低廉等优势。超声成像可获得肌肉收缩实时动态图像，被广泛应用于研究骨骼肌形态学参数，骨骼肌超声成像为临床观察和诊断肌肉组织特性提供了重要工具。

骨骼肌由肌腹和肌腱两部分构成，肌腹由肌纤维构成，肌腱由肌束末端致密结缔组织汇聚形成，肌腹通过肌腱牵引骨骼从而产生运动。肌束附着在腱膜一侧，解剖学上把靠近皮肤的腱膜叫浅层腱膜，靠近骨骼的腱膜叫深层腱膜。如图1所示，肌肉厚度是肌肉深层腱膜和浅层腱膜间的距离；羽状角是深层腱膜和肌束方向形成的锐角；肌肉长度是肌束与深层腱膜和浅层腱膜交点间的距离。

图1 骨骼肌形态结构示意图Fig.1 Simplified structure of musculoskeletal

1 动态骨骼肌超声图像处理常用算法

早期骨骼肌超声都是获取静态超声图像，处理方法通常都是手动处理。随着技术的发展，更多研究都是采集肌肉运动过程中的动态超声图像，用手动方法处理不但耗时，而且结果不稳定。因此各种图像处理方法被提出。本节对常用动态骨骼肌超声图像处理算法进行简要概述，各种算法应用将在第2节具体介绍。

1.1 直线检测算法

在超声图像中肌纤维呈现近似直线结构，因此各种直线检测算法被用于肌纤维检测。Hough变换和Radon变换是图像处理中常用的直线检测算法。

1.1.1 Hough变换

霍夫变换(Hough Transform)是图像处理中从图像中识别几何形状的基本方法之一。由Hough P于1962年提出，最初只用于二值图像直线检测，后来扩展到多种形状检测。Hough变换基本原理在于通过改变曲线表达形式，将原始图像空间给定曲线变为参数空间的一个点。这样就把原图像中给定曲线检测问题转化为寻找参数空间中峰值的问题。Hough变换检测肌纤维结构一般分为四步：①对图像做二值化和边缘检测，作为Hough变换预处理步骤；②对边缘检测后的图像做Hough变换，返回Hough变换矩阵、角度和半径用于后续直线提取；③设定好需要检测Hough变换矩阵峰值数（即直线条数），用Hough变换峰值检测函数找出符合条件的峰值；④在找出的峰值中按预先设定的直线检测条件，如最小直线间距、直线长度，用Hough变换直线检测函数找出符合条件的直线，即为检测的肌纤维位置。

1.1.2 Radon变换

Radon变换（Radon Transform）常用于二维或三维数据重建，相较于Hough变换Radon变换可以直接处理灰度图像。Radon变换是对一个平面f (x, y)，在极坐标系中沿着不同的方向投影，对每一条投影线计算平面f (x, y)线积分，得到的积分结果就是f (x, y)的Radon变换。Radon变换基本原理在于计算图像在各个方向的积分，即亮度累加，这样就把直线检测转化为在变换区域对亮点、暗点的检测。Radon变换检测肌纤维结构一般分为四步：①把肌肉图像转化为灰度图像、设置Radon变换角度范围，作为Radon变换预处理步骤；②对灰度图像进行Radon变换，得到Radon变换矩阵和半径；③设定好需要检测直线条数，用峰值检测函数找出符合条件的峰值；④将峰值点坐标由极坐标系转化到直角坐标系，从而得到图像中肌束位置。

1.2 特征点/特征区域追踪算法

超声图像中有的肌纤维有一定弧度呈现出非直线形状，这时需要手动在第一帧超声图像上选取特征点/特征区域，然后运用相关算法辨别它在下一帧的位置，从而计算出每帧超声图像中肌肉参数。区域匹配算法和光流法是常用于特征点/特征区域追踪的算法。

1.2.1 区域匹配算法

区域匹配算法通过跟踪特定区域位移变化提取出组织运动信息，用于追踪目标组织中平移或接近于平移运动，如肌肉腱膜移动。区域匹配算法追踪肌肉运动一般分为三步：①在第一帧超声图像中选择目标追踪区域；②根据某种匹配准则，在下一帧图像中搜索与目标区域最相似的区域；③根据追踪结果计算目标区域移动，依此类推实现目标区域动态追踪。常用的匹配准则包括：互相关、快速互相关、平方差、绝对差分求和等，选用不同准则得到的追踪结果也不相同。互相关算法需要在全局范围内计算两帧图像的互相关函数，运算量较大；快速互相关算法通过灵活选择搜索方法可以降低计算量；绝对差分求和算法直接计算两个检测区域的绝对差分和，可以提高目标区域检测速度。

1.2.2 光流法

光流法是利用图像序列中像素在时间域上的变化、相邻帧之间物体亮度相关性来找到上一帧跟当前帧之间存在的对应关系，计算相邻帧之间物体运动信息的一种方法。光流法大多是建立在Horn-Schunck算法和Lucas-Kanade算法基础之上。肌肉收缩具有时间连贯性，采用光流法可以有效捕捉肌束连续运动[3]。光流法追踪肌肉运动方法和区域匹配法类似也分三步：①假设目标区域运动时保持亮度不变，这是光流法最基本的约束条件；②根据亮度约束条件和全局约束条件、局部约束条件等附加约束条件求解出目标区域在后一帧的位置信息；③用第二步的方法追踪当前帧和后一帧的位置关系，实现目标区域动态追踪。相对于区域匹配算法，光流法计算量相对较大，但是鲁棒性强、可用于追踪平移之外的其它肌肉形变。

2 骨骼肌超声算法研究

现有对肌肉参数的研究主要集中在肌肉羽状角、长度和厚度动态测量，本节分别介绍这三个参数不同提取方法以及对这些方法进行简要比较。

2.1 羽状角

测量肌肉羽状角主要在于确定深层腱膜位置和肌束方向。ZHOU等[4]提出迭代Hough变换（Revoting Hough Transform, RVHT）检测肌束方向，这种方法通过将超声边缘图像映射到Hough空间，找出图像中最明显的直线，然后移除这条直线设定线宽范围内所有特征点，再用移除这些特征点后的图像进行下一次Hough变换检测肌束方向，重复上述过程直至低于自定义阈值。这种方法可以检测出图像中每一个可能肌束方向，但检测误差会因迭代过程累加。深浅腱膜在肌肉超声图像中是最明显也是最容易检测出来的，ZHOU等[5]利用肌肉的这一结构特性改进了迭代Hough变换。这种方法核心思想是对高斯滤波和二值化后的图像用迭代Hough变换检测出深浅腱膜方向，再在检测出的深浅腱膜间选择一个合适的感兴趣区域(Region of Interest, ROI)用于肌束方向检测。这种方法缩小了Hough变换区域，减小了计算量。Hough变换检测直线特性要经过判断检测峰值是否符合要求的投票步骤，增加了计算复杂程度，为了尽可能避免这一过程LING等[6]提出了一个肌束方向检测新框架，这种方法只对有较小间隔肌纤维部分进行迭代Hough变换，达到了降低计算量的目的。

Hough变换对斑点噪声敏感，极大依赖于边缘检测算子的性能，只能对二值图像处理，Radon变换受斑点噪声影响小，不必进行边缘检测还可以直接对灰度图像进行处理。RANA等[7]提出用Radon变换检测肌束方向。由于Radon变换是沿不同方向和角度对整个图像区域投影，那些长度较短的肌纤维投影值较小，不容易被检测出来，ZHAO等[8]利用肌肉结构中肌纤维方向角度差异小这一特性提出了局部Radon变换，其主要方法是对肌肉图像在接近肌纤维方向进行投影。这样不但可以检测出长度较短的肌束，而且极大缩小了Radon变换角度范围减少了计算量。温慧莹等[9]提出了基于局部Radon变换和卡尔曼滤波的超声图像肌束方向跟踪方法，用基于约束互信息的自由变换算法和局部Radon变换实现了ROI区域自动选取和肌束方向检测。

以上这些方法是利用超声图像直线特征来提取肌束方向，存在几个缺点：①肌束通常是由多段组成，各段尺寸明显小于图像总尺寸，长度较小的肌纤维不容易被检测出来；②图像上肌束像素灰度在不同位置分布不均匀，肌束直线所对应Radon空间极值点强度被削弱；③标准Radon变换是整个图像平面内像素点灰度值的积分，直线方向上像素点灰度值能够产生极值，同时投影方向积分路径较长也能产生极值。虽然局部Radon变换通过限制肌束位置和方向来增强准确性，但这些方法都是检测肌束直线特征，并未从根本上克服以上问题，检测性能存在局限。

LI等[10]利用了肌束并行排列的纹理特征提出了归一化Radon变换方法，该方法原理是利用归一化Radon变换对超声图像进行投影，当投影方向与肌束方向重合时，归一化Radon变换结果在不同截距上变化达到最大，再在所有Radon方向利用分段-方差(Segmented-Variance)法求取不同截距上统计变化量，变化最大的方向即为肌束方向，这种方法很好避免了由于积分路径不同对测量结果的影响。CHEN等[11]创新性提出了频域Radon变换，先对超声图像进行傅里叶变换，然后再对频域进行Radon变换，最终找到肌束方向。与检测单一肌束方向相比，这两种方法可以检测出某一区域所有肌束的平均方向，提高了算法鲁棒性与准确性。

2.2 长度

肌肉长度是肌束与深层腱膜和浅层腱膜交点间的距离，因此很多长度测量方法都利用肌肉结构这一特点，直接计算两点间的距离。这些方法通常先选取一帧肌肉超声图像为基准帧，然后在肌束与腱膜交点处选取两个特征点/特征区域，计算出两个特征点/特征区域的距离即为肌束长度，最后用特征点/特征区域追踪算法追踪相关特征点/特征区域，实现肌束长度动态追踪。

2005年LORAM等[12]提出块匹配算法追踪腓肠肌和比目鱼肌运动，它的核心步骤是用互相关算法寻找相关性最大区域。这种方法对于肌肉运动范围较小的超声图像序列有明显检测效果，DARBY等[13]提出了Kanade-Lucas-Tomasi特征追踪算法，这种方法和上一种方法相比用平方差算法代替互相关算法寻找相关性最大区域，当肌肉长度变化较小时两种方法效果相似，但当肌肉长度变化较大时这种方法比互相关算法表现更好。CHUANG等[14]提出基于多核块匹配的光流算法追踪腱膜运动，其主要方法是将ROI区域分成很多小块，对每小块用块匹配算法找近似匹配块，然后组合所有小块匹配结果即为整个区域的匹配结果。这种方法相对于前两种方法在寻找ROI区域方面准确性有了提高。

基于块匹配的斑点追踪算法只对肌肉平行移动有较好追踪效果，而肌肉运动往往是不规则的有时会发生形变。光流法可以检测物体不规则运动，CRONIN等[15]提出了基于Lucas-Kanade光流算法的肌束长度自动追踪方法，这种方法主要是用光流算法追踪ROI区域从而确定肌束长度。FARRIS等[16]利用光流算法开发的UltraTrack软件可以用来测量肌束长度动态变化，对肌肉超声数据进行离线分析。

以上这些方法是通过追踪端点/ROI的运动来评估肌肉长度变化，但肌肉是由很多有序排列的肌纤维构成的，这些方法只能计算出很小区域内肌肉参数。ZHOU等[17]提出了一种更宏观的肌肉参数计算方法，从整幅超声图像入手，先对超声图像用Gabor变换，增强超声图像中肌肉线性结构，然后用Hough变换检测腱膜方向，标准Radon变换检测肌束方向，然后用定义求肌肉长度。这种方法是对整帧超声图像做变换，扩大了肌肉参数检测有效面积，可以提高检测准确性，但是超声设备所成肌肉图像中还会夹杂着与肌束方向不同的血管， ZHOU等[18]进一步提出了在肌束方向敏感段用光流算法评估肌束长度的方法。利用肌束方向粘结性分割出超声图像中有相同肌纤维方向的区域即方向敏感段，再从方向敏感段提取腱膜和肌束信息，从而获得肌肉参数。

2.3 厚度

一般测量肌肉厚度有两种方法，一种是基于点的追踪方法，这种方法和肌束长度追踪方法类似，在图像基准帧中选取两个关键点，然后在连续帧用特征区域/特征点追踪方法定位关键点，最后计算出两点间的距离即为肌束厚度。另一种是基于肌肉腱膜轮廓的追踪方法，这种方法是先检测出深浅腱膜然后测量两腱膜间距离。

CARESIO等[19]提出了肌束厚度自动提取方法，它对高斯滤波和导数增强肌肉线性结构后的图像用标准Hough变换检测深浅腱膜方向，然后计算腱膜上特征点间距离即为肌肉厚度。ZHENG[20-21]等提出了用互相关算法连续测量肌束厚度，这种方法是追踪手动选取在深浅腱膜上的矩形区域，然后测量矩形区域中水平中心线间距离即为肌肉厚度。但是这种方法需要手动选择ROI，重复性和可靠性不能得到很好保证。ZHENG等[22]进一步提出了一种高效基于压缩追踪算法的Coarse-to-Fine方法，这种方法涉及多次滤波、采样和判断，检测效果有一定提升，但计算量大并且由于浅层和深层腱膜运动不一致，两个ROI水平位置在运动中不能保证对齐造成测量结果不准。

为克服上述ZHENG等方法的问题LI等[23-24]创新性提出了一种基于光流的超声图像肌肉厚度测量方法。这种方法，根据肌肉腱膜与测量位置之间的几何关系计算肌肉厚度，分别在第一帧图像深浅腱膜上人工标记一个ROI，为了保证ROI对齐，在ROI肌束上再各标记两个点，对每帧超声图像均采用光流法跟踪ROI，利用两点间短线动态跟踪肌束。最后在超声图像一定像素位置引入一条垂线，垂线与两短横线间交点间的距离即为肌肉厚度。

上述几种方法是基于点或者较小范围追踪腱膜的方法，只能获得腱膜上几个点或者两段直线间的距离，这些方法依赖于操作者经验，获得的结果具有一定偶然性，不能完全反映肌肉厚度信息，因此有人提出了基于轮廓的追踪方法。HAN等[25]假定超声图像中深浅腱膜是主要肌束特征，用迭代Hough变换定位腱膜，深浅腱膜上所有点间距离的平均值即为肌肉厚度。LING等[26]提出了一种基于局部和全局强度拟合模型（Local and Global Intensity Fitting, LGIF）的肌束厚度追踪算法，这种方法首先人工在超声图像的深浅腱膜各选一点作为特征点，然后用区域生长方法，在特征点临近像素范围内寻找与特征点性质相似的点，并与之合并形成新的生长区域，这样不断向外扩张直至得到肌束深浅腱膜轮廓，然后计算腱膜间厚度平均值。与基于点的追踪方法相比基于轮廓的追踪方法准确性和鲁棒性更高，得到的肌肉厚度更具代表性，能更好反映肌肉厚度变化。

3 总结与展望

本文简述了骨骼肌超声图像处理算法方面的研究进展。Hough变换处理二值图像只需要处理前景或者背景像素，计算量小、速度较快，但由于二值图像不连续性，难以对一幅图像进行恰当二值分割，所以在一般情况下Radon变换要比Hough变换更精准。区域匹配方法用于追踪特征点/特征区域原理简单，但只对肌肉平移变化敏感，光流算法追踪超声图像中特征点/特征区域有很好的效果，结果可靠，但其计算量大，对计算设备配置要求较高。这些肌肉参数提取和追踪算法选择在实际应用中要结合超声图像自身特点、设备性能以及需要得到的参数等多方面综合考虑。肌肉参数计算结果是否准确，不仅和超声图像处理算法有关，还与超声探头、超声射频信号重建质量有很大关系。不同探头放置方法获得的超声图像不同，怎样放置超声探头才能获得最利于图像处理的超声图像，到目前为止还没有标准可供参考。骨骼肌长度偏长，正常成年人股肌长度为8～12 cm，而一般线阵超声探头视野只有4～6 cm，不能完全获取肌束长度图像，宽景成像技术通过移动探头可以行成更大视野超声图像。骨骼肌形态参数获取也不局限于骨骼肌超声图像，在超声射频信号中提取肌肉参数是一种新思路[27]。近几年深度学习发展迅速，用深度学习的方法可以更准确地分割出肌束方向[28]，深度学习为动态骨骼肌超声图像处理提供了一种新的研究方法。在图像处理领域现有肌肉图像离线处理算法已经比较成熟，也取得了丰硕研究成果。但到目前为止这种技术还没有应用到临床上，将骨骼肌图像处理算法实现到超声设备上，研发肌肉形态参数检测专用超声设备有广阔应用前景[29]。