一种基于正逆向多体动力学的个体下颌运动虚拟仿真方法和装置

1.本技术属于计算机仿真领域，特别涉及一种基于正逆向多体动力学的个体下颌运动虚拟仿真方法和装置。


背景技术：

2.虚拟仿真下颌运动在颞下颌关节疾病诊断、功能评价、牙科材料实验、颌面关节疾病治疗、颌面修复重建、人工关节以及仿生材料等方面均具有重要作用。颞下颌关节是颌面部的重要关节，其与周围肌肉共同决定下颌运动。下颌运动过程中关节、韧带、肌腱、颌骨、牙齿等的力与运动参数虽然可以从离体实验中获取，但与真实人体仍存在差异。由于功能解剖的复杂性及伦理所限，目前并不能基于在体实验虚拟重现完整精确的下颌运动。人体中执行下颌运动功能的下颌骨、咀嚼肌、关节、韧带等组织可以视作一个精密的多体系统，运用多体系统动力学方法研究下颌运动已经进入广大学者的视野。现有技术主要包括以下几类仿真模型：
3.(1)个体化的下颌运动虚拟显示模型仅作为下颌运动的动态显示，本质上并不包含力与运动的关系，无法作为动力学研究模型、仿生材料平台等应用。
4.(2)下颌运动与力的三维有限元分析是针对某一运动状态时下颌关节受力的分析，本质上是对被动受力的分析，没有主动力和力驱动的运动，并且，该方法运算慢、时间长，不能满足实时运动需求。
5.(3)肌肉正向驱动的下颌运动模型是使用骨骼肌的力学模型，设定包括粘弹性等力学参数，通过肌肉的力与长度的关系，以及肌肉骨骼的连接方式，驱动下颌骨的运动，该方案所得运动结果的准确性受限于肌肉物理模型的设定，物理模型无法真实模拟不同个体人的肌肉力学和运动。
6.(4)逆向运动轨迹控制的下颌运动模型，是以设定的运动轨迹为目标，迭代求解实现运动轨迹的控制系统动力学参数，主要用于咀嚼机器人控制。其求解结果往往与人体真实肌肉活动相去甚远，往往用于实验室材料分析等用途，难以应用于直接与人体疾病和康复有关的需求。


技术实现要素：

7.为解决现有技术中的问题，本发明提供了一种基于正逆向多体动力学的个体下颌运动虚拟仿真方法和装置，所述方法结合个性化正向动力学参数以及个性化逆向动力学参数，在虚拟仿真过程中考虑舌骨的位置变化，采用正逆双向仿真结合的方式对患者的下颌运动进行虚拟仿真，该虚拟仿真最大程度地预测患者术后的运动形态。
8.本技术的目的在于提供以下几个方面：
9.第一方面，本技术提供一种基于正逆向多体动力学的个体下颌运动虚拟仿真方法，所述方法包括：
10.获取个性化正向动力学参数；
11.通过逆向动力学方法获取个性化逆向动力学参数；
12.根据所述个性化正向动力学参数以及所述个性化逆向动力学参数采用正向动力学方法生成动力学仿真数据。
13.在一种可实现的方式中，所述获取个性化正向动力学参数包括：
14.采集表面肌电数值；
15.对所述表面肌电数值进行预处理生成正向动力学参数。
16.进一步地，对所述表面肌电数值进行预处理包括：
17.将采集到的表面肌电数据进行翻正处理；
18.翻正处理后的表面肌电进行滤波；
19.获取各时间滑窗内滤波后数据均值的最大值；
20.根据所述各时间滑窗内滤波后数据均值的最大值对表面肌数据进行归一化处理。
21.在一种可实现的方式中，所述获取个性化逆向动力学参数包括：
22.获取下颌运动学参数，所述下颌运动学参数包括舌骨运动学参数；
23.获取骨骼肌肉几何参数；
24.根据所述下颌运动学参数以及所述骨骼肌肉几何参数生成逆向动力学参数。
25.进一步地，所述获取下颌运动学参数包括：
26.采集下颌运动轨迹；
27.根据所述下颌运动轨迹计算质心运动轨迹r以及旋转矩阵r；
28.根据所述质心运动轨迹r以及旋转矩阵r生成下颌骨上任意点在实验室坐标系内的运动轨迹。
29.进一步地，所述获取骨骼肌肉几何参数包括：
30.获取口颌肌群几何参数；
31.获取骨骼几何参数。
32.可选地，所述获取骨骼几何参数包括：
33.获取cbct数据；
34.根据cbct数据重建颌面部骨骼解剖模型；
35.根据所述颌面部骨骼解剖模型生成骨骼几何参数。
36.可选地，所述获取口颌肌群几何参数包括：
37.建立咀嚼肌解剖模型；
38.建立骨骼肌解剖模型；
39.根据所述咀嚼肌解剖模型与所述骨骼肌解剖模型生成口颌肌群解剖模型；
40.根据所述口颌肌群解剖模型获取骨骼肌肉解剖数据。
41.进一步地，所述建立咀嚼肌解剖模型包括：
42.获取颌面部骨骼解剖模型；
43.根据所述颌面部骨骼解剖模型建立咀嚼肌解剖模型。
44.可选地，所述根据所述咀嚼肌解剖模型与所述骨骼肌解剖模型生成口颌肌群解剖模型包括：
45.确定所述咀嚼肌解剖模型中肌肉附着点；
46.确定所述骨骼肌解剖模型中肌肉截面轴线；
47.连接所述肌肉附着点与所述肌肉截面轴线。
48.在一种可实现的方式中，所述根据所述个性化正向动力学参数以及所述个性化逆向动力学参数生成动力学仿真数据包括：
49.建立个性化肌骨系统虚拟仿真模型；
50.获取所述个性化正向动力学参数以及所述个性化逆向动力学参数；
51.根据所述个性化正向动力学参数以及所述个性化逆向动力学参数利用所述个性化肌骨系统虚拟仿真模型生成动力学仿真数据。
52.第二方面，本技术还提供一种基于正逆向多体动力学的个体下颌运动虚拟仿真装置，所述装置包括：
53.个性化正向动力学参数，用于获取个性化正向动力学参数；
54.个性化逆向动力学参数获取模块，用于通过逆向动力学方法获取个性化逆向动力学参数；
55.动力学仿真数据生成模块，用于根据所述个性化正向动力学参数以及所述个性化逆向动力学参数采用正向动力学方法生成动力学仿真数据。
56.在一种可实现的方式中，所述个性化正向动力学参数获取模块具体用于：
57.采集表面肌电数值；
58.对所述表面肌电数值进行预处理生成正向动力学参数。
59.进一步地，对所述表面肌电数值进行预处理包括：
60.将采集到的表面肌电数据进行翻正处理；
61.翻正处理后的表面肌电进行滤波；
62.获取各时间滑窗内滤波后数据均值的最大值；
63.根据所述各时间滑窗内滤波后数据均值的最大值对表面肌数据进行归一化处理。
64.在一种可实现的方式中，所述个性化逆向动力学参数获取模块具体用于：
65.获取下颌运动学参数，所述下颌运动学参数包括舌骨运动学参数；
66.获取骨骼肌肉几何参数；
67.根据所述下颌运动学参数以及所述骨骼肌肉几何参数生成逆向动力学参数。
68.进一步地，所述获取下颌运动学参数包括：
69.采集下颌运动轨迹；
70.根据所述下颌运动轨迹计算质心运动轨迹r以及旋转矩阵r；
71.根据所述质心运动轨迹r以及旋转矩阵r生成下颌骨上任意点在实验室坐标系内的运动轨迹。
72.进一步地，所述获取骨骼肌肉几何参数包括：
73.获取口颌肌群几何参数；
74.获取骨骼几何参数。
75.可选地，所述获取骨骼几何参数包括：
76.获取cbct数据；
77.根据cbct数据重建颌面部骨骼解剖模型；
78.根据所述颌面部骨骼解剖模型生成骨骼几何参数。
79.可选地，所述获取口颌肌群几何参数包括：
80.建立咀嚼肌解剖模型；
81.建立骨骼肌解剖模型；
82.根据所述咀嚼肌解剖模型与所述骨骼肌解剖模型生成口颌肌群解剖模型；
83.根据所述口颌肌群解剖模型获取骨骼肌肉解剖数据。
84.进一步地，所述建立咀嚼肌解剖模型包括：
85.获取颌面部骨骼解剖模型；
86.根据所述颌面部骨骼解剖模型建立咀嚼肌解剖模型。
87.可选地，所述根据所述咀嚼肌解剖模型与所述骨骼肌解剖模型生成口颌肌群解剖模型包括：
88.确定所述咀嚼肌解剖模型中肌肉附着点；
89.确定所述骨骼肌解剖模型中肌肉截面轴线；
90.连接所述肌肉附着点与所述肌肉截面轴线。
91.在一种可实现的方式中，所述根据所述个性化正向动力学参数以及所述个性化逆向动力学参数生成动力学仿真数据包括：
92.建立个性化肌骨系统虚拟仿真模型；
93.获取所述个性化正向动力学参数以及所述个性化逆向动力学参数；
94.根据所述个性化正向动力学参数以及所述个性化逆向动力学参数利用所述个性化肌骨系统虚拟仿真模型生成动力学仿真数据。
95.第三方面，本技术还提供一种基于正逆向多体动力学的个体下颌运动虚拟仿真的程序，所用程序用于执行时实现上述第一方面所述基于正逆向多体动力学的个体下颌运动虚拟仿真方法的步骤。
96.第四方面，一种计算机可读存储介质，其上存储有计算机指令，该指令被处理器执行时实现上述第一方面所述基于正逆向多体动力学的个体下颌运动虚拟仿真方法的步骤。
97.第五方面，一种检测设备，所述检测设备包括：至少一个处理器；以及与所述至少一个处理器通信连接的存储器；其中，所述存储器存储有可被所述一个处理器执行的指令，所述指令被所述至少一个处理器执行，以使所述至少一个处理器执行上述第一方面所述基于正逆向多体动力学的个体下颌运动虚拟仿真方法。
98.与现有技术相比，本技术提供的方法以个性化逆向动力学参数为基础，结合直接采集于个体的咀嚼肌肌电活动以及下颌运动轨迹等个性化正向动力学参数，采用正逆向相结合的方法，对肌骨动力关系进行修正和优化，获得更接近于真实的个体化仿真。
附图说明
99.图1示出本技术提供一种基于正逆向多体动力学的个体下颌运动虚拟仿真方法的流程图。
具体实施方式
100.这里将详细地对示例性实施例进行说明，其示例表示在附图中。下面的描述涉及附图时，除非另有表示，不同附图中的相同数字表示相同或相似的要素。以下示例性实施例
中所描述的实施方式并不代表与本发明相一致的所有实施方式。相反，它们仅是与如所附权利要求书中所详述的、本发明的一些方面相一致方法的例子。
101.下面通过具体的实施例对本技术提供的基于正逆向多体动力学的个体下颌运动虚拟仿真方法及装置进行详细阐述。
102.首先，对本方案的使用场景作简要介绍。
103.传统的正向个体下颌运动虚拟仿真方法所获得的仿真数据仅仅颅骨模型为个性化的，而其运动轨迹则不是个性化的，即，使用统一的运动学参数嵌入于不同的颅骨模型中，然而，这种仿真模拟具有较大的失真，仅可作为示意演示，对于用于严谨的医疗设计的基础，例如，在下颌运动仿真基础上，设计手术和术后康复方案，其仿真效果显然是不足的。
104.基于多体系统动力学方法建立的肌肉—骨骼系统(简称肌骨系统)模型，为模拟颞下颌关节运动提供了数值仿真工具。近年来，随着计算机技术的进步，用于骨肌系统建模的软件平台被不断开发，目前已形成opensim、anybody、lifemod等肌肉—骨骼系统仿真软件，用于膝关节、脊柱等人体大型关节的运动仿真中。人体中执行下颌运动功能的下颌骨、咀嚼肌、关节、韧带等组织就可以视作一个精密的多体系统，运用多体动力学方法研究下颌运动已经进入广大学者的视野。
105.图1示出本技术提供一种基于正逆向多体动力学的个体下颌运动虚拟仿真方法的流程图，如图1所示，所述方法包括以下步骤s100至步骤s300：
106.步骤s100，获取个性化正向动力学参数。
107.在本实例中，所述获取个性化正向动力学参数包括步骤s101至步骤s102：
108.步骤s101，采集表面肌电数值。
109.在本实例中，采集表面肌电数值可以使用现有技术中任意一种可行的方案。由于面颌部的肌肉数量多，因此，可以结合多种采集手段，以采集尽量多的表面肌电数值。可以理解的是，如果采集方案可同时采集其它信息，在本步骤中也可同时采集，而不必加以阻断，以减少实际操作的复杂度，提高处理效率。
110.例如，可以同步采集表面肌电数值与本实例中步骤s211需要使用的下颌运动轨迹。具体地，可以使用电子面弓(如，winjaw,zebris medical gmbh,isny,germany)采集患者下颌运动，该电子面弓依据超声感应原理采集六自由度刚体运动信息。采集过程首先需要将上颌颌板与患者的上颌牙列通过硅橡胶(奥凯贝特，ak-y-1，奥凯迪伟，中国)固定，该上颌颌板设置有3个标志点；随后使用自凝树脂(luxatemp automix plus,dmg,hamburg,germany)将下颌颌叉固定在被试的下牙列唇面上，同时保证其不影响被试的下颌运动。用75％酒精棉球擦拭被试者咬肌和颞肌肌束最突出处的皮肤表面，粘贴配套的zebris双极电极片并连接肌电测量模块。
111.使用与电子面弓配套的winjaw+软件(v1.4.8；zebris medical gmbh)进行运动采集，记录上颌颌板与下颌颌叉初始位置。患者按照实验人员口令进行前伸、大张、左侧方、右侧方等下颌边缘运动，以及咀嚼等下颌功能运动。每个动作从牙尖交错位起始，最后动作结束返回牙尖交错位，每个动作重复3次，每次大约持续5秒钟。在这一过程中，zebris表面肌电仪同步记录双侧咬肌和颞肌肌电变化。将记录的上颌颌板标志点运动轨迹保存为xml文件，咀嚼肌表面肌电保存为csv文件。
112.此外，患者还需要采集最大咬合工况下的静态相应数值，亦保存为xml文件，咀嚼
肌表面肌电亦保存为csv文件。
113.步骤s102，对所述表面肌电数值进行预处理生成正向动力学参数。
114.在本实例中，对所述表面肌电数值进行预处理生成正向动力学参数可以使用现有技术中任意一种对表面肌电数值进行预处理的方法，实现对直接采集所得的表面肌电数值进行归一化处理，例其能够应用于后续步骤中，以下以一种实现方式为例加以说明，在本实例中，本步骤具体可以包括步骤s121至步骤s124：
115.步骤s121，将采集到的表面肌电数据进行翻正处理。
116.在matlab环境中，将颞肌与咬肌等面颌肌肉的表面肌电翻正，以利于后续数据处理，将电信号映射为肌肉激活程度。
117.步骤s122，翻正处理后的表面肌电进行滤波。
118.在本实例中，通过四阶butterworth低通滤波对步骤s121所得结果进行滤波处理，从而良好降噪，得到未归一化的解析肌电。
119.可以理解的是，还可以采用现有技术中其它滤波方案进行滤波处理。
120.步骤s123，获取各时间滑窗内滤波后数据均值的最大值。
121.在本实例中，选取合适长度的时间滑窗，获取滤波后的表面肌电在该长度的时间滑窗内的均值，即，该滑窗内的平均肌电，所述时间滑窗的长度可根据实际情况而具体设定，例如，可以根据信号长度而设定。
122.对于同一块肌肉，从所有滑窗内平均肌电中选取最大值即为最大自主收缩(maximal voluntary contraction，mvc)。可以理解的是，每块肌肉仅对应一个mvc值。
123.步骤s124，根据所述各时间滑窗内滤波后数据均值的最大值对表面肌电数据进行归一化处理。
124.在本实例中，以步骤s123获得的mvc值对步骤s122获得的滤波后表面肌电进行归一化，得到介于0至1之间的解析表面肌电时间序列ui(t)，1≤i≤n
mus
，n
mus
表示咬合肌群的肌束数目，以映射肌肉的激活程度。
[0125][0126]
其中，emg表示肌电图(electromyography)。
[0127]
在本步骤中，在进行归一化处理过程中所用解析参考肌电，可以通过以下方式获得：对在最大咬合工况下测得的同一肌肉静态表面肌电信号做如步骤s121至步骤s123同样处理。
[0128]
步骤s200，通过逆向动力学方法获取个性化逆向动力学参数。
[0129]
在本实例中，所述个性化逆向动力学参数为通过逆向动力学方法获取的一系列参数，这些参数可进一步应用于后续正向动力学模型中，从而获得具有个性化特点的动力学仿真数据，以提高仿真结果与真实情形的相似程度。
[0130]
具体地，本步骤可以包括以下步骤s201至步骤s203：
[0131]
步骤s201，获取下颌运动学参数，所述下颌运动学参数包括舌骨运动学参数。
[0132]
在本实例中，所述获取下颌运动学参数包括步骤s211至步骤s213：
[0133]
步骤s211，采集下颌运动轨迹。
[0134]
在本实例中，本步骤所需要采集的下颌运动轨迹可以步骤s101中同步采集，具体
方法可参见步骤s101，在此不再赘述。
[0135]
可以理解的是，下颌运动轨迹为随时间变化体坐标的位置。
[0136]
步骤s212，根据所述下颌运动轨迹计算质心运动轨迹r以及旋转矩阵r。
[0137]
将ct的dicom文件导入数字化软件，在所述数字化软件中分割上颌骨运动轨迹与下颌骨运动轨迹，并且，依据分割后的内容分别重建得到上颌骨三维模型(stl文件)与下颌骨三维模型(stl文件)。
[0138]
在上颌骨三维模型中选取牙列上的牙尖、边缘嵴等解剖标志，在上述标志点配准的基础上，实现上颌骨三维模型(stl文件)与口内扫描stl文件的全局配准。其中，所述口内扫描stl文件可用口内扫描仪(trios,intraoral scanner,3shape,denmark)扫描得到，并且，所述口内扫描stl文件记录有上颌颌叉—上颌牙列相对位置关系。
[0139]
通过同样的方式，将口内扫描stl文件与标准上颌颌板stl文件进行配准，所述标准上颌板stl文件预存于系统中，从而使上颌骨三维模型与标准上颌颌板三维模型实现配准。通过本技术提供的方案，将原本无共同部分的两个模型实现配准。
[0140]
进一步在，在ct与标准上颌颌板完成配准后得到一个带有轨迹标志点的颌骨模型。
[0141]
在本实例中，所有配准均可选用现有技术中任意一种配准对齐的方法，例如，可以采用多个标志点对齐方法、三维形态对齐方法或者曲面对齐配准方法等。
[0142]
下面示例性说明配准基本转换原理：
[0143]
假设存在两个相关联的点集dj和mj，其中，j的取值范围为1到n，n表示静态三维模型中点的数量，那么这两个点集的关系如下式(2)表示：
[0144]dj
＝rmj+t+vjꢀꢀꢀꢀꢀꢀꢀ
式(2)
[0145]
其中，r是一个标准的3
×
3旋转矩阵，t是一个三维平移向量，vj是干扰向量，求解最优变换[r^,t^]映射集合mj在dj通常需要最小化一个最小二乘误差准则，所述最小二乘误差准则可如下式(3)表示:
[0146][0147]
在点集中的点不存在错误对应的情况下，这个最小二乘方程的最优解就是成立的。
[0148]
进一步地，假设下颌处在闭口位时，上颌骨与下颌骨的局部坐标系重合。导入步骤s101获得的xml轨迹文件，根据奇异值分解方法，求解下颌运动过程中，下颌骨质心运动轨迹r与刚体旋转矩阵r’，并将旋转矩阵转换为基于euler角定义的下颌姿态角度α，β，γ。
[0149]
步骤s213，根据所述质心运动轨迹r以及旋转矩阵r生成下颌骨上任意点在实验室坐标系内的运动轨迹。
[0150]
本步骤可采用现有技术中任意一种方案，根据下颌骨质心运动轨迹r与刚体旋转矩阵r，通过刚体运动学变换得到下颌骨上任意点在实验室坐标系内的运动轨迹。
[0151]
步骤s202，获取骨骼肌肉几何参数。
[0152]
进一步地，所述获取骨骼肌肉几何参数包括步骤s221和s222：
[0153]
步骤s221，获取骨骼几何参数，包括步骤s2211至步骤s2213。
[0154]
步骤s2211，获取cbct数据。
[0155]
在本实例中，可对被试者进行cbct拍摄(newtom vg，newtom，italy)，拍摄范围下自舌骨，上至关节窝。
[0156]
拍摄参数可根据需要设定，例如，视野15cm*15cm，层厚0.3cm，扫描过程中被试保持牙尖交错位咬合。
[0157]
步骤s2212，根据cbct数据重建颌面部骨骼解剖模型。
[0158]
在本实例中，本步骤可以使用现有技术中任意一种根据cbct数据重建颌面部骨骼解剖模型的方法，例如，可以使用以下方法：
[0159]
将步骤s2211获得的cbct-dicom图像导入mimics17.0(mimics，v17.0.0.436，materialise，belgium)软件中，使用阈值分割技术，分割出1250hu～4095hu的ct图像，擦除寰枢关节以下的颈椎图像，分别根据解剖标志在剩余ct图像上分离出上颌骨、下颌骨、舌骨的影像，并按照软件中high quality选项进行3d重建，重建所得骨骼解剖模型后分别以stl格式文件输出。
[0160]
进一步，将得到的骨骼解剖stl模型导入到insides多体动力学软件。在所述软件中，定义颅骨、下颌骨与舌骨三个刚体，对于每个刚体来讲，质量m与转动惯量矩阵i均可根据ammr模型库中的标准口颌模型线性缩放得到，缩放比例scale根据标准口颌模型与实测左、右两侧颞下颌关节中心距离相比得到。
[0161]
步骤s2213，根据所述颌面部骨骼解剖模型生成骨骼几何参数。
[0162]
在本实例中，本步骤可采用现有技术任意一种方式进行。
[0163]
步骤s222，获取口颌肌群几何参数。
[0164]
可选地，本步骤可以包括步骤s2221至步骤s2224：
[0165]
步骤s2221，建立咀嚼肌解剖模型。
[0166]
可选地，所述建立咀嚼肌解剖模型包括：
[0167]
步骤s22211，获取颌面部骨骼解剖模型；
[0168]
在本实例中，可以使用前述步骤已经建立的颌面部骨骼解剖模型，具体地，在步骤s2212中所建立的颌面部骨骼解剖模型。
[0169]
步骤s22212，根据所述颌面部骨骼解剖模型建立咀嚼肌解剖模型。
[0170]
在本实例中，在所述颌面部骨骼解剖模型中确定各块肌肉的起止附着点，并在所述颌面部骨骼解剖模型中插入咀嚼肌。
[0171]
具体地，所插入的咀嚼肌包括：颞肌的前、中、后束，翼外肌上层和翼外肌下层、咬肌的深层和浅层、翼内肌、二腹肌前束、下颌舌骨肌的前束和后束、颏舌骨肌。
[0172]
进一步地，各块肌肉的起止附着点可按现有技术关于肌肉解剖图谱进行定位，例如，按照gray和sicher主编的解剖图谱进行精确定位，并由具有解剖学背景知识的专业人员在颌面部骨骼几何stl模型上点出。
[0173]
步骤s2222，建立骨骼肌解剖模型。
[0174]
本实例中，本步骤可采用现有技术任意一种方法建立骨骼肌解剖模型，例如，可以使用guo等(multibody syst dyn(2020)49:315
–
336)提出的柔性索有限单元建模骨骼肌，在该方案中，每一有限单元包含2个节点(1，2)，每个节点的运动基于部分拉格朗日—欧拉描述，从而得到的单元广义坐标取为其中，ra是节点在实验室坐标系的空间坐标，pa是凝聚肌肉质量的物质坐标，a为节点数，取值为1或者2。
[0175]
步骤s223，根据所述咀嚼肌解剖模型与所述骨骼肌解剖模型生成口颌肌群解剖模型。
[0176]
在本实例中，本步骤可以包括：
[0177]
步骤s2231，确定所述咀嚼肌解剖模型中肌肉附着点。
[0178]
在本实例中，本步骤可以由前述步骤s22212一并确定。
[0179]
步骤s2232，确定所述骨骼肌解剖模型中肌肉截面轴线。
[0180]
在本步骤中，连接步骤s2231确定的附着点，其中，索单元的路径定义为对应肌肉截面轴线。
[0181]
可选地，肌肉内力可以建模为hill模型，其中，每束肌肉的最大主动肌力由标准模型取值根据scale
1.5
缩放得到，其中，所述标准模型可以为预存模型，所述包含肌肉主动收缩产生的主动肌力f
ce
，以及肌纤维被动伸长产生的被动肌力f
pe
。
[0182]
在本步骤中，所述颞下颌关节的几何模型建模为髁突与关节面间柔性接触，所述颞下颌关节接触模型设为flores建立的软物质接触函数，其中，法向接触力并忽略关节面之间的切向摩擦力。
[0183]
其中，δ表示关节面间的嵌入深度，表示关节面间的嵌入速度，k表示法向接触刚度，n表示接触力的非线性指数，χ由关节面接触的恢复系数确定。
[0184]
在本实例中，可将所述髁突接触几何简化为椭球，特别地，所述椭球几何参数是对真实髁突几何做最小二乘拟合得到；颞下颌关节的接触面由真实关节面拟合成样条扫掠面得到。
[0185]
此外，在上、下颌中切牙与第一磨牙间建立刚性hertz接触模型，避免咬合过程中出现上下牙列嵌入。
[0186]
步骤s2233，连接所述肌肉附着点与所述肌肉截面轴线。
[0187]
在本实例中，本步骤可以手动实现，也可以利用计算机自动实现。
[0188]
步骤s224，根据所述口颌肌群解剖模型获取骨骼肌肉解剖数据。
[0189]
在本实例中，所述骨骼肌肉解剖数据包括骨骼上多个特征点所在位置以及各块肌肉起始点和终止点位置。
[0190]
在本实例中，本步骤可以采用现有技术中任意一种方式由所述口颌肌群解剖模型获取相应骨骼肌肉解剖数据。
[0191]
步骤s203，根据所述下颌运动学参数以及所述骨骼肌肉几何参数生成逆向动力学参数。
[0192]
在本实例中，本步骤可以具体为获取下颌骨质心的运动轨迹与姿态角度的时间序列。
[0193]
在本实例中，所述二腹肌等降颌肌群的肌纤维长度由舌骨与下颌骨时变位置共同确定，特别地，本实例中考虑舌骨牵拉引起的降颌肌群长度变化，该变化不可忽略。
[0194]
为获取个性化舌骨运动轨迹，本实例对测得的闭口位、最大张口位舌骨相对下颌骨位置线进行性插值，从而获得下颌骨姿态角度对应的舌骨位置。
[0195]
在已知上述骨骼运动的条件下，逆向求解降颌肌群在下颌骨与舌骨附着点之间的距离，作为进一步正向求解的目标肌纤维长度λi，1≤i≤n
mus
，其中，n
mus
表示各模型中的肌
纤维个数。
[0196]
步骤s300，根据所述个性化正向动力学参数以及所述个性化逆向动力学参数采用正向动力学方法生成动力学仿真数据。
[0197]
在本实例中，考虑在大张口动作中头颈部在矢状面后倾θ
skull
。具体地，可以假设在最大张口时，头颈部后伸15
°
，可以理解的是，该角度还可设置为其它角度，使仿真结果更真实情况。
[0198]
本实例中，所述根据所述个性化正向动力学参数以及所述逆向动力学参数生成动力学仿真数据包括步骤301至步骤s303：
[0199]
步骤s301，建立个性化肌骨系统虚拟仿真模型。
[0200]
在本实例中，所述个性化肌骨系统虚拟仿真模型可以使用现有技术中任意一种可用的个性化肌骨系统虚拟仿真模型，可以理解的是，相关参数可以进行优化。
[0201]
优选地，所述个性化肌骨系统的虚拟仿真模型的控制方程可写为如下式(4)和式(5)所示的微分—代数方程，基于bdf方法进行数值积分：
[0202][0203]
c(q，t)＝0
ꢀꢀꢀ
式(5)
[0204]
步骤s302，获取所述个性化正向动力学参数以及所述个性化逆向动力学参数。
[0205]
在本实例中，设定肌纤维仿真长度li为降颌肌群的正向求解结果，以其与步骤s203中由逆向仿真得到的肌纤维目标长度λi的差距e
λ
(t)＝l
i-λi为控制量，建立确定降颌肌激活度模型的控制方程可以为如下式(6)所示的比例—微分控制器：
[0206][0207]
其中，表示调节控制器比例系数，表示微分系数，所述二者使仿真模型达到最大张口度时，降颌肌群的肌肉激活度a
λ
(t)最大。
[0208]
进一步地，任意一条肌束的解析肌电信号可通由下式(7)表示的一阶肌肉激活动力学方程确定：
[0209][0210]
其中，表示升颌肌群的肌肉激活度，可由解析表面肌电结果确定。
[0211]
τa表示包含归一化解析肌电与肌肉激活度的非线性函数，特别地，可由下式(8)表示：
[0212][0213]
其中，取值范围为0～1。
[0214]
步骤s303，根据所述个性化正向动力学参数以及所述个性化逆向动力学参数利用所述个性化肌骨系统虚拟仿真模型生成动力学仿真数据。
[0215]
在本实例中，将口颌肌群的激活度a
λ
(t)代入到式(6)表示的柔性肌肉有限单元模型中求解肌肉内力，从而获得正向驱动下颌骨运动数据，实现患者个性化肌骨系统虚拟仿
真模型的正逆向耦合动力学仿真。
[0216]
在本实例中，在获得动力学仿真数据后，可以进一步利用这些数据在坐标系中动态展现。具体地，转换坐标系为原始ct坐标系，选择下颌模型上的下切牙点、双侧髁突点，输出三个点在原始ct坐标系中随时间的轨迹坐标列表，以stl或obj格式输出虚拟模型上颌骨、下颌骨、舌骨的三维物体。
[0217]
在可以动态显示三维物体的软件中，通过数据接口，输入上述虚拟模型三维物体，确定下颌骨上的下切牙点、双侧髁突点，通过数据接口，输入上述下颌模型上三个点的轨迹坐标。
[0218]
对于每一个时间点，将下颌骨模型上确定的三点与轨迹坐标采用最小距离法对齐，获得下颌骨的位置，将所有时间点下颌骨位置随时间变化，获得下颌运动轨迹的显示，虚拟模型可以输出附着于颌骨上的多组咀嚼肌肉随着下颌骨的运动而激活的程度，可以通过颜色深浅反映激活程度，也可通过时间轴窗口实时查看激活数值。
[0219]
本技术提供的方案以肌骨动力物理参数为基础，结合可采集的下颌运动轨迹和表层咀嚼肌肌电，逆向求解降颌肌群长度等深层肌肉的物理参数，采用正逆向耦合动力学方法进行仿真，利用正逆向结合的肌肉激活度控制器对肌骨动力关系进行修正和优化，获得更接近个体真实下颌运动轨迹和相关咀嚼肌的动力学特征的多体动力学仿真效果，甚至是个体化仿真，并对舌骨运动位置的虚拟仿真。
[0220]
具体地，本技术提供的方案依据下颌运动捕捉、表面肌电采集、cbct测量等方法，创建患者个性化肌骨系统生物动力学模型，解决现有虚拟仿真模型无法结合受试者真实测试数据的问题。基于柔性多体动力学方法，实现口颌骨骼肌质量沿肌纤维分布，骨骼肌模型包含主、被动肌力模型，解决口颌运动过程中，降颌肌群、翼内肌等深层肌群肌肉激活情况的计算问题。
[0221]
以上结合具体实施方式和范例性实例对本技术进行了详细说明，不过这些说明并不能理解为对本技术的限制。本领域技术人员理解，在不偏离本技术精神和范围的情况下，可以对本技术技术方案及其实施方式进行多种等价替换、修饰或改进，这些均落入本技术的范围内。本技术的保护范围以所附权利要求为准。