面向椎体成形术的骨水泥注射过程可视化仿真

Visual Simulation of Bone Cement Injection Process for Vertebroplasty

摘要

在骨科经皮椎体成形手术中，骨水泥注射核心参数（如注射量、压力及角度）的确定仍依赖经验判断，而这些参数对防止渗漏和并发症至关重要.为此，本文提出一种集成多相非牛顿流体模型的先进仿真框架.该框架采用基于Voronoi图的仿生骨表征方法，用于预测骨小梁结构内的水泥流动特性，并支持不同工况下水泥弥散情况的高保真可视化分析.实验结果表明，本仿真框架在术前规划优化、手术培训辅助及渗漏风险评估方面展现出潜在应用价值：术前可辅助医护人员预演注射过程以优化操作方案；在培训场景中，能模拟不同参数下的水泥弥散效果，为新手医师熟悉操作逻辑提供参考；同时，其对渗漏风险的预测功能，也有望为提升临床手术安全性提供技术层面的支持.

Abstract

Determining key injection parameters of bone cement， such as volume， pressure and angle， in percutaneous vertebroplasty still largely relies on empirical judgment， and they have a critical role in preventing leakage and related complications. An advanced simulation framework that integrated a multiphase non-Newtonian fluid model was introduced. The framework employed a Voronoi-based biomimetic bone representation to predict cement flow behavior within trabecular structures， enabling high-fidelity visualization and analysis of cement dispersion under varying operating conditions. Experimental results demonstrate the framework’s potential value in preoperative planning optimization， surgical training assistance， and leakage risk assessment. For preoperative use， it can assist clinicians in rehearsing injection procedures and refining operation strategies； in training scenarios， it can simulate cement dispersion patterns under different parameter settings， thereby helping novice surgeons understand procedural logic； furthermore， its capability to predict leakage risks holds promise for providing technical support to improve clinical surgical safety.

Graphical abstract

关键词

医学可视化 / 骨填充模拟 / 多相非牛顿流体建模 / 骨水泥

Key words

medical visualization / bone filling simulation / multiphase non-Newtonian fluid modeling / bone cement

全球人口老龄化导致骨质疏松症及骨折发病率上升，给医疗系统带来沉重负担，目前影响约2亿人.骨水泥增强技术，特别是经皮椎体成形术等微创手术^［1-2］，对于恢复椎体结构完整性具有关键作用（见图1）.

然而，骨水泥的制备与注射仍面临技术挑战.临床医师通常依据指南和个人经验确定混合比例及注射参数，这易导致水泥分布欠佳，并引发渗漏、肺栓塞及神经功能缺损等并发症^［3-4］.多组骨水泥填充案例直观呈现了这些风险（见图2）.传统培训方法成本高昂且耗时，亟须开发高效的手术模拟工具.

针对上述挑战，本研究提出一种骨水泥注射仿真框架，整合了骨骼建模、多相非牛顿流体模拟和高保真渲染技术.本文方法采用基于体素和粒子采样的骨骼建模技术，结合基于构象张量的隐式混合模型进行流体模拟，并通过真实感渲染实现直观可视化.与现有方法^［5］相比，本方法在模型精度、工作流程标准化和材料色彩映射等方面均有显著提升，为手术规划和培训提供了一种经济高效的解决方案.

1 相关工作简介

1.1 骨水泥材料研究

骨水泥在骨科关节置换及骨质疏松治疗中得到应用.自首次将其应用于股骨假体固定以来^［6］，聚甲基丙烯酸甲酯（PMMA）已成为应用最广泛的骨水泥材料.该材料由PMMA粉末与甲基丙烯酸甲酯单体组成，混合后发生聚合反应.PMMA具有良好的生物相容性，并可通过化学修饰调整其性能参数.

磷酸钙骨水泥（CPC）^［1］​是另一种常用生物材料，其主要成分为羟基磷灰石（HAP）.由于与骨矿物质成分相似，CPC表现出优异的生物相容性和生物活性，可逐渐被吸收并由新生骨组织替代.该材料已广泛应用于骨再生领域，包括颌面缺损修复^［7］及生物活性复合材料^［8］等方向.

1.2 计算机模拟技术

作为植入人体的器械类产品，不可能进行实体实验，因此，基于物理原理开展的不同条件下骨水泥注射过程（如注射角度、速度及成分）仿真技术成为临床专业人员和手术培训的重要替代方案.

流体模拟作为计算机图形学的重要领域，目前已发展出多种求解器来实现对压力、黏度、边界和细节的处理功能.拉格朗日方法（如光滑粒子流体动力学（SPH）^［9］和基于位置的动力学（PBD）^［10］）提供了稳定的压力求解方案，包括无散SPH（DFSPH）^［11］和基于位置的流体（PBF）^［12］.欧拉方法则包含格子玻尔兹曼方法（LBM）^［13］和有限体积法（FVM）^［14］等基于网格的求解框架.

现有研究主要集中于单相牛顿流体模拟，已实现高精度和视觉逼真的效果^［15-17］.相比之下，多相流模拟研究仍较为有限，主要集中在多相流框架构建^［18-19］、气泡动力学^［20］和相变模拟^［21］等方面.非牛顿流体的黏度会随剪切速率非线性变化，包括剪切稀化和剪切增稠2种类型.骨水泥作为一种具有剪切稀化特性的多相非牛顿流体，其模拟需要特定的技术方法，而目前该领域的研究仍显不足.

2 仿真与可视化框架

在骨水泥手术模拟中，必须同时考虑松质骨建模和骨水泥流体模拟算法.下面将详细介绍本文所提的研究方法.

2.1 骨建模

椎体建模方法通过分别表征皮质骨与松质骨结构，将模型转化为基于粒子的格式以实现流固耦合，并通过整合Voronoi纹理^［22］、行进立方体算法^［23］及泊松圆盘采样^［24］生成均匀粒子分布用于模拟，构建流程如图3所示.

2.2 基于SPH方法的流体模拟简介

SPH方法将连续流体介质离散化为拉格朗日粒子群，作为一种无网格计算框架，实现了空间场及其微分算子的离散化近似.在该方法中，流体密度、质量、速度、压力等物理场变量均以粒子为载体进行定义，从而实现对流体运动的完全拉格朗日描述.

SPH方法通过核函数加权的邻近粒子贡献来计算每个粒子在下一时间步的状态演化（见图3）.具体而言，物理场A_i（在粒子i处采样）可表示为

A_i=\sum _{j\in N\left(i\right)}\frac{m_j}{{\rho }_j}{A}_j{W}_{ij}.(1)

其中：N\left(i\right)是影响粒子i的邻域；m是粒子质量；\rho \mathrm{是}粒子密度；W是类高斯核函数（例如3次样条核函数）；A_j是目标物理场，如速度场、应力场等.

进一步地，利用SPH标准离散化方法估计物理场A的梯度、散度和拉普拉斯算子，公式如下：

\left.\begin{array}{l}\nabla A_i=\sum _{j\in N\left(i\right)}\frac{m_j}{{\rho }_j}{A}_j\otimes \nabla W_{ij},\\ \nabla \cdot A_i=\sum _{j\in N\left(i\right)}\frac{m_j}{{\rho }_j}{A}_j\nabla W_{ij},\\ {\nabla }^{\mathrm{2}}{A}_i=\sum _{j\in N\left(i\right)}\frac{m_j}{{\rho }_j}{A}_j{\nabla }^{\mathrm{2}}{W}_{ij}.\end{array}\right\}(2)

其中：\mathit{a}\otimes \mathit{b}=\mathit{a}{\mathit{b}}^{\mathrm{T}}；\nabla²为拉普拉斯算子；\nabla为那勃勒算子.

2.3 改进的隐式混合模型相态传递算法

隐式混合模型是混合模型的一种变体，能够实现多相流体的通用模拟，通过一系列参数描述其流动特性.相浓度由相体积分数表示.对于粒子i的所有k个相，其相体积分数{\phi }_{i,k}之和满足归一化条件：

\sum _k{\phi }_{i,k}=\mathrm{1}.(3)

本文及后续内容采用A_{i,k}表示粒子i中与相k相关的物理场A.混合物速度场{\mathit{v}}_{i,\mathrm{m}\mathrm{i}\mathrm{x}}通过相速度场{\mathit{v}}_{i,k}重构得到，具体表达式为

{\mathit{v}}_{i,\mathrm{m}\mathrm{i}\mathrm{x}}=\sum _k{\phi }_{i,k}{\mathit{v}}_{i,k}.(4)

混合粒子的密度{\rho }_{i,\mathrm{m}\mathrm{i}\mathrm{x}}的计算公式为

{\rho }_{i,\mathrm{m}\mathrm{i}\mathrm{x}}=\sum _k{\phi }_{i,k}{\rho }_k^{\mathrm{0}}.(5)

其中，{\rho }_k^{\mathrm{0}}为第k相的静止密度.

在隐式混合模型^［25］中，相级的压力、黏性力及外力的计算公式如下：

\begin{array}{c}\frac{D{\mathit{v}}_{i,k}^p}{Dt}=\frac{{\mathit{M}}_{i,\mathrm{m}\mathrm{i}\mathrm{x}}^p}{{\rho }_{i,\mathrm{m}\mathrm{i}\mathrm{x}}}\left(C_{\mathrm{d}}+\left(\mathrm{1}-C_{\mathrm{d}}\right)\frac{{\rho }_{i,\mathrm{m}\mathrm{i}\mathrm{x}}}{{\rho }_k^{\mathrm{0}}}\right),\end{array}(6)

\frac{D{\mathit{v}}_{i,k}^{\mathrm{e}\mathrm{x}\mathrm{t}}}{Dt}=\mathit{g}.(7)

其中：\frac{D}{Dt}\mathrm{表示}\mathrm{物质}\mathrm{导数};\frac{D\mathit{v}}{Dt}表示由不同力产生的加速度，上标p和\mathrm{e}\mathrm{x}\mathrm{t}分别代表压力和外力； M 表示动量源项；C_{\mathrm{d}}∈［0，1］是隐式混合模型的调节参数，用于控制各相对混合体的影响程度； g 表示重力加速度.

经典相态传递由漂移速度{\mathit{v}}_{i,k}^{\mathrm{d}}和扩散系数定义，其中漂移速度表达式为

{\mathit{v}}_{i,k}^{\mathrm{d}}={\mathit{v}}_{i,k}-{\mathit{v}}_{i,\mathrm{m}\mathrm{i}\mathrm{x}}.(8)

相体积分数的变化由以下2个因素共同决定：

\left.\begin{array}{l}\frac{D{\phi }_{i,k}}{Dt}=-\sum _{j\in N\left(i\right)}{V}_{\mathrm{0}}\left({\phi }_{i,k}{\mathit{v}}_{i,k}^{\mathrm{d}}+{\phi }_{j,k}{\mathit{v}}_{j,k}^{\mathrm{d}}\right)\nabla \cdot W_{ij},\\ {\nabla }^{\mathrm{2}}{\phi }_{i,k}=C_{\mathrm{f}}\sum _{j\in N\left(i\right)}\left({\phi }_{i,k}-{\phi }_{j,k}\right)\frac{{\mathit{x}}_{ij}\cdot \nabla W_{ij}}{\Vert {\mathit{x}}_{ij}{\Vert }^{\mathrm{2}}+ϵ}.\end{array}\right\}(9)

其中：C_{\mathrm{f}}为扩散系数；V_{\mathrm{0}}为粒子的静止体积；\mathit{x}为粒子位置，{\mathit{x}}_{ij}={\mathit{x}}_i-{\mathit{x}}_j；ϵ为小的正则化常数（例如ϵ≃0.001）.在原隐式混合模型中，C_{\mathrm{d}}∈［0，1］是用于调节相漂移效率的因子.C_{\mathrm{d}}越大，相态传递越困难.在原始方法^［25］中C_{\mathrm{d}}为常数，假设影响恒定，这对骨水泥或玉米淀粉水溶液等反应性溶液并不准确.因此，需要一种机制来计算溶质浓度对相态传递的影响.

为简化分子间相互作用的复杂计算，本文提出一种键合效应网络，基于溶质浓度动态调整C_{\mathrm{d}}.为每种混合物设置基础值C_{\mathrm{d}}^{\mathrm{0}}，并采用SPH方法估算当前动态C_{\mathrm{d}}.

针对骨水泥应用，本文聚焦于两相流体，其中{\phi }_{k\mathrm{1}}代表液相体积分数，{\phi }_{k\mathrm{2}}代表聚合物相体积分数.

C_{\mathrm{d}}=C_{\mathrm{d}}^{\mathrm{0}}+\left(\mathrm{1}-C_{\mathrm{d}}^{\mathrm{0}}\right)\sum _{j\ne i}{V}_{\mathrm{0}}{\phi }_{j,k_{\mathrm{2}}}{W}_{ij}.(10)

当粒子周围溶质浓度较高时，C_{\mathrm{d}}→1，这将阻碍相态传递以模拟多相之间的强融合效应.

2.4 构象张量方法拓展

构象张量可用于描述牛顿流体与非牛顿流体中的物质分布特性^［26］.对于三维空间中的构象张量 Q （3\times3矩阵），其经典更新公式为

\frac{D\mathit{Q}}{Dt}=\mathit{Q}\nabla \mathit{v}+{(\nabla \mathit{v})}^{\mathrm{T}}\mathit{Q}-\frac{\mathrm{1}}{\tau }\left(\mathit{Q}-\mathit{I}\right).(11)

其中：\tau表示用于描述流体黏弹性的弛豫时间；\mathit{I}表示构象张量的单位矩阵，代表物体的未变形状态.该模型无法模拟剪切稀化流体，因此本文采用扩展方程：

\frac{D\mathit{Q}}{Dt}=\mathit{Q}\nabla \mathit{v}+{(\nabla \mathit{v})}^{\mathrm{T}}\mathit{Q}-\frac{\mathrm{1}}{\tau }\left(\mathit{Q}-\mathit{I}\right)-\gamma \left(\mathit{Q}-\mathit{I}\right)\mathit{Q}.(12)

其中:\gamma \left(\mathit{Q}-\mathit{I}\right)\mathit{Q}为可模拟剪切稀化效应的非线性项，是改进模型的关键组成部分；\gamma∈［0，1］为稀化因子.\gamma值越大，稀化效应越显著.基于构象张量的混合应力{\sigma }_{i,\mathrm{m}\mathrm{i}\mathrm{x}}定义为

{\sigma }_{i,\mathrm{m}\mathrm{i}\mathrm{x}}=c{\eta }_{\mathrm{s}}\left({\mathit{Q}}_i-\mathit{I}\right).(13)

其中:c表示溶液的聚合物浓度，在本文的多相框架中等于{\phi }_{k\mathrm{2}}；{\eta }_{\mathrm{s}}表示溶液的黏度.根据相与混合层之间的关系{\sigma }_{\mathrm{m}\mathrm{i}\mathrm{x}}=\sum _k{\sigma }_k，可以计算相级应力{\sigma }_{i,k}：

{\sigma }_{i,k}={\phi }_{i,k}{\sigma }_{i,\mathrm{m}\mathrm{i}\mathrm{x}}.(14)

通过这种方式，相速度可根据混合物应力进行更新，其表达式为

\frac{D{\mathit{v}}_{i,k}^{\mathrm{v}\mathrm{i}\mathrm{s}\mathrm{c}}}{Dt}=\frac{\nabla \cdot {\sigma }_{i,k}}{{\phi }_{i,k}{\rho }_k^{\mathrm{0}}}.(15)

其中，上标visc代表黏性力.

算法1（表1）展示了本文提出的端到端模拟与可视化方法流程，本文称之为IMM-CT算法.

3 实验与分析

本文通过模拟骨水泥的各种调制和混合/注射过程来测试所提的方法，设计了2组实验来验证所提模型的优势.图4分别展示了所有的场景设置和模型，所用模型全部基于真实医疗器材建模.图5为临床所使用的骨水泥注射螺钉.

调制混合实验.本文通过不同黏度参数配置开展流体混合实验，并选取其中2组典型结果进行展示.根据骨水泥制备所用2种物质的特性，将流体配置为溶剂与溶质两部分，实际可视化结果见图6.所有实验场景中溶剂黏度均设为0.01 Pa⋅s，而溶质黏度在图6a中设为5 Pa⋅s，溶质与溶剂的黏度比为500，在图6b中设为10 Pa⋅s，溶质与溶剂的黏度比为1 000.实验结果表明，相较于现有单相与多相流体求解器，本方法展现出更高的稳定性.本文在图7和表2中提供了各种求解器的性能对比数据与实时计算耗时，时间步长随黏度比增大而减小，IMM-CT在各个区间均具优势，曲线表明本求解器在更广的流体黏度范围内均具有优越性能.

椎骨填充仿真实验.本文通过椎骨填充仿真验证本方法.骨骼建模流程如图3所示，在获取多孔松质骨与致密皮质骨后，应用粒子采样配合临床螺钉模型（图4）构建对应模型.本文开发了粒子发射器模拟螺钉上的6个注射端口，可自定义喷射方向、速度及出口尺寸.图8展示2种黏度下的骨水泥注射可视化结果：1） 高黏度组（上）与低黏度组（下），每组包含2个视角的RGB（red-green-blue）渲染图像.骨水泥显示为绿色，骨组织为半透明状，图像按时间序列排列，同列表示相同模拟时间轴；2） 特别提供时间t_{\mathrm{5}}的填充状态细节表征，图中红框标注区域显示流体边界，放大视图展示流体具体形态，多视角呈现帮助医师观察骨水泥流动状态与流体边界，以评估手术方案的可行性.本方法能清晰呈现填充过程与流体边界，辅助医师选择最佳注射策略，同时在注射场景配置上具有高度灵活性，包括螺钉定位、出口方向、速度、黏度等参数均可调.

椎骨填充对比实验.为验证本研究所提方法的可靠性，本文以相同物理参数配置的骨水泥为研究对象，通过对比2种数值模拟方法的仿真结果与真实场景下骨水泥注射后的实际分布情况，开展有效性验证实验.其中，PBF（position-based fluids）^［12］方法是流体力学领域中一种高效且通用的流体模拟算法.为保证对比实验的一致性，本研究将 PBF 方法的模拟流体黏度设置为与骨水泥实际黏度一致，随后通过2种方法分别观测骨水泥仿真注射过程中的流体边界形态演变特征（图9）.实验结果表明，在仿真注射后期，PBF方法因在黏度计算模型与流变行为表征方面存在固有缺陷，其模拟的流体边界出现显著偏差，甚至表现出非物理性渗漏现象；而本研究提出的方法可始终维持稳定的模拟过程，且其输出的流体边界形态与真实注射场景下的骨水泥边界特征吻合度更高，进一步证实了本方法在骨水泥仿真注射中的可靠性.

4 结 语

本研究提出了一种骨水泥注射仿真的计算框架，整合了先进的骨建模技术、多相非牛顿流体求解器和高保真渲染系统.该方法能够精确重建皮质骨与松质骨结构，同时模拟不同黏度下骨水泥的混合与流动过程.通过真实感色彩映射和流体表面重建技术增强可视化效果，为术前规划和手术培训提供有力支持.实验结果表明，相较于传统单相流及早期多相流方法，本框架在稳定性和计算效率方面具有显著优势，可实现实时仿真以优化注射策略，同时降低手术风险与成本.

未来研究将重点拓展两方面工作：一是深化骨水泥动态行为模拟，在现有流动特性仿真基础上，进一步纳入固化过程的时间依赖性变化，并引入拉伸力、毛细力等多物理场作用，以更精准地复现骨水泥在体内复杂环境中的材料特性与行为规律；二是推进与临床医师的协作研究，通过合作明确临床实际需求（如不同手术场景下的参数阈值、操作偏好等），为系统优化提供方向支撑.

参考文献

[1] Saha S， Pal S. Mechanical properties of bone cement： a review［J］. Journal of Biomedical Materials Research， 1984， 18（4）： 435-462.

[2] Duffy R K， Shafritz A B. Bone cement［J］. Journal of Hand Surgery， 2011， 36（6）： 1086-1088.

[3] Clarius M， Heisel C， Breusch S J. Pulmonary embolism in cemented total hip arthroplasty［M］//The Well-Cemented Total Hip Arthroplasty. Berlin/Heidelberg： Springer-Verlag， 2006： 320-331.

[4] Moussazadeh N， Rubin D G， McLaughlin L， et al. Short-segment percutaneous pedicle screw fixation with cement augmentation for tumor-induced spinal instability［J］. The Spine Journal， 2015， 15（7）： 1609-1617.

[5] Shen L， Zhang Y L， Frey S， et al. Visual simulation of bone cement blending and dynamic flow［C］//2024 IEEE International Conference on Bioinformatics and Biomedicine （BIBM）. Lisbon， 2025： 6226-6233.

[6] Charnley J. Anchorage of the femoral head prosthesis to the shaft of the femur［J］. The Journal of Bone & Joint Surgery British Volume， 1960， 42（1）： 28-30.

[7] Tan B W， Tang Q， Zhong Y J， et al. Biomaterial-based strategies for maxillofacial tumour therapy and bone defect regeneration［J］. International Journal of Oral Science， 2021， 13： 9-32.

[8] Jeong J， Kim J H， Shim J H， et al. Bioactive calcium phosphate materials and applications in bone regeneration［J］. Biomaterials Research， 2019， 23（1）： s40824-018-0149-3.

[9] Wang X K， Xu Y R， Liu S N， et al. Physics-based fluid simulation in computer graphics： survey， research trends， and challenges［J］. Computational Visual Media， 2024， 10（5）： 803-858.

[10] Macklin M， Müller M， Chentanez N. XPBD： position-based simulation of compliant constrained dynamics［C］//Proceedings of the 9th International Conference on Motion in Games. Burlingame， California， 2016： 49-54.

[11] Bender J， Koschier D. Divergence-free smoothed particle hydrodynamics［C］//Proceedings of the 14th ACM SIGGRAPH / Eurographics Symposium on Computer Animation. Los Angeles， 2015： 147-155.

[12] Macklin M， Müller M. Position based fluids［J］. ACM Transactions on Graphics， 2013， 32（4）： 104.1-104.12.

[13] Ma Y H， Xiao X Y， Li W， et al. Hybrid LBM-FVM solver for two-phase flow simulation［J］. Journal of Computational Physics， 2024， 506： 112920.

[14] Liu C， Hu C H. An adaptive multi-moment FVM approach for incompressible flows［J］. Journal of Computational Physics， 2018， 359： 239-262.

[15] Panuelos J， Goldade R， Grinspun E， et al. PolyStokes： a polynomial model reduction method for viscous fluid simulation［J］. ACM Transactions on Graphics， 2023， 42（4）： 124.1-124.13.

[16] Liu S S， He X W， Wang W C， et al. Adapted SIMPLE algorithm for incompressible SPH fluids with a broad range viscosity［J］. IEEE Transactions on Visualization and Computer Graphics， 2022， 28（9）： 3168-3179.

[17] Su H Z， Xue T， Han C， et al. A unified second-order accurate in time MPM formulation for simulating viscoelastic liquids with phase change［J］. ACM Transactions on Graphics， 2021， 40（4）： 119.1-119.18.

[18] Ren B， Li C F， Yan X， et al. Multiple-fluid SPH simulation using a mixture model［J］. ACM Transactions on Graphics， 2014， 33（5）： 171.1-171.11.

[19] Yang T， Chang J， Ren B， et al. Fast multiple-fluid simulation using Helmholtz free energy［J］. ACM Transactions on Graphics， 2015， 34（6）： 233.1-233.11.

[20] Li W， Ma Y H， Liu X P， et al. Efficient kinetic simulation of two-phase flows［J］. ACM Transactions on Graphics， 2022， 41（4）： 114.1-114.17.

[21] Yan X， Jiang Y T， Li C F， et al. Multiphase SPH simulation for interactive fluids and solids［J］. ACM Transactions on Graphics， 2016， 35（4）： 79.1-79.11.

[22] Li J L， Yang Y B， Sun Z W， et al. Integrated evaluation of biomechanical and biological properties of the biomimetic structural bone scaffold： biomechanics， simulation analysis， and osteogenesis［J］. Materials Today Bio， 2024， 24： 100934.

[23] Newman T S， Yi H. A survey of the marching cubes algorithm［J］. Computers & Graphics， 2006， 30（5）： 854-879.

[24] Bridson R. Fast Poisson disk sampling in arbitrary dimensions［C］// International Conference on Computer Graphics and Interactive Techniques. San Diego，2007： 22.

[25] Xu Y R， Wang X K， Wang J M， et al. An implicitly stable mixture model for dynamic multi-fluid simulations［C］//SIGGRAPH Asia 2023 Conference Papers. Sydney NSW，2023：263697252.

[26] Bejan A. The constructal law of organization in nature： tree-shaped flows and body size［J］. Journal of Experimental Biology， 2005， 208（9）： 1677-1686.

基金资助

国家自然科学基金资助项目(62306032)