随着汽车安全领域的快速发展，乘员约束系统有了显著改进，交通事故伤亡率大幅下降。然而，下肢损伤在交通事故损伤中仍然占有很高的比重^[1]，乘员约束系统对下肢保护效果的提升也不如头和胸部^[2]。虽然下肢损伤往往不是致命伤，但会使人丧失运动与劳动能力，降低生活质量，给家庭和社会带来沉重负担，Kuppa等^[1]提出用“因损伤造成的生命丢失”(life-years lost to injury, LLI)估算下肢损伤造成的生活质量下降程度，在其研究的事故中，简易损伤评级(abbreviated injury scale, AIS)为2+的下肢损伤带来的年均LLI为1.28×10⁵年，这一数值比由头部、颈部和胸部AIS 2+损伤造成的年均LLI之和还要高，对乘员下肢损伤风险和损伤机理的研究应该引起更多的重视。研究表明，人群中的个体特征差异对损伤风险有显著影响，与中等尺寸中年男性相比，肥胖和老年人群在交通事故中的伤亡风险更高^[3-4]。

有限元仿真分析是开展相关研究的有效工具并得到广泛应用，但现有的各类人体有限元模型却存在诸多不足之处。已有多个研究基于计算机断层扫描(computed tomography, CT)或者磁共振成像(magnetic resonance imaging, MRI)数据建立了精细的下肢骨有限元模型^[5-7]，但这些模型都是基于某一个体的扫描数据，未能反映人群个体特征差异对骨骼形状的影响，模型的密质骨厚度也无法更改，而已有研究表明长骨密质骨厚度对其力学性能有显著影响，且该厚度分布与年龄、身体质量指数(body mass index, BMI)和性别等因素有关^[8-9]。Klein等^[10]基于98例美国人CT扫描数据建立了考虑个体特征差异和密质骨厚度的股骨参数化几何模型，但其所用的密质骨厚度估算方法为全局固定阈值法，该方法对所选用阈值非常敏感，且无法对密质较薄的区域给出准确估计。Bryan等^[11]和Park等^[12]基于CT扫描数据，应用非刚性表面配准法和基于径向基函数的网格变换方法建立了一系列个体股骨模型，但模型没有考虑密质骨厚度的差异，也未能建立参数化股骨模型。胫骨参数化模型还未见于文献。THUMS模型^[13]和GHBMC模型^[14]是具有代表性的全身人体模型，但其开发思路与碰撞假人类似，只包括高大男性、中等尺寸男性和矮小女性3种尺寸，然而，这两种人体模型非常适用于作为建立参数化模型的基准模型。值得注意的是，上面提到的参数化模型都是基于西方人体数据建立，目前还没有使用中国人体数据建立的参数化模型。

在本研究中，以男性中等尺寸GHBMC复杂模型的股骨和胫骨作为基准模型，首次基于国内95例临床CT扫描数据建立了以性别、年龄、BMI和长骨长度(身高)为预测因子的股骨和胫骨参数化几何和密质骨厚度模型，并讨论了个体特征对骨骼几何形状和密质骨厚度分布的影响。本研究所建立的参数化模型可以用于组装建立完整的参数化人体有限元模型，也可用于对比分析东西方人群下肢几何形状和密质骨厚度分布的差异；通过赋予材料属性并施加边界条件，还可以用于对比分析东西方人群下肢的损伤风险和损伤机理，具有一定研究应用价值。

1 方法

本研究基本流程如下：通过CT扫描数据得到目标个体股骨和胫骨的几何外表面，通过网格变换和投影将基准有限元模型外表面与目标个体几何表面贴合，应用局部阈值法计算密质骨厚度，应用主成分和多元回归分析建立参数化几何模型和密质骨厚度模型，最终通过内部网格重建得到参数化模型。

1.1 股骨和胫骨几何的提取与标志点识别

本研究使用的人体CT扫描数据来源于国内多家医院，这些数据的获取和使用流程经过清华大学伦理委员会审核批准，项目编号为20160011。经过筛选，共有95例CT扫描数据作为研究对象，CT扫描的层厚为0.625 mm至1.25 mm，每层有512×512个像素，像素尺寸介于0.623 mm和1.079 mm之间。扫描对象中包含59例男性和36例女性个体，个体年龄范围为16~83岁，BMI范围为17.31~32.05 kg/m²，年龄和BMI之间不存在显著相关性。研究对象的个体特征分布如图 1所示。

股骨和胫骨几何模型的提取和重建在Mimics (19.0，Materialise NV)软件中完成，通过设定合适的阈值可以将每一层CT扫描切面中的骨骼部分提取出来，再通过三维重建和一系列操作得到光滑无孔洞的股骨和胫骨的几何外表面。将这些几何表面导入到Rhinoceros 3D (5.0，Robert McNeel & Associates)软件中进行标志点标记。在股骨的远心端和近心端共需标记23个标记点，在胫骨则需标记16个标记点，这些标记点为解剖学特征点和易于识别区分的特征点，例如股骨大转子最顶端顶点和髁间窝的中点。对于骨干区域，计算提取出来的每一层CT扫描切面的几何中心点，这些点组成了股骨和胫骨的解剖轴，在解剖轴上，通过线性插值的方法得到额外11个均匀分布在轴上的标记点。骨干的上标志点没有放置在几何表面，这样做可以避免在后续网格变换的过程中出现网格扭曲。所有标志点的识别均由同一人完成，以尽量保持标志点的一致性。图 4中“标志点识别”部分为股骨和胫骨的标志点。

1.2 网格变换和投影方法

本研究选用男性中等尺寸GHBMC模型的股骨和胫骨作为有限元基准模型，在基准模型的外表面，同样标记了与节1.1中描述一致的标志点，通过这些标志点建立了从有限元基准模型外表面到目标个体几何表面的映射关系。网格变换采用基于径向基函数的空间插值方法，该方法在多个研究中已经得到广泛应用^{[10, 15-16]}。在变换过程中，基准模型的节点坐标按照标志点之间的对应关系变换到目标个体的位置。该过程并不能保证所有节点落在目标个体的几何表面上，因还需要将变换后的基准模型向目标个体几何投影，使基准模型更加贴合目标个体的几何表面。本研究提出一种基于法向向量的投影方法，首先，计算经网格变换后的基准模型的每个单元的面法向向量(需将所有单元改为三角形单元)，对任一节点V，按单元面积加权计算所有包含该节点的单元的面法线向量的均值作为该节点的法向向量。

$ {\mathit{\boldsymbol{N}}_V} = \frac{{\sum\limits_{i = 1, \cdots , n} {{\mathit{\boldsymbol{N}}_{{F_i}}} \cdot {S_i}} }}{{\sum\limits_{i = 1, \cdots , n} {{S_i}} }}. $

(1)

其中：N_V为节点V的法向向量，F_i为包含该节点的第i个单元，N_Fi为单元F_i的面法向向量，S_i为单元F_i的面积。

然后在目标个体的几何表面上计算每个三角形面片的几何中心，与给定节点V法向向量N_V距离最近的几何中心点所在的面片即为该节点投影的目标几何面。最终，将该节点沿目标几何面的面法向向量投影到该几何面上即可。

该投影方法如图 2所示，下方网格为变换后的基准模型的局部(黄色)，上方面片为目标几何表面的局部(绿色和浅蓝色)，以点O为几何中心的三角形面即为节点V的投影目标面(浅蓝色)，N_V和N_O分别为点V和面O的法向向量，点V′为投影结果。从图 5中“网格变换和投影部分”中股骨的投影结果可以看出，投影之后的节点和单元(红色)能够非常好的与目标几何面(绿色)贴合。整个网格变换和投影过程由代码驱动在MATLAB (R2016a，Mathworks)中自动完成，处理一个个体的股骨和胫骨模型在本地计算机(Intel i5 @ 3.5 GHz, 16 GB RAM)上共需要2 min。

1.3 密质骨厚度估计方法

在前期研究中，李沛雨等^[17]提出了一种基于临床CT扫描的局部阈值法估计肋骨某一点处的密质骨厚度。本研究将该方法推广至下肢骨，并提出针对薄密质区域(如股骨头)和厚密质区域(如长骨骨干)应在计算公式中分别采用不同的密质骨系数计算厚度，该系数通过尸体试验测量标定，并得到了验证。以股骨CT扫描为例，简述该计算方法，如图 3所示。首先通过计算得到过测量点O′沿其法向向量的线段MN的亨氏单位(Hounsfield unit, HU)值—位置曲线及HU导数值—位置曲线，然后确定图中B、C点的位置和C点的HU值HU_C。

$ {\rm{DH}}{{\rm{U}}_{{\rm{thres}}}} = {\mu _{{D_{{\rm{HU}}}}}} + 4.27{\sigma _{{D_{{\rm{HU}}}}}} $

(2)

其中：${\mu _{{D_{{\rm{HU}}}}}}$和${\sigma _{{D_{{\rm{HU}}}}}}$分别为HU导数值—位置曲线MX段的均值和标准差，DHU_thres为特征阈值。B点为HU导数值—位置曲线上第一个取值等于或大于DHU_thres的点，C点为HU值—位置曲线上与B点位置相同的点。A点为HU值—位置曲线上的最大值点，其HU值为HU_A，那么，适用于该测量点的局部阈值HU_cor的计算公式为

$ {\rm{H}}{{\rm{U}}_{{\rm{cor}}}} = ({\rm{H}}{{\rm{U}}_A} - {\rm{H}}{{\rm{U}}_C}) \times {K_{{\rm{cor}}}} + {\rm{H}}{{\rm{U}}_C}. $

(3)

其中：K_cor为密质骨计算系数，在薄密质区域，K_cor=0.860，在厚密质区域，K_cor=0.522。最后，HU值—位置曲线上DE点之间的距离即为该测量点的密质骨厚度。节1.2中的基准模型经过网格变换和投影之后，外表面的每一个节点就是密质骨厚度的测量点。

1.4 参数化几何模型和厚度模型的建立

每个目标个体在空间中的位置并不一致，经过网格变换和投影之后得到的股骨和胫骨表面需要使用广义Procrustes方法在空间中进行刚性校准对齐，这是一种没有尺寸缩放和形状改变的对齐方法^{[10, 16]}。在建立关于个体特征参数的多元回归模型之前，需要应用主成分分析^{[10, 15-16]}对对齐后的数据进行降维处理，主成分分析可以将一组可能存在相关性的变量正交变换为一组线性不相关的变量，即主成分。选定合适数量的主成分，可以在一个新的降维空间内重新描述原始数据，同时尽可能保留原始数据间的信息。以主成分为因变量，选定合适的自变量，就可以对其进行回归分析。本研究分别建立了男性和女性的几何和密质骨厚度分布模型，回归分析中的自变量为年龄、BMI和股骨(胫骨)长度。股骨长度定义为股骨头中心到髁间窝中点的距离，胫骨长度定义为胫骨粗隆向前最凸起点到内踝内测点的距离，以上距离都可以通过标志点计算得到。

给定个体特征参数，结合回归分析和主成分分析结果，就可以生成相应的预测模型(几何形状和密质骨厚度分布)。

$ P = \bar G + \mathit{\boldsymbol{C}} \cdot \left( {\mathit{\boldsymbol{M}} \cdot F'} \right) $

(4)

其中：P为预测模型，$\bar G$为平均模型，C为主成分系数矩阵，M为回归参数矩阵，F=[年龄BMI骨长1]为个体特征参数矩阵。

通过式(4)得到的模型包含股骨和胫骨的几何外表面和表面上每一节点处的密质骨厚度信息，还需对股骨和胫骨的内表面重建以得到完整的几何模型。在GHBMC模型中，股骨和胫骨骨干区域为带有骨髓腔的中空结构，密质骨厚度通过多层六面体单元表征，因此，需要根据预测得到的密质骨厚度对体单元尺寸进行调整，本节以股骨为例介绍该过程。在基准模型的骨干区域，最外层和最内层的节点和单元并没有一一对应关系，无法直接调整两层节点之间的距离实现密质骨厚度的变化。应用节1.2中的投影方法，将基准模型的内表面节点投影到外表面，构建从内层节点到内层投影后的节点之间的“伪内外对应关系”来解决这一问题(如图 4“密质骨厚度调整”所示)。原外层表面节点处的密质骨厚度可由式(4)计算得到，利用内层投影后的节点与外层节点位置的对应关系，使用基于径向基函数的插值方法估计内层投影后的节点位置处的密质骨厚度值。通过调整原内层节点位置，就可以得到想要的厚度值，在调整内层的节点位置之后，以其为内边界条件，外表面为外边界条件，使用节1.2中的网格变换方法，将基准模型的所有其他节点再进行一次变换，就可以得到与基准模型具有相同节点和单元数的股骨模型。

骨干以外的区域主要为松质骨，在GHBMC模型中以四面体单元表征，这些区域的密质骨由覆盖在松质骨外表面的壳单元表征，通过修改相应关键字，可以直接将厚度值赋予壳单元的每个节点，这里不再赘述。图 4概括总结了开发建立股骨和胫骨参数化模型的整个过程。

2 实验结果

每一个个体经网格变换和投影得到(简称“变换得到”)的几何表面和经CT影像重建得到(简称“重建得到”)的几何表面之间的差异为^[18]

$ \begin{array}{l} {\rm{D}}\left( {X, Y} \right) = \frac{1}{2}\left[ {\frac{1}{U}\sum\limits_{x \in X} {\mathop {\min d}\limits_{y \in Y} } } \right.\left( {x, y} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\left. {\frac{1}{V}\sum\limits_{y \in Y} {\mathop {\min d}\limits_{x \in X} \left( {x, y} \right)} } \right]. \end{array} $

(5)

其中：x和y分别为表面X和Y中的各个顶(节)点，d(x, y)为两点之间的Euler距离，U和V分表面X和Y的总顶(节)点数量。所有股骨样本经变换得到和经重建得到的几何表面之间的平均差异为1.02 mm，胫骨样本的平均差异为1.35 mm。结果表明通过变换得到的几何表面具有很高的精度，可以很好的反映个体特征差异。

对于几何形状模型，男性和女性的股骨模型分别选取了56和34个主成分，胫骨模型分别选取了50和33个主成分，这些主成分能够描述原数据中99%以上的差异。对于密质骨厚度模型，男性和女性的股骨和胫骨模型都分别选取了56和34个主成分，这些主成分能够描述原数据中99%以上的差异。每个回归模型的决定系数R²为

$ {R^2} = 1 - \frac{{{\rm{RSS}}}}{{{\rm{TSS}}}}. $

(6)

其中：RSS为预测模型和观测模型的均方差(对应节点的坐标或者厚度值)，TSS为平均模型和观测模型的均方差(对应节点的坐标或者厚度值)，每个模型的R²值列于表 1中。对于股骨，男性和女性几何预测模型和观测模型节点位置的平均绝对误差分别为2.93和2.56 mm，密质骨厚度预测模型和观测模型的平均绝对误差分别为0.56和0.58 mm。对于胫骨，以上误差分别为3.00、2.45、0.60和0.50 mm。

回归模型中，对前5个主成分的各个预测变量(即自变量)的p值分别列于表 2和3中。结果表明：骨长对股骨和胫骨的几何模型的影响最为显著(如图 4“参数化模型”所示)，年龄和BMI对女性股骨和胫骨的几何模型具有显著影响。年龄和BMI对股骨厚度模型具有显著影响，只有年龄对胫骨厚度模型具有显著影响。骨长对所有厚度模型都没有显著影响，因此将其从厚度模型的预测变量中移除。以股骨为例，选取骨干区域5个位置的截面可以进一步验证上述部分论述，如图 5所示，截面旁标记的数字为该截面的平均密质骨厚度值(mm)。

参数化模型与基准模型具有相同数量的节点和单元，网格质量方面，基准模型体单元的Jacobian行列式值最小为0.3，参数化模型体单元的Jacobian行列式值一般在0.2以上，当个体特征参数取值与均值偏差过大而位于两端时，模型偏差会变大，对网格质量造成一定影响，经过多次验证，生成的所有模型的Jacobian行列式值都大于0.15。在现阶段，本研究建立的参数化模型的材料属性与基准模型保持一致。

3 讨论

本文几何模型中，没有选择身高作为回归分析自变量的原因主要有两点，一是本研究的最终目的是建立完整的人体参数化有限元模型，骨长可以通过其他的人体外表皮参数化模型计算得到^[19]，使用骨长作为自变量，可以确保本研究预测得到的股骨和胫骨模型能够与人体外表皮模型兼容；二是骨长与身高成显著正相关关系，骨长可以通过身高预测得到，同时也可以很容易地将模型中的骨长用身高替换重新进行回归分析。

股骨和胫骨几何模型的R²较高，而股骨和胫骨厚度模型的R²值较低，说明模型中的现有预测变量未能很好的预测密质骨厚度的分布，这表明需要增加其他个体特征参数来提高模型的预测能力。例如，个体饮食运动习惯和骨质疏松等情况也可能会影响密质骨厚度分布。研究结果表明骨长对几何模型的影响最为显著，年龄和BMI对女性股骨和胫骨几何模型的影响较为显著，股骨部分的结论与Klein等^[10]的研究结果一致。本研究表明：BMI对股骨密质骨厚度分布有显著影响，这与Klein等^[10]对美国人群的研究结果吻合。从图 5中可以看出，来自BMI较高个体的股骨具有更厚的骨密质。本研究还表明年龄对股骨和胫骨厚度模型都起到显著影响作用，从图 4“密质骨厚度分布”中的分布云图可以看出，老年人的密质骨厚度要显著低于年轻人，这与Clarke^[20]的研究结论一致，该研究指出，健康成年人随着年龄的增长，骨髓腔会变大而密质骨会变薄。目前文献中尚没有关于胫骨模型的讨论，但本研究表明BMI对胫骨密质骨厚度的影响并不显著。

本研究中，对于男性和女性个体分别建立了模型，并没有将“性别”作为回归模型的因变量，这主要是为了提高参数化模型的准确性。然而，个体特征对模型的影响在不同性别之间却存在不同之处，因此，性别因素对模型的影响还需要进一步研究，这需要建立更为复杂的模型，例如线性混合模型。女性样本数量相对较少，年龄分布也有一定缺失，未来还需要补充更多的女性样本。

现阶段参数化模型的材料属性来自于基准模型，材料属性可以依据CT数据和个体特征进行修改，Park等^[21]和Helgason等^[22]介绍了HU值、骨密度和弹性模量之间的关系，Dokko等^[23]研究了年龄对骨骼弹性模量和屈服强度的影响。未来工作包括考虑个体特征对材料属性的影响，对模型进行广泛验证，并将其应用到有限元仿真分析中。

4 结论

本文提出了一套建立人体股骨和胫骨参数化模型的流程和方法，首次使用中国人体CT扫描数据建立了反映个体特征差异的股骨和胫骨参数化几何模型和密质骨厚度分布模型，讨论了个体特征对股骨和胫骨几何和密质骨厚度分布的影响。这些参数化模型是开发完整人体参数化有限元模型的基础，并可广泛应用于下肢损伤风险和损伤机理的研究。