2019年，全球骨关节炎发病率逐年上升，已超过7%，人数超过5亿^[1]，其中3.44亿处于中重度阶段，亟需康复治疗改善^[2]. 康复治疗须根据个体差异设计个性化的治疗方案，对人体工学产品定制的要求较高. 针对个性化的运动关节康复需求，本研究1）通过精确重建关节运动姿态，优化外骨骼关节机构，在运动域中减小瞬心轨迹误差，实现变瞬心仿生运动，提升人机协同性能. 2）通过应力分布模拟和载荷特性分析，在体积和应力约束下，优化外骨骼零件的材料分布和结构设计，实现外骨骼的定制设计.

识别和了解人体关节角度与姿态有助于评估关节康复情况，为运动能力改善提供基准，并在治疗方案制定中发挥关键作用. Lalwala等^[3]研究动态负荷的下肢运动学，Vianello等^[4]根据姿势和末端执行器速度预测关节运动，Mousse等^[5]提出提高跌倒检测准确性的姿态识别算法，Takano等^[6]利用2D相机图像估计人体姿态. 明确人体在运动过程中的运动模式和骨骼关节活动机制，能够为外骨骼设计提供必要的数据支持. Simon等^[7]进行肩-髋-膝关节角度分析，设计出增强背部支撑的外骨骼；Lerner等^[8]制造出关节完全伸展的外骨骼；Barrutia等^[9]使用机械关节模拟器分析行走中的矢状面运动学特征；陈栋等^[10]基于步态重建和关节运动信息设计柔性踝关节外骨骼机器人，通过优化模态参数实现人机一体化. Missiroli等^[11]设计出提供关节屈曲伸展重力支撑的外骨骼，使肌肉活动减少了32%；De Groof等^[12]设计出具有伸展扭矩功能的外骨骼，旨在防止关节伸展不足和过度屈曲；Perry等^[13]在可变加载的条件下研究外骨骼装置对关节屈曲性能的影响.

获取准确的关节运动信息对推动康复设备发展至关重要. 常见的关节运动学检测方法包括皮肤标记和运动捕捉. Hachaj等^[14]通过关节运动捕捉计算加速度，Dobos等^[15]比较了单摄像头无标记和3D光学标记技术在运动学数据采集中的准确度，Ziegler等^[16]利用3D标记位置识别下肢几何和步态轨迹. 摄像机和惯性传感器可捕捉实时数据，但成本高且操作复杂^[17]. 关节运动学与配准技术密切相关，配准技术能将不同数据源的信息对齐，为运动学分析提供准确输入. Saadat等^[18]提出在X光透视图像中插值预测3D位置并微调的方法，通过从粗到精的配准在提高配准速度的同时保持了位置预测的准确性；Matsuki等^[19]利用图像配准研究胫骨相对于股骨的平移和旋转. 配准须消耗大量计算资源，且数据质量对精度影响较大. 分析关节运动过程中的生物力学特性有助于疾病预防^[20]. Thienkarochanakul等^[21]发现膝关节软骨层和半月板内侧区域承受较高应力，易受损；Sidhu等^[22]在模拟器上测量运动过程中关节肌肉和组织松弛度；Ng等^[23]采用有限元方法分析了胫骨近端骨移植物的生物力学；Park等^[24]发现行走过程中膝关节的内侧软骨接触压力高于外侧.

在前期工作^[25-28]基础上，针对支持问题重构与方案演化的人体工学产品设计难题，本研究提出基于锥规划的运动关节康复外骨骼设计方法. 1）以膝关节为例，通过医学图像重建微观骨骼结构，利用锥规划迭代计算得到形貌特征匹配关系，重建关节运动姿态. 2）分析膝关节瞬时旋转中心的运动轨迹，优化外骨骼结构设计，通过反向双摇杆模拟变瞬心运动，在运动域中提高轨迹匹配性. 3）通过应力分布模拟，在体积和应力约束下，优化材料分布和结构形态，完成康复外骨骼定制设计.

基于锥规划的运动关节姿态重建方法利用半正定约束降低计算复杂度，提高重建的计算效率，并通过形貌特征描述符进行匹配计算，提升姿态重建的精度. 锥规划的目标是在某些锥约束下最小化线性目标函数，典型特点是变量必须落在某个锥形区域内. 锥规划的数学公式为

(1)$ \left. \begin{gathered} \min {\text{ }}{\boldsymbol{Hx}}. \\ {\text{s}}{\text{.t}}{\text{. }}{\boldsymbol{A}}x = {\boldsymbol{b}},\;{\boldsymbol{x}} \in K. \\ \end{gathered} \right\} $

式中：$ {\boldsymbol{x}} $为规划向量，$ H $为目标函数系数矩阵，$ {\boldsymbol{A}} $为约束条件系数矩阵，$ {\boldsymbol{b}} $为约束条件常数向量，$ K $为特定锥. 当约束条件为半正定矩阵时，形成半正定锥规划. 令$ {\boldsymbol{y}} = {\boldsymbol{x}}{{\boldsymbol{x}}^{\text{T}}} $，式(1)可转换为

(2)$ \left. \begin{gathered} \min {\text{ tr}}\;({\boldsymbol{Hy}}{\text{)}}. \\ {\text{s}}{\text{.t}}{\text{. }}{S_k}({\boldsymbol{y}}) \geqslant 0. \\ \end{gathered} \right\} $

式中：$ {S_k}( \cdot ) $为矩阵的顺序主子式，如果所有顺序主子式均非负，则矩阵$ {\boldsymbol{y}} $为半正定矩阵；$ {\text{tr}}\;( \cdot ) $为矩阵的迹，为矩阵的特征值之和. 特征值通过求解特征方程得到.

(3)$ \det\; ({\boldsymbol{y}} - \lambda {\boldsymbol{E}}) = 0. $

式中：$ \det\; ( \cdot ) $为矩阵的行列式，$ {\boldsymbol{E}} $为单位矩阵，$ \lambda $为特征值；如果所有特征值均非负，则$ {\boldsymbol{y}} $为半正定矩阵. 在锥规划的基础上，根据运动关节的形貌和咬合特征，分别选取源关节咬合曲面$ C $和目标关节咬合曲面$ M $，计算$ C $的包围盒并进行给定体素网格数量$ {N_v} $划分.

(4)$ \left. \begin{gathered} {N_x} = \left\lfloor {\frac{{{x_{\max }} - {x_{\min }}}}{{{N_v}}}} \right\rfloor +1, \\ {N_y} = \left\lfloor {\frac{{{y_{\max }} - {y_{\min }}}}{{{N_v}}}} \right\rfloor +1, \\ {N_z} = \left\lfloor {\frac{{{z_{\max }} - {z_{\min }}}}{{{N_v}}}} \right\rfloor +1. \\ \end{gathered} \right\} $

式中：$\left\lfloor {\cdot} \right\rfloor $表示向下取整，$ {x_{\min }} $、$ {x_{\max }} $、$ {y_{\min }} $、$ {y_{\max }} $、$ {z_{\min }} $和$ {z_{\max }} $分别对应轴向的最大值和最小值，$ {N_x} $、$ {N_y} $和$ {N_z} $分别为沿X、Y、Z向划分的体素个数. 更新包围盒中的所有体素，如果体素包含$ C $的部分，则用体素中心表示该部分. 所有体素中心形成源咬合曲面$ {P_{{C}}} $. 更新$ M $的包围盒的体素，形成目标咬合曲面$ {Q_{{M}}} $. 对$ {P_{{C}}} $和$ {Q_{{M}}} $进行咬合特征检测，生成源特征集$ {H_{{P}}} $和目标特征集$ {H_{{Q}}} $. 针对任意特征$ p $及其正方体邻域，计算$ p $对应的局部特征矩阵$ {{\boldsymbol{V}}_{{p}}} $，

(5)$ {{\boldsymbol{V}}_{{p}}} = \frac{1}{{\displaystyle\sum\nolimits_{i = 1}^{{n_{{p}}}} {{\omega _i}} }}\sum\nolimits_{i = 1}^{{n_{{p}}}} {{\omega _i}} \left[ {\begin{array}{*{20}{c}}{{q_{i,x}}{p_x}}&{{q_{i,x}}{p_y}}&{{q_{i,x}}{p_z}}\\{{q_{i,y}}{p_x}}&{{q_{i,y}}{p_y}}&{{q_{i,y}}{p_z}}\\{{q_{i,z}}{p_x}}&{{q_{i,z}}{p_y}}&{{q_{i,z}}{p_z}}\end{array}} \right]. $

式中：$ {\omega _i} $为$ p $的权重，$ {q_i} $为$ p $正方体邻域包含的某个特征，$ {n_{{p}}} $为$ p $正方体邻域中包含的特征数量， q_i,x、 q_i,y、q_i,z分别为q_i的x、y、z坐标，p_x、p_y、p_z分别为p的x、y、z坐标. 通过局部体素网格描述符对$ {H_{{P}}} $和$ {H_{{Q}}} $分别进行描述，生成源咬合曲面特征描述符$ {F_{{P}}} $和目标合曲面特征描述符$ {F_{{Q}}} $，流程如下.

1）在$ {H_{\text{P}}} $中选取任意特征$ p $为中心，以给定边长$ {l_{\text{E}}} $构建正方体邻域，计算对应的局部特征矩阵$ {{\boldsymbol{V}}_{{p}}} $. 对$ {{\boldsymbol{V}}_{{p}}} $进行特征分解，以3个特征向量为轴建立局部参考坐标系.

2）在局部参考坐标系中，根据给定的体素尺寸沿X、Y、Z向划分体素. 如果体素内包含特征，则该体素标记为1，否则标记为0.

3）连接所有体素标记，形成局部体素网格描述符.

4）在$ {H_{{P}}} $中进行全域更新，生成$ {F_{{P}}} $.

5）在$ {H_{{Q}}} $中进行全域更新，生成$ {F_{{Q}}} $.

计算$ {F_{{P}}} $和$ {F_{{Q}}} $的描述符距离比$ {R_{{\mathrm{d}}}} $，

(6)$ {R_{\text{d}}} = \frac{{\left| {{F_{Pi}} - {F_1}} \right|}}{{\left| {{F_{Pi}} - {F_2}} \right|}}. $

式中：$ {F_{Pi}} $为第i个源咬合曲面特征描述符，$ {F_1} $、$ {F_2} $为$ {F_{{Q}}} $中距离$ {F_{Pi}} $最近的2个目标合曲面特征描述符. 给定距离比阈值$ \alpha $，通过比较$ {R_{\mathrm{d}}} $判断$ {F_{{P}}} $与$ {F_{{Q}}} $的匹配关系. 如果$ {R_{\text{d}}} \leqslant \alpha $，则认为$ {F_{{P}}} $与$ {F_{{Q}}} $匹配；否则，认为不匹配. 随机选取n个符合匹配关系的特征，形成源曲面$ P $和目标曲面$ Q $. 设定$ a,b \in \{ 1,n\} $为$ P $中的特征索引，在$ P $中，任意选取曲面特征$ {P_a} $和$ {P_b} $；设定$ c,d \in \{ 1,n\} $为目标曲面$ Q $中的特征索引，在$ Q $中，任意选取曲面特征$ {Q_c} $和$ {Q_d} $. 计算$ {P_a} $和$ {P_b} $、$ {Q_c} $和$ {Q_d} $形成的线段长度差$ {l_{\text{d}}} $，

(7)$ {l_{\text{d}}} = \delta ({P_{{a}}} - {P_{{b}}}) - \delta({P_{{c}}} - {P_{{d}}}). $

式中：$ \delta({P_{{a}}} - {P_{{b}}}) $为$ {P_a} $和$ {P_b} $之间的欧氏距离，$ \delta({P_{{c}}} - {P_{{d}}}) $为$ {Q_c} $和$ {Q_d} $之间的欧氏距离. 计算$ {P_a} $和$ {P_b} $、$ {Q_c} $和$ {Q_d} $形成的夹角$ \beta $，

(8)$ \beta = \arccos\; ({{\boldsymbol{V}}_1} - {{\boldsymbol{V}}_2}) - \arccos\; ({{\boldsymbol{V}}_3} - {{\boldsymbol{V}}_4}). $

式中：$ {{\boldsymbol{V}}_1} $、$ {{\boldsymbol{V}}_2} $、$ {{\boldsymbol{V}}_3} $、$ {{\boldsymbol{V}}_4} $别为与$ {P_a} $、$ {P_b} $、$ {Q_c} $和$ {Q_d} $对应的$ {{\boldsymbol{V}}_{{p}}} $. 给定长度差阈值$ {l_{\text{s}}} $和角度阈值$ {\beta _{\text{s}}} $，计算$ {P_a} $和$ {P_b} $、$ {Q_c} $和$ {Q_d} $形成的得分$ {A_{{\text{ac,bd}}}} $，

(9)$ {A}_{\text{ac,bd}}=\left\{ \begin{array}{*{20}{l}}\mathrm{exp}\;(-\left|{l}_{\text{d}}\right|),& \left|{l}_{\text{d}}\right| < {l}_{\text{S}}且\beta < {\beta }_{\text{S}}; \\0,& 其他. \end{array} \right.$

对$ P $和$ Q $上的特征组成的线段进行广度优先搜索，动态更新不同线段对应的$ {A_{{\text{ac,bd}}}} $，组成得分矩阵$ {{\boldsymbol{M}}_{{A}}} $，

(10)$ {{\boldsymbol{M}}_{{A}}} = \left[ {\begin{array}{*{20}{c}} {{A_{1,1}}}& \cdots &{{A_{1,{n^2}}}} \\ \vdots & {} & \vdots \\ {{A_{{n^2},1}}}& \cdots &{{A_{{n^2},{n^2}}}} \end{array}} \right]. $

其中$ {\mathrm{ac}} = a+(c - 1)n $为$ {{\boldsymbol{M}}_{{A}}} $的行索引，$ {\mathrm{bd}} = b+(d - 1)n $为$ {{\boldsymbol{M}}_{{A}}} $的列索引. 对$ {{\boldsymbol{M}}_{{A}}} $进行锥规划计算，得到最大得分矩阵对应的匹配关系. 锥规划过程表示为

(11)$ \left. \begin{gathered} \max {\text{ }}{\boldsymbol{X}}{{\boldsymbol{M}}_A}{\boldsymbol{X}}. \\ {\text{s}}{\text{.t}}{\text{. }}\;\sum\nolimits_i^{{n_i}} {{X_{i,j}} \leqslant 1} ,\;\;\;\; i \in \{ 1,2, \cdots ,{n_i}\} ; \\ \qquad\sum\nolimits_j^{{n_j}} {{X_{i,j}} \leqslant 1,\;\;\; j \in \{ 1,2, \cdots ,{n_j}\} } . \\ \end{gathered} \right\} $

令$ {\boldsymbol{Y}} = {\boldsymbol{X}}{{\boldsymbol{X}}^{\text{T}}} $，则可转换成半正定锥规划问题：

(12)$ \left. \begin{gathered} \max {\text{ tr}}\;({{\boldsymbol{M}}_A}{\boldsymbol{Y}}). \\ {\text{s}}{\text{.t}}{\text{. }}\;{S_k}({\boldsymbol{Y}}) \geqslant 0; \\ \qquad\sum\nolimits_i^{{n_{\text{Y}}}} {{Y_{i}} \leqslant 1} . \\ \end{gathered} \right\} $

式中：$ {n_{\text{Y}}} $为Y中元素的个数. 在$ {{\boldsymbol{M}}_{{A}}} $值最大的匹配关系下，计算符合匹配关系的$ P $和$ Q $的半正定匹配函数. 通过半正定匹配函数将$ Q $重建为目标曲面$ Q' $，计算$ P $和$ Q' $的均方根误差RMSE，

(13)$ {\text{RMSE = }}\sqrt {\frac{1}{N}\sum\nolimits_i^N {\mathop {\min }\limits_j } {{({P_i} - Q_j^{'})}^2}} . $

式中：N为P上的特征数量，$ {P_i} $为P上第i个特征，$ Q_j^{'} $为$ Q' $上的第j个特征点. RMSE越小，代表重建的精度越高. 给定最大迭代次数，对半正定锥规划匹配函数进行迭代计算，迭代过程中将最小均方根误差对应的半正定锥规划重建匹配函数作为最优解，重建屈曲过程中的人体骨骼关节运动学姿态.

基于载荷传递和力线的拓扑优化能够准确地反映复杂载荷条件下的力学行为，为结构设计提供更有效的解决方案. 在三维流形模型中，表面上任意点的坐标可由张量积函数表示：

(14)$ t({\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}}) = \sum\nolimits_{i = 1}^{{n_{\text{R}}}} {\sum\nolimits_{j = 1}^{{m_{\text{R}}}} {\sum\nolimits_{k = 1}^{{l_{\text{R}}}} {R_{i,j,k}^{{\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}}}{T_{i,j,k}}} } } . $

式中：$ t({\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}}) $为表面上的任意点坐标，$ {T_{i,j,k}} $为控制点的坐标，$ R_{i,j,k}^{{\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}}} $为非均匀有理B样条的基函数，$ {\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}} $为三维流形模型的3个维度向量，$ i,j,k $分别为$ {\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}} $上控制点的索引，$ {n_{\text{R}}} $、$ {m_{\text{R}}} $、$ {l_{\text{R}}} $分别为$ {\boldsymbol{u}},{\boldsymbol{v}},{\boldsymbol{w}} $上控制点的数量. 对于线弹性问题，单元的刚度矩阵为

(15)$ {{\boldsymbol{K}}_{\text{c}}} = \int_\varOmega {{{\boldsymbol{B}}^{\text{T}}}{\boldsymbol{DB}}} {\text{d}}\varOmega {\text{.}} $

式中：$ \varOmega $为单元域，B为应变-位移矩阵，D为弹性矩阵. 对于各向同性的弹性体材料，

(16)$ {\boldsymbol{D}} = \frac{E}{{1 - {\mu ^2}}}\left[ {\begin{array}{*{20}{c}} 1&\mu &0 \\ \mu &1&0 \\ 0&0&{\dfrac{{1 - \mu }}{2}} \end{array}} \right]{\text{.}} $

式中：$ E $为材料的弹性模量，$ \mu $为泊松比. 全局刚度矩阵K由$ {{\boldsymbol{K}}_{\text{c}}} $构造得到，

(17)$ {\boldsymbol{K}} = \sum\nolimits_{e = 1}^{{n_{\text{e}}}} {{\boldsymbol{T}}_e^{\text{T}}{{\boldsymbol{K}}_{\text{c}}}{{\boldsymbol{T}}_e}} . $

式中：$ {{\boldsymbol{T}}_e} $为单元节点自由度与全局自由度之间的映射矩阵，e为单元编号，$ {n_e} $为单元数量. 位移通过求解线性方程得到，

(18)$ {\boldsymbol{KU}} = {\boldsymbol{F}}. $

式中：$ {\boldsymbol{K}} $为位移矩阵，$ {\boldsymbol{F}} $为作用的外力矩阵. 结合材料本构关系，通过有限元法进一步计算应力$ \sigma $.

(19)$ \sigma = {\boldsymbol{DBU}}. $

对于三维流形模型，主应力通过求解应力的特征方程得到，

(20)$ \det \;(\sigma - {\lambda _\sigma }{\boldsymbol{E}}) = 0. $

式中：$ {\lambda _\sigma } $为$ \sigma $的特征值，即主应力. 主应力形成的力线显示载荷传递的特征. 在三维流形模型的每个控制点上引入密度值$ \rho _{i,j,k}^{'} $，通过Shepard插值函数进行光滑处理. 光滑后的控制点密度$ {\rho _{i,j,k}} $：

(21)$ {\rho _{i,j,k}} = \sum\nolimits_{i = 1}^{{n_{\text{R}}}} {\sum\nolimits_{j = 1}^{{m_{\text{R}}}} {\sum\nolimits_{k = 1}^{{l_{\text{R}}}} {{S_{\text{T}}}\rho _{i,j,k}^{'}} } } . $

式中：$ {S_{\text{T}}} $为Shepard插值函数在控制点处的值. 以最小化应变能为目标，在体积和应力约束下，建立优化模型，确定设计结构中的材料分布：

(22)$ \left. \begin{gathered} \min {\text{ }}{J_{\text{e}}} = \frac{1}{2}{\int_\varOmega {\boldsymbol{U}} ^{\text{T}}}{\boldsymbol{DU}}{\mathrm{d}}\varOmega . \\ {\mathrm{s.t.}}{\text{ }}0 \leqslant {\rho _{i,j,k}} \leqslant 1,\; \sigma \leqslant {\sigma _{\text{s}}}, \; \frac{V}{{{V_0}}} \leqslant {V_{\text{s}}}. \\ \end{gathered} \right\} $

式中：$ {J_{\text{e}}} $为结构应变能，$ {\sigma _{\text{s}}} $为给定的应力阈值，V为优化后结构的体积，$ {V_0} $为优化前结构的体积，$ {V_{\text{s}}} $为给定的体积约束. 对$ {\rho _{i,j,k}} $进行Heaviside投影，形成增材制造约束，确保边界清晰并隐式控制最小尺寸.

膝关节是人体最大、最复杂的关节之一，在支撑和运动中作用关键. 以左膝关节为例，进行屈曲运动姿态的重建. 如图1所示，根据CT图像进行三维重建，重建的股骨体积为115 935.90 mm³，重建的胫骨体积为73 515.11 mm³. 股骨截面面积S和截面面积差分S_d随分层层数R的变化曲线如图2所示. S最大值为3 001.01 mm²（第2 840层），最小值为0.07 mm²（第12 500层），平均值为1 377.81 mm²，标准差为803.51 mm². S_d最大值为3.14 mm²（第1 109层），最小值为−1.11 mm²（第12 066层），平均值为−8.80×10⁻⁶ mm²，标准差为0.69 mm². 第6 923~11 471层相邻截面的面积差分变化平缓，为股骨骨干区域；第1~6 923层相邻截面的面积差分变化剧烈，为股骨髁区域. 计算胫骨截面面积和截面面积差分，识别胫骨平台和胫骨骨干区域. 根据运动关节的形貌和咬合特征，通过锥规划迭代计算，对运动关节的屈曲姿态进行精确重建. 定义膝关节的屈曲角度为矢状面上股骨与胫骨解剖轴线之间的角度. 不同方法（表面匹配^[29]、分级匹配^[30]、锥规划重建）重建的左膝关节运动姿态均方根误差迭代曲线如图3所示. 图中，I为迭代次数，$ \theta $为屈曲角度. 可以看出，均方根误差均随迭代次数的增加而减小. 不同重建方法的均方根误差对比如表1所示，其中RMSE_MAX、RMSE_MIN、$ \overline {{\text{RMSE}}} $分别为RMSE的最大值、最小值和平均值.

根据屈曲运动特征，对股骨和胫骨的咬合曲面进行分割. 如图4所示为矢状面上左胫股关节的锥规划重建运动姿态. 膝关节运动是旋转副与移动副的耦合运动，当关节进行屈曲运动时，其旋转中心的变化轨迹为J形曲线，不同屈曲角度下的瞬时旋转中心参数如表2所示，其中$ {d}_{{\mathrm{c}}} $为瞬时旋转中心到股骨解剖轴距离，$ {c}_{{\mathrm{c}}} $为瞬时旋转中心曲率. 可以看出，随着膝关节屈曲角度的增大，瞬时旋转中心与股骨解剖轴的距离逐渐增大，瞬心曲率呈现先增后减的趋势. 除均方根误差外，平均绝对误差MAE和最大误差ME也是常见的动态重建的效果评价指标. 平均绝对误差衡量整体重建精度，计算式为

(23)$ {\text{MAE}} = \frac{1}{N}\sum\nolimits_i^N {\mathop {\min }\limits_j } \left| {{P_i} - Q_j^{'}} \right|. $

最大误差用于检测极端重建失误，计算式为

(24)$ {\text{ME}} = \mathop {\max }\limits_i \left( {\mathop {\min }\limits_j \left| {{P_i} - Q_j^{'}} \right|} \right). $

如表3所示为不同方法的运动左膝关节姿态重建精度对比. 可以看出，相比表面匹配方法，锥规划重建方法的均方根误差降低了42.76%，平均绝对误差降低了44.85%，最大误差降低了28.26%；相比分层匹配方法，锥规划重建方法的均方根误差降低了10.95%，平均绝对误差降低了12.29%，最大误差降低了6.05%.

参考人体尺寸数据和膝关节屈曲重建姿态，基于人体工学，在运动域中对左膝关节的康复外骨骼结构进行匹配性定制设计. 如图5所示，反向双摇杆变瞬心机构位于大腿组件和小腿组件之间，在运动域中实现外骨骼运动关节的瞬时旋转中心J形曲线运动，提高人机协同运动的匹配性. 采用正运动学方法，结合封闭矢量方程与三心定理，建立如图6所示的左膝关节外骨骼变瞬心机构运动学模型，计算外骨骼变化的瞬时旋转中心. 其中BC杆与大腿固定，AD杆与小腿相连，AD杆和BC杆的交点P为反向双摇杆变瞬心机构的瞬时旋转中心. 反向双摇杆变瞬心机构的矢量封闭方程式为

(25)$ \left. \begin{gathered} {L_1}\cos {\alpha _1}+{L_4}\cos {\alpha _4} = {L_2}\cos {\alpha _2}+{L_3}\cos {\alpha _3}, \\ {L_1}\sin {\alpha _1}+{L_4}\sin {\alpha _4} = {L_2}\sin {\alpha _2}+{L_3}\sin {\alpha _3}. \\ \end{gathered} \right\} $

根据三心定理，机构瞬时旋转中心P表示为

(26)$ \left[ {\begin{array}{*{20}{c}} {{x_{{P}}}} \\ {{y_{{P}}}} \end{array}} \right] = {\left[ {\begin{array}{*{20}{c}} {{y_{{C}}} - {y_{{A}}}}&{{x_{{A}}} - {x_{{C}}}} \\ {{y_{{B}}} - {y_{{D}}}}&{{x_D} - {x_B}} \end{array}} \right]^{ - 1}}\left[ {\begin{array}{*{20}{c}} {{x_{{A}}}{y_{{C}}} - {x_{{C}}}{y_{{A}}}} \\ {{x_{{D}}}{y_{{B}}} - {x_B}{y_{{D}}}} \end{array}} \right]. $

式中：$ ({x_{{P}}},{y_{{P}}}) $为机构瞬时旋转中心P的坐标. 瞬心的运动轨迹误差采用欧几里得距离误差:

(27)$ {E_{{P}}} = \sum\nolimits_{i = 1}^{{n_{\text{E}}}} {[{{(x_{{P}}^i - x_0^i)}^2}+{{(y_{{P}}^i - y_0^i)}^2}]} . $

式中：$ {E_{{P}}} $为关节瞬时旋转中心的运动轨迹误差，$ {n_{\text{E}}} $为屈曲运动过程中关节瞬时旋转中心的数量，$ (x_{{P}}^i,y_{{P}}^i) $为外骨骼关节运动过程中第i个瞬时旋转中心的坐标，$ (x_0^i,y_0^i) $为人体关节运动过程中第i个瞬时旋转中心的坐标. 以最小化瞬心的运动轨迹误差为目标，对反向双摇杆变瞬心机构进行优化设计，杆长约束（单位：mm）为

(28)$ \left. \begin{gathered} {\text{min }}{E_p}. \\ {\text{s}}{\text{.t}}{\text{. }}20 \leqslant {L_1} \leqslant 80, \; 20 \leqslant {L_2} \leqslant 80, \\ \quad\;\; 20 \leqslant {L_3} \leqslant 80, \; 20 \leqslant {L_4} \leqslant 80. \\ \end{gathered} \right\} $

不同方案的变瞬心机构尺寸如表4所示. 在矢状面上，按照方案三设计的外骨骼瞬时旋转中心轨迹与人体关节屈曲运动的瞬时旋转中心轨迹对比如图7所示，其中x为X向的位移，z为Z向的位移. 膝关节与变瞬心机构瞬时旋转中心的RMSE=0.90 mm、MAE=0.82 mm、ME=1.61 mm. 经过尺寸优化后，机构能够精确拟合人体关节瞬时旋转中心轨迹，提高关节运动的匹配性，实现高精度的人机协同变瞬心运动.

膝关节在负重活动中承受的力和接触力随屈曲角度的增大而增大^[31-32]，以$ \theta $=116.35°的康复外骨骼小腿组件为例，进行有限元分析，验证该小腿组件安全性与可靠性. 外骨骼材料选用铝合金，E=72 GPa，μ=0.33，密度为2 720 kg/m^3[33]. 负载主要来自体质量（约70 kg），在外骨骼与大腿和小腿接触处施加固定约束. 基于赫兹接触理论，通过有限元法计算得到的等效应力$ \sigma $及载荷传递分布如图8所示. 以最小化康复外骨骼小腿组件的应变能为优化目标，引入体积和应力约束，沿力线和载荷传递方向对材料分布进行优化设计，实现结构轻质性设计. 康复外骨骼左小腿组件在拓扑优化过程中的应变能$ {J_{\text{e}}} $随迭代次数I的变化曲线如图9所示. 不同体积约束下拓扑优化迭代过程中结构的应变能最大值和最小值如表5所示，其中J_MAX、J_MIN分别为$ {J_{\text{e}}} $的最大值和最小值. 在拓扑优化过程中，应变能越小，结构的刚度越大，变形越小，能够更好地承受外部载荷. 设计的康复外骨骼可以根据不同个体的需求，在运动轨迹、结构尺寸和材料分布上进行定制化设计，提供个性化和多样化的设计方案. 康复外骨骼左小腿组件样品可在数字光处理打印设备上使用光敏树脂打印而成，如图10所示. 打印精度为55 μm，分层厚度为100 μm.

（1）提出基于锥规划重建的运动关节康复外骨骼定制设计方法. 针对个性化的运动关节康复需求，基于锥规划重建对人体工学的康复外骨骼进行定制设计，有效避免了局部最优解的产生，提高了运动姿态重建的精度. 相比分层匹配方法，锥规划重建方法的均方根误差降低了10.95%，平均绝对误差降低了12.29%，最大误差降低了6.05%. 通过锥规划重建关节的运动姿态定制化设计康复外骨骼，能够提高康复外骨骼运动舒适性.

（2）设计考虑下肢关节变瞬心运动特性的定制康复外骨骼结构. 通过锥规划精确重建个体膝关节的运动姿态变化，获得瞬时旋转中心的精确运动轨迹，定制化设计康复外骨骼，有效匹配人体关节的变瞬心J形曲线运动. 设计的外骨骼反向双摇杆变瞬心机构与人体运动关节的瞬时旋转中心均方根误差为0.90 mm，平均绝对误差为0.82 mm，最大误差为1.61 mm. 在运动域中，定制的康复外骨骼可精确拟合人体膝关节瞬时旋转中心轨迹，提高与个体关节运动轨迹的匹配性.

（3）提出基于载荷传递的结构拓扑优化方法. 通过重建关节运动姿态，利用有限元精确分析应力分布、力线和载荷传递情况，以应变能最小为目标进行拓扑优化，在体积和应力约束下优化结构的材料分布. 提出的优化方法可以适应个性化需求的结构定制设计，实现不同条件下的轻质性设计.

今后，将在研究中进一步探讨不同材料和结构设计对外骨骼承载性的影响. 未来将在结构域、性能域进行综合分析，结合多域性能需求，进一步优化康复外骨骼的设计.