基于有限元和机器学习的脊柱近端交界角的手术优化方法与流程

本发明涉及医学图像处理技术领域，尤其涉及基于有限元和机器学习的脊柱近端交界角的手术优化方法。

背景技术：

随着新设备设计的出现，用于治疗脊柱畸形的脊柱内固定装置不断发展，为外科手术提供了许多可供选择的矫形方法。这些改变引领了许多手术方案的革新，包括近端和远端融合节段的选择，种植物类型例及内固定装置材料的选择。在最新公布的专利中大多都是针对手术器械的改进，例如中国专利cn110916732a一种用于脊柱手术拉钩的保护套，cn110338897a一种骨科脊柱侧弯手术用可调弯折板。但在脊柱矫形这类复杂手术过程中每个组成部分的生物力学效应，以及这些变量如何通过生物力学效应而改变病人的预后情况，尚不完全清楚。在脊柱矫形领域，计算机辅助设计的最新进展迅速改变了骨科手术的面貌，随着有价值的模型应用，患者的临床结果得到了显著改善。对于后路脊柱融合术，通过植入物和连接杆纠正脊柱在3个维度上的异常曲度。在虚拟矫形手术中，脊柱近端交界角度(pja)会重新和手术过程中的椎体去旋转匹配，每个pja的改变都伴随着脊椎的移位和旋转，都可以通过有限元精确的模拟出来。假设术中pja选择不当，导致术后上肢体重在最上端固定椎体(upperinstrumentedvertebrae,uiv)表面应力分布不均匀，极可能导致近段交界性后凸(proximaljunctionalkyphosis,pjk)的发生，对该术后并发症采用目前广泛公认的定义：最上端固定椎体的下终板与其上位2个椎体的上终板所形成的cobb角(即pja)满足大于10°并且比术前测量值增大10°及以上)。“biomechanicalriskfactorsforproximaljunctionalkyphosis”(spine第39卷第8期2014年1月)一文中通过有限元生物力学研究得出保留更多的近端后段韧带，在最上端固定椎体使用横突钩和过渡杆，降低连接杆的杆曲度将有助于降低术后近端交界性后凸的发生。但该研究着力于从统计学的角度描述多种手术方案可能造成的差异，并未针对每一位患者做出个性化的手术矫形方案推荐，且在脊柱曲度的具体矫形数值方面也未给出的优化建议。类似的生物力学研究最大的局限性在于需要手动选择纳入变量的数值或类别范围，过程相当繁琐而且效率低下，还无法做到给定范围区间进行自动的数值优化选择。

技术实现要素：

本发明旨在克服现有研究方法的不足，为避免术后近段交界性后凸发生，结合有限元受力分析和机器学习高效模拟的优势，为患者提供一种自动化、个性化且精准的手术方案优化方法。

为达到上述目的，本发明是采用以下技术方案实现的：

基于有限元和机器学习的脊柱近端交界角的手术优化方法，包括以下步骤：

步骤s1：收集患者信息，所述信息包括患者的性别、年龄、体重、术前ct、术前的全脊柱x光以及术后的全脊柱x光；

步骤s2：对所述术前ct进行图像分割；

步骤s3：从患者术前的全脊柱x光以及术后的全脊柱x光上分别测量得到患者的pja角度α；

步骤s4：计算患者在直立情况下最上端固定椎体头端的椎间盘术前术后所受到的压力f；

步骤s5：通过软件hypermesh进行静力分析，得到由手术过程中患者pja角度变化所产生的患者在直立状态下椎间盘上表面的应力分布改变数据；

步骤s6：根据患者信息以及步骤s3中所得到的应力分布改变数据来建立机器学习模型g(·)；

步骤s7：通过libsvm软件包来训练所述机器学习模型g(·)；

步骤s8：通过最小化理想输出结果与机器学习模型g(·)所得到的预测结果之间的绝对差异，得出术后pja的最优解

优选的，在步骤s2中，对术前ct进行图像分割是指将患者术前ct导入mimics软件并对最上端固定椎体头端椎间盘的纤维环和髓核部分进行自动的分割。

优选的，在步骤s4中，所述压力f的计算方法包括以下步骤：

步骤s401：根据uiv的位置以及患者的体重计算出其对应的椎间盘所受到的重力g；

步骤s402：对重力g进行垂直于椎间盘上表面的力的分解，从而得到椎间盘受到的压力f。

进一步的，在步骤s402中，压力f的计算式为：

f＝g*cos(α)。

优选的，在步骤s5中，所述静力分析包括以下步骤：

步骤s501：对纤维环和髓核进行网格的划分，选择0.3mm的三角形网格；

步骤s502：分别为纤维环和髓核赋予不同的材料特征，其弹性模量分别为3.0mpa和1.0mpa,其泊松比分别为0.45和0.3；

步骤s503：限制椎间盘底面的在xyz方向上的平移和旋转自由度，再对椎间盘表面施加对应的集中力即压力f，得到椎间盘的所有单元的应力分布区域；

步骤s504：将应力分布区域划分为共8个区域，分别是右前髓核区、右前纤维环区，右后髓核区、右后纤维环区，左前髓核区、左前纤维环区，左后髓核区以及左后纤维环区；

步骤s505：通过将对应区域术前、术后的压力f平均值或最大值进行相减，得到由手术过程中pja角度变化导致患者在直立状态椎间盘上表面的应力分布改变数据，所述应力分布改变数据为16个对应区域的应力改变值。

进一步的，在步骤s6中，将患者的年龄、性别以及16个应力改变值共18个变量作为输入变量，以术后两年患者是否发生脊柱近端交界性后凸作为二元的结局变量y，建立机器学习模型g(·)；所述结局变量y为1表示呈阳性，为0表示呈阴性，结局变量y通过下式表示：

y＝g[δσ(gi,αi,e,v)¹，δσ(gi,αi,e,v)²,…,δσ(gi,αi，e,v)ⁿ，agei，genderi,w]

其中w代表模型参数，i表示其为第i个患者，δσ(·)代表的是由有限元受力分析计算得出的八个区域所对应的16个应力改变值，每一个应力改变值由患者的重力g、pja角度α、弹性模量e以及泊松比v决定；age和gender表示患者的年龄和性别。

进一步的，在步骤s7中，机器学习模型g(·)以下列函数作为损失函数来优化模型参数：

δσ(gi,αi,e,v)²,…,δσ(gi,αi，e,v)ⁿ，agei，genderi,w||l2norm

其中是模型的最优参数，yi是指该患者的真实结局变量。

进一步的，在步骤s8中，将理想输出结果设定为0，即理想输出结果设定为呈阴性，则术后pja的最优解可用下列表达式表示：δσ(gi,αi,e,v)²,…,δσ(gi,αi，e,v)ⁿ，agei，genderi,

优选的，在步骤s5中，通过软件hypermesh进行静力分析。

优选的，在步骤s7中，通过libsvm软件包来训练所述机器学习模型g(·)。

本发明具有以下有益效果：

本发明通过有限元得出了不同的脊柱近端交界角度与椎间盘应力之间的对应关系，并且分区域展示分析的同时，还利用机器学习模型合理地模拟了患者术后两年并发症发生的概率与年龄、性别以及术前术后应力变化值等变量之间的非线性关系，除可用于预测脊柱侧弯患者术后情况外，还可以根据患者的术前信息自动化地计算出手术过程中最优的近段交界角度值以便减少术后并发症的发生。

附图说明

图1为通过mimics半自动分割并通过hypermesh进行网格划分后的最上端固定椎体头端椎间盘的纤维环和髓核；

图2为最上端固定椎体头端椎间盘上表面所受压力的分解示意图；

图3为hypermesh软件中的压力施加以及边界条件设定的示意图；

图4为椎间盘受力的区域分割图。

具体实施方式

为了使本发明的目的、技术方案及优点更加清楚明白，以下结合附图，对本发明进行进一步详细说明。

基于有限元和机器学习的脊柱近端交界角的手术优化方法，包括以下步骤：

步骤s1：收集患者信息，所述信息包括患者的性别、年龄、体重、术前ct、术前的全脊柱x光以及术后的全脊柱x光。

步骤s2：对所述术前ct进行图像分割。

在步骤s2中，对术前ct进行图像分割是指将患者术前ct导入mimics软件并对最上端固定椎体头端椎间盘的纤维环和髓核部分进行自动的分割。如图1所示，图1表示最上端固定椎体头端的椎间盘的分割效果图，其中橙色区域代表髓核，深红色区域代表纤维环。

步骤s3：从患者术前的全脊柱x光以及术后的全脊柱x光上分别测量得到患者的pja角度α。

步骤s4：计算患者在直立情况下最上端固定椎体头端的椎间盘术前术后所受到的压力f。

在步骤s4中，所述压力f的计算方法包括以下步骤：

步骤s401：根据uiv的位置以及患者的体重计算出其对应的椎间盘所受到的重力g；步骤s402：对重力g进行垂直于椎间盘上表面的力的分解，从而得到椎间盘受到的压力f。在步骤s402中，压力f的计算式为：f＝g*cos(α)。如图2所示，x光上蓝色色块代表脊柱最上端固定椎体，红色色块代表最上端固定椎体头端一个椎体，绿色色块代表最上端固定椎体头端两个椎体。最下方白色粗线为最上端固定椎体的下终板平面，最上方白色粗线为其头端两个椎体的上终板平面，这两条粗线延长线所形成的夹角即为近端交界角(pja,α)。黄色线段表示重力(g)的分解，f为施加于椎间盘表面的压力，由此可以得出不同的α在直立状态下产生的压力值f。

步骤s5：进行静力分析，得到由手术过程中患者pja角度变化所产生的患者在直立状态下椎间盘上表面的应力分布改变数据。

在具体实施中，通过软件hypermesh进行静力分析。在步骤s5中，具体地，所述静力分析包括以下步骤：

步骤s501：对纤维环和髓核进行网格的划分，选择0.3mm的三角形网格；步骤s502：分别为纤维环和髓核赋予不同的材料特征，其弹性模量分别为3.0mpa和1.0mpa,其泊松比分别为0.45和0.3；步骤s503：限制椎间盘底面的在xyz方向上的平移和旋转自由度，再对椎间盘表面施加对应的集中力即压力f，得到椎间盘的所有单元的应力分布区域，如图3所示，图3为在hypermesh中对椎间盘进行压力f的施加并产生的受力分布图；步骤s504：将应力分布区域划分为共8个区域，分别是右前髓核区、右前纤维环区，右后髓核区、右后纤维环区，左前髓核区、左前纤维环区，左后髓核区以及左后纤维环区，如图4所示，其中f代表纤维环，n代表髓核，a代表前方，p代表后方，l代表左，r代表右，以便进行术前术后对应区域压力值变化的计算；步骤s505：通过将对应区域术前、术后的压力f平均值或最大值进行相减，得到由手术过程中pja角度变化导致患者在直立状态椎间盘上表面的应力分布改变数据，所述应力分布改变数据为16个对应区域的应力改变值。

步骤s6：根据患者信息以及步骤s3中所得到的应力分布改变数据来建立机器学习模型g(·)。

在步骤s6中，将患者的年龄、性别以及16个应力改变值共18个变量作为输入变量，以术后两年患者是否发生脊柱近端交界性后凸作为二元的结局变量y，建立机器学习模型g(·)；所述结局变量y为1表示呈阳性，为0表示呈阴性，结局变量y通过下式表示：

y＝g[δσ(gi,αi,e,v)¹，δσ(gi,αi,e,v)²,…,δσ(gi,αi，e,v)ⁿ，agei，genderi,w]

其中w代表模型参数，i表示其为第i个患者，δσ(·)代表的是由有限元受力分析计算得出的八个区域所对应的16个应力改变值，每一个应力改变值由患者的重力g、pja角度α、弹性模量e以及泊松比v决定；age和gender表示患者的年龄和性别。

步骤s7：训练所述机器学习模型g(·)。在具体实施中通过libsvm软件包来训练所述机器学习模型g(·)。

在步骤s7中，机器学习模型g(·)以下列函数作为损失函数来优化模型参数：

δσ(gi,αi,e,v)²,…,δσ(gi,αi，e,v)ⁿ，agei，genderi,w||l2norm

其中是模型的最优参数，yi是指该患者的真实结局变量。

步骤s8：通过最小化理想输出结果与机器学习模型g(·)所得到的预测结果之间的绝对差异，得出术后pja的最优解

在步骤s8中，将理想输出结果设定为0，即理想输出结果设定为呈阴性，则术后pja的最优解可用下列表达式表示：

δσ(gi,αi,e,v)²,…,δσ(gi,αi，e,v)ⁿ，agei，genderi,

当然，本发明还可有其它多种实施例，在不背离本发明精神及其实质的情况下，熟悉本领域的技术人员可根据本发明作出各种相应的改变和变形，但这些相应的改变和变形都应属于本发明所附的权利要求的保护范围。