基于深度学习的髋关节应力分布算法研究
刘远平，宋昱锴，张小燕，刘贤强

【摘 要】针对髋关节软骨的应力分布算法研究问题，设计了一个基于深度学习模型来代替有限元分析。该深度学习模型分为无监督学习模块和有监督学习模块，首先使用无监督学习模块对髋关节的软骨和股骨进行形状编码；之后实现对应力分布数据的编码与解码，使得应力数据能够与神经网络相结合；然后通过监督学习，利用编码好的应力数据进行监督，使神经网络学习得到一个从髋关节软骨和股骨的形状码到应力分布的应力码的映射关系；最终得到一个拟合的深度学习模型。此模型能够在一定程度上模拟有限元分析方法，但是由于其平均绝对误差和归一化平均绝对误差比较大，所以还不能完全替代有限元分析方法。在此基础上，进一步探索了新模型在特征利用上的局限，并提出了改进的方向。
【关键词】髋关节软骨；深度学习；应力分布算法；FEA 替代算法

引用格式  刘远平, 宋昱锴, 张小燕, 刘贤强.基于深度学习的髋关节应力分布算法研究. 智能科学与技术学报[J], 2019, 1(3):  260-268

Research on hip joint stress distribution algorithms
based on deep learning
LIU Yuanping, SONG Yukai, ZHANG Xiaoyan, LIU Xianqiang

Abstract Aiming at the problem of the stress distribution algorithm of hip cartilage, a deep learning model to replace the finite element analysis (FEA) was proposed. This deep learning model was divided into unsupervised learning module and supervised learning module. Firstly, an unsupervised learning module was adopted to encode the shape of hip cartilage and femur. Then the coding and decoding of stress distribution implement was implemented so that stress data can be combined with the neural network. Next a supervised learning module supervised by the stress data was used, and the model uses neural networks to learn a mapping relationship from the shape code of the hip cartilage and femur to the stress code of the stress distribution. Finally, a fitted deep learning model was obtained. This deep learning model can simulate the FEA method to a certain extent. But the mean absolute error and the normalized mean absolute error are still larger than that of the FEA method, so the FEA method cannot be completely replaced by our deep learning model.
Meanwhile, the limitations of the deep learning model in the use of input features were studied, and a direction to improve the performance of the model was proposed.
Key words hip cartilage, deep learning, stress distribution algorithm, FEA surrogate algorithm

Citation LIU Yuanping.Research on hip joint stress distribution algorithms based on deep learning. Chinese Journal of Intelligent Science and Technology[J], 2019, 1(3):  260-268

1 引言

应力分析在工业领域有广泛的应用，包括航空航天、汽车、建筑以及生物医学工程空间等^{[1,2,3,4,5,6,7,8,9]}。但是当前应力分析的方法主要是以有限元分析（finite element analysis，FEA）方法为代表的较为传统的分析方法。虽然这些分析方法能取得较理想的计算结果，但是随着模型和应用的复杂性发生变化，利用这些方法模拟计算并得到结果可能需要几分钟^[10]到几天^[11]的时间。在医学领域，当需要快速诊断或治疗时间敏感时，运用这些方法变得不切实际。髋关节的应力分析同样遇到了这个问题。为了解决这个问题，笔者使用深度学习的方法对其进行模拟计算，以期得到一个保证精度的深度学习模型，替代已有的有限元分析方法，并将模拟计算的时间压缩到可接受的范围内。

对于髋关节的应力分析，有限元分析方法是许多课题的选择^{[12,13,14,15]}，其主要特点是：能够将连续体离散为有限单元进行数值求解，具有极大的灵活性和通用性。而使用深度学习与神经网络进行分析的项目少之又少。可以看到的是，参考文献[16]使用神经网络的方法代替了有限元分析方法，他们进行的是人胸主动脉的应力分析，这与髋关节的应力分析很相似。首先对人胸主动脉进行编码，使其能够输入神经网络中进行模拟计算，再通过神经网络的无监督学习方法和有监督学习方法进行学习，最终得到深度学习模型的参数，使得人胸主动脉的应力分析时间从半小时缩减为1 s^[16]，这给予了笔者极大的启发。笔者也希望通过深度学习的方法，对髋关节的应力进行分析。

本文主要研究并实现基于深度学习模型的髋关节软骨应力分布计算。本文实现的深度学习模型分为两部分：第一部分使用无监督模块对髋关节软骨和股骨进行形状编码，并对应力数据进行编码和解码，使得现有的应力数据能够编码为神经网络的输出形式，也能将神经网络的输出解码为应力数据；第二部分进行有监督学习，利用神经网络学习得到一个从髋关节软骨和股骨的形状码到应力分布的应力码的映射关系。

完成上述工作之后，可以得到一个髋关节软骨和股骨到软骨应力分布的深度学习模型。此模型能够在一定程度上模拟 FEA 对髋关节软骨应力分布的分析效果，但是却不能完全替换FEA。最后，本文提出了一个改进深度学习模型的方向，即需要重新评估有限元分析的给力过程，通过深度学习网络对其进行描述。

2 实验数据与分析

本文实验数据中的髋关节软骨和股骨均为未受力时的状态，还需要说明的一点是：实验数据中的应力数据是通过摆动股骨形状到具体的位置，并给定一个适当的力，再通过有限元分析方法进行分析得到的。本节将对实验数据进行相应的分析与预处理。

模型的输入数据共有794个样本，每个样本可以分成4个部分，即股骨（femur）的形状（如图1黄色部分）、髋关节软骨（cartilage）的形状（如图1紫色部分）、股骨对软骨的作用力（force）的大小以及软骨的受力情况（stress，如图2 所示）。股骨和软骨的描述方式都是点集，股骨由37 630个点组成，每个点包含x轴、y轴和z轴3个坐标的值，因此股骨数据的维数为37 630×3=112 890。软骨由2 645个点组成，每个点也包含x轴、y轴和z轴3个坐标的值，共有维数2 645×3=7 935。作用力的大小可使用三维空间中的一个向量来表示，包含x轴的分量F_x、y轴的分量F_y和z轴的分量F_z这3个数据。软骨的受力情况是指每个软骨上的点都有一个力（是一个数值，即应力），因此软骨的受力情况共有2 645个数据。4个部分加起来可以得到每个样本数据共112 890+7 935+3+2 645=123 473维。

图1

新窗口打开| 下载原图ZIP| 生成PPT

图1   包含股骨、髋关节软骨和骨盆的解剖模型示意

图2

图2   软骨示意

2.1 数据的降维分析

因为股骨和软骨的维数较高，而且样本数非常少，如果对原始数据进行训练时非常容易出现过拟合的现象，因此考虑对其进行降维，使用较少的维数对其进行表示。主成分分析（principal component analysis，PCA）被广泛用于统计形状分析^[17]，本文使用PCA对股骨和软骨进行形状分析。

对于股骨数据，PCA之后得到特征值，特征值保留整数部分，按照从大到小的顺序排列之后，取前18个特征值，见表1。

表1   股骨数据进行PCA之后得到的特征值

从表1可以看到，排在前面的特征值数量级远大于后面的特征值。特征值越大，表示此特征方向上的信息越多。故经过 PCA 降维后，取前 9 个特征值即可在误差较小的情况下表示原始数据。

而对于软骨数据，因为数据生成的过程就是假定软骨不发生改变的，所以每个样本的软骨数据都是一致的，不必进行PCA。

对于应力数据，因为软骨上的每一个点都有值，所以应力数据共有 794个样本，每个样本有2 645个值。可以进行奇异值分解（singular value decomposition，SVD）分析。得到的奇异值共有794个，取前20个并保留整数部分，见表2。从表2 可以看出，应力数据的奇异值下降得很慢，前 20 个数据只占了总量的 94.03%，前40个数据占了96.95%，要达到99.00%需要97个数据，达到99.50%需要140个数据，而达到99.90%需要245个数据。可见对应力数据进行SVD不是很理想，但是依然可以使用SVD的方法进行低秩估计。

表2   应力数据进行SVD之后得到的奇异值

2.2 数据的聚类分析

使用 k-means 对股骨数据进行聚类，5 类结果如图3所示。由图3可知，除某类数据只有少数几个样本外，股骨数据可大致分为4类。

图3

新窗口打开| 下载原图ZIP| 生成PPT

图3   股骨数据聚类示意

对股骨数据做完聚类之后，使用PCA和t分布随机邻域嵌入（t-distributed stochastic neighbor embedding，t-SNE）对应力数据进行降维。将应力数据降到三维之后，再把股骨的聚类标签加到应力样本中，得到的颜色效果如图4所示。笔者发现，上述 4 类数据（绿色、蓝色、紫色和红色）的大致分布聚集得比较好。但也有少数几个样本脱离了原来的类别，笔者暂且判定这是正常现象，因为决定应力值的不只是股骨数据，还有施加给股骨数据的力。

2.3 软骨的可视化分析与预处理

软骨数据只有一个样本，将软骨数据进行可视化，可以得到如图5所示的效果。

图4

新窗口打开| 下载原图ZIP| 生成PPT

图4   使用股骨数据聚类的标签（颜色）对应力数据的聚类

进一步地，每次只可视化 24×6=144 个点，把这144个点分为6组，每组24个数据（作为一段），各组用不同的颜色表示，做5次可视化后，得到如图6（a）~图6（e）所示的结果，将前5次可视化的结果叠加，得到第6次的结果，如图6（f）所示。从图6（f）可以看出，软骨的数据排列得非常不规则，对这样的数据进行卷积计算，难以从局部进行特征提取。

图5

新窗口打开| 下载原图ZIP| 生成PPT

图5   软骨数据的可视化

对软骨数据进行 PCA 处理，找到它在三维空间中的主成分向量，如图7（a）所示。

图6

新窗口打开| 下载原图ZIP| 生成PPT

图6   每24个点作为一段得到软骨形状分布

适当平移中心位置之后（如图7（b）所示），将第2个向量绕第3个向量旋转360°，每2.25°划分一块区域，可将软骨划分成160块区域。由于有些区域数据点较少或者完全没有数据点，最后只筛选数据点较多的100个区域，得到一个近似的软骨模型。笔者希望这个软骨模型的每个区域都有相同的数据点数，这里选择24个点。

图7

新窗口打开| 下载原图ZIP| 生成PPT

图7   软骨示意

对数据点数多于24 个的区域进行去点操作：计算数据点与中心点的距离，并根据这个距离对数据点进行排序，再计算相邻两点之间的距离，得到一个点与两个邻居的距离之和，每次去掉距离之和最小的那个数据点，直到这个区域的数据点数为24为止。

对于数据点数不够 24 个的区域，按照以下方法进行插值，得到新的数据点：计算数据点到中心点的距离，并根据这个距离对数据点进行排序，再计算相邻两点之间的距离，得到一个点与前一个邻居的距离，每次在与前一个邻居的距离最大的那个点前面插入一个点，并将这两个点的应力值取平均之后赋给新的点。

经过上述操作后，得到一个比较平整、规则的24×100的软骨模型，如图8所示。最后，在拟合得到了当前软骨模型的应力值，要恢复原始模型的应力值的情况下，如果原模型的某个数据点缺失了应力值，则找到离它最近的有应力值的点，使缺失应力值的点与该点一致。

图8

图8   平整、规则的24×100的软骨模型示意

3 基于深度学习的髋关节应力分析模型

本节介绍深度学习模型的主要实现细节，包括形状编码的实现、神经网络非线性映射的实现以及应力编码和解码的实现，并给出了两个性能评价指标，对深度学习模型的性能进行评价。

3.1 形状编码器

形状编码指的是通过低维的数据表示形状的方法。采用 PCA 对股骨数据进行降维表示。对于每一个股骨数据X，有：

3.2 神经网络层

图9

新窗口打开| 下载原图ZIP| 生成PPT

图9   神经网络层

3.3 应力的编码与解码

髋关节软骨上的应力分布可以由一组标量值表示（即编码）和重构（即解码）。编码和解码的过程可以通过具有多层线性单元的双向神经网络在概念上描述。应力编码只在训练阶段使用，而应力解码是深度学习模型的一部分。为了便于处理，将每个软骨样本分成20 个区域，每个区域由相同的数量点（即每个区域120个节点）组成，使用低秩近似（low rank approximation，LRA）^[19]进行应力编码和解码。

在解码过程中，整个软骨的应力码β_n将会转换为软骨各区域的应力码。由式（3）给出：

3.4 模型性能

把第3.1~3.3节中的模型搭建完成之后，随机选择 90%的样本数据进行训练，以获取模型的参数，剩下的 10%的样本数据作为测试集，以得到模型的性能指标。

模型的目的是得到一个近似 FEA 结果的应力分布，在训练中使用FEA结果作为每一个样本的标签，因此以模型在测试集上的预测结果与FEA结果的误差来评价模型性能。在这里，采取两个性能指标来评估模型的预测应力分布的准确度：平均绝对误差（mean absolute error，MAE）和归一化平均绝对误差（normalized mean absolute error，NMAE）。对于每个形状，MAE定义为：

上述模型经过训练和测试之后得到的应力值的MAE和NMAE见表3。

表3   模型的性能指标

对于具体的样本来说，如111号样本，作为测试集中的一个样本，它的MAE为0.11 MPa，NMAE为5.69%，可视化如图10所示。

图10

图10   111号样本的应力分布对比

同样为测试集中的样本的 542 号样本的MAE为0.05 MPa，NMAE为2.26%，可视化如图11所示。

图11

图11   542号样本的应力分布对比

对比深度学习模型得到的应力结果与 FEA 得到的应力结果可知，深度学习方法能够在一定程度上模拟FEA。但是由于其MAE和NMAE还比较大，这样的模型不能完全替代FEA。

4 模型改进方向

使用有限元分析方法时，为了得到应力值，需要给股骨施加一个力，现在分析一下深度学习模型有没有考虑到这个力的作用。在前文实现的模型中，输入数据由两部分组成，一部分是股骨和软骨的形状数据，另一部分是股骨对软骨的作用力大小。在之前的模型中，这两部分只是在神经网络输入的时候简单地叠加到了一起。对于每一个样本，股骨和软骨降维后有9维，再加上这个3维的力，每个样本在输入神经网络的时候是12维。这种简单叠加在一起的方式没有考虑量纲的问题，避震对其做出了一些调整。在这样的情况下，将力的大小从模型的输入中去掉，只将 9 维的股骨和软骨的形状坐标值输入神经网络中，如图12 所示，得到的结果为：新模型的MAE是0.09±0.16 MPa， NMAE是4.11%±3.16%，见表4。新模型的效果更好一些。

图12

新窗口打开| 下载原图ZIP| 生成PPT

图12   简化后的深度学习模型

表4   简化后的模型的性能指标

5 结束语

本文实现了一个基于深度学习的髋关节软骨的应力分布预测模型，该模型在分析软骨应力分布时，平均绝对误差为0.11± 0.22 MPa，归一化平均绝对误差为5.07% ± 4.47%，能够在一定程度上模拟有限元分析方法，但是由于其平均绝对误差和归一化平均绝对误差比较大，故这个模型还不能完全替代有限元分析方法。

除此之外，本文还进一步探索了深度学习模型的局限性^[19,20]，比如上述模型不能充分利用现有的力的信息对髋关节软骨的应力分布进行分析。笔者希望在以后的研究中解决这个问题，使深度学习模型在髋关节软骨应力的计算上有更好的效果。

参考文献

[1] AURICCHIO F, CONTI M, MORGANTI S, et al. Simulation of transcatheter aortic valve implantation: a patient-specific finite element approach[J]. Computer Methods in Biomechanics and Biomedical Engineering, 2014, 17(12): 1347-1357.
[2] CAPELLI C, BOSI G M, CERRI E, et al. Patient-specific simulations of transcatheter aortic valve stent implantation[J]. Medical & Biological Engineering & Computing, 2012, 50(2): 183-192.
[3] DWYER H A, MATTHEWS P B, AZADANI A N, et al. Migration forces of transcatheter aortic valves in patients with noncalcific aortic insufficiency[J]. The Journal of Thoracic and Cardiovascular Surgery, 2009, 138(5): 1227-1233.
[4] MORGANTI S, CONTI M, AIELLO M, et al. Simulation of transcatheter aortic valve implantation through patient-specific finite element analysis: two clinical cases[J]. Journal of Biomechanics, 2014, 47(11): 2547-2555.
[5] WANG Q, KODALI S, PRIMIANO C, et al. Simulations of transcatheter aortic valve implantation: implications for aortic root rupture[J]. Biomechanics and Modeling in Mechanobiology, 2015, 14(1): 29-38.
[6] WANG Q, SIROIS E, SUN W, et al. Patient-specific modeling of biomechanical interaction in transcatheter aortic valve deployment[J]. Journal of Biomechanics, 2012, 45(11): 1965-1971.
[7] MORGANTI S, BRAMBILLA N, PETRONIO A S, et al. Prediction of patient-specific post-operative outcomes of TAVI procedure: the impact of the positioning strategy on valve performance[J]. Journal of Biomechanics, 2015: 49(12): 2513-2519.
[8] DE JAEGERE P, DE SANTIS G, RODRIGUEZOLIVARES R, et al. Patient-specific computer modeling to predict aortic regurgitation after  transcatheter aortic valve replacement[J]. Jacc-cardiovascular Interventions, 2016, 9(5): 508-512.
[9] SUN W, MARTIN C, PHAM T, et al. Computational modeling of cardiac valve function and intervention[J]. Annual Review of Biomedical Engineering, 2014, 16(1): 53-76.
[10] LIANG L, LIU M, MARTIN C, et al. A machine learning approach to investigate the relationship between shape features and numerically predicted risk of ascending aortic aneurysm[J]. Biomechanics & Modeling in Mechanobiology, 2017, 16(5): 1-15.
[11] WANG Q, PRIMIANO C, MCKAY R G, et al. CT image-based engineering analysis of transcatheter aortic valve replacement[J]. Jacc-cardiovascular Imaging, 2014, 7(5): 526-528.
[12] BUENO-PALOMEQUE F L, CORTÉS-RODRÍGUEZ C J, GARCÍASARMIENTO C D, et al. A comparative case study of the pre and post operative stress distribution in a dysplastic hip joint[C]//Proceedings of the 11th International Symposium, Computer Methods in Biomechanics and Biomedical Engineering, April 3 - 7, 2013, Salt Lake City, Utah, USA. [S.l.:s.n.], 2013
[13] HELLWIG F L, TONG J, HUSSELL J G. Hip joint degeneration due to cam impingement: a finite element analysis[J]. Computer Methods in Biomechanics & Biomedical Engineering, 2016, 19(1): 41-48.
[14] 王正军, 孙立明. ANSYS 18.0 有限元分析基础与实例教程[M]. 北京: 机械工业出版社, 2018.
WANG Z J, SUN L M. ANSYS 18.0 finite element analysis basics and examples tutorial[M]. Beijing: China Machine Press, 2018.
[15] ZAMEER S, MOHAMMED M A H. Evaluation of contact stress distribution of hip joint model using finite element method[J]. Esrsa Publications, 2014, 3(6): 1603-1610.
[16] LIANG L, LIU M, MARTIN C, et al. A deep learning approach to estimate stress distribution: a fast and accurate surrogate of finite-element analysis[J]. Journal of the Royal Society Interface, 2018, 15(138).
[17] HEIMANN T, MEINZER H. Statistical shape models for 3D medical image segmentation: a review[J]. Medical Image Analysis, 2009, 13(4): 543-563.
[18] DUGAS C, BENGIO Y, BELISLE F, et al. Incorporating second-order functional knowledge for better option pricing[C]//Neural Information Processing Systems, November 27-December 2, 2000, Denver, USA. Boston: MIT Press, 2000: 472-478.
[19] 张钹. 人工智能进入后深度学习时代[J]. 智能科学与技术学报,2019, 1(1): 4-6.
ZHANG B. Artificial intelligence is entering the post deep-learning era[J]. Chinese Journal of Intelligent Science and Technology, 2019, 1(1): 4-6.
[20] 郑南宁. 人工智能新时代[J]. 智能科学与技术学报, 2019, 1(1): 1-3.
ZHENG N N. The new era of artificial intelligence[J]. Chinese Journal of Intelligent Science and Technology, 2019, 1(1): 1-3.

作者简介