【www.zhangdahai.com--其他范文】
楼狄凯 ,张丹 ,梁华庚
(1.浙江工业大学信息工程学院,浙江杭州 310014;2.华中科技大学同济医学院附属协和医院,湖北武汉 430022)
穿刺手术是现代医学中使用较为普遍的一种微创手术,具有对患者创伤小,伤口恢复快的优点,广泛应用于病灶组织采样、药剂注射、近距离放化疗等领域,是医学诊断、治疗的重要手段.传统刚性针由于组织变形、针尖受力不均匀等带来的较大误差,在穿刺过程中难以对针尖位置进行有效的控制,容易造成穿刺失败[1].穿刺任务往往需要多次重复的进行才能完成,这会延长手术时间,增加患者痛苦,影响治疗效果.且在一些进针过程中需要避开神经、血管骨骼等障碍物,传统刚性针很难满足这个要求.针对传统穿刺针的缺点,Webster III等人[1]提出了柔性针的概念.柔性针采用弹性较好的合金材料,针尖采用非对称斜角的形式,通过组织施加在针尖的非对称力使针尖轨迹发生偏转,通过控制针尖斜面的朝向控制针尖的偏转方向,从而实现对针尖位置的控制.针尖位置的准确控制对成功穿刺至关重要.Webster III等人[1]首先提出了柔性针的概念,根据柔性针的运动形式建立了柔性针的“自行车模型”.Kallem等人[2]将图像引导与柔性针穿刺系统相结合,在Webster III等人提出的“自行车模型”基础上引入了图像引导系统,通过基于观测器的反馈控制来实现对柔性针针尖位置的控制,并给出“自行车模型”的状态空间模型.Park 和Kim等人[3]在“自行车模型”的基础上将柔性针模型简化为“独轮车模型”.Zhao等人[4]以实际轨迹和误差分析为依据,对“独轮车模型”进行改进,提出了带返程独轮车模型.Abayazid等人[5]以“独轮车模型”为基础提出了一个包含轨迹切割角的柔性针运动学模型,同时结合针与组织之间的相互作用,提出了基于力学的柔性针模型.
需要指出的是,柔性针在实际穿刺的过程中,针尖的位姿往往无法完全获取.Benam和Talebi等人[6]提出了一种高增益观测器对针尖状态进行估计,并将估计状态用于反馈控制率设计,通过高增益观测器和输出反馈控制可以提高柔性针穿刺的精度.Huo等人[7]将无迹卡尔曼滤波算法(unscent Kalman filter,UKF)应用到针尖位姿的估计中,并提出了一种基于针尖可达范围的控制方法对柔性针穿刺过程进行控制.Favaro等人[8]提出了一种使用基于2-D运动学的扩展卡尔曼算法(extend Kalman filter,EKF)估计针尖的姿态,并通过实验验证了该算法的性能.Kaya等人[9]提出了一种基于Gabor滤波器的图像处理算法对针尖进行估计,并用卡尔曼滤波器降低估计噪声,实现针尖的实时跟踪估计.Agarwal等人[10]在柔性针穿刺过程中将图像处理技术与卡尔曼滤波算法相结合来提高算法的估计精度.Yan等人[11]提出了一种结合改进的压缩跟踪算法和自适应卡尔曼滤波器的跟踪算法,以实现在2D超声环境下的柔性针针尖实时跟踪.Mignon等人[12]人提出了一种新的3D超声柔性针检测方法,利用一种新的基于力学的预测模型和卡尔曼滤波器提高已有算法的鲁棒性减少检测噪声和反误检测.Lapouge等人[13]在“独轮车模型”的基础上提出了一种多速率无迹卡尔曼滤波器应用于柔性针针尖的转向,实现来自3D图像、机器人传感器和术前弹性成像测量的异步数据融合.Favaro和Secoli等人[14]提出了一种可操纵针全位姿估计方法,该方法在建立了一种新型的PBN三维运动学模型,并通过扩展卡尔曼滤波来预测针在插入过程中的完全位姿.Zhao等人[15]在“自行车模型”的基础上,分别采用自适应卡尔曼滤波算法(adaptive Kalman filter,AKF)、无迹滤波算法和自适应无迹卡尔曼滤波算法(adaptive unscent Kalman filter,AUKF)3种算法对柔性针的针尖位置进行估计,并比较了3种算法的差异.然而在实际穿测过程中,由于活体组织本身的复杂性使得穿刺过程存在不确定性[1],同时由于呼吸作用等影响,超声成像等量测设备工作环境产生变化,导致目标观测噪声的方差实时变化[16].现有的估计算法并未考虑柔性针实际穿刺过程中的不确定性问题和超声成像等设备存在的量测噪声统计特性不准确性问题,这也是本文的研究动机.
本文针对柔性针穿刺过程存在的不确定性问题以及超声成像等设备存在的量测噪声统计特征不准确性问题,提出了一种带有噪声估计器的自适应奇异值分解无迹卡尔曼滤波算法(adaptive singular value decomposition unscented Kalman fliter,ASVD–UKF).该算法采用自适应因子实时修正动力学模型误差,通过奇异值分解抑制系统状态协方差矩阵的负定性,利用Sage-Husa估计器在线估计系统噪声的统计特性,提高滤波精度,减小系统状态的估计误差.并且将新算法应用于带有不定性的柔性针穿刺模型进行计算仿真,结果证明,新算法能有效改善滤波性能,提高系统状态的估计精度.
Webster III等人首先提出了柔性针穿刺过程的自行车模型[1].当斜面的柔性针插入组织时,在针尖不对称力的作用下,针尖轨迹将遵循一个曲率为κ的圆弧轨迹.如果针的材料属性预期穿刺组织的属性相匹配,针体的轨迹将几乎完全遵循针尖的轨迹,即针尖的轨迹可以完全代表针体轨迹.
如图1所示,定义世界坐标系A和体坐标系B,C,带斜面的柔性针由两个针基的控制输入驱动u1和u2,即针的进给和自转.此系统为非完整约束系统,针尖在体坐标系B中的速度可以表示为
图1 运动学自行车模型[1]Fig.1 The kinematic bicycle model[1]
其中:u1为柔性针穿刺速度,使针在体坐标系z轴上的移动,同时由于针尖斜面的不对称作用力产生针体关于x轴的旋转;u2为柔性针针基旋转速度,使针体在体坐标系z轴方向上的旋转;V1,V2为对应的控制向量,
其中:ei为单位向量,κ为穿刺曲率.定义一组广义坐标q=[x y z α β γ]T,用于表示针尖的姿态,其中x,y,z表示针尖在坐标系A中的位置.α,β,γ分别为绕坐标系A中x,y,z轴的角度.针尖在世界坐标系中的速度可以表示为
其中RT为坐标系A,B之间的旋转矩阵,针尖的运动学模型可以写成如下形式:
当针尖稳定至y–z平面时,系统状态[y z α]对其他状态[x β γ]的动态过程没有影响,柔性针模型可以降阶为下述三维非线性模型:
值得注意的是,当针体插入活体组织时,由于组织本身并不是均一介质,具有不定性,这有可能会使针尖轨迹的曲率发生变化而并非保持在一个常数[1].本文记曲率的变化为∆κ,实际的曲率为
记s=[s1s2s3]T=[x β γ],设针基插入速度恒定为u1=1,柔性针运动模型可以表示为
其中s1为针尖的x轴坐标,可通过超声成像设备直接量测.实际穿刺过程中,由于穿刺环境复杂性和传感器噪声,量测值将受到量测噪声影响而偏离实际真值.
穿刺模型为单输入,单输出的仿射非线性系统,对非线性模型进行反馈线性化处理.记z=[z1z2z3]T,
其中:∆a1=∆κ·κ,∆a1=∆κ/κ.
设计线性反馈控制率如下:
反馈线性化形式的状态方程为
其中:
本文针对柔性针穿刺过程存在的不确定性问题以及超声成像等设备存在的量测噪声统计特征不准确性问题,提出了一种带有噪声估计器的自适应奇异值分解无迹卡尔曼滤波算法.该算法采用自适应因子实时修正动力学模型误差,通过奇异值分解抑制系统状态协方差矩阵的负定性,利用Sage-Husa估计器在线估计未知系统噪声的统计特性,减小状态估计误差.
3.1 奇异值分解
由于迭代过程中的截断误差、模型扰动等因素,误差协方差矩阵Pk容易失去对称正定性,致使算法无法进行cholesky分解[17].本文采用奇异值分解(singular value decomposition,SVD)替代cholesky分解进行Sigma采样.奇异值分解是一种具有良好数值稳定性的矩阵分解算法,其特点是不需要被分解矩阵是对称正定矩阵[17].具体如下:
假设P∈Rm×n(m≥n),则矩阵P的SVD分解可以表示为
其中:U∈Rm×m,V∈Rn×n分别为A的左、右奇异向量,Λ∈Rm×n,S=diag{s1,s2,···,sr},s1≥s2≥···≥sr为P的奇异值.
3.2 自适应因子
考虑如下非线性系统:
其中:xk为k时刻的系统状态向量,ωk为系统噪声,yk为观测向量,φk为量测噪声,f(·)为状态模型,h(·)为观测模型.过程噪声和测量噪声的统计特性分别为ωk ∼N(0,Q),φk+1∼N(0,R).以预测残差为变量,构造状态模型的误差判别统计量及自适应因子.预测残差表示为
其中:αk表示自适应因子,c为检验常数,通常情况下1 对于穿刺过程中观测噪声的不定性,用Sage-Husa估计器对进行估计.对于时变噪声,应强调新近数据的作用,对于陈旧数据的作用应该渐渐遗忘和消失,不同于算术平均,应对和式中的每项乘以不同的加权系数,采用渐消记忆指数加权法实现,次优Sage-Husa估计器与卡尔曼滤波器结合,可以得到改进的适用于带时变噪声特性的线性系统的自适应滤波器[16].可由时变噪声估计器递推获得 其中dk−1=(1−b)/(1−bk),b为遗忘因子,通常b的范围为0.95 2)基于奇异值分解构建Sigma 点对协方差矩阵Pk−1进行奇异值分解 式中:Uk−1,Vk−1列向量分别为方差矩阵Pk−1的左、右奇异向量,Λ为方差矩阵Pk−1的奇异值矩阵.经过SVD分解得到的Sigma点矩阵如下: 3)时间更新 4)量测更新 注1近年来,卡尔曼滤波器在许多工程系统中得到了实际应用.例如,Julier等人[18]提出了基于无迹变换的卡尔曼滤波算法,利用无迹变换估计系统状态的均值和协方差,能够较好的处理非线性系统的滤波问题.与扩展卡尔曼滤波算法相比,无迹滤波算法不需要计算Jacobian矩阵,计算更为简单,估计精度更高于扩展卡尔曼滤波算法[19].但由于无迹滤波算法可能存在协方差阵的负定性,致使滤波病态问题的出现.Wang等人[17]提出了一种自适应奇异值分解无迹滤波算法应用于车辆质心侧偏角估计,通过用奇异值分解代替cholesky分解,抑制系统状态,协方差矩阵负定性问题,同时通过自适应因子调节动力学模型误差,提高算法的跟踪能力.但由于柔性针穿刺环境的复杂性,会导致观测噪声的方差发生实时变化.卡尔曼滤波器在噪声环境变化时会出现滤波发散[16].Shi等人[20]将Sage-Husa估计器与无迹滤波算法结合,提出一种自适应的无迹滤波算法,在线估计噪声的统计特性,减小状态估计误差.本文将奇异值分解无迹卡尔曼滤波算法与Sage-Husa估计器相结合,应用于柔性针的穿刺状态估计.相较于现有的UKF,AUKF算法,本文提出的自适应的SVD–UKF算法能够利用自适应因子实时修正穿刺过程中曲率不定性带来的估计偏差,同时用SVD分解替代cholesky分解,抑制穿刺估计时滤波病态问题的出现,通过Sage-Husa估计器减小由于超声成像系统等设备带来的观测噪声干扰,改善滤波性能,提高系统状态的估计精度. 本文采用MATLAB软件进行仿真,通过比较针尖位姿的估计与真实位姿验证观测器的估计精度与鲁棒性.为了评估曲率不定性对观测器估计精度的影响,仿真模拟了在不同的曲率偏差下无迹卡尔曼滤波算法对位姿估计的情况.设定系统初始状态为 采样时间∆T=0.05 s,轨迹曲率κ=1/122.设定曲率偏差在±2%范围内随机波动,UKF估计曲线如图2所示. 图2 ±2%曲率偏差UKF估计曲线Fig.2 The UKF estimation with±2%curvature deviation 图2中红色虚线和蓝色虚线分别表示没有曲率偏差和存在±2%曲率偏差情况下估计器的跟踪曲线.当曲率偏差在±2%范围内随机波动时,观测曲线难以精确跟踪实际的轨迹曲线,观测器估计精度下降.为了精确描述曲率不定性对估计精度的影响,仿真分别在曲率偏差为±2%和±5%的情况下进行模拟50次,得到如表1所示各情况的数据统计结果. 表1 不同曲率偏差下的UKF估计误差Table 1 The UKF estimation errors under different curvature deviations 从表1可以看出,当曲率偏差存在±2%的随机波动时,UKF算法的平均估计误差从0.0538 mm增大到0.1008 mm,估计误差增大87.36%.当曲率偏差以±5%波动时,估计平均误差增大至为0.1242 mm,均方根误差增大至0.207.由于曲率不定性带来的偏差将会降低观测器的估计精度. 为了评估噪声统计特征不准确的情况对滤波算法估计精度的影响,仿真分别模拟了量测噪声在±5%和±10%随机偏差情况下,UKF滤波算法的估计情况. 图3中红蓝绿三色虚线分别表示在没有噪声统计特征偏差情况下UKF估计器的跟踪曲线和在±5%噪声偏差情况下以及在±10%噪声偏差情况下估计器的跟踪曲线.从图3可以看出,当噪声统计特征存在偏差时,无迹滤波算法难以精确跟踪真实曲线,噪声偏差越大,跟踪曲线波动越大.为了定量描述噪声统计特征偏差估计精度的影响,仿真分别在偏差为±5%和±10%的情况下模拟50次,统计结果如表2所示. 图3 噪声偏差下UKF 估计曲线Fig.3 The UKF estimation with measurement noise deviation 表2 不同噪声偏差下的UKF估计误差Table 2 The UKF estimation errors under different measurement noise deviations 从表2可以看出,当噪声存在±5%的偏差随机波动时,UKF观测器的平均估计误差从0.0946 mm增大到0.1422 mm.当噪声偏差以±10%波动时,估计平均误差为0.1738 mm,均方根误差增大至0.206.量测噪声统计特征的不准确将直接影响滤波算法的估计精度. 将自适应SVD–UKF算法应用到柔性针穿刺系统,选择Sigma 采样参数α=1,β=1,η=0,ρ=1.22,自适应参数b=2,c=2.0.设定系统噪声统计参数,R=0.02以及量测噪声统计特征存在±5%随机偏差.将估计效果与现有的UKF算法和AUKF算法[15]进行仿真对比,结果如图4–7所示. 图4–6分别为UKF估计器、AUKF估计器和ASVD–UKF估计器的状态估计结果.在曲率存在±5%的随机偏差且噪声统计特征不完全准确的情况下,UKF算法估计曲线存在较大的估计偏差,无法精确跟踪真实曲线.AUKF算法较UKF 算法能更好的跟踪真值,但依旧无法修正由于曲率偏差带来的估计误差.自适应SVD–UKF算法较AUKF算法能更快的收敛,更好地跟踪实际曲线. 图4 UKF滤波算法跟踪曲线Fig.4 The estimated result of UKF filter 图5 AUKF滤波算法跟踪曲线Fig.5 The estimated result of AUKF filter 图6 ASVD–UKF滤波算法跟踪曲线Fig.6 The estimated result of ASVD–UKF estimator 图7表示各个滤波算法跟踪实际状态的误差曲线.自适应SVD–UKF从误差曲线上相较于UKF算法和AUKF算法,误差波动更小且更趋近于真值,能较好的跟踪真实轨迹.为了准确地比较3种滤波器算法,本文采用均方根误差的形式描述了滤波器的预测精度.每个算法模拟50次,取平均值.各个滤波算法估计状态与真实系统最大误差、平均误差、均方根误差如表3所示. 表3 3种估计算法的估计误差Table 3 The estimation errors of three filters 图7 3种滤波算法估计误差Fig.7 The estimation errors of three filters 由仿真结果可知,在相同初始条件下,ASVDUKF算法的平均误差和均方误差明显小于无迹滤波算法,与AUKF算法算法相比,平均误差减小了51.96%,均方误差减小了82.46%.自适应奇异值分解UKF算法通过自适应因子修正由于曲率偏差带来的估计误差,并利用Sage-Husa估计器在线估计系统噪声的统计特性,有效提高了算法的估计精度. 以下讨论奇异值分解和cholesky分解在无迹滤波算法中的区别.基于cholesky分解的无迹滤波算法在进行Sigma采样近似时,容易出现协方差矩阵非正定的情况,致使病态问题出现[17].病态问题的出现将致使滤波算法无法对状态进行估计跟踪.在不出现病态问题的情况下,分别使用SVDKF算法和UKF算法对系统状态进行跟踪估计.在不考虑模型不定性和噪声偏差的前提下,两算法跟踪情况如图8所示. 图8 SVDKF与UKF跟踪曲线Fig.8 The estimated result of UKF and SVDKF 从仿真结果,可以看出SVDKF算法相较于UKF算法,跟踪曲线更加平滑,能更好的吻合真实状态.对跟踪效果进行数据统计,得到如表4所示结果. 表4 UKF算法与SVDKF算法估计误差Table 4 The estimation errors of UKF and SVDKF 从统计结果,可以看出SVDKF算法相较于UKF算法具有更小的估计误差.以SVD分解代替cholesky分解的滤波算法具有更好的数值稳定性和鲁棒性. 柔性针在实际穿刺过程中,由于生物组织的非均匀性,会造成穿刺曲率发生变化,导致柔性针穿刺过程存在参数不确定性问题.本文针对该问题建立带有曲率不定性∆κ的柔性针穿刺的状态空间模型.同时针对柔性针穿刺过程存在的不定性问题以及超声成像设备等量测噪声统计特征不准确性情况,提出了一种带有噪声估计器的自适应奇异值分解无迹卡尔曼滤波算法.特别地,采用自适应因子实时修正动力学模型误差,抑制模型中误差不定性带来的影响,用奇异值分解替代cholesky分解,抑制系统状态协方差矩阵的负定性.同时利用Sage-Husa估计器在线估计噪声的统计特性,减小由于设备存在的量测噪声统计特征不准确情况带来的误差.最后,将新的自适应SVDUKF算法应用于带有曲率不定性的柔性针穿刺模型进行计算仿真,结果表明,新算法能有效改善滤波性能,提高穿刺状态的估计精度.3.3 次优Sage-Husa估计器
3.4 ASVD–UKF算法步骤
