基于N线模型的超声探头自动标定方法

李文骏 丁 辉 王广志

(清华大学医学院生物医学工程系，北京 100084)

引言

超声成像技术具有实时、廉价、无电离辐射等优点，在临床诊断和引导下的介入手术等方面得到了广泛应用[1-2]。在超声引导下的介入治疗、超声融合成像和徒手三维超声成像中[1-3]，为了跟踪探头的空间位置姿态，常常在超声探头上固定光学或电磁位置传感器，通过一系列的空间坐标变换，将超声图像和术前影像、器械位姿等统一到同一个坐标系中。笔者在本文中所讨论的超声探头空间标定，是指确定其中位置传感器坐标系到超声图像坐标系之间的变换关系的过程，这对保证超声引导介入、融合成像或三维重建的系统精度有重要的意义，标定过程中的坐标变换关系如图1所示。通过模型中的位置确定N形标记线，在超声成像中切割出对应的标志点，以及在模型与超声探头上定位传感器的相对位置，可以形成闭环的空间变换关系，其中Ts←i就是所需要标定的变换关系。

图1 标定过程中的坐标变换关系Fig.1 Coordinate systems and transforms in calibration

国内外的研究者已经提出了多种超声探头空间标定方法[4]，其中基于N线模型的方法是应用最为广泛的方法之一[5-6]。这种模型在水箱中设置多组由绷直的细线构成的N形目标，每个N形目标经超声成像形成3个亮斑。拾取图像上各组亮斑的坐标并确定其对应的N线编号后，就可以利用模型的结构参数解得最终的标定变换。

目前，拾取亮斑的坐标、判别哪些亮斑属于同一组N线、确定N线组的编号等工作，是靠人工完成的，这种人工识别和点选亮斑标志点的方法存在一些问题。首先，识别和确定N线编号要求操作者了解N线模型的设计信息，对医护人员或普通工人的实际操作有较大难度;而由于模型中N形目标的数量较多，即使是有经验的实验人员也难免出现错误判断N线编号的情况。其次，人工点选超声图像上的亮点坐标带有主观性，会引入随机误差。再次，为获得较好的标定精度，往往需要获得多帧图像上的几十甚至几百个N形目标点，人工操作耗时较长，这也更增加了发生编号错误的概率。

为解决上述问题，研究者提出了若干基于N线模型的自动或半自动超声标定方法。Lindseth等[7]和宋等[8]人工在亮点附近指定初步位置，在其邻域内找到亮度脊，并利用几何约束进行修正。这种方法减小了选点时的主观误差，但并未解决标志点自动编号的问题。Hsu等在模具表面加装特制的薄膜，通过检测在超声图像上形成的直线，实现了自动标定，但增加了模具设计和制作的复杂度[9]。Chen等则把N线组的数目减少到2组，降低了N线识别的难度，但需要拍摄更多帧的图像[10]。

针对以上问题，笔者在不改变N线模型设计和实验流程的条件下，提出一种自动化的超声探头标定处理方法，提高了标定的效率和精度。

1 方法

1.1 标定原理

1.1.1 N线模型目标点重建

用于超声探头标定的N线模型由多层尼龙线安装在水槽中组成。这些尼龙靶线构成一系列被编号的N线组，每个N线组包括两条平行线和一条斜线。N线组的成像过程可以抽象为超声图像平面与靶线相交。如图2(a)所示，3条靶线与超声成像平面相交形成3个交点，即图中的E、G、F。其中，点F在模型坐标系下的三维坐标可以利用超声图像和模型结构信息求得，这个过程被称为N形目标点的“重建”。从超声图像中拾取标志点E、F、G的图像坐标，计算线段EF与EG长度的比值，有

由于三角形△BEF和△CGF相似，得

对已知编号的N线组，式(2)中B点和C点的坐标可由N线模型的设计信息确定，由此即可重建出F点在模型坐标系下的三维坐标。

图2 N形目标点重建。(a)坐标重建原理;(b)实际成像光斑Fig.2 Target point reconstruction of N-fiducials.(a)Schematic;(b)Sample image

在实际图像中，由于超声成像平面有一定的厚度，且超声探头横向聚焦并非绝对理想，因此靶线所成的标志点呈现亮斑的形式，如图2(b)所示。采用人工拾取标志点坐标，一方面难以确定亮斑中心的精确位置，另一方面易受背景噪声等因素影响而出现错误。在超声成像视场中N线组很多的情况下，人工所记录的编号也可能出现错误，这些都可以通过自动化手段加以改进。

1.1.2 空间变换关系

标定实验涉及图1所示的几个坐标系，包括世界坐标系{w}、模型坐标系{m}、传感器坐标系{s}和图像坐标系{i}。在标定中，各坐标系之间均考虑为三维刚性变换，包含3个旋转参数和3个位移参数，共6个自由度。利用齐次坐标，刚性变换可以用矩阵乘法表示。用Τa←b表示从坐标系{b}到坐标系{a}的变换，用xa代表点在{a}系下的坐标，则有

超声探头标定的最终目标是求解图像坐标系到传感器坐标系的变换，即Ts←i。在实际标定过程中，由三维定位系统的输出数据可得Ts←w。在每组实验开始时，先用三维定位探针点选出模型外壁和上缘的定位孔，并将小孔在世界坐标系中的位置和在模型坐标系中的位置配准，可以求解出变换Τw←m。对每个完整可见的N线组，用本文1.1.1节的目标点重建算法，求得F点在模型坐标系下的坐标xm，并由式(4)变换到传感器系，即xs。另外，由预先标定好的超声像素尺寸，可以直接从超声图像中求得标志点在图像坐标系下的坐标xi。由3个以上的(xs-xi)坐标对，可求解标定变换Ts←i。在最小二乘意义下，即为

至此，标定问题就抽象为一个已知配对关系的点对刚性配准问题，可利用Horn等提出的解析算法[1-2]来求解。

1.2 自动标定算法

从图像中，自动提取亮斑中心准确的坐标，并确定其对应的N线组的编号，是实现基于N线模型的超声探头自动标定方法的关键。自动提取亮斑的亮度中心，能消除主观点选坐标带来的随机误差，提高标志点坐标的精确度。另外，进一步通过自动算法确定各标志点对应的靶线编号，能有效防止实验者在确定N线编号时偶然出错的情况。

然而，实现N线组的自动编号需要克服一些困难。为了减小在超声成像平面法线方向上的标定误差，应使探头在不同角度下获得的图像均包含数量较多的N线目标，并且其分布应尽量遍布整个超声图像范围，因此在模型中通常要设置较多的N线组[9]。所设计的模型包含67条靶线(30个N线组，部分靶线由两个相邻的N线组共用)，每帧图像中有30～50个亮斑标志点，对应8～16个N线组。也就是说，每帧图像都只能观察到模型中的一部分靶线，而这部分靶线的编号和位置都不能事先知道。此外，图像上出现的亮斑粘连、缺失，以及由于超声在容器壁多次反射等原因造成的伪影，也是干扰N线标志点识别的因素。针对上述N线模型及其超声图像的特点，设计了一套自动标定算法，其流程如图3所示。

图3 自动标定算法流程Fig.3 Flow chart of the automatic calibration method

1)预处理。预处理步骤的目的是初步提取图像中的标志点亮斑。首先，在本实验的条件下，以96作为灰度阈值，对图像进行二值化。由于超声成像平面的厚度效应和超声横向聚焦不理想等原因，标志点亮斑常沿着超声束的切向扩散，呈横向扁平状;而容器壁的反射像(如图4(a)右下部的高亮线段)、超声波多次反射相干形成的条带状伪影(见图4(a)下部)和一般的图像噪声，则没有上述横向扩散的特征。因此，遍历二值化后图像的各个连通域，将竖直方向上延展长度超过30像素的连通域删除。用调整后的二值图作掩膜，提取出可能的标志点区域图像。预处理后的图像如图4(b)所示。

2)计算标志点预测位置。这一步是利用已知信息来预测靶线标志点在图像中的位置，而这些预测值将作为后面搜索步骤的初始位置。该过程用到了N线模型的设计形状尺寸信息，由此得到的标志点预测位置与模型靶线是一一对应的，且编号已知，这为后续步骤中自动确定N线组编号奠定了基础。

图4 自动标定算法各步骤处理结果。(a)原始图像;(b)预处理结果;(c)标志点预测位置;(d)邻域搜索获得的亮斑中心位置;(e)含编号的N线组识别结果Fig.4 Results of the automatic calibration algorithm.(a)Original image;(b)Results of preprocessing;(c)Predicted positions of N-fiducials;(d)Results of neighborhood search that represent refined positions of N-fiducials;(e)Results of N-fiducial identification

处理第一帧图像和后续其他图像的方法有所不同。采集第一帧图像时，让探头对准N线模型的中间位置，并保持超声成像平面与水模壁基本平行。由此所得的图像标志点的分布样式能被预先确定，从而能给出标志点的预测位置。处理第二帧及之后的各帧图像时，之前已经由前面的图像解得标定变换Ts←i的一个初始解;利用这个初始解，可以由三维定位传感器的读数来预测超声成像平面的位置和姿态，求出这个虚拟截面与模型靶线的交点，即为标志点的预测位置。为了数学上的方便，将模型靶线两端的控制点坐标P1和P2转换到图像坐标系下，有:

图5 线段被x-O-y平面所截Fig.5 Line segment intersecting the x-O-y plane

在式(6)中，Τs←w由三维定位系统的实时读数直接得到，Τw←m则在每组实验开始时配准得到。如图5所示，在图像坐标系下，超声成像平面就是x-O-y平面。设靶线端点坐标分别为P1(x1，y1，z1)和P2(x2，y2，z2)，则模型靶线的直线方程可写为

即为该靶线与超声成像平面交点的坐标。注意，交点必须在两端点之间，其限制条件为

计算每条靶线与虚拟超声平面的交点，并排除不满足限制条件的情况，就得到各标志点坐标的预测值，如图4(c)所示。

3)邻域搜索和N线编组。在各预测位置的一定大小的邻域内搜索亮度超过阈值的点。如果不存在，说明该点没有被成像，则该组N线无法用于标定，自动加以排除;如果存在，则找到邻域内面积最大的亮斑，并求出其亮度中心，认为是该标志点的坐标，如图4(d)所示。

通过以上处理，所搜得的标志点对应的靶线编号就是已知的，能方便地确定哪些标志点属于同一个N线组。为了避免搜索到图像中的噪声或伪影，或标错靶线编号，处理中进一步剔除N线组三点共线性太差(连线夹角大于3°)以及其中某两点重合或接近重合(距离小于10像素)的情况。对符合要求的N线标志点，标记其N线编号，并记录下对应的k值供重建时使用，如图4(e)所示。

4)重建和解算。至此，获得了一帧图像中所有能够找到的N线组编号及对应的k值。参照文中“标定原理”部分的方法，可以得到N线目标点在传感器坐标系和图像坐标系下的坐标，即增加了式(5)中新的(xs-xi)坐标对。利用Horn算法[11-12]，可以更新最终的标定变换Τs←i。

一帧处理完毕后，判断是否已处理完全部的标定图像。若否，读取下一帧图像，再次更新标定变换，直到所有图像处理完毕并输出结果。

1.3 实验验证

为了验证上述自动标定算法的性能，使用自行设计的N线模型进行了标定实验验证。模型用有机玻璃制成，外尺寸为28 cm×11 cm×24 cm;模型内部分布了7层深度为4～16 cm的靶线，由直径0.14 mm的尼龙线穿成并绷紧，每层有4～5个N线组;在模型外壁和上缘分布有多个定位孔，定位孔及穿线孔的加工精度均达到0.1 mm;在模型内壁和底面放置吸声材料，以减小反射噪声。

笔者采用的超声图像采集系统包括:超声成像仪(DC-6，迈瑞医疗)，配临床常用的3.5 MHz腹部探头;交流脉冲电磁三维定位系统(Aurora，Northern Digital Inc.);计算机(xw8400，HP)。在计算机端，用视频采集卡捕捉超声图像，并进行数字化。与图像同步的超声探头三维定位数据，通过串口发送到计算机。实验中，超声成像的深度为17.5 cm。标定时的传声介质为室温下去除气泡的纯净水。由于标定的结果将应用于人体组织的成像和测量，因此使用超声仪自带的液性声速补偿功能，可减小扫描线方向距离测量的误差。

标定实验共进行6组，每组包含7帧图像。每组的第1帧图像均在探头正直地对准水槽中央位置时采集，并使水模第4层N线组位于图像的中央位置，其余6帧图像对应的探头位置和角度不固定。

为了通过实验考察标志点自动识别的成功率，验证自动识别的可靠性，笔者从精密度(precision)和准确度(accuracy)两方面来验证新算法的精确度。精密度是指算法在多次实验中所得结果的一致程度，即重复性;准确度则是指测量结果与真值之间的误差大小。

采用一个直观评价精密度的方法，就是考察标定变换矩阵各元素在多次实验中的标准差[7]。但是，由于各研究中实际的标定变换互不相同，这种方法并不适用于标定方法的横向比较。以往广泛采用的评价方法是将超声图像中的某个固定点(如图像中心和4个角点)变换到传感器坐标系下，统计在多次实验所得的标定变换下点云的散布情况。记某点在超声图像坐标系下的坐标为xi，用各次实验所得的变换将其变换到传感器坐标系下的坐标均值为，则标定精密度可定义为

常用目标配准误差(target registration error，TRE)作为标定算法准确度的评价指标[5，7，9]。每组实验分别抽出一帧图像，将该帧图像内的N形目标点作为测试用的注册目标。利用其余各帧图像(训练集)得到标定变换矩阵，计算测试帧中各目标点由所得的变换转换到位置传感器系下的坐标，与由模型结构信息重建所得坐标间的欧氏距离的方均根，即

将自动方法与人工方法进行对比，给出了自动标定算法在识别正确率、标定精度和运行速度方面的结果。

3 结果

3.1 标志点自动识别情况

由于实际标定所采集的超声图像存在噪声、模糊和对比度不均匀等图像退化，因此人工点选和自动搜索都有未被提取的目标。对实验用的6组共42帧图像，通过自动识别或人工点选的方法，总共确定了461个N形目标，其中被自动算法成功识别的有414个，占总数的89.8%。N形目标识别结果的分类统计如表1所示。对所有自动识别得到的N形目标点进行了检查，未发现误识别(即把噪声误认为靶线标志点)或误编号的情况。自动算法未成功识别的目标点占总数的10.2%，数量可以接受，一般是由于标志点附近有较强的反射回声伪影引起的。比较自动识别和人工点选的目标点坐标，两者平均距离为1.25像素。以上结果说明，标志点自动识别是可靠的。

表1 自动和人工选点结果分类统计Tab.1 Classification of automatically and manually identified N-fiducials

3.2 精确度验证结果

如前所述，精确度包含精密度(precision)和准确度(accuracy)两个方面。精密度是指算法在多次实验中所得结果的一致程度，即重复性;准确度则是指测量结果与真值之间的误差大小。

根据式(11)，将采集的6组实验(每组包含7帧图像)由自动标定方法所得的精密度结果与传统方法(人工点选N形标志点，经坐标重建后利用解析算法，求解标定变换[5])的进行比较，如表 2所示。

表2 标定精密度对比Tab.2 Comparison of calibration precision

另外，式(12)所得出的准确度结果(见表3)表明，新算法与传统算法相比，目标配准误差显著减小，标定准确度有所提高。

表3 目标注册误差比较Tab.3 Comparison of target registration error(TRE)

随着参与计算的图像数量的增加，标定精确度逐渐得到改善，并最终收敛到一个较小的值上。以第4组实验为例，给出目标注册误差随帧数的关系曲线，如图6所示。采用新算法得到的TRE随帧数增加呈指数曲线下降，且收敛速度比传统算法更快。从第6帧到第7帧，TRE下降的比例为1.5%，表明已经基本收敛。

3.3 运行速度

本研究的自动标定算法全部在Matlab R2011a平台上实现。在搭载Intel®CoreTMi5-2450M CPU@2.50 GHz、4 GB内存的计算机上运行，处理每帧图像所需的时间平均为0.067 s，最大为0.13 s。相比人工逐个点选标志点，大大提高了效率。

4 讨论

图6 目标注册误差随帧数的关系曲线Fig.6 Target registration error(TRE)in terms of number of frame

笔者查询了近年来国内外利用N线模型进行超声探头标定且被标探头参数与本方法相近的文献，本方法所得标定精确度与相关文献比较的情况如表4所示。可见，本方法获得的精密度和准确度与国内文献比较均有较为明显的提高;在设备相对落后的条件下，超过或接近了国外文献报道的水平。

表4 本方法标定精确度与近年国内外文献的比较Tab.4 Comparison of calibration precision and accuracy to the recent literature

笔者进一步分析了自动识别标志点提高标定精确度的原因。对由自动识别或人工点选得到的所有标志点进行分类统计，如表1所示。其中，被自动方法选出但在人工点选时被略去的标志点约占13.0%，主要分布在超声图像的上部和边缘位置。这部分标志点一般是亮斑亮度偏暗、面积偏小，人工点选时难以辨识，因此只好被放弃。首先，自动方法利用了这些标志点，使N形目标点更均匀地覆盖到超声图像的全部范围，这是提高最终标定精确度的主要原因之一。其次，自动方法通过计算标志点亮斑的中心位置，消除了人工点选标志点过程中引入的随机位置误差。对自动和人工方法均选出的标志点(占总数的76.8%)，由两种方法所确定的坐标的平均距离为1.25像素，这反映了自动方法在消除随机误差上的贡献。再次，自动方法所得标志点的共面性比人工点选更优。标志点是模型靶线被超声成像平面截得的交点，理想情况下每帧图像所得的标志点应当是共面的。但由于误差的存在，实际重建的标志点并不完全共面;研究表明，共面性的丢失是影响最终标定精度的重要因素[13]。以第4组实验为例，比较了在各帧图像中由自动和人工方法分别得到的标志点坐标与其拟合平面间的平均距离和最大距离，如表5所示。可见，自动方法所得标志点的共面性有显著改善，从而提高了最终的标定精度。

表5 自动和人工方法所得标志点坐标与其拟合平面间的距离比较Tab.5 Statistics of distance from reconstructed fiducials to the fitting plane

笔者通过本研究提出了自动标定方法，也提高了超声探头标定的效率。采用传统的人工点选方法，实验者点选处理一帧图像的时间一般需要几分钟。而如本文3.2节所述，自动标定算法每处理一帧的时间平均只需0.67 s，因此大大节约了人工点选标志点所需的时间。

5 结论

超声探头标定是确定超声图像与三维定位系统传感器之间的空间变换关系的过程，是超声引导介入手术、超声融合成像、徒手三维超声等应用中的基础技术环节。面向目前广泛使用的基于N线模型的标定方法，本研究对传统的人工点选标志点的方法进行了改进，提出了一种新的自动化方法。通过实验验证，证实新方法能提高超声探头标定的精度和效率。

[1]Xu S，Krucker J，Turkbey B，et al.Real-time MRI-TRUS fusion for guidance of targeted prostate biopsies[J].Computer Aided Surgery，2008，13(5):255-264.

[2]Solberg O，Lango T，Tangen G，et al.Navigated ultrasound in laparoscopic surgery[J].Minimally Invasive Therapy＆ Aided Technologies，2009，18(1):36-53.

[3]Prager R，Rohling R，Gee A，et al.Rapid calibration for 3-D freehand ultrasound systems[J].Ultrasound in Medicme ＆Biology，1998，24(6):855-869.

[4]Mercier L，Lango T，Lindseth F，et al.A review of calibration techniques for freehand 3-D ultrasound systems[J].Ultrasound in Medicine＆ Biology，2005，31(2):143-165.

[5]Pagoulatos N，Haynor D，Kim Y.A fast calibration method for 3-D tracking of ultrasound images using a spatial localizer[J].Ultrasound in Medicine＆ Biology，2001，27(9):1219-1229.

[6]罗杨宇，徐静，鲁通，等.基于磁定位器的手动三维超声图像标定[J].中国生物医学工程学报，2008，27(2):250-254.

[7]Lindseth F，Tangen G，Lango T，et al.Probe calibration for freehand 3-D ultrasound[J].Ultrasound in Medicine＆ Biology，2003，29:1607-1623.

[8]宋章军，徐静，陈恳.一种快速、自动的三维超声标定方法[J]. 北京生物医学工程，200928(6):601-605，613.

[9]Hsu PW，Prager R，Gee A，Treece G，et al.Real-time freehand 3D ultrasound calibration.Ultrasound in Medicine ＆ Biology，2008，34(2):239-251.

[10]Chen TK，Thurston A，Ellis R，et al.A real-time freehand ultrasound calibration system with automatic accuracy feedback and control[J].Ultrasound in Medicine ＆ Biology，2009，35(1):79-93.

[11]Horn BKP.Closed-form solution of absolute orientation using unit quaternion[J].J Opt Soc Am A，1987，4(4):629-642.

[12]Horn BKP，Hilden HM，Negahdaripour S.Closed-form solution of absolute orientation using orthogonal matrices[J].J Opt Soc Am A ，1988，5(7):1127-1135.

[13]朱立人，李文骏，丁辉，等.一种提高N线模型探头标定精度的解算方法[J].中国生物医学工程学报，2012，31(3):337-343.