A BONE MECHANICAL STUDY INTEGRATED BY DEEP-LEARNING METHOD AND MECHANICAL MODELING
-
摘要: 骨缺损是骨科临床常见且复杂的疾患, 根据患者体内缺损区骨组织力学性能, 设计生物力学性能匹配的个性化骨假体, 有望提升临床骨缺损诊疗的水平. 然而, 当前的个性化骨缺损诊疗, 在体内骨组织微观结构分析、非均质各向异性力学行为表征和建模等方面存在诸多问题, 难以实现生物力学性能的适配, 导致骨重建效果不佳. 针对上述问题, 提出了一种融合数据驱动与力学建模的骨缺损重建方法, 以实现临床条件下骨组织力学性能的准确表征. 首先, 以羊股骨远端为对象, 提出了基于临床CT影像的多神经网络模型, 通过建立低分辨率临床CT下的宏观骨密度分布与micro-CT下松质骨微结构形态特征的映射关系, 能够直接通过临床CT对体内骨组织非均匀骨密度分布和结构张量等组织形态学参数进行准确预测. 其次, 建立了基于非均匀骨密度和结构张量的松质骨各向异性本构模型和实验表征方法. 通过贝叶斯反演识别本构模型参数, 修正了实验中由于材料主方向与加载方向偏离引入的系统误差. 实验结果验证了所建立本构模型与参数反演方法的准确性, 并揭示了不同部位松质骨力学行为与微结构生长方向的关系. 文章通过深度学习与力学建模融合的骨力学性能研究, 解决了临床医学影像下松质骨微观结构分析的难题, 为个性化骨假体的设计奠定了基础.Abstract: Bone defects are common and complex conditions in orthopedic clinics. Designing personalized bone implants with biomechanical properties that match the mechanical properties of the bone tissue in the patient's defect area holds great promise for desired bone defect reconstruction. However, the current design of personalized bone implants faces numerous challenges in the microstructural analysis of bone tissues in vivo, the characterization and modeling of heterogeneous anisotropic mechanical behavior, making it difficult to achieve mechanical property matching, which results in suboptimal bone reconstruction outcomes. To address these issues, this paper establishes an integrated approach combining data-driven and mechanical modeling for the mechanical theory, computation, and experimental methods of bone defect reconstruction, enabling accurate characterization of the mechanical properties of bone tissue under clinical conditions. Firstly, the distal femoral bone tissues of sheep are adopted as the experimental subject, a data-driven model for accurately predicting the morphological parameters of cancellous bone tissue under clinical CT imaging was proposed. A multi-neural network model combining high-resolution micro-CT and clinical CT was established. By correlating the macroscopic bone density distribution from low-resolution clinical CT and the microstructural morphological characteristics of cancellous bone from high-resolution micro-CT, a mapping relationship between bone density distribution and microstructural characteristics was established, which realized the accurate prediction of morphological parameters such as heterogeneous bone density distribution and fabric tensor of in vivo bone tissue using clinical CT. Furthermore, an anisotropic constitutive model and a Bayesian calibrated experimental method for cancellous bone based on heterogeneous bone density and fabric tensor were developed, which revealed the relationship between the mechanical behavior of cancellous bone at different locations and the growth direction of its microstructure. Combining a Bayesian-based method for identifying constitutive model parameters, the systematic errors introduced by the deviation between the principal material direction and the loading direction in cancellous bone experiments were corrected. The accuracy of the established constitutive model and parameter identification method was validated through experiments, addressing the challenges of microstructural analysis of cancellous bone under clinical medical imaging, and lays the foundation for the design of personalized bone implants.
-
引 言
骨缺损是因组织病变、创伤或手术导致的骨质缺失. 在临床上, 骨科骨缺损重建多集中在脊柱、髋关节、股骨远端以及胫骨平台等松质骨丰富的区域. 图1所示为股骨一侧的骨干部位缺损以及股骨远端缺损. 相关研究表明[1], 除体内生物化学环境、血供和细胞活动等生物因素外, 合适的力学刺激是诱导骨重建与再生的条件之一, 过小的力难以有效刺激成骨, 而过大的力学刺激则往往造成病理性骨生长, 难以形成具备承载力的新生骨组织. 骨缺损的临床治疗需要通过内、外固定以及骨植入物, 使自然骨与固定装置或植入物产生相互作用, 从而诱导骨应力的产生, 使其处于正常成骨的应变阈值内, 有效促进骨重建. 因此, 骨缺损的精准诊疗依赖于对患者骨组织和骨结构力学行为的分析. 其中, 骨组织的弹性模量、疲劳强度和断裂韧性等力学因素对骨组织力学性能的表征以及骨缺损的诊疗具有重要作用.
![]()
图 1 临床影像下典型的骨缺损形式及其数字建模Figure 1. Typical form of bone defect and its digital modeling under clincial imaging骨是典型的多层级结构[2], 宏观上可以按照孔隙率与组织形态差异分为外部包裹的密质皮质骨及内部的海绵状松质骨. 皮质骨由较为致密的成熟骨单位组成Haversian系统[3-4], 提供刚度和强度性能, 孔隙率为3% ~ 5%, 模量为10 ~ 20 GPa, 强度为110 ~ 220 MPa. 松质骨为典型的小梁骨组成的多孔点阵结构, 小梁骨的特征尺寸为100 ~ 200 μm, 孔隙率为50% ~ 90%, 等效模量为0.02 ~ 3.7 GPa, 提供血供和动态负载的顺应性调节等功能. 在宏观尺度下, 皮质骨表现为横观各向同性, 松质骨表现为正交各向异性. 对于松质骨, 其等效力学性质受多种因素的影响, 具有较大的分散性, 例如解剖部位、体积分数和材料方向的差异性. 医学上对骨结构的检测, 常用的影像检测手段有普通CT (computed tomography)、定量CT (quantitative CT, qCT)、外周高分辨率CT (HR-p-QCT)、显微CT (micro-CT)、双能X射线吸收DEXA (dual energy X-ray absorbtiometry)以及磁共振成像MRI (magnatic resonance imaging). 目前, 临床常用的低分辨率普通CT和定量CT等由于分辨率较低, 难以直接捕捉骨组织的微观结构特征. 而以高分辨率的micro-CT等影像技术虽然能够捕捉清晰的骨组织微结构, 但因其过大的辐射量, 难以用于临床大段骨的检查.
近年来, 骨缺损重建成为力学和医学领域的研究热点, 同时也面临诸多挑战性难题: 首先, 骨组织的力学非均质性和结构各向异性规律复杂; 其次, 临床用于骨缺损诊断的CT影像分辨率较低, 难以在体内环境下建立骨微观组织微观结构与其力学性能的关系. 这些难题构成了骨缺损临床治疗的瓶颈. 特别是对于需要人工骨移植的骨缺损患者, 由于难以在临床条件下精准表征缺损区骨组织本身的力学性质, 当前骨科植入物未能充分考虑患者骨缺损部位的个体差异, 不同患者的治疗效果差异性较大. 因此, 当前骨缺损的精准诊疗亟需针对当前医学检查技术, 建立合适的力学理论模型、计算方法和实验方法, 定量描述缺损区的骨结构及其功能. 为解决上述骨缺损个性化治疗亟待解决的关键问题, 就需要深入认知骨组织的力学性质与微结构形式, 找到无创、快速及精准地预测松质骨组织力学性能的方法.
当前, 机器学习在工程和科学的各个领域均显示出巨大的应用前景: 基于物理信息的深度学习网络PINN (physics-informed neural network)可以用于快速求解偏微分方程[5-8], 深度神经网络DNN (deep neural network)可以作为代理模型[9-10], 实现高效地优化和设计任务[11-12]. 特别是随着以深度学习方法为代表的人工智能技术在医学以及力学领域的应用不断深入[12-14], 通过数据驱动与力学模型结合的方法, 将有望解决当前临床低分辨率医学影像以及骨组织生物力学表征等领域的关键性问题[15], 实现满足临床治疗需求的缺损区骨组织生物力学功能的分析研究.
综上, 为解决临床骨缺损诊疗中松质骨组织力学性能表征的问题, 本文提出micro-CT与临床CT影像结合的数据驱动方法, 分析不同位置松质骨微结构形态学特征, 在临床CT影像数据下精确表征与松质骨非均质各向异性力学性能相关的关键微结构特征与参数. 此外, 本文发展临床CT下基于骨密度与结构张量的松质骨非均质各向异性参数化本构模型, 并建立多轴加载实验方法与参数识别算法. 结合贝叶斯推断的模型参数识别方法, 实现了在三维应力状态下对模型参数的反演与不确定性分析, 从而提高参数识别的准确性, 并揭示不同部位松质骨微结构特征与力学性能之间的关系.
1. 深度学习CT影像表征松质骨力学特征
1.1 基于临床CT影像的微结构特征分布重建
首先, 本文选用成年山羊股骨远端展开相关研究, 如图2所示. 挑选了20只大于1周岁的成年山羊作为实验动物, 截取外形与受力环境较为复杂的股骨远端进行相关研究. 分别通过临床定量CT以及micro-CT扫描设备对股骨远端标本进行高精度的数字重建. 临床定量CT成像空间分辨率为0.6 mm, micro-CT成像空间分辨率为0.06 mm. 显然, 定量CT下松质骨小梁的网络结构模糊, 难以分辨清晰的骨小梁结构; 作为对比, 图2(c)所示的micro-CT (空间分辨率0.06 mm)可以清晰地重现骨小梁微结构. 此处, 由于需要使两种分辨率下的体积胞元VE (volume element)一一对应, 因此需要对异源设备下重建的数字影像进行配准, 如图2(d)所示.
传统方法通常依赖高分辨率的医学影像如micro-CT或HR-MRI等方式获取骨组织的微结构形态. 由于整体结构外形复杂, 松质骨的研究往往采用体积胞元的方式展开. 利用micro-CT影像进行数字重建获取的基于高分辨率松质骨体积胞元, 通过图形学分析获取对应的形态学参数, 通过有限元计算获取等效各向异性力学性质, 并通过三维卷积神经网络获取其微结构特征规律.
针对临床条件下的松质骨力学性能预测, 本文的方法建立了低分辨率临床CT与微结构特征规律的关系, 并通过该特征规律实现各向异性等效力学性质与形态学参数. 该方法需要克服两个难点: 首先面对微结构信息确实的临床低分辨率LR (low-resolution) CT影像, 能否从中找到准确的微结构特征规律; 其次, 这些微结构特征规律能否与各向异性等效力学性质和形态学参数建立联系.
基于上述分析, 本文通过数据驱动方法, 一方面将松质骨的非均匀骨密度和微结构取向分布等与松质骨非均质各向异性力学性能表征强相关的形态学特征作为预测目标; 另一方面通过基于高分辨率CT (micro-CT)影像建立高精度有限元模型, 计算基于RVE的各向异性弹性矩阵标签数据库. 将数据驱动方法提供的关键松质骨微结构形态学特征, 与能反映骨组织非均质性、各向异性以及生长方向性的各向异性等效力学性质相结合, 形成数据驱动与力学模型融合的多尺度建模方法.
具体地, 本文通过micro-CT中提取准确的微结构特征规律, 并以此为标签, 建立低分辨率CT影像到准确的微结构特征规律之间的映射关系. 同时, 引入两次上采样模型对低分辨率CT影像进行处理, 提升临床低分辨率CT影像到目标性能的预测精度. 具体的建模与分析流程如图3所示.
![]()
下面采用${{{\boldsymbol{P}}}^H}$表示待预测的宏观等效力学性质或形态学参数, 并以松质骨等效力学性能为例对本模型进行分析讨论. 假设骨小梁为各向同性线弹性材料, ${{{\boldsymbol{P}}}^H}$可表示为
$$ {{{\boldsymbol{P}}}^H} = {{{\boldsymbol{P}}}^H}\left\{ {\lambda ,\mu ,{{{\boldsymbol{\varPsi}} }_{SFI}}\left[ {{{\boldsymbol{K}}}\left( {{{\text{χ}}} } \right)} \right]} \right\} $$ (1) 其中, $\lambda $和$\mu $为拉梅常数, ${{\boldsymbol{K}}}\left( {{\text{χ}} } \right)$表示松质骨的微结构分布, ${{\text{χ}} }$表示医学影像的采样点分布. ${{{\boldsymbol{\varPsi}} }_{SFI}}\left[ {{{\boldsymbol{K}}}\left( {{\text{χ}} } \right)} \right]$则表示形状特征函数, 用以描述${{\boldsymbol{K}}}\left( {{\text{χ}} } \right)$中对${{{\boldsymbol{P}}}^H}$有贡献的形状特征函数.
在高分辨率的CT影像下, 材料域${{\boldsymbol{K}}}$内的采样点为${{{\text{χ}} }_{HR}}$. 对应地, 在临床低分辨率影像下, 采样点为${{{\text{χ}} }_{LR}}$. 假设${{\boldsymbol{K}}}\left( {{{{\text{χ}} }_{HR}}} \right)$能够充分捕捉可准确预测松质骨的组织形态学参数和各向异性力学性能, 则可通过图中所描述的VAE[16-19]的模型重建${{\boldsymbol{K}}}\left( {{{{\text{χ}} }_{HR}}} \right)$并获取符合标准正态分布的编码特征向量${{\boldsymbol{f}}}_{HR}^{{\mathrm{encoded}}}$. 一般地, 对于训练好的VAE模型, 编码特征向量${{\boldsymbol{f}}}_{HR}^{{\mathrm{encoded}}}$能够包含松质骨复杂空间结构与材料分布的所有信息. 此处, 表示重建误差的损失函数为
$$ Los{s_{{\mathrm{recons}}}} = \frac{1}{N}\sum\limits_{i = 1}^N {\left\| {{{\left[ {{{\boldsymbol{K}}}\left( {{{{\text{χ}} }_{HR}}} \right)} \right]}_i} - \left[ {{{\boldsymbol{K}}}\left( {{{{\text{χ}} }_{HR}}} \right)} \right]_i^{{\mathrm{recons}}}} \right\|} _2^2 $$ (2) 其中, $ \left[ {{{\boldsymbol{K}}}\left( {{{{\text{χ}} }_{HR}}} \right)} \right]_i^{{\mathrm{recons}}} $为第$i$个样本的重建样本.
KL散度[10,20-22]部分则通过特征向量${{\boldsymbol{f}}}_{HR}^{{\mathrm{encoded}}}$的均值和方差描述, 即
$$\begin{split} &Los{s_{KL - {\mathrm{divergence}}}} = - 0.5 \sum\limits_{i = 1}^N \left[ 1 + \lg \left( {Var{{\left[ {{\boldsymbol{f}}_{HR}^{{\mathrm{encoded}}}} \right]}_i}} \right) -\right.\\ &\qquad \left.{E^2}{\left[ {{\boldsymbol{f}}_{HR}^{{\mathrm{encoded}}}} \right]_i} - Var{\left[ {{\boldsymbol{f}}_{HR}^{{\mathrm{encoded}}}} \right]_i} \right]\end{split} $$ (3) 因此, 结合式 (2)和式 (3), 图3中第1步通过micro-CT重建的高精度松质骨模型提取松质骨微结构特征所使用的3D VAE模型损失函数可表示为$ Los{s_{VAE}} = Los{s_{{\mathrm{reconstruction}}}} + Los{s_{KL - {\mathrm{divergence}}}} $. 对于每一个样本, 在临床低分辨率下对应的样本为${{\boldsymbol{K}}'}\left( {{{{\text{χ}} }_{LR}}} \right)$. 考虑到${{\boldsymbol{K}}'}\left( {{{{\text{χ}} }_{LR}}} \right)$中的采样点过少, 这里先通过空间三次样条插值, 增加采样点的数量至${{{\text{χ}} '}_{LR}}$. 进而, 通过超分辨率上采样模型FSRCNN, 是将低分辨率下的像素矩阵${{\boldsymbol{K}}'}\left( {{{{{\text{χ}} '}}_{LR}}} \right)$扩充至$ {{\boldsymbol{K}}''}\left( {{{{{\text{χ}} '}}_{LR}}} \right) $. 采用两步上采样模型的目的是通过${{\boldsymbol{K}}'}\left( {{{{\text{χ}} }_{LR}}} \right)$, 以${{\boldsymbol{f}}}_{HR}^{{\mathrm{encoded}}}$为标签, 实现对微结构特征的精准预测. 具体地, 定义$Los{s_{{\mathrm{up}} - {\mathrm{sampling}}}}$为超分辨率模型重建的损失函数, $Los{s_{{\mathrm{feature}}}}$为以特征向量预测的RMSE误差的损失函数. 因此, 第2步基于临床CT的两步上采样并提取微结构特征向量的损失函数表示为
$$\left.\begin{split} & Los{s_{{\text{two-step}}\;{\text{model}}}} = Los{s_{{\mathrm{up}} - {\mathrm{sampling}}}} + Los{s_{{\mathrm{feature}}}} \\ & Los{s_{{\mathrm{up}} - {\mathrm{sampling}}}} = \frac{1}{N}\sum\limits_{i = 1}^N {\left\| {{{\left[ {{{\boldsymbol{K}}'}\left( {{{{{\text{χ}} '}}_{LR}}} \right)} \right]}_i} - {{\left[ {{{\boldsymbol{K}}''}\left( {{{{{\text{χ}} '}}_{LR}}} \right)} \right]}_i}} \right\|} _2^2 \\ & Los{s_{{\mathrm{feature}}}} = \frac{1}{N}\sum\limits_{i = 1}^N {\left\| {{{\left[ {{{\boldsymbol{f}}}_{HR}^{{\mathrm{encoded}}}} \right]}_i} - {{\left[ {{{{\boldsymbol{f}}}_{LR}}} \right]}_i}} \right\|_2^2} \end{split}\right\} $$ (4) 由此, 临床低分辨率CT直接预测准确的松质骨微结构特征向量, 就可以进一步表示为利用 对目标等效力学性质 进行预测. 此处, 可通过全连接的神经网络直接实现, 该神经网络的损失函数为
$$ Los{s_{{\mathrm{predict}}}} = \frac{1}{N}\sum\limits_{i = 1}^N {\left\| {{{\left[ {{{{\boldsymbol{P}}}^H}} \right]}_i} - \left[ {{{{\boldsymbol{P}}}^H}} \right]_i^{{\mathrm{predict}}}} \right\|} _2^2 $$ (5) 基于上述深度学习模型, 建立临床CT和从micro-CT中编码得到的微结构特征向量, 可实现临床CT松质骨组织形态学参数和力学性能的准确预测. 本文所建立的3D VAE模型结构中, 编码器使用$2 \times 6$层卷积层, 采用ReLU为激活函数, 输出为1024维的特征向量. 对应地, 解码器使用$2 \times 6$层反卷积层, 采用ReLU为激活函数, 卷积核均采用(3, 3, 3), 步幅为1, 填充为1. 超分辨率特征提取模型由超分辨率模块和特征抽取模块组成, 超分辨率模块由8层卷积层和1层反卷积层组成, 输出为$100 \times 100 \times 100$的像素矩阵. 特征抽取模块采用5层卷积层, 输出为1024维的特征向量, 采用ReLU为激活函数, 卷积核均采用(3, 3, 3) 步幅为1, 填充为1. 微结构特征向量为输入的多层感知机预测模型的隐藏层, 由5层全连接层(512, 256, 256, 128, 64)组成, 可用于预测各向异性力学性能、组织形态学参数和结构张量.
1.2 松质骨各向异性力学参数
松质骨各向异性力学参数的计算, 需要考虑弹性矩阵${{\boldsymbol{C}}}$以及材料主方向的计算.
研究发现松质骨的宏观等效力学性质基本满足正交各向异性力学行为[15]. 对于正交各向异性弹性矩阵的计算, 独立参数共有9个, 需要3个单轴压缩、3个纯剪切共计6个独立的加载工况.
此处, 骨小梁基体材料假设为各向同性线弹性材料, 杨氏模量为12 GPa, 泊松比为0.3. 对于待计算的均匀化弹性矩阵${\boldsymbol{C}}_{ijkl}^H$[23–25], 即
$$ {\boldsymbol{C}}_{ijkl}^H = \frac{1}{{\left| \varOmega \right|}}\int_\varOmega {{C_{pqrs}}\left( {\varepsilon _{pq}^{0\left( {ij} \right)} - \varepsilon _{pq}^{\left( {ij} \right)}} \right)\left( {\varepsilon _{rs}^{0\left( {kl} \right)} - \varepsilon _{rs}^{\left( {kl} \right)}} \right){{\mathrm{d}}}\varOmega } $$ (6) 其中, ${C_{pqrs}}$为骨小梁的弹性常数, $\varepsilon _{pq}^{0\left( {ij} \right)}$为施加在RVE上的宏观均匀应变场, $\varepsilon _{pq}^{\left( {ij} \right)}$为内部骨小梁的应变, 可以表示为
$$ \varepsilon _{pq}^{\left( {ij} \right)} = {\varepsilon _{pq}}\left( {{{{\boldsymbol{X}}}^{ij}}} \right) = \frac{1}{2}\left( {{{\boldsymbol{X}}}_{p,q}^{ij} + {{\boldsymbol{X}}}_{q,p}^{ij}} \right) $$ (7) 其中, ${{{\boldsymbol{X}}}^{ij}}$为与虚位移${{{\boldsymbol{X}}}^*}$对应的位移场, 定义为
$$ \int_\varOmega {{C_{ijkl}}{\varepsilon _{ij}}\left[ {{{\left( {{\boldsymbol{X}}} \right)}^*}} \right]{\varepsilon _{pq}}\left( {{{{\boldsymbol{X}}}^{kl}}} \right){{\mathrm{d}}}\varOmega } = \int_\varOmega {{C_{ijkl}}{\varepsilon _{ij}}\left[ {{{\left( {{\boldsymbol{X}}} \right)}^*}} \right]\varepsilon _{pq}^{0\left( {kl} \right)}{{\mathrm{d}}}\varOmega } $$ (8) 此外, 松质骨宏观等效的弹性力学行为可描述为$ {{{\boldsymbol{\sigma}} }^{{\mathrm{macro}}}} = {{{\boldsymbol{C}}}^H}:{{{\boldsymbol{\varepsilon}} }^{{\mathrm{macro}}}} $. 那么, 假设此时有一个沿方向${{\boldsymbol{d}}}\left( {{{\boldsymbol{d}}} \in {\mathbb{R}^2}} \right)$载荷[26-27], 沿该方向的等效杨氏模量可表示为
$$ {{\boldsymbol{E}}}\left( {{\boldsymbol{d}}} \right) = \frac{1}{{\left( {{{\boldsymbol{d}}} \otimes {{\boldsymbol{d}}}} \right):{{\left( {{{{\boldsymbol{C}}}^H}} \right)}^{ - 1}}:\left( {{{\boldsymbol{d}}} \otimes {{\boldsymbol{d}}}} \right)}} $$ (9) 进而, 材料主方向可通过结构张量${{\boldsymbol{M}}}$进行计算, 具体可通过平均截距长度法MIL (mean intercept length)[28-32]以及Minkowski张量[31,33-34]. MIL定义为扫掠方向${{\boldsymbol{n}}}$上射线经过的材料域长度总和, 即下式. 本文结构张量${{\boldsymbol{M}}}$使用MIL方法获得, 即
$$\left.\begin{split} & {{\boldsymbol{M}} = }\sum\limits_{i = 1}^3 {{\lambda _i}{{{\boldsymbol{M}}}_i}} = \sum\limits_{i = 1}^3 {{\lambda _i}\left( {{{{\boldsymbol{m}}}_i} \otimes {{{\boldsymbol{m}}}_i}} \right)} \\ & 1/{L^2}\left( {{\boldsymbol{n}}} \right) = {{\boldsymbol{n}}} \cdot {{\boldsymbol{M}}} \cdot {{\boldsymbol{n}}} \end{split} \right\}$$ (10) 其中, ${\lambda _i}\left( {i = 1,2,3} \right)$和${{{\boldsymbol{m}}}_i}\left( {i = 1,2,3} \right)$为${{\boldsymbol{M}}}$的特征值和特征向量. 特征值${\lambda _1}$, ${\lambda _2}$和${\lambda _3}$以及特征向量${{{\boldsymbol{m}}}_1}$,${{{\boldsymbol{m}}}_2}$和${{{\boldsymbol{m}}}_3}$分别对应于松质骨的3个材料主方向.
由于生物组织本身的个体差异和明显的离散性, 通过MIL得到的结构张量${{\boldsymbol{M}}}$可能难以满足正交各向异性材料的特点. 为使结构张量${{\boldsymbol{M}}}$符合正交各向异性材料的特点, 本文采用数值优化方法[35-36]调整结构张量${{\boldsymbol{M}}}$. 首先, 对MIL得到的结构张量${{\boldsymbol{M}}}$进行特征值分解$ {{\boldsymbol{M}}}{\boldsymbol{ = }}{{\boldsymbol{V}} }\Lambda{{{\boldsymbol{V}}}^{\mathrm{T}}} $, 其中${{\boldsymbol{M}}}$是特征向量矩阵, 是对角特征值矩阵. 将特征值分解的结果作为初始值, 并构建目标函数$J\left( {{\boldsymbol{M}}} \right)$评估结构张量的特征值和特征向量与理想情况的偏差, 即
$$ J\left( {{\boldsymbol{M}}} \right) = \sum\limits_{i = 1}^3 {{{\left( {\frac{{{{\boldsymbol{v}}}_i^{\mathrm{T}}{{\boldsymbol{M}}}{{{\boldsymbol{v}}}_i}}}{{{\lambda _i}}} - 1} \right)}^2}} $$ (11) 通过BFGS优化算法对式(11)所示的目标函数进行最优化. 当$J\left( {{\boldsymbol{M}}} \right)$达到最小值时, 结构张量${{\boldsymbol{M}}}$的特征值和特征向量取得最符合正交各向异性材料特性的值. 该方式能够有效处理由于生物组织的个体差异和离散性引起的问题, 更准确地确定了松质骨材料的主方向和各向异性特性.
1.3 预测模型验证与分析
本文所建立训练样本集中共有10000例, 表1对所使用训练样本的分布范围进行了统计.
表 1 训练样本集的4个组织形态学参数分布范围Table 1. The distribution range of four histomorphological parameters of the training sample setHistomorphological parameters Distribution range BV/TV $0.48 \pm 0.21$ DA $0.51 \pm 0.13$ Conn.D $4.55 \pm 1.76$ BS/TV $3.93 \pm 0.49$ 对于组织形态学参数[30,37-38]的计算, 这里选取体积分数BV/TV (bone volume/tissue volume)、各向异性程度DA (degree of anisotropy)、连通度Conn.D (connectivity degree)以及比表面积BS/TV (bone surface/tissue volume)作为形态学参数, 用以评估所有样本的范围.
图4(a) ~ 图4(d)分别为BV/TV、DA、Conn.D以及BS/TV的预测结果, 预测值与真实值的相关系数分别为0.9981, 0.9829, 0.9711以及0.9812, 均方根误差分别为0.0092, 0.011, 0.153以及0.172 mm−1. 上述结果表明, 通过低分辨率CT影像重构的微结构特征向量${{{\boldsymbol{f}}}_{LR}}$准确地预测了与松质骨等效力学性能和内部微结构分布表征密切相关的组织形态学参数, 即非均匀的骨密度分布能够反映松质骨微结构形态学特征. 表1给出了所有样本的形态学参数分布范围. 该统计结果说明, 本文所使用的山羊股骨远端松质骨样本拥有足够大的覆盖范围.
此处共使用3000例样本组织形态学参数的预测结果进行验证, 预测结果如图4所示.
弹性刚度矩阵以及结构张量的预测评估, 则通过下式所示的标量误差进行分析, 即
$$ \left.\begin{split} & {\delta _M} = \sqrt {\frac{{\left( {{{\boldsymbol{M}}}^* - {{\boldsymbol{M}}}} \right) \cdot \left( {{{\boldsymbol{M}}}^* - {{\boldsymbol{M}}}} \right)}}{{{{\boldsymbol{M}}} \cdot {{\boldsymbol{M}}}}}} \\ & {\delta _C} = \sqrt {\frac{{\left( {{{\boldsymbol{C}}}^* - {{\boldsymbol{C}}}} \right):\left( {{{\boldsymbol{C}}}^* - {{\boldsymbol{C}}}} \right)}}{{{{\boldsymbol{C}}}:{{\boldsymbol{C}}}}}} \end{split}\right\} $$ (12) 对于弹性矩阵预测准确性的评估, 则通过相关系数评估预测结果. 如图5所示, 对于结构张量, 预测的平均误差为1.63%, 最大误差为5.73%, 说明该模型能够准确预测结构张量.
![]()
图 5 正交各向异性弹性刚度矩阵预测结果Figure 5. The prediction results of orthotropic elastic stiffness matrix通过预测模型得到的弹性矩阵与均匀化分析得到的弹性矩阵的相关系数分别为0.9973, 0.9929, 0.9921, 0.9888, 0.9849, 0.9722, 0.9645, 0.9342和0.9566. 结果表明, 该模型能够有效建立从重构的微结构特征到弹性矩阵以及结构张量的映射关系. 因此, 即使在临床低分辨率CT影像的限制下, 本方法仍能从中提取足够多的有效信息, 重构松质骨的微结构特征, 并最终利用这些微结构特征实现对各向异性力学性能的准确预测.
此外, 为验证本模型的整体预测精度, 此处选择3个松质骨立方体模型, 比较基于micro-CT图像的有限元方法和基于低分辨率临床CT影像的预测模型给出弹性矩阵${{\boldsymbol{C}}}$和结构张量${{\boldsymbol{M}}}$的预测结果. 这里, ${{\boldsymbol{C}}}$和${{\boldsymbol{M}}}$则取非0的独立项, 可简化为
$$\left. \begin{split} & {{\boldsymbol{C}}} = \left[ {\begin{array}{*{20}{c}} {{C_{11}}}&{{C_{22}}}&{{C_{33}}} \\ {{C_{44}}}&{{C_{55}}}&{{C_{66}}} \\ {{C_{12}}}&{{C_{13}}}&{{C_{23}}} \end{array}} \right]\;\; \\ & {{\boldsymbol{M}}} = \left[ {{{{\boldsymbol{M}}}_1},{{{\boldsymbol{M}}}_2},{{{\boldsymbol{M}}}_3}} \right] = \left[ {\begin{array}{*{20}{c}} {{M_{11}}}&{{M_{12}}}&{{M_{13}}} \\ {{M_{21}}}&{{M_{22}}}&{{M_{23}}} \\ {{M_{31}}}&{{M_{32}}}&{{M_{33}}} \end{array}} \right] \end{split}\right\} $$ (13) 由图6可知, 预测的立方体体积分数分别为0.5677, 0.4995和0.4887. 棕色立方体模型均使用60 μm空间分辨率的micro-CT影像, 以骨小梁微结构形态表征. 相应的灰色立方体提取0.6 mm空间分辨率的低分辨率临床CT影像. 此处需要注意的是, 传统图形学分析得到的$ {{\boldsymbol{M}}} $不为单位向量, 模长不能保持为1. 为方便对比并保持其清晰的物理解释, 这里需要对预测的特征向量进行归一化.
![]()
图 6 Micro-CT有限元计算结果与临床低分辨率CT预测结果的对比Figure 6. Comparison of finite element calculation results of micro-CT and clinical low resolution CT prediction results与基于micro-CT影像的有限元法和传统图形学分析方法计算结构张量的结果相比, 本文所建立的基于深度学习的预测方法能够充分挖掘低分辨率临床CT影像中的松质骨微结构信息. 该方法成功预测的各向异性弹性力学性能和结构张量与通过micro-CT影像建立的高精度数字模型计算结果基本一致, 从而进一步验证了所建立模型的优势.
2. 临床CT影像下松质骨各向异性力学性能的表征方法
2.1 股骨远端松质骨力学实验
松质骨非均质各向异性力学性能精准实验表征是长期以来的难题. 在实验室条件下, 松质骨样本量少, 试样难以制备, 无法对单一试件进行重复多次的材料实验. 同时, 基于解剖经验制备的松质骨试样忽略了依赖于微结构空间分布的材料主方向和加载方向的差异所产生的系统误差, 试样在假设的“单轴”压缩实验下处于多轴应力状态, 影响了松质骨材料力学行为的准确表征.
针对股骨远端形状复杂且硬质的特性, 以及内部松质骨切割成骨块试样后难以定位的问题, 样本被切割成尺寸为$5.6{\text{ mm}} \times 5.6{\text{ mm}} \times 5.6{\text{ mm}}$的立方体试样, 如图7所示. 需要指出的是, 在试件加工过程中, 试样均按照解剖轴进行切割. 因此, 几乎所有试样的材料主方向都不与试样的棱边平行或垂直. 在这种情况下, 试样几乎全部处于多轴应力状态. 实验使用最大载荷为5 kN的材料试验机, 对骨立方试样进行准静态压缩试验, 加载速率设置为1 mm/min.
首先, 基于解剖学经验对试样进行切割加工, 将解剖面法向定义为加载坐标系主轴方向$\left( {{{{\boldsymbol{d}}}_1},{{{\boldsymbol{d}}}_2},{{{\boldsymbol{d}}}_3}} \right)$. 同时, 通过计算每个试样的结构张量$ {{\boldsymbol{M}}} $, 可以获得正确的材料主轴$\left( {{{{\boldsymbol{e}}}_1},{{{\boldsymbol{e}}}_2},{{{\boldsymbol{e}}}_3}} \right)$, 这里, 将材料主轴作为参考坐标, 以加载面法向用解剖学主轴与材料主轴的欧拉角$\left( {\theta ,\varphi ,\psi } \right)$表示两者的差异[39].
如图8所示, 图中的S1 ~ S11为11个松质骨独立试样的编号, 每个标号$S \cdot \left( {\theta ,\varphi ,\psi } \right)$则是解剖学主轴和材料主轴的欧拉角. 参考Rincón-Kohli等[40]的股骨颈松质骨的单轴多轴压缩实验, 图8所示为材料失效前的应力-应变曲线, 不同试样的力学响应可分为两种类型: 试件S1, S3, S4和S5屈服后保持较长的平台应力状态; 其余试件S2, S6和S7 ~ S11屈服后由于损伤软化, 应力明显下降.
![]()
图 8 不同应力状态下的松质骨名义应力-名义应变曲线Figure 8. Nominal stress-nominal strain curves of cancellous bone under different stress states通过结合各试样的解剖学主轴与材料主轴的欧拉角, 可以发现欧拉角可以定性地表示试样所处的应力状态及材料主方向与加载方向的偏离程度. 当解剖学主轴接近材料主轴时, 试样在压缩屈服后表现出较长的平台应力段; 而当欧拉角较大, 即解剖学主轴远离材料主轴时, 试样在屈服后表现出明显的软化行为. 这两种现象分别与文献中单轴和多轴加载的应力应变曲线一致. 结果表明, 大多数试样在测试时处于多轴加载应力状态. 因此, 不能再按照单轴加载进行计算, 而需要考虑其真实的多轴应力状态, 以确保得到的各向异性材料力学性能的准确性.
2.2 基于贝叶斯推断的材料参数反演框架
由于松质骨呈现出的正交各向异性力学行为与其微结构密切相关, 可使用基于骨密度和结构张量的正交各向异性模型描述其在弹性阶段的响应. 本文采用Zysset-Crunier模型[36,41-42], 即
$$ \begin{split} &\boldsymbol{C}(\bar{\rho}, \boldsymbol{M}) =\sum_{\substack{i=1}}^3\left(\lambda_0+2 \mu_0\right) \bar{\rho}^k m_i^{2 l} \boldsymbol{M}_i \otimes \boldsymbol{M}_i + \\ &\qquad\sum_{\substack{i, j=1 \\ i\ne j}}^3 \lambda'_0 \bar{\rho}^k m_i^l m_j^l \boldsymbol{M}_i \otimes \boldsymbol{M}_j +\\ &\qquad\sum_{\substack{i, j=1 \\ i\ne j}}^3 2 \mu_0 \bar{\rho}^k m_i^l m_j^l \boldsymbol{M}_i \overline{\underline{\otimes }} \boldsymbol{M}_j \end{split} $$ (14) 其中, $\rho $和$ {{\boldsymbol{M}}} $基于松质骨的微结构分布, 与具体的材料性质无关. ${\lambda _0},$ ${\mu _0},$ ${\lambda '_0},$ $k$和$l$为材料参数. 另外, $\rho $为相对密度, 取值范围在0 ~ 1, ${m_i}$为结构张量的特征值. $\overline{\underline{\otimes }} $表示张量积.
式(14)所示的原始Zysset-Crunier模型需要micro-CT等临床难以使用的高精度CT影像获得的微结构特征可通过上节所建立的预测模型直接获取, 但相对密度难以反映骨小梁本身的密度非均匀性, 使得该模型在临床中难以使用, 且模型的泛化能力差. 为解决上述问题, 本文通过临床定量CT获取真实的骨密度并引入Zysset-Crunier模型, 反映骨结构在不同层级上的非均匀性. 此外, $\bar \rho $和$ {{\boldsymbol{M}}} $均可通过上节所建立的预测模型直接获取, 使该模型可在临床影像下直接使用.
基于改进后的Zysset-Crunier模型, 结合实验数据可发现, 所有试样均处于多轴加载的工况, 无法通过单一的应变状态确认独立的材料参数. 因此, 需要考虑基于复杂应力状态下的本构模型材料参数识别方法.
针对本构模型参数的辨识, 可建立多参数识别的贝叶斯反演模型[43]. 首先, 实验中采集到的力信号可通过上述所涉及到的结构参数和等效材料参数表示, 即
$$ {F_r} = f\left( {{\boldsymbol{\varTheta }};\rho ,{\boldsymbol{M}},\bar u\left( {\boldsymbol{\varepsilon }} \right),{\boldsymbol{a}}} \right) $$ (15) 其中, ${{\boldsymbol{\varTheta}} = }\left[ {{\lambda _0},{\mu _0},{{\lambda '}_0},k,l} \right]$为本构模型的材料参数, 试样尺寸参数为${{\boldsymbol{a}}} = \left[ {{a_x},{a_y},{a_z}} \right]$, 压缩实验中所施加的位移载荷为$\bar u$, 所对应的传感器采集到的反力为${F_r}$.
实际上, 在一次独立的压缩测试中, 仅${F_r}$和$\bar u\left( {{\boldsymbol{\varepsilon}} } \right)$可以通过实验直接测量, 压缩位移$\bar u\left( {{\boldsymbol{\varepsilon}} } \right)$所对应的真实应变状态则无法测量. 在此情况下, 通过已知的${F_r}$和$\bar u$以及结构参数$\rho $和$ {{\boldsymbol{M}}} $, 求解材料参数${{\boldsymbol{\varTheta}} }$将成为不适定的反问题, 适合建立贝叶斯参数反演模型[43-44], 将原本参数求解的确定性问题转化为根据实验结果求解参数取值的概率分布问题, 使模型参数的求解成为适定问题.
本文建立的基于贝叶斯推断的松质骨参数识别整体框架如图9所示.
![]()
图 9 基于贝叶斯推断的材料反演整体框架Figure 9. The whole framework of material parameters identification based on Bayesian inference首先, 对羊的股骨远端松质骨试样进行临床定量CT医学影像扫描, 完成股骨远端整体模型的数字建模, 并制备标准的立方体压缩试样. 通过结构张量$ {{\boldsymbol{M}}} $确定材料主方向后, 对试样进行压缩实验, 记录不同加载状态下试样的应力-应变响应曲线.
其次, 利用修正的Zysset-Crunier模型, 对该试样在压缩载荷下的应力-应变曲线进行数值计算, 形成多材料参数的训练集, 建立正问题代理模型, 快速完成不同材料参数组合下试样应力-应变响应曲线的预测, 实现高效、准确的三维问题计算. 同时, 建立松质骨材料参数反演的贝叶斯推断模型, 建立各材料参数的后验概率分布模型, 结合实验数据, 对各材料参数进行推断.
最终, 基于所获得的材料参数, 表征松质骨的非均质性与各向异性, 准确评估松质骨整体力学性能.
具体地, 每个试件的实验结果均被考虑为一次对公共参数空间的观测, 其结果通过${{{\boldsymbol{y}}}_{\exp }}$记录. 正问题代理模型的功能是能够根据不同的结构参数和材料参数组合, 在相同的加载条件下, 给出对应的反馈值$F_r^{{\mathrm{surrogate}}}.$ ${{{\boldsymbol{y}}}_{\exp }}$和$F_r^{{\mathrm{surrogate}}}$的一致性作为判断参数是否合适的验证指标. 具体地, 贝叶斯推断参数的过程是一个多次迭代的过程, 在第$i$次迭代时, 贝叶斯模型将当前迭代给出的参数${{{\boldsymbol{\varTheta}} }_i}$传递到正问题代理模型, 代理模型给出对应的反力值${\left( {F_r^{{\mathrm{surrogate}}}} \right)^i}$并传递回贝叶斯模型, 通过判断第$i$次反力曲线和实验曲线${{{\boldsymbol{y}}}_{\exp }}$的差别, 进行参数更新, 得到更新后的参数${{{\boldsymbol{\varTheta}} }_{i + 1}}$传至代理模型, 直至代理模型得到的反力值和实验值保持一致. 通过统计有效次迭代过程中, 每次迭代的参数值以及对应的${\left( {F_r^{{\mathrm{surrogate}}}} \right)^i}$和${{{\boldsymbol{y}}}_{\exp }}$差别的大小, 可得到各参数的后验概率分布, 即
$$ p\left( {{\boldsymbol{\varTheta }},{\boldsymbol{h}}|{{\boldsymbol{y}}_{\exp }}} \right) \propto p\left( {{{\boldsymbol{y}}_{\exp }}|{\boldsymbol{\varTheta }}} \right)p\left( {{\boldsymbol{\varTheta }}|{\boldsymbol{h}}} \right)p\left( {\boldsymbol{h}} \right) $$ (16) 其中, $p\left( {{{\boldsymbol{y}}_{\exp }}|{\boldsymbol{\varTheta }}} \right)$是${{\boldsymbol{\varTheta}} }$参数条件下${{{\boldsymbol{y}}}_{\exp }}$的先验概率分布, $p\left( {{\boldsymbol{\varTheta }}|{\boldsymbol{h}}} \right)$是超参数${{\boldsymbol{h}}}$下的${{\boldsymbol{\varTheta}} }$的先验概率分布, $p\left( {\boldsymbol{h}} \right)$则为超参数${{\boldsymbol{h}}}$的先验概率分布.
评估实验响应与数值模拟响应差异需要采用合适的灵敏度函数, 评估实验响应与数值模拟响应两个向量的差距. 这里采用2-范数的平方$\left\| {{{{\boldsymbol{y}}}_{\exp }} - \psi \left( {{\boldsymbol{\varTheta }}} \right)} \right\|_2^2$进行评估. 为保证采样方式的一致性, 这里将各试样破坏前的应力-应变曲线等间隔取100个数据点. 因此, 对应于式(16)的似然函数具体形式为
$$ p\left( {{{\boldsymbol{y}}_{\exp }}|{\boldsymbol{\varTheta }}} \right) = \exp \left[ { - \frac{1}{2}{{\boldsymbol{\varGamma }}^{ - \tfrac{1}{2}}}\sum\limits_{i = 1}^N {\sum\limits_{k = 1}^K {{{\left( {F_{k,i}^{\exp } - F_{k,i}^{{\mathrm{sim}}}} \right)}^2}} } } \right] $$ (17) 其中, $F_{k,i}^{\exp }$和$F_{k,i}^{{\mathrm{sim}}}$分别为等间隔采样的实验数据和仿真数据, $i$表示此时的数据为第$i$个试样的数据, $k$表示此时的数据点为第$k$个数据点的数据.
考虑$p\left( {{\boldsymbol{\varTheta }}|{\boldsymbol{h}}} \right)$待求解参数空间包含5个参数, 适合建立5层贝叶斯模型进行计算, 即
$$ p\left( {{\boldsymbol{\varTheta }}|{\boldsymbol{h}}} \right) = \mathop \prod \limits_{i = 1}^N \mathop \prod \limits_{j = 1}^5 {p_i}\left( {{{\varTheta }_{j,i}}|{\theta _j},\sigma _j^2} \right) $$ (18) 其中, $\left[ {\left( {{\theta _1},{\sigma _1}} \right),\left( {{\theta _2},{\sigma _2}} \right),\left( {{\theta _3},{\sigma _3}} \right),\left( {{\theta _4},{\sigma _4}} \right),\left( {{\theta _5},{\sigma _5}} \right)} \right]$表示5个待定参数高斯分布的均值和标准差, 也就是贝叶斯模型需要确定的超参数. $ {{\varTheta }_{j,i}} $为第$ i $个试样的第$ j $层模型的参数.
对应地, 超参数的先验分布模型为
$$ p\left( {\boldsymbol{h}} \right) = \mathop \prod \limits_{i = 1}^5 {p_i}\left( {{\theta _i}|\hat \theta ,{{\hat \sigma }^2}} \right){p_i}\left( {{\sigma _i}|\hat \theta ,{{\hat \sigma }^2}} \right) $$ (19) 其中, $\hat \theta $和$\hat \sigma $是参数$ \theta $和$\sigma $的高斯分布的均值和期望.
基于超参数中的均值${\theta _1} \sim {\theta _5}$, 均值的不确定性行为为模型参数的合理估计提供了依据, 可用于模拟实验响应. 此外, ${\theta _1} \sim {\theta _5}$和${\sigma _1} \sim {\sigma _5}$解释了模型参数${{\boldsymbol{\varTheta}} }$的分布, 这些参数决定了由噪声数据表征的${{\boldsymbol{\varTheta}} }$的可信区间.
本模型采用Hamiltonian Monte Carlo (HMC)方法[45], 利用随机抽样过程探索后验分布的搜索空间. 此处, 建立5条马尔可夫链, 用于估计超参数的分布. 每条链的采样数为5500个, 为了减少初始采样阶段的干扰, 前500个采样结果被废弃. 在MCMC方法对后验概率空间探索和采样的过程中, 代理模型的计算开销是巨大的. 通常需要根据采样生成的参数组合, 带入正问题模型获取对应的响应结果$y$. 有限元分析能够获取可靠的计算结果, 但本问题面临多个参数的反演, 参数组合多, 大量计算有限元的成本高, 为完成一个多链MCMC模型的计算, 往往需要上万次的有限元计算. 因此, 将贝叶斯方法识别推广至三维模型的参数识别就需要高效率、鲁棒性强的计算模型.
此处, 采用基于深度神经网络DNN的快速预测模型取代基于有限元的三维模型正问题计算, 通过有限数量的数值计算结果, 利用参数空间${{\boldsymbol{\theta}} }$以及结构形态学参数$\rho $和${{\boldsymbol{M}}}$, 训练参数空间到应力-应变空间的映射关系. 并将训练好的正问题快速预测模型作为计算模块, 用以提供快速准确的响应预测值${{f}^{{\mathrm{surrogate}}}}$, 显著降低计算开销. 此处, 正问题代理模型使用6层隐藏层的MLP神经网络(256, 256, 256, 256, 128, 64), 激活函数采用Leaky ReLU, 并通过PyTorch搭建神经网络框架
2.3 参数反演结果的验证与分析
通过2.2节所述计算步骤, 可获得材料参数${{\boldsymbol{\theta}} }$的高斯分布, 如图10所示.
根据采样分布, 可以拟合出高斯分布中每个模型参数对应的超参数组合. 5个弹性本构模型参数的期望值是基于当前实验观测的最优估计, 可用于估计松质骨的各向异性力学性能. 在当前超参数设置下, 基于实验结果所获得的参数最佳估计, 即${\lambda _0} = 4089{\text{ MPa}}$, ${\mu _0} = 2102{\text{ MPa}}$, ${\lambda '_0} = 1633{\text{ MPa}}$, $k = 3.36$和$l = 1.19$, 可用于精准预测松质骨的弹性力学响应. 对于5个参数$\left[ {{\lambda _0},{\mu _0},{{\lambda '}_0},k,l} \right]$, 以$\sigma /\theta $表示标准差$\sigma $相对于均值$\theta $的大小, 5个参数的$\sigma /\theta $分别为2.06%, 3.43%, 5.57%, 11.01%和13.03%.
基于本方法进行校正后的本构模型参数和使用传统试样方法拟合的本构模型参数[46]的对比结果, 如图11所示.
![]()
图 11 基于贝叶斯反演的松质骨力学性能表征与传统方法的对比Figure 11. The comparison between proposed calibrated method and traditional method based on Bayesian inversion针对最常用的弹性响应阶段, 利用所确定的本构模型参数建立松质骨试样的均匀化有限元模型, 施加与实验相同的载荷边界, 获取力-位移响应信息. 提取模拟结果以及实验结果计算弹性响应阶段的结构刚度${K^{{\mathrm{FEM}}}}$以及${K^{{\mathrm{Exp}}}}$, 统计分析两者的相对误差以及相关系数. 从图11(a)发现, 对用于贝叶斯推断计算的240个试样进行统计分析, 相关性系数为0.957, 均方差RMSE = 15.00 N/mm, 最大误差为15.87%, 平均误差为7.98%. 同时, 为验证所获取参数的普遍性, 即该参数能否用于新的松质骨材料力学行为分析, 这里使用来自于全新个体且未参与贝叶斯参数校正的36个试样作为验证组, 进行模拟与实验结果的对比分析, 结果如图11(b)所示. 结果显示验证组试样所计算得到的相关性系数为0.968, 均方差为28.84 N/mm, 最大误差为16.55%, 平均误差为8.35%.
因此, 基于参数高斯分布的期望值选取的最优参数组合能够对松质骨进行高精度的力学行为分析. 所识别的参数能够反映同类型个体同年龄段相同部位的松质骨材料力学性能, 进一步说明了参数的合理性.
另外, 通过实验组和验证组的结果对比, 可以发现两组试样的最大误差和平均误差基本一致, 均方差相差较小, 说明基于统计规律所测定的模型参数能够描述绝大多数松质骨的力学性能. 对于个别预测误差较大的个体, 说明和该个体类似的试样数量较少, 未能充分体现这些个体对统计规律的影响. 这也说明, 通过增加试样的数量, 本方法能够进一步提高参数的精度, 最终实现对不同个体不同部位松质骨力学性能的预测.
基于本方法进行校正后的本构模型参数, 可进一步对完整松质骨骨组织的非均质各向异性的力学行为进行更为准确的描述. 图12展示了股骨远端松质骨在主方向上的模量空间分布. 结果表明, 股骨远端松质骨区域沿3个材料主方向的主方向模量分布呈现出明显的非均匀性和轴向优势[47-48].
![]()
图 12 基于修正模型的股骨远端松质骨力学性能Figure 12. The mechanical properties of distal femur based on proposed calibrated model如图12(a)所示, 定义骨干方向为${{{\boldsymbol{n}}}_1}$, 垂直髁线方向为${{{\boldsymbol{n}}}_2}$, 沿髁线方向为${{{\boldsymbol{n}}}_3}$. 参考Allen等[49]和图2(a)对羊股骨远端骨肌系统的解剖学特征, 健康羊的股骨远端主要承受上部骨干沿${{{\boldsymbol{n}}}_1}$方向的压力, 以及胫骨平台沿${{{\boldsymbol{n}}}_2}$方向对股骨后髁施加的作用力. 进一步地, 考虑骨组织往往随所承受载荷大小以及方向进行生长改建的规律, 股骨远端骨组织各向异性力学性能在骨干方向为${{{\boldsymbol{n}}}_1}$和垂直髁线方向为${{{\boldsymbol{n}}}_2}$上表现出明显的轴向优势.
图12(b)和图12(c)所示的股骨远端主方向模量幅值与主方向分布表明, 3组主方向模量${E_1}$, ${E_2}$和${E_3}$在股骨后髁的幅值高于前髁, 股骨后髁为羊股骨远端与胫骨平台的主要受力区域. 在股骨后髁区域, ${E_1}$主要沿${{{\boldsymbol{n}}}_2}$方向分布, ${E_2}$沿${{{\boldsymbol{n}}}_1}$方向分布, 且这两个方向的模量远高于主要沿${{{\boldsymbol{n}}}_3}$方向的模量${E_3}$.
综上所述, 基于本部分建立的参数识别模型所获得的围关节骨组织微结构非均质性与各向异性建模, 能够充分反映松质骨材料力学性能对整个骨结构力学响应的影响. 从整体骨结构的角度来看, 该模型反映了松质骨各向异性力学性能分布的轴向优势, 并与解剖学特征相一致. 进一步验证了所得到的松质骨力学性能符合羊股骨远端的受力特性, 揭示了复杂围关节区域的非均匀力学性能.
为进一步验证所获得修正模型对骨结构整体力学响应的作用, 本文建立了包含完整股骨远端的有限元模型.
如图13(a)所示, 该模型包含完整的股骨远端(完整股骨髁和一半的股骨骨干). 松质骨的各向异性力学性能分别由本文所得到的修正后的材料参数和传统的单轴方法得到的材料参数提供. 股骨干顶部截面施加固支边界条件. 胫骨平台采用刚性圆柱模型, 圆柱底部施加垂直于底面的125 N的均布压力(采用羊重力的${1 \mathord{\left/ {\vphantom {1 4}} \right. } 4}$). 然后, 通过股骨远端与刚性圆柱的接触行为, 将相互作用力施加到股骨远端后髁位置. 此处, 该模型只考虑在重力作用下股骨远端与胫骨平台的相互作用力, 因此两者之间的法向接触行为定义为硬接触, 切向行为定义为摩擦系数为0.2的摩擦接触.
![]()
图 13 股骨远端CT影像、有限元建模与结果分析Figure 13. The micro-CT images, FEM modeling and results of distal femur数值计算由ABAQUS/Standard模块实现. 基于骨密度和结构张量的松质骨各向异性本构模型通过ABAQUS UMAT子程序实现. 利用USDFLD和UFIELD两个子程序定义整个股骨远端模型骨密度的非均匀分布和结构张量特征值. 利用所建立的松质骨组织形态学参数的预测模型, 计算结构张量的特征向量, 并通过ORIENT子程序来定义局部材料3个主方向.
首先, 如图13(b)所示, 本文比较了修正方法与传统方法(采用模量-密度幂律关系)的Mises应力场. 可以发现, 利用本文材料参数建立的模型Mises应力值明显高于传统方法. 在线弹性材料假设下, 两个模型具有相同的几何形状和加载状态, 在弹性阶段, 两种方法的Mises应力场的分布情况相似.
其次, 考虑到股骨远端在日常生活中承受压缩负荷, 本研究选取同一矢状面截面, 对比分析两个模型最小主压应力的分布差异. 图13(a)中, 蓝色区域表示BV/TV大于0.2的解剖位置, 红色区域表示BV/TV小于0.2的解剖位置. 体积分数高的解剖部位骨小梁更多, 意味着该区域承载更大的力, 而体积分数低的解剖部位骨组织往往为离散的骨小梁, 该区域几乎不能受力. 为评估所建立有限元模型的合理性和可靠性, 本文分析了压应力分布和BV/TV分布的一致性. 图13(b)中黑色虚线区域表示小应力区域(最小主应力幅值小于0.1 MPa). 在本文材料参数建立的模型中, 黑色虚线内椭圆区域的位置与低BV/TV区域基本一致, 而高BV/TV区域则具有更高的应力. 不同区域的主应力方向与骨小梁方向基本一致. 考虑到材料分布的非均质性及空间微结构变化引起的材料主方向变化, 本文所建立的模型显示应力分布与空间微结构分布基本一致. 该结果也表明, 松质骨的非均匀生长和分布与自然状态下整个骨组织内部的非均匀应力分布密切相关.
作为对比, 采用传统方法所得到参数建立的有限元模型, 局部主应力方向并不总是与小梁方向一致. 例如, 在股骨髁前方, 局部主应力方向与小梁方向无明显相关性. 其次, 可以观察到黑色虚线中矩形区域的位置与低BV/TV区域的位置并不一致. 主应力的幅值分布与BV/TV的分布也不一致. 低BV/TV区域仍承受较大的压应力.
综上所述, 准确的非均质各向异性模型参数对整个股骨远端模型的力学响应分析能够产生显著的影响.
3. 结 论
随着精准医疗概念在骨缺损治疗领域的不断发展, 骨缺损患者的个性化治疗已经成为相关领域的研究热点. 然而, 当前临床条件下, 无创、快速及精准的个性化骨力学分析仍旧是挑战性问题. 本文旨在建立融合数据驱动与力学建模的骨缺损重建理论、计算和实验方法, 以解决体内骨组织微观结构分析、非均质各向异性力学行为表征和建模方面的诸多问题, 实现临床条件下骨组织力学性能的准确表征.
首先, 针对当前骨缺损临床诊疗中临床骨肌影像系统无法准确评估松质骨非均质各向异性力学性能的问题, 本文提出基于数据驱动与力学模型融合的建模方法, 通过micro-CT与临床CT影像关联的深度学习模型, 预测骨组织的形态学特征与结构力学性能. 使用羊股骨远端松质骨CT影像数据集验证了该模型的表征能力, 实现了仅通过临床CT影像对非均匀骨组织形态学参数的精准预测. 本方法降低了对临床CT影像分辨率的要求, 使其应用于高精度的骨组织形态学参数与力学行为表征成为可能, 解决了骨科临床中难以表征骨组织结构非均质性和各向异性的难题.
其次, 结合临床CT下对松质骨微结构特征规律表征的数据驱动模型, 本文发展了基于非均匀骨密度与结构张量的松质骨各向异性参数化本构模型. 通过贝叶斯修正传统松质骨实验方法中因材料主方向与加载方向偏差引入的系统误差. 通过羊股骨远端松质骨实验验证所提出本构模型与参数反演方法的准确性, 揭示不同部位松质骨力学行为与微结构生长方向的关系. 此外, 基于有限元数据集建立的DNN代理模型的使用显著提高了MCMC方法的迭代效率, 实现了三维问题下的贝叶斯多参数反演辨识, 并通过实验和计算表征了松质骨各向异性力学性能, 揭示了骨力学行为与微结构生长方向的内在关系.
当前临床骨缺损诊疗中, 个性化人工骨结构的设计需要首先考虑骨的力学性能. 然而, 现有方法缺乏基于体内环境的精确力学设计目标, 未充分考虑骨结构各向异性和力学非均质性, 导致假体结构的力学性能与自然骨不匹配, 进而引起应力屏蔽效应和假体-自然骨系统的不稳定. 本研究结合松质骨, 建立了临床条件下骨组织力学性能精准表征的理论、实验和计算方法, 有望更准确地实现人体骨力学性能的临床表征, 指导骨缺损的临床治疗.
此外, 本文所提出的骨力学性能研究方法并未针对特定的骨组织结构, 因此具有较强的可扩展性与普适性. 通过建立包含脊柱、髋关节、股骨远端及胫骨平台等常见的松质骨区域的影像数据集, 可进一步训练所建立的神经网络模型, 使其精准预测不同部位的骨组织的组织形态学参数. 结合生物力学实验, 可获得不同部位的骨组织的本构模型参数, 精准预测骨组织各向异性等效力学性能.
目前, 本研究中尚存在部分问题有待未来进一步开展: 首先, 本文工作主要围绕微观组织复杂、非均质性与各向异性较强的松质骨假体设计开展, 尚未对力学性能较为均匀的皮质骨进行分析. 其次, 目前的研究材料均来自成羊股骨髁的新鲜尸体骨, 尚未评估不同部位松质骨是否具有相同的等效力学参数. 未来工作中将尝试引入不同部位的骨材料, 评估其力学性能差异. 最后, 考虑到动物和人体骨组织在骨小梁微观结构分布和特征尺寸上的区别, 未来研究将采用人体尸体骨, 研究实验动物和人体骨组织结构的尺度率关系, 从而将研究成果应用于临床, 评估患者骨缺损区的力学性能.
-
![]()
图 1 临床影像下典型的骨缺损形式及其数字建模
Figure 1. Typical form of bone defect and its digital modeling under clincial imaging
![]()
![]()
图 5 正交各向异性弹性刚度矩阵预测结果
Figure 5. The prediction results of orthotropic elastic stiffness matrix
![]()
图 6 Micro-CT有限元计算结果与临床低分辨率CT预测结果的对比
Figure 6. Comparison of finite element calculation results of micro-CT and clinical low resolution CT prediction results
![]()
图 8 不同应力状态下的松质骨名义应力-名义应变曲线
Figure 8. Nominal stress-nominal strain curves of cancellous bone under different stress states
![]()
图 9 基于贝叶斯推断的材料反演整体框架
Figure 9. The whole framework of material parameters identification based on Bayesian inference
![]()
图 11 基于贝叶斯反演的松质骨力学性能表征与传统方法的对比
Figure 11. The comparison between proposed calibrated method and traditional method based on Bayesian inversion
![]()
图 12 基于修正模型的股骨远端松质骨力学性能
Figure 12. The mechanical properties of distal femur based on proposed calibrated model
![]()
图 13 股骨远端CT影像、有限元建模与结果分析
Figure 13. The micro-CT images, FEM modeling and results of distal femur
表 1 训练样本集的4个组织形态学参数分布范围
Table 1 The distribution range of four histomorphological parameters of the training sample set
Histomorphological parameters Distribution range BV/TV $0.48 \pm 0.21$ DA $0.51 \pm 0.13$ Conn.D $4.55 \pm 1.76$ BS/TV $3.93 \pm 0.49$ 
下载: 导出CSV
-
[1] 乔爱科. 生物力学的两个核心问题. 力学与实践, 2020, 42(1): 119-123 (Qiao Aike. Two core issues of biomechanics. Mechanics in Engineering, 2020, 42(1): 119-123 (in Chinese) Qiao Aike. Two core issues of biomechanics. Mechanics in Engineering, 2020, 42(1): 119-123 (in Chinese)
[2] Koons GL, Diba M, Mikos AG. Materials design for bone-tissue engineering. Nature Reviews Materials, 2020, 5(8): 584-603 doi: 10.1038/s41578-020-0204-2
[3] Zhang M, Lin R, Wang X, et al. 3D printing of Haversian bone-mimicking scaffolds for multicellular delivery in bone regeneration. Science Advances, 2020, 6(12): eaaz6725 doi: 10.1126/sciadv.aaz6725
[4] Malachanne E, Dureisseix D, Cañadas P, et al. Experimental and numerical identification of cortical bone permeability. Journal of Biomechanics, 2008, 41(3): 721-725 doi: 10.1016/j.jbiomech.2007.09.028
[5] Lu L, Jin P, Pang G, et al. Learning nonlinear operators via DeepONet based on the universal approximation theorem of operators. Nature Machine Intelligence, 2021, 3(3): 218-229 doi: 10.1038/s42256-021-00302-5
[6] Liu X, Tian S, Tao F, et al. A review of artificial neural networks in the constitutive modeling of composite materials. Composites Part B: Engineering, 2021, 224: 109152 doi: 10.1016/j.compositesb.2021.109152
[7] Yan W, Lin S, Kafka OL, et al. Data-driven multi-scale multi-physics models to derive process–structure–property relationships for additive manufacturing. Computational Mechanics, 2018, 61(5): 521-541 doi: 10.1007/s00466-018-1539-z
[8] Liu B, Kovachki N, Li Z, et al. A learning-based multiscale method and its application to inelastic impact problems. Journal of the Mechanics and Physics of Solids, 2022, 158: 104668 doi: 10.1016/j.jmps.2021.104668
[9] Shukla K, Jagtap AD, Blackshire JL, et al. A physics-informed neural network for quantifying the microstructural properties of polycrystalline nickel using ultrasound data: A promising approach for solving inverse problems. IEEE Signal Processing Magazine, 2021, 39(1): 68-77
[10] Lookman T, Balachandran PV, Xue D, et al. Active learning in materials science with emphasis on adaptive sampling using uncertainties for targeted design. NPJ Computational Materials, 2019, 5(1): 21
[11] 李想, 严子铭, 柳占立. 机器学习与计算力学的结合及应用初探. 科学通报, 2019, 64(7): 635-648 (Li Xiang, Yan Ziming, Liu Zhanli. Combination and application of machine learning and computational mechanics. Chinese Science Bulletin, 2019, 64(7): 635-648 (in Chinese) Li Xiang, Yan Ziming, Liu Zhanli. Combination and application of machine learning and computational mechanics. Chinese Science Bulletin, 2019, 64(7): 635-648 (in Chinese)
[12] 李想, 严子铭, 柳占立等. 基于仿真和数据驱动的先进结构材料设计. 力学进展, 2021, 51(1): 82-105 (Li Xiang, Yan Ziming, Liu Zhanli, et al. Advanced structural material design based on simulation and data-driven method. Advances in Mechanics, 2021, 51(1): 82-105 (in Chinese) doi: 10.6052/1000-0992-20-012 Li Xiang, Yan Ziming, Liu Zhanli, et al. Advanced structural material design based on simulation and data-driven method. Advances in Mechanics, 2021, 51(1): 82-105 (in Chinese) doi: 10.6052/1000-0992-20-012
[13] Li X, Liu Z, Cui S, et al. Predicting the effective mechanical property of heterogeneous materials by image based modeling and deep learning. Computer Methods in Applied Mechanics and Engineering, 2019, 347: 735-753 doi: 10.1016/j.cma.2019.01.005
[14] Li X, Ning S, Liu Z, et al. Designing phononic crystal with anticipated band gap through a deep learning based data-driven method. Computer Methods in Applied Mechanics and Engineering, 2020, 361: 112737 doi: 10.1016/j.cma.2019.112737
[15] 庄茁, 严子铭, 姚凯丽等. 固体力学跨尺度计算若干问题研究. 计算力学学报, 2024, 41(1): 40-46 (Zhuang Zhuo, Yan Ziming, Yao Kaili, et al. Several problem studies in solid mechanics by solid mechanics by spanning the scale computation analysis. Chinese Journal of Computational Mechanics, 2024, 41(1): 40-46 (in Chinese) doi: 10.7511/jslx20230909003 Zhuang Zhuo, Yan Ziming, Yao Kaili, et al. Several problem studies in solid mechanics by solid mechanics by spanning the scale computation analysis. Chinese Journal of Computational Mechanics, 2024, 41(1): 40-46 (in Chinese) doi: 10.7511/jslx20230909003
[16] Wu J, Zhang C, Xue T, et al. Learning a probabilistic latent space of object shapes via 3D generative-adversarial modeling//Proceedings of the 30th International Conference on Neural Information Processing Systems, 2016: 82-90
[17] Wang L, Chan YC, Ahmed F, et al. Deep generative modeling for mechanistic-based learning and design of metamaterial systems. Computer Methods in Applied Mechanics and Engineering, 2020, 372: 113377 doi: 10.1016/j.cma.2020.113377
[18] Dan Y, Zhao Y, Li X, et al. Generative adversarial networks (GAN) based efficient sampling of chemical composition space for inverse design of inorganic materials. NPJ Computational Materials, 2020, 6(1): 84 doi: 10.1038/s41524-020-00352-0
[19] Sundar S, Sundararaghavan V. Database development and exploration of process–microstructure relationships using variational autoencoders. Materials Today Communications, 2020, 25: 101201 doi: 10.1016/j.mtcomm.2020.101201
[20] Fuglede B, Topsoe F. Jensen-Shannon divergence and Hilbert space embedding BT-IEEE International Symposium on Information Theory, 2004, 31
[21] Jin X, Wu L, Li X, et al. Predicting aesthetic score distribution through cumulative jensen-shannon divergence. arXiv, 2017: 77-84
[22] Yang L, Meng X, Karniadakis GE. B-PINNs: Bayesian physics-informed neural networks for forward and inverse PDE problems with noisy data. Journal of Computational Physics, 2021, 425: 109913 doi: 10.1016/j.jcp.2020.109913
[23] Li Q, Xu R, Wu Q, et al. Topology optimization design of quasi-periodic cellular structures based on erode–dilate operators. Computer Methods in Applied Mechanics and Engineering, 2021, 377: 113720 doi: 10.1016/j.cma.2021.113720
[24] Xu S, Shen J, Zhou S, et al. Design of lattice structures with controlled anisotropy. Materials and Design, 2016, 93: 443-447 doi: 10.1016/j.matdes.2016.01.007
[25] Dong G, Tang Y, Zhao YF. A 149 line homogenization code for three-dimensional cellular materials written in MATLAB. Journal of Engineering Materials and Technology, 2019, 141(1): 011005 doi: 10.1115/1.4040555
[26] Zheng L, Kumar S, Kochmann DM. Data-driven topology optimization of spinodoid metamaterials with seamlessly tunable anisotropy. Computer Methods in Applied Mechanics and Engineering, 2021, 383: 113894 doi: 10.1016/j.cma.2021.113894
[27] Kumar S, Tan S, Zheng L, et al. Inverse-designed spinodoid metamaterials. NPJ Computational Materials, 2020, 6(1): 73 doi: 10.1038/s41524-020-0341-6
[28] Mittra E, Rubin C, Qin YX. Interrelationship of trabecular mechanical and microstructural properties in sheep trabecular bone. Journal of Biomechanics, 2005, 38(6): 1229-1237 doi: 10.1016/j.jbiomech.2004.06.007
[29] Maquer G, Musy SN, Wandel J, et al. Bone volume fraction and fabric anisotropy are better determinants of trabecular bone stiffness than other morphological variables. Journal of Bone and Mineral Research, 2015, 30(6): 1000-1008 doi: 10.1002/jbmr.2437
[30] Xiao P, Zhang T, Dong XN, et al. Prediction of trabecular bone architectural features by deep learning models using simulated DXA images. Bone Reports, 2020, 13: 100295 doi: 10.1016/j.bonr.2020.100295
[31] Callens SJP, Tourolle né Betts DC, Müller R, et al. The local and global geometry of trabecular bone. Acta Biomaterialia, 2021, 130: 343-361 doi: 10.1016/j.actbio.2021.06.013
[32] Steiner L, Synek A, Pahr DH. Bone reports comparison of different microct-based morphology assessment tools using human trabecular bone. Bone Reports, 2020, 12: 100261 doi: 10.1016/j.bonr.2020.100261
[33] Mickel W, Kapfer SC, Schröder-Turk GE, et al. Shortcomings of the bond orientational order parameters for the analysis of disordered particulate matter. Journal of Chemical Physics, 2013, 138(4): 044501
[34] Schröder-Turk GE, Mickel W, Kapfer SC, et al. Minkowski tensor shape analysis of cellular, granular and porous structures. Advanced Materials, 2011, 23(22-23): 2535-2553 doi: 10.1002/adma.201100562
[35] Zysset PK, Curnier A. An alternative model for anisotropic elasticity based on fabric tensors. Mechanics of Materials, 1995, 21(4): 243-250 doi: 10.1016/0167-6636(95)00018-6
[36] Zysset PK. A review of morphology-elasticity relationships in human trabecular bone: Theories and experiments. Journal of Biomechanics, 2003, 36(10): 1469-1485 doi: 10.1016/S0021-9290(03)00128-3
[37] Kirby M, Morshed AH, Gomez J, et al. Three-dimensional rendering of trabecular bone microarchitecture using a probabilistic approach. Biomechanics and Modeling in Mechanobiology, 2020, 19(4): 1263-1281 doi: 10.1007/s10237-020-01286-8
[38] Yan Z, Hu Y, Li X, et al. Data-driven based characterization of anisotropic mechanical properties for cancellous bone from low-resolution CT images. IEEE Transactions on Biomedical Engineering, 2024, 71(2): 689-699 doi: 10.1109/TBME.2023.3315846
[39] Yan Z, Hu Y, Shi H, et al. Experimentally characterizing the spatially varying anisotropic mechanical property of cancellous bone via a Bayesian calibration method. Journal of the Mechanical Behavior of Biomedical Materials, 2023, 138: 105643 doi: 10.1016/j.jmbbm.2022.105643
[40] Rincón-Kohli L, Zysset PK. Multi-axial mechanical properties of human trabecular bone. Biomechanics and Modeling in Mechanobiology, 2009, 8(3): 195-208 doi: 10.1007/s10237-008-0128-z
[41] Ovesy M, Voumard B, Zysset P. A nonlinear homogenized finite element analysis of the primary stability of the bone–implant interface. Biomechanics and Modeling in Mechanobiology, 2018, 17(5): 1471-1480 doi: 10.1007/s10237-018-1038-3
[42] Charlebois M, Jirásek M, Zysset PK. A nonlocal constitutive model for trabecular bone softening in compression. Biomechanics and Modeling in Mechanobiology, 2010, 9(5): 597-611 doi: 10.1007/s10237-010-0200-3
[43] Meng X, Yang L, Mao Z, et al. Learning Functional Priors and Posteriors from Data and Physics. Journal of Computational Physics, 2022, 457: 111073 doi: 10.1016/j.jcp.2022.111073
[44] Nguyen T, Francom DC, Luscher DJ, et al. Bayesian calibration of a physics-based crystal plasticity and damage model. Journal of the Mechanics and Physics of Solids, 2021, 153: 104284
[45] Brunton SL, Proctor JL, Kutz JN, et al. Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences of the United States of America, 2016, 113(15): 3932-3937
[46] Munford MJ, Ng KCG, Jeffers JRT. Mapping the multi-directional mechanical properties of bone in the proximal tibia. Advanced Functional Materials, 2020, 30(46): 2004323 doi: 10.1002/adfm.202004323
[47] Liu P, Liang X, Li Z, et al. Decoupled effects of bone mass, microarchitecture and tissue property on the mechanical deterioration of osteoporotic bones. Composites Part B: Engineering, 2019, 177: 107436 doi: 10.1016/j.compositesb.2019.107436
[48] Li Z, Liu P, Yuan Y, et al. Loss of longitudinal superiority marks the microarchitecture deterioration of osteoporotic cancellous bones. Biomechanics and Modeling in Mechanobiology, 2021, 20(5): 2013-2030 doi: 10.1007/s10237-021-01491-z
[49] Allen MJ, Houlton JE, Adams SB, et al. The surgical anatomy of the stifle joint in sheep. Veterinary Surgery, 1998, 27: 596-605 doi: 10.1111/j.1532-950X.1998.tb00536.x
计量
- 文章访问数: 217
- HTML全文浏览量: 81
- PDF下载量: 119
