一种基于时空相关性的心脏表面目标点运动预测方法

摘要

本发明公开了一种基于时空相关性的心脏表面目标点运动预测方法，利用心脏表面目标点运动轨迹的时间自相关性，以及与周围辅助点之间的空间相关性，建立预测目标点的线性模型GLM，再基于当前GLM模型下获得预测值，从而解决了长时间跨度预测时，由于时间自相关性减弱导致的预测误差急剧增大的问题。

说明书

技术领域

本发明属于运动预测技术领域，更为具体地讲，涉及一种基于时空相关性 的心脏表面目标点运动预测方法。

背景技术

近年来，机器人技术越来越多的用于微创手术领域，用于减轻患者痛苦， 降低手术医生的工作强度，提高手术精度和降低手术难度。世界各地的研究机 构都在积极开展机器人辅助外科手术技术的研究。而很多先进的机器人辅助手 术技术的临床应用，都需要建立在对手术器官表面目标点的精准定位之上。

在现有技术中，机器人辅助心脏外科手术中的“心跳同步”技术，手术中需要 实时跟踪心脏表面目标点的运动，并主动控制手术器械与其同步运动，从而为 医生提供一个虚拟稳定的操控环境，使得手术医生在进行切割或缝合等精准手 术操控时，无需再手动的克服快速心跳运动的干扰。为此，人们采用了不同的 传感器系统对心脏表面目标点的运动进行测量，如基于立体内窥镜的视觉测量 系统、基于超声波探测器的测量系统和基于激光探测器的测量系统等。

然而，在实际临床应用中，想要实时、准确和稳定的测量和跟踪快速跳动 的心脏表面目标点的运动并不容易。手术器械的遮挡、器官表面软组织流血、 电刀切割产生的烟雾等手术过程中的各种动态干扰都可能中断测量系统对目标 点的测量，使其暂时无法获得测量值。为了弥补测量空白并在干扰消失后重启 测量系统，就需要对目标点的运动进行精准的预测。另外，从提高机械手臂控 制精度的角度，也需要运动预测技术。由于机械手臂和手术器械自身的质量和 控制环节的时滞，要驱动和控制其与快速运动的心脏表面目标点同步运动时， 仅仅依靠常规的反馈控制无法实现，需要借助预测控制技术，提前预测目标点 的运动，给出相应的控制量。

现有的运动预测方法，多从心跳运动的准周期性出发，基于目标点运动轨 迹在时间上的自相关性，利用最近的若干历史测量值通过建立预测模型，对当 前或未来的目标点进行预测。这类方法在进行长时间跨度的预测时，由于长时 间未获得有效的测量值，随着时间的推移，较早时刻获取的历史测量值与目标 点当前位置之间的相关性急剧减弱，形成了利用预测值的预测值的…预测值来 进行预测的局面，预测误差会不断的累计和放大。另一方面，心跳运动的频率 和幅度等特性是随时间变化的，目标点历史测量值与当前的运动之间的相关性 也会随时间变化，因此，在长时间跨度的预测中，这种基于时间自相关性的预 测方法会引入较大的预测误差。

发明内容

本发明的目的在于克服现有技术的不足，提供一种基于时空相关性的心脏 表面目标点运动预测方法，利用目标点与其周围辅助点在空间上的相关性，解 决长时间跨度预测，同时降低预测误差急剧增大的问题。

为实现上述发明目的，本发明一种基于时空相关性的心脏表面目标点运动 预测方法，其特征在于，包括以下步骤：

(1)、建立并初始化线性预测模型GLM

(1.1)、建立GLM模型：利用心脏表面目标点p的N个历史测量值和M个 辅助点建立L＝M+N+1阶的GLM模型，用方程表示为：

$\hat{p}\left(k\right)=Q(k-1)w(k-1)$

其中，为三维列向量，表示心脏表面目标点p在k时刻的三维空间坐标 预测值；Q(k-1)是3×L维的模型设计矩阵，可表示为：

它是由心脏表面目标点p在k时刻之前的N个历史测量值： p(k-N),p(k-N+1),...,p(k-1)和M个辅助点在k-1时刻的测量值 h₁(k-1),h₂(k-1),...,h_M(k-1)，以及一个元素全为1的三维列向量组成；w(k-1)是k-1 时刻的模型参数，是由L个权值系数组成的列向量；

(1.2)、初始化GLM模型：令k＝1，将k＝1之前的N个历史测量值和0时刻 的M个辅助点测量值都初始化为0向量，即：

p(1-N)＝p(2-N)＝…＝p(0)＝0

h₁(0)＝h₂(0)＝…＝h_M(0)＝0

并将0时刻的模型参数初始化为0向量，即w(0)＝0，将方差矩阵初始化为 V(0)＝10000I_L×L，其中I_L×L表示L×L维的单位矩阵；

(2)、判断k时刻心脏表面目标点p是否测量成功，若测量系统提供测量值， 则测量成功，获得该时刻心脏表面目标点p的测量值p(k)，然后执行步骤(3)； 反之测量失败，则执行步骤(4)；

(3)、更新GLM模型参数

基于迭代最小二乘滤波(RLS)原理，利用当前获取的心脏表面目标点p的 测量值p(k)，更新模型参数w(k)及其方差矩阵V(k)，待更新完毕后，跳入步骤(5)；

(4)、基于当前GLM模型获得心脏表面目标点p的预测值，具体如下：

$\hat{p}\left(k\right)=Q(k-1)w(k-1)$

用预测值代替测量值，即令然后执行步骤(5)；

(5)、更新GLM模型设计矩阵

利用当前时刻的心脏表面目标点p的测量值p(k)和M个辅助点的测量值 h₁(k-1),h₂(k-1),...,h_M(k-1)更新模型设计矩阵，得到

(6)、当前时刻值k加1，即：k＝k+1，再返回步骤(2)，进入下一时刻k+1 的处理流程。

本发明的发明目的是这样实现的：

本发明基于时空相关性的心脏表面目标点运动预测方法，利用心脏表面目 标点运动轨迹的时间自相关性，以及与周围辅助点之间的空间相关性，建立预 测目标点的线性模型GLM，再基于当前GLM模型下获得预测值，从而解决了 长时间跨度预测时，由于时间自相关性减弱导致的预测误差急剧增大的问题。

同时，本发明基于时空相关性的心脏表面目标点运动预测方法还具有以下 有益效果：

(1)、本发明充分利用心脏表面目标点与其周围辅助点之间的空间相关性， 提高预测精度和鲁棒性。现有技术只考虑目标点运动在时间上的自相关性，利 用目标点的历史测量值预测未来运动。这类方法在进行长时间连续多步的预测 时，预测误差会随时间推移急剧增大。本发明方法不仅利用目标点在时间上的 自相关性，还利用了目标点和周围辅助点的运动在空间上存在的相关性进行预 测，因而在进行长时间大跨度的预测时，仍能保持极高的预测精度。另一方面， 由于目标点与辅助点之间存在天然的物理联系，而这种物理联系几乎不随时间 变化；因而可以获得更准确的预测结果。

(2)、本发明方法利用了简单的一般线性模型(GLM)对上述时空相关性 进行建模，所设计的递归流程，在测量系统成功获取测量值时，能够实时在线 的更新模型参数，确保预测模型的准确，运算的复杂度低，实时性好，可以满 足手术过程中对心脏表面目标点的实时跟踪和预测。

附图说明

图1是本发明基于时空相关性的心脏表面目标点运动预测方法流程图；

图2是本发明实施例中心脏表面目标点和辅助点的位置示意图；

图3是本发明实施例中心脏表面目标点测量失败时的运动预测结果示意图。

具体实施方式

下面结合附图对本发明的具体实施方式进行描述，以便本领域的技术人员更 好地理解本发明。需要特别提醒注意的是，在以下的描述中，当已知功能和设 计的详细描述也许会淡化本发明的主要内容时，这些描述在这里将被忽略。

实施例

本发明是利用心脏表面目标点运动的时间自相关性和与周围辅助点的空间 相关性，预测目标点未来的运动，预测的结果可用于当测量系统测量失败时， 对目标点位置做出准确估计，或用对手术机器人机械臂的预测控制。

具体的说是通过在目标点测量值与其历史测量值以及辅助点测量值之间建 立GLM模型的方式实现。测量系统同时对目标点和若干个辅助点进行测量，外 科手术中，目标点由医生根据病人病情确定，其位置无法随意选择；而辅助点 则可根据需要，选择具有明显特征(容易被测量)且不易被手术中的动态因素 影响的自然特征点或放置了人工测量标识物的点。当无法获取目标点测量值时， 仍可利用辅助点的最新测量值，对目标点进行准确预测，解决现有方法在长时 间预测时误差急剧增大的问题。

在本实施例中，测量目标点运动的测量系统可以是基于立体内窥镜的视觉 测量系统、基于超声波探测的测量系统或基于激光和视觉的测量系统等。无论 何种测量系统，在其能对目标点进行直接测量时，本发明所述的方法会在每个 测量采样时刻，利用最新获取的目标点和辅助点测量值在线的更新GLM模型的 参数和设计矩阵，确保GLM模型可以准确描述目标点与其历史测量值以及辅助 点测量值之间时空相关性；当测量系统无法获取测量值时，则通过GLM模型预 测目标点的位置，弥补测量系统的测量空白。因此，对于任何测量系统因任何 原因(如手术器械遮挡目标点、外部噪声干扰等)短时无法获取目标点位置时， 都可采用本方法进行预测，获得精确的目标点位置预测值。

下面结合附图和实施例对本发明作进一步说明。

图1是本发明基于时空相关性的心脏表面目标点运动预测方法流程图。

在本实施例中，利用本发明所述的方法对目标点近4个时刻的历史测量值 和周围3个辅助点，建立8阶GLM模型进行预测。

下面结合图1所示的流程图，对预测的过程进行详细说明，具体包括如下 步骤：

S1、建立并初始化GLM预测模型(k＝1)，包括如下两个子步骤：

S1.1、建立GLM模型：利用心脏表面目标点p的4个历史测量值和3个辅 助点，建立8阶GLM预测模型，用方程表示为：

$\hat{p}\left(k\right)=Q(k-1)w(k-1)$

其中，为三维列向量，表示心脏表面目标点p在k时刻的三维空间坐标 预测值；Q(k-1)是3×8维的模型设计矩阵，可表示为：

Q(k-1)＝[p(k-4)p(k-3)p(k-2)p(k-1)h₁(k-1)h₂(k-1)h₃(k-1)1]

其中，p(k-4),p(k-3),p(k-2),p(k-1)是心脏表面目标点p在k之前的4个历史 测量值，h₁(k-1),h₂(k-1),h₃(k-1),h₄(k-1)是4个辅助点在k-1时刻的测量值；w(k-1) 是k-1时刻的模型参数，是由8个权值系数组成的列向量。

S1.2、初始化GLM模型：令k＝1，对k＝1时刻之前的心脏表面目标点p和辅 助点测量值、模型参数及其方差矩阵进行初始化：

p(-3)＝p(-2)＝p(-1)＝p(0)＝0

h₁(0)＝h₂(0)＝h₃(0)＝0

w(0)＝0

V(0)＝10000I_8×8

其中，I_8×8表示8×8维的单位矩阵。

S2、判断k时刻心脏表面目标点p是否测量成功，若测量系统提供测量值， 则测量成功，获得该时刻心脏表面目标点p的测量值p(k)，然后执行步骤S3；反 之测量失败，则执行步骤S4。

S3、更新GLM模型参数

基于迭代最小二乘滤波(RLS)原理，利用当前获取的心脏表面目标点p的 测量值p(k)，更新模型参数w(k)及其方差矩阵V(k)，具体包括如下步骤：

S3.1、计算模型当前的预测误差和计算增益矩阵：

预测误差：$e=p\left(k\right)-\hat{p}\left(k\right)$

增益矩阵：K＝V(k-1)Q^T(k-1)[λI_3×3+Q(k-1)V(k-1)Q^T(k-1)]

其中遗忘因子λ＝0.98，I_3×3表示3×3维的单位矩阵；

S3.2、利用预测误差和增益矩阵更新模型参数及其方差矩阵

模型参数：w(k)＝w(k-1)+K·e

方差矩阵：V(k)＝λ^-1V(k-1)-λ^-1KQ(k-1)V(k-1)

待模型参数更新完毕，执行步骤S5。

S4、基于当前GLM模型获得心脏表面目标点p的预测值，具体如下：

$\hat{p}\left(k\right)=Q(k-1)w(k-1)$

用预测值代替测量值，即令然后执行步骤S5。

S5、更新GLM模型设计矩阵

利用当前时刻的心脏表面目标点p的测量值p(k)和4个辅助点测量值 h₁(k),h₂(k),h₃(k),h₄(k)更新模型设计矩阵，得到

Q(k)＝[p(k-3)p(k-2)p(k-1)p(k)h₁(k)h₂(k)h₃(k)1]

然后执行步骤S6。

S6、将当前时刻值k加1，即：k＝k+1，再返回步骤S2，进入下一时刻k+1的 处理流程。

图2是本发明实施例中心脏表面目标点和辅助点的位置示意图。

在本实施例中，测量系统采用基于立体内窥镜的视觉测量系统，对心脏表 面1个目标点和周围3个辅助点进行测量，3个辅助点选取了特征较明显的像素 点，相较目标点更容易测量。实际应用中，对于辅助点也可以采用完全相同的 预测方法(此时，原目标点和其它2个辅助点则变为辅助点)。这种处理可以确 保当辅助点在某一时刻无测量时，由其预测值替代，对目标点的GLM预测模型 进行有效更新。

图3是本发明实施例中心脏表面目标点测量失败时的运动预测结果示意图。

在本实施例中，目标点在第908帧～987帧，即k＝908～k＝987期间，因手术 器械的遮挡，无法获得测量值，而3个辅助点未受遮挡，仍可正常测量。图3 中显示了测量系统获得的目标点和辅助点的运动曲线，并在目标点运动曲线中 显示了其在第908帧～987帧受遮挡期间，利用本发明方法得到的预测值的曲线。 在本实施例中，根据图1中的处理流程：

在第908帧，测量系统无法提供目标点测量值时，经步骤S2判断后，执 行步骤S4，基于当前GLM模型进行预测，获得目标点预测值，具体方程为：

$\hat{p}\left(908\right)=Q\left(907\right)w\left(907\right)$

其中，Q(907)和w(907)是在第907帧正常更新后的模型设计矩阵和模型参数。 接下来，在执行步骤S5时，由于第908帧无测量值，故令利用 目标点的预测值和辅助点的测量值h₁(908),h₂(908),h₃(908),h₄(908)更新设计矩阵，得 到：

Q(908)＝[p(905)p(906)p(907)p(908)h₁(908)h₂(908)h₃(908)1]

而模型参数未更新，仍为前一帧的值，即：w(908)＝w(907)，其方差矩阵仍有： V(908)＝V(907)。

在第909帧，与第908帧一样，在步骤4，目标点的经GLM模型方程预 测得利用预测值(令)和辅助点测量值 h₁(909),h₂(909),h₃(909),h₄(909)对模型参数矩阵更新，得到Q(909)，而模型参数仍未更 新，有w(909)＝w(908)和V(909)＝V(908)。

直到，第987帧，测量系统恢复，重新获得目标点测量值p(987)，此时，在 步骤S3中，恢复对模型参数矩阵的正常更新。

尽管上面对本发明说明性的具体实施方式进行了描述，以便于本技术领域 的技术人员理解本发明，但应该清楚，本发明不限于具体实施方式的范围，对 本技术领域的普通技术人员来讲，只要各种变化在所附的权利要求限定和确定 的本发明的精神和范围内，这些变化是显而易见的，一切利用本发明构思的发 明创造均在保护之列。