在体器官组织生物力学参数估计方法

1.本发明涉及一种在体器官组织生物力学参数估计方法。


背景技术：

2.人体器官的生物力学模型在虚拟手术仿真、软组织运动计算等方面有重要意义，然而人体器官不仅有着复杂的内部结构，还处于多器官环绕的复杂内部环境中，使得内脏的生物力学模型的估计存在着较大不确定性。
3.随着现代医学和计算机技术的发展，术中导航技术成为在医学领域的重要领域之一，其中对内脏器官形变的精准模拟是术中导航技术成功的关键技术之一，对器官形变的模拟依赖于个性化的生物力学模型，包括器官复杂的边界条件和个性化的弹性生物力学参性，生物力学模型的准确构建与否关系到术中器官形变跟踪精度。
4.目前国内外已经有许多人提出了生物力学建模的方法。在生物力学参数估计方面，现有的超声弹性成像和核磁共振弹性成像技术都可以直接用于测量人体器官的相对弹性分布，但是目前这些技术的精度不能保证，主要用于临床中人体组织的定性分析，例如评估肝脏组织纤维化程度。逆有限元方法是近年来计算生物力学中人体组织建模的主要探索方向，这些方法通过测量一定应力下的物体形变，然后基于应力-应变关系进行反向求解得到材料参数，hajhashemkhani等为了估计在部分边界条件已知时的材料力学参数，基于采集到的部分位移信息和相应的预测信息，基于高斯-牛顿算法完成对目标的逆有限元分析，但是该方法中考虑的边界条件仍相对简单。近年来，部分学者提出结合贝叶斯推断来减少由不确定性因素造成的估计误差，peterlik等使用了降阶无迹卡尔曼滤波器(roukf)对目标的生物力学参数进行估计，以提升模拟过程中的鲁棒性，nikolaev等则将roukf应用到了边界条件突变情况下生物力学参数的估计过程，以此证明了roukf对边界条件变化情况也有很好的适应效果，但是这些方法需要一个精确的系统模型，否则系统模型本身的不确定性会导致这些方法结果发生偏离甚至不收敛。


技术实现要素：

5.有鉴于此，有必要提供一种在体器官组织生物力学参数估计方法，能够通过其内部的自适应阶段来减小生物力学参数估计过程中的不确定性误差。
6.本发明提供一种在体器官组织生物力学参数估计方法，该方法包括如下步骤：a.进行在体器官组织几何建模，得到在体器官组织的四面体网格结构；b.根据得到的在体器官组织的四面体网格结构，进行在体器官组织生物力学模型建模；c.根据所述在体器官组织生物力学模型，进行在体器官组织边界条件估计；d.使用强跟踪卡尔曼滤波器，在当前边界条件估计结果下对在体器官组织的生物力学参数进行估计。
7.具体地，所述的步骤a具体包括：
8.首先将在体器官组织几何模型由4d mri数据进行重建获得在体器官组织重建表面，再将在体器官组织重建表面进行网格划分得到在体器官组织的四面体网格结构。
9.具体地，所述的步骤b具体包括：
10.将在体器官组织边界条件等效为一组在体器官组织周围解剖组织碰撞接触而施加到在体器官组织表面的等效外力λ，再使用常用弹性参数杨氏模量和泊松率作为表示在体器官组织弹性材料属性的参数；所述在体器官组织的状态平衡方程为：
11.f(q
t
)+f+j
t
λ＝0，
12.其中，f(q
t
)表示外部力，f表示内部力，λ和j分别是要求解的边界等效外部力和外力约束的雅克比矩阵。
13.具体地，所述的步骤c具体包括：
14.首先对不同形变状态下的在体器官组织模型进行非刚性配准，获取两个形变状态下在体器官组织几何模型表面点的对应关系；
15.然后计算在体器官组织在没有边界状态下的自由形变，此时在体器官组织位移信息在对应点上进行计算：
16.δ＝dot(q
s-q
t
，nc)，
17.qs和q
t
分别是无边界条件下和有边界条件下在体器官组织表面对应点的位移，nc是指向在体器官组织各点形变方量的单位向量；
18.最后通过求解约束问题对边界条件λ进行估计：
19.wλ+δ＝0，
20.其中，w＝jcj
t
为常见的delassus算子，c表示在体器官组织弹性属性的顺应性矩阵；最终基于高斯-赛德尔算法求解得到估计的边界条件λ。
21.具体地，所述的步骤d具体包括：
22.使用强跟踪卡尔曼滤波器在当前边界条件估计结果下，分三个阶段对在体器官组织的生物力学参数进行估计；所述三个阶段包括：预测阶段、自适应阶段以及更新阶段。
23.具体地，所述的预测阶段包括：
24.使用高斯分布表示包含在体器官组织生物力学参数和在体器官组织初始形状信息的状态向量基于初始状态预测当前帧在体器官组织形变后的状态以及相应的协方差矩阵
[0025][0026][0027]
其中，是初始状态高斯分布组中的第i(i＝1，2，...，2n+1)个sigma组，和分别表示均值权重和协方差权重，qk表示高斯的噪声。
[0028]
具体地，所述的自适应阶段包括：为了尽可能的减少系统模型中不确定性误差的影响，自适应阶段结合在体器官组织历史形变的观测信息yk对协方差矩阵不断更新，直到模型的不确定性误差(即马氏距离)少于特定的阈值，的更新式子为：
[0029][0030]
其中，是观测数据的创新向量，m和分别表示历史窗口大小及其对应权重，rk表是观测数据中的高斯白噪声。
[0031]
具体地，所述的更新阶段包括：结合观测到的当前帧在体器官组织部分表面位移信息，利用计算的卡尔曼增益kk对预测的在体器官组织状态结果进行修改：
[0032][0033][0034]
其中是观测数据的协方差矩阵。
[0035]
本技术的有益效果包括：
[0036]
其一，本技术所处理的边界条件更复杂，通过使用等效外力代表边界条件，更完整地反应在体器官组织周围解剖组织对在体器官组织的影响，同时又不会使要处理的问题过于复杂，更全面地反应了在体器官组织表面所受到的约束情况，同时本技术对其它场景也更具泛化性。
[0037]
其二，本技术所使用的强跟踪卡尔曼滤波器能够更好地处理生物力学参数估计过程中的不确定性误差，强跟踪卡尔曼滤波器的自适应阶段通过结合历史信息，将系统中的不确定性误差减小到规定阈值内，使估计的生物力学参数更加精确。
附图说明
[0038]
图1为本发明在体器官组织生物力学参数估计方法的流程图。
具体实施方式
[0039]
下面结合附图及具体实施例对本发明作进一步详细的说明。
[0040]
参阅图1所示，是本发明在体器官组织生物力学参数估计方法较佳实施例的作业流程图。
[0041]
步骤s1，进行在体器官组织几何建模，得到在体器官组织的四面体网格结构。具体而言：
[0042]
在本实施例中，首先将在体器官组织几何模型由4d mri数据进行重建获得在体器官组织重建表面，再将在体器官组织重建表面进行网格划分得到在体器官组织的四面体网格结构。
[0043]
步骤s2，根据得到的在体器官组织的四面体网格结构，进行在体器官组织生物力学模型建模。具体而言：
[0044]
首先将在体器官组织边界条件等效为一组在体器官组织周围解剖组织碰撞接触
而施加到在体器官组织表面的等效外力λ，再使用常用弹性参数杨氏模量和泊松率作为表示在体器官组织弹性材料属性的参数。所述在体器官组织的状态平衡方程为：
[0045]
f(q
t
)+f+j
t
λ＝0，
[0046]
其中，f(q
t
)表示外部力，f表示内部力，λ和j分别是要求解的边界等效外部力和外力约束的雅克比矩阵。
[0047]
步骤s3，根据所述在体器官组织生物力学模型，进行在体器官组织边界条件估计。具体而言：
[0048]
首先对不同形变状态下的在体器官组织模型进行非刚性配准，获取两个形变状态下在体器官组织几何模型表面点的对应关系。
[0049]
然后计算在体器官组织在没有边界状态下的自由形变，此时在体器官组织位移信息在对应点上进行计算：
[0050]
δ＝dot(q
s-q
t
，nc)，
[0051]qs
和q
t
分别是无边界条件下和有边界条件下在体器官组织表面对应点的位移，nc是指向在体器官组织各点形变方量的单位向量。
[0052]
最后通过求解约束问题对边界条件λ进行估计：
[0053]
wλ+δ＝0，
[0054]
其中，w＝jcj
t
为常见的delassus算子，c表示在体器官组织弹性属性的顺应性矩阵。最终基于高斯-赛德尔算法求解得到估计的边界条件λ。
[0055]
步骤s4，使用强跟踪卡尔曼滤波器，在当前边界条件估计结果下对在体器官组织的生物力学参数进行估计。具体而言：
[0056]
使用强跟踪卡尔曼滤波器在当前边界条件估计结果下对在体器官组织的生物力学参数进行估计，主要分三个阶段：
[0057]
首先是预测阶段：使用高斯分布表示包含在体器官组织生物力学参数和在体器官组织初始形状信息的状态向量基于初始状态预测当前帧在体器官组织形变后的状态以及相应的协方差矩阵
[0058][0059][0060]
其中，是初始状态高斯分布组中的第i(i＝1，2，...，2n+1)个sigma组。和分别表示均值权重和协方差权重，qk表示高斯的噪声。
[0061]
然后是自适应阶段：为了尽可能的减少系统模型中不确定性误差的影响，自适应阶段结合在体器官组织历史形变的观测信息yk对协方差矩阵不断更新，直到模型的不确定性误差(即马氏距离)少于特定的阈值，的更新式子为：
[0062][0063]
其中，是观测数据的创新向量，m和分别表示历史窗口大小及其对应权重，rk表是观测数据中的高斯白噪声。
[0064]
最后是更新阶段：结合观测到的当前帧在体器官组织部分表面位移信息，利用计算的卡尔曼增益kk对预测的在体器官组织状态结果进行修改：
[0065][0066][0067]
其中是观测数据的协方差矩阵。
[0068]
本技术考虑了复杂边界条件对在体器官组织形变的影响，并采用改进的强跟踪卡尔曼滤波器，通过其内部的自适应阶段来减小生物力学参数估计过程中的不确定性误差。
[0069]
虽然本发明参照当前的较佳实施方式进行了描述，但本领域的技术人员应能理解，上述较佳实施方式仅用来说明本发明，并非用来限定本发明的保护范围，任何在本发明的精神和原则范围之内，所做的任何修饰、等效替换、改进等，均应包含在本发明的权利保护范围之内。