基于单位四元数的手术导航注册方法*

陶玲 王跃 钱志余

(南京航空航天大学自动化学院，江苏南京210016)

手术导航系统的定位精度是脑立体定向导航手术中的关键问题，决定着手术治疗和预后效果［1］.系统定位精度一般受医学影像质量、定位仪系统误差、注册误差和术中脑组织移位产生的误差等因素的影响.目前，医学影像误差通常不超过1 mm，定位仪系统误差在0.5 mm左右，注册误差一般为2～3mm，而脑组织移位产生的误差对整个系统精度的影响最大，一般为前三者的2～3倍.显然，注册误差和脑组织移位产生的误差是影响系统导航精确度的主要因素［2］.针对脑组织移位引起的误差问题，文献［3］中提出了基于近红外光谱(NIRS)的神经导航穿刺路径纠错方法，并开发了一套基于NIRS的脑立体定向手术导航系统软件平台，能实现病灶靶点识别、实时导航与纠错，这里主要讨论如何提高注册精度的问题，以期解决目前手术导航系统中注册误差较大的缺点.文中通过将单位四元数中的刚体姿态变换理论应用于手术导航系统的电磁定位导航技术中求解最优变换矩阵，来达到明显降低注册误差的目的.

1 电磁定位技术

手术导航系统采用的是基于电磁定位的导航技术.电磁定位子系统在系统中主要起两方面的作用［4］:首先，手术注册时获取手术标记点在实时手术空间中的位置信息，与三维模型相结合，获得从实时手术空间到计算机空间的空间变换关系;然后，根据电磁定位跟踪子系统实时获得手术器械空间方位信息，引导医生完成手术导航治疗.

电磁定位系统主要由磁场发生器、接收传感器以及信号处理器3个主要部分构成.磁场发生器主要由平行的3个线圈组成，系统向3个线圈依次通电，但任意时刻只有一个线圈工作，这样会在空间形成一个交变的电磁场B，并可知此电磁场的强度和分布.传感器一般安放在跟踪目标上，当磁场发生变化时，传感器接收线圈内会产生感应电流，信号处理器通过检测接收到的感应电流的大小，采用简化的磁偶极子模型来获得跟踪目标相对于磁场发生器的空间方位信息［5］.根据磁偶极子的理论，其共轴耦合系数是共面耦合系数的两倍，而且若共轴耦合系数大于零，则两个线圈同时旋转180°后的共面耦合系数小于零，用一个对角阵可表示为D=diag(2，-1，-1).

电磁定位系统与跟踪目标的相对空间位置关系如图1所示.定义磁场发生器的3个轴为参考坐标系，跟踪目标相对磁场发生器的空间位置信息可由二者相对距离r、方位角α以及仰角β描述，跟踪目标相对磁场发生器的空间姿态信息可用滚动角φ、偏航角φ以及俯仰角θ描述，由以上3个空间位置参数和3个空间姿态参数就可以实时表示跟踪目标的空间位置信息，这6个参数也被称为跟踪目标的6 个自由度［6］.

图1 电磁定位系统与跟踪目标的相对空间位置关系Fig.1 Space relation between electromagnetic localization system and tracking target

根据法拉第电磁感应定律可得近距离电磁跟踪系统的模型［7］:

式中:Y为接收矩阵，U为天线输出电压矩阵，M为位置转换矩阵，S为姿态转换矩阵，且Y、U、S、M均为3阶方阵;G为电磁定位系统参数.

刚体的姿态是刚体相对观察者的几何角度关系的简称.最一般的姿态描述方法是矩阵法，但矩阵法有6个约束方程、9个变量，计算复杂，不利于工程运用.欧拉角法也是一种常用的刚体姿态表示法，但在解算参数时会遇到奇异点，导致解算错误.用单位四元数法描述刚体姿态可以克服以上缺点，并且计算方便，不存在奇异点，速度快且稳定性好［8］.

2 基于单位四元数的姿态描述方法

四元数是一种用于研究物体空间运动状态的数学工具.四元数通过4个参数来描述刚体的旋转运动，其表达式为Q=q0+q1i+q2j+q3k，其共轭复数为Q*=q0-q1i-q2j-q3k，由复数知识可知QQ*=，若 Q*Q=QQ*=1，则 Q 为单位四元数［9］.

将式(1)的系统参数忽略后，有N=YU-1=SM，N 为正交矩阵，则 NTN=3X+I，其中，X=MTSTSM，I为3阶单位矩阵，又由ST=S-1，有

式中，E为 3阶基本矩阵，利用四元数法改写式(2)，有:

设p'=pλ为假设的四元数初始值，λ=λ0+λ1i+λ2j+λ3k 为误差四元数.其中，λ0、λ1、λ2和λ3为常数.令，则，当逐渐逼近于0时，λ0相应地逼近于1，则p'越逼近真实值p.

用p'代替p并对式(3)进行变换得

式中右侧最后一项很小，可忽略.根据误差四元数的性质有λ1=0，则等式右边只有两个未知数λ2和λ3，而等式左边是和X有关的已知量，这样由式(4)可以解得λ2、λ3，又由误差四元数定义以及λ1、λ2和λ3可求出λ0，最后便可解得误差四元数λ.由误差四元数λ可以求出新的跟踪四元数p的初始值p'，依次迭代，当λ满足设定的精度条件时，得到的四元数p就是目标的位置参数.

采用迭代算法进行运算，即将前一次获得的跟踪目标的空间位置数据迭代入当前的运算中，可获得表示此时跟踪目标的空间位置参数的四元数值.四元数与跟踪目标空间姿态参数有如下的数学关系［9］:

利用式(5)求解得到的滚动角φ、偏航角φ及俯仰角θ值，便可以推出当前跟踪目标的空间姿态转换矩阵S.求解时，通常迭代运算6次后就可达到符合系统正常工作时的精度要求，运算便捷且便于计算机编程实现.

3 手术导航系统注册内容

脑立体定向手术导航系统的患者注册过程实际上就是实际物理空间(H)、计算机空间(I)和实时手术空间(P)的空间映射过程.

3.1 从物理空间到计算机空间的转换

注册前需在病人脑部标记至少4个不能共面且分布均匀的医学定位标记点［10］，将4个标记点记为F0、F1、F2和 F3，建立物理空间 H，定义 4 个定位标记点在 H 中的坐标分别为:F0(0，0，0)、F1(1，0，0)、F2(0，1，0)和 F3(0，0，1)，记为 FH.

对患者的头部进行CT或MRI扫描，将扫描结果三维可视化，从三维模型中提取出4个标记点，建立对应物理空间H的计算机空间I.根据这4个标记点在三维模型中的位置信息，可得到其在计算机空间 I中的坐标为:F'0(x0，y0，z0)、F'1(x1，y1，z1)、F'2(x2，y2，z2)以及 F'3(x3，y3，z3)，记为 FI.

依据医学断层序列影像数据的三维可视化理论，显然从物理空间H到计算机空间I的空间变换是刚性的［11］，因此可以用如下的齐次转换矩阵T来实现两个空间的转换:

3.2 从物理空间到实时手术空间的转换

以电磁定位仪中心为原点，建立实时手术空间，显然从物理空间到实时手术空间的变换为刚体变换.设齐次转换矩阵RPH为从物理空间H到实时手术空间P的旋转变换矩阵，RPT为物理空间H在实时手术空间P下的平移变换矩阵，

为定位标记点在实时手术空间P下的位置坐标，

为定位标记点在物理空间H下的坐标，则从物理空间到实时手术空间的转换关系为

在电磁定位系统中，由电磁定位仪跟踪患者头部的定位标记点可获得7个元素al(l=0，1，…，6)，al实时表示手术空间相对物理空间的刚体变换，其中，前4个元素表示单位四元数，后3个元素表示位移量，则平移变换矩阵为

根据电磁定位仪工作原理可知a=a0+a1i+a2j+a3k，则对应的旋转矩阵为

将式(8)、(9)代入式(7)，得到从物理空间到手术空间的转换关系如下:

3.3 手术器械在3个空间中的坐标转换

手术开始前，通过定义手术器械中的电磁传感器的位置信息，建立手术器械坐标系FC.在实时手术空间坐标系中，通过电磁定位测得手术器械的旋转矩阵为RPH1、平移矩阵为RPT1，其中测得单位四元数中的7个元素为bl(l=0，1，…，6).设手术器械在实时手术空间下的坐标为FP1，在物理空间下的坐标为FH1，在计算机空间下的坐标为FI1，令

由于 FP1=AFH1，FP1=BFC，则有 FH1=A-1BFC，将该关系式代入式(7)，可获得手术器械在计算机空间的位置坐标为FI1=TA-1BFC.

4 手术导航系统注册算法及其过程

令电磁定位仪采集到的标记点在实时手术空间的坐标为{Pn}，影像三维可视化后得到的标记点在计算机空间的坐标为{In}，n=1，2，…，N，N 为标记点总个数.设W为平移变换矩阵，V为旋转变换矩阵，s为伸缩因子，则当{Pn}与{In}的变换矩阵S满足条件时，当前的变换矩阵SP-I就是{Pn}与{In}的最优变换矩阵，采用单位四元数和齐次坐标理论［12］可求出SP-I，具体算法过程如下.

(1)求出标记点组在实时手术空间和计算机空间下的中心坐标，其中Po为标记点组在实时手术空间中的中心坐标，Io为标记点组在计算机空间下的中心坐标，

(2)将标记点组在两个空间下的坐标转变为关于中心坐标的相对坐标值，其中P'n为第n个标记点在实时手术空间中的相对坐标值，I'n为第n个标记点在计算机空间下的相对坐标值，

(3)根据矩阵论和四元数理论，可用q(q0，q1，q2，q3)表示最优旋转变换矩阵 V［9］:

(5)求解最优平移变换矩阵:

(6)由下式可解出最优变换矩阵:

利用该最优变换矩阵，可解出脑立体定向手术导航系统中从实时手术空间到计算机空间的空间映射关系，完成手术过程中手术器械的注册过程.

5 注册结果分析

5.1 注册点数量对手术导航系统精度的影响

取一个长方体的纸盒子，将其6个面均以坐标纸封上，坐标纸的最小单位为2 mm，以此纸盒作为该实验中脑部的模型，其示意图如图2所示，图中每一方格表示最小单位为2 mm.相对于实物模型，方形纸盒模型操作简便，更利于实验中数据的采集，减小实验中的误差.

图2 实验模型示意图Fig.2 Sketch of experiment model

选取模型某一个表面，在该面内选定边长为80mm的正方形.注册时，注册点对应于医学标记点，为便于分析注册点个数变化所带来的影响，实验中将取点范围限制在一个表面的正方形中.当注册点为4个时，可以分别取某一个面上的4个顶点，得到软件平台上4个点在计算机空间的坐标值;当注册点为5个时，在4个注册点包围的面上再任取一点，获取计算机空间坐标;当注册点为6个时，在同样的面上再增加一点，依此类推，直至取到15个点.利用文中提出的注册算法对注册点进行注册，记录实验中注册点数量改变时所对应的系统定位误差.每个注册点误差值的定义就是计算机空间的坐标值和实际手术空间的坐标值的差值的绝对值，注册点在实际手术空间的坐标就是实验模型上注册点的实际坐标值，所开发的导航软件平台上则会给出相应注册点在计算机空间的坐标值，导航软件在Windows XP操作系统下使用Visual C++6.0开发.结果如表1所示.

表1 注册点数量与定位误差Table 1 Numbers of registration points and lo calization error

从表1可以看出，系统的注册误差均在2mm以下，显然本研究所提方法的注册误差有了明显的下降，注册精度得到了提高;同时可以看出，当注册点数量越多时，系统的注册误差会逐渐变小.当注册点个数多于7时，系统注册误差减小幅度较大.但实际进行手术时，注册点越多，医生进行手术注册的工作量越大，并且发生错误注册的概率也会增大，因此注册点个数也不应太多.结合临床应用，通常进行手术导航时选取7～12个点为最佳.

5.2 注册点分布范围对手术导航系统精度的影响

依然利用图2的实验模型，制作边长分别为8、20、40、60、80、100 和 120 mm 的 7 个正立方体，正立方体共有8个顶点，在每个实验中，将正立方体的1个顶点定义为坐标系原点，其余7个顶点为一组，作为手术注册的点集，此时该正立方体范围就是注册点所包围的范围.每进行一次实验，则增加正立方体的边长并依次取7组点集进行手术注册.8、20、40、60、80、100 和120mm 对应的导航定位误差分别为0.47、1.04、1.27、1.40、1.59、1.71 和 1.89mm.显然，系统的注册误差同样均在2 mm以下，并且随模型正立方体边长的增加系统注册误差也增大，即在注册点个数相同的情况下，注册点所包围的范围边长越长，系统注册误差也越大.因此注册点的选取是越紧密围绕病灶区域越好.

综合多次实验结果发现:保证手术器械尖端一直都在注册点有效作用范围内是提高手术导航的注册精度的必要条件，通常需要用7～12个注册点围绕病人头部病灶部位，注册点需排放对称、密度均匀，同时包裹范围适中.如果将注册点包裹患者整个头部，则会造成导航精度误差较大，目标偏移严重，影响手术效果;如果注册点包裹范围较小，未完全包裹病灶，则手术时手术器械尖端可能超出注册点包裹的范围，也会导致导航精度误差较大.应用文中的注册算法做了大量的模型试验，注册误差均在2mm以下，注册速度均可达到实时，并且稳定性和重复性好.

6 结语

手术注册过程是脑神经外科手术的重要环节，注册精度是导航手术治疗效果的关键.文中重点分析了基于电磁定位的脑立体定向手术导航系统中的实际物理空间、实时手术空间、计算机空间的建立以及三者之间的映射变换;在此基础上，提出了根据单位四元数和矩阵轮理论求解从实时手术空间到计算机空间的最优转换矩阵新算法，减少了注册误差，加快了迭代速度，提高了系统的准确性、实时性和稳定性;最后通过实验数据分析验证了所提算法可以有效地提高手术导航系统的定位精度，为临床应用提供帮助.目前系统平台和注册算法还需要在大量动物和临床实验中不断修正和完善，整个系统的稳定性和软件平台的操作方便性还需进一步提高.

［1］Bjartmarz H，Rehncrona S.Comparison of accuracy and precision between frame-based and frameless stereotactic navigation for deep brain stimulation electrode implantation ［J］.Stereotactic and Functional Neurosurgery，2007，85(5):235-242.

［2］DHaese P F，Pallavaram S，Konrad P E，et al.Clinical accuracy of a customized stereotactic platform for deep brain stimulation after accounting for brain shift［J］.Stereotactic and Functional Neurosurgery，2010，88(2):81-87.

［3］朱晓芬，陶玲，钱志余.基于近红外光谱的神经导航穿刺路径纠错方法研究［J］.生物医学工程学杂志，2010，27(4):887-892.Zhu Xiao-fen，Tao Ling，Qian Zhi-yu.Study on navigation path correction in neuro-navigation based on near infrared spectroscopy ［J］.Journal of Biomedical Engineering，2010，27(4):887-892.

［4］田和强，吴冬梅，王继虎，等.基于电磁定位的手术导航探针可视化与实时跟踪技术［J］.机器人，2011，33(1):59-65.Tian He-qiang，Wu Dong-mei，Wang Ji-hu，et al.Visualization and real-time tracking technologies of the probe used in surgical navigation based on electromagnetic positioning［J］.Robot，2011，33(1):59-65.

［5］李雪丽，王毅.电磁定位系统中的定位算法研究［J］.微计算机信息，2009，25(7):298-299.Li Xue-li，Wang Yi.The quaternion solution in the remote object position and orientation determining system ［J］.Microcomputer Information，2009，25(7):298-299.

［6］杨杰，邹锦，李跃波.六自由度电磁跟踪定位系统的设计研究 ［J］.自动化与仪表，2007，22(5):34-37.Yang Jie，Zou Jin，Li Yue-bo.Research on six-stance magnetic position and orientation tracking system［J］.Automation ＆ Instrumentation，2007，22(5):34-37.

［7］Chao H，Meng M Q H，Mandal M.A linear algorithm for tracing magnets position and orientation by using 3-axis magnetic sensors［J］.IEEE Transactions on Magnetics，2007，43(12):4096-4101.

［8］Paperno E，Sasada I，Leonovich E.A new method for magnetic position and orientation tracking［J］.IEEE Transactions on Magnetics，2001，27(4):1928-1940.

［9］孙冬梅，田增山，韩令军.捷联惯导系统中四元素法求解姿态角仿真模拟［J］.弹箭与制导学报，2009，29(1):51-53.Sun Dong-mei，Tian Zeng-shan，Han Ling-jun.Simulation on quaternion calculate attitude angle of the strap down inertial navigation system ［J］.Journal of Projectiles，Rockets，Missiles and Guidance，2009，29(1):51-53.

［10］Galloway R L.The process and development of imageguided procedures［J］.Annual Review Biomedical Engineering，2001，3(1):83-87.

［11］Casap N，Wexler A，Tarazi E.Application of a surgical navigation system for implant surgery in a deficient alveolar ridge postexcision of an odontogenic myxoma［J］.Journal of Oral and Maxillofacial Surgery，2005，63(7):982-988.

［12］Berthold K P.Closed-form solution of absolute orientation using unit quaternions ［J］.Journal of the Optical Society of America，1987，4(4):629-642.