三通路反射控制的人体姿势平衡预测仿真

Three-Pathway Reflex Control for Predictive Simulation of Human Postural Balance

摘要: 反馈控制模型已被用于人类姿势控制的运动模拟研究。本研究基于肌肉骨骼建模与控制，改进了一种三通路反馈控制器，用于预测不同感觉运动时延与扰动幅度下的人体姿势平衡响应。基于改编的gait10dof18musc模型，在100、150和200 ms三种反馈延迟和0 mm~80 mm的不同幅度下对扰动平衡进行了模拟。结果显示，较长时延导致更高的峰值肌肉激活；随着扰动幅度增大，踝与膝关节的活动范围增加而骨盆与髋关节活动范围呈下降趋势；多块关键肌在不同时延间的峰值激活存在显著差异；此外，相较原控制器，改进的反馈控制器在更广的时延范围和扰动幅度内维持了更小的CoM与CoP偏差，四对主要的拮抗肌共收缩指数明显降低。这些发现表明，该控制器在更广时延和扰动范围内鲁棒性更好，在扰动跟踪与姿态稳定方面更为有效。

Abstract: Feedback-control models have been widely used to simulate human postural control. Building on musculoskeletal modeling and control, this study improves a three-path feedback controller to predict human postural responses under varying sensorimotor delays and perturbation amplitudes. Using a modified two-dimensional gait10dof18musc model, we simulated perturbed balance at feedback delays of 100, 150, and 200 ms and amplitudes of 0 mm~80 mm. The results show that longer delays lead to higher peak muscle activations; as perturbation amplitude increases, the ranges of motion at the ankle and knee increase, whereas those of the pelvis and hip decrease; and several key muscles exhibit significant between-delay differences in peak activation. Moreover, relative to the original controller, the improved controller maintains smaller CoM-CoP deviations over a broader range of delays and amplitudes and markedly reduces the co-contraction index for four principal antagonist pairs. These findings indicate that the improved controller achieves greater robustness across delays and perturbations and is more effective in perturbation tracking and postural stabilization.

1. 引言

平衡控制是一个复杂的过程，人体平衡并非意识控制的功能，而是由外部扰动引发的大脑皮层反应与中枢神经系统多种反射机制共同调节的复杂过程，这些机制包括前庭系统、本体感觉系统和视觉系统[1]。这些系统向大脑传递人体位置与运动的感官信息，并协同多种反射和控制机制以维持平衡。维持平衡的基本运动学约束是将质心(Center of Mass, CoM)保持在支撑基底(Base of Support, BoS)范围之内：当CoM超出BoS，机体必须通过快速的感觉–运动反馈将CoM拉回，否则可能会导致跌倒；平衡的另一项关键指标是压力中心(Center of Pressure, CoP)，CoP的运动直接受到CoM的影响[2]。

在人体平衡研究中，倒立摆模型为理解全身运动学与动力学提供了简洁而有效的理论框架[3] [4]。此类基于倒立摆的模型已广泛用于研究人体在矢状面的平衡与步态动力学，成为研究的常用起点。然而，该类模型通常不考虑肌肉和神经肌肉控制。最优控制将姿势调节表述为跨时域最小化代价并结合状态估计，长时域目标与规范性解释强，但它往往依赖较强模型假设与对代价函数的先验设定，映射到肌肉层面时需额外处理冗余与时延，计算代价随非线性与约束复杂度显著上升[5]。模型预测控制(Model Predictive Control, MPC)以滚动优化与显式约束见长，便于前瞻与安全域处理，在人类平衡研究中已用于抗扰与轨迹跟踪。其局限在于对模型精度与扰动预测较敏感，实际应用中常需模型降阶，肌骨冗余与神经传导时延的显式融合并不直接[6]。因此，肌肉骨骼(Musculoskeletal, MSK)建模被用来研究人体平衡神经运动控制。神经肌肉骨骼模型有助于研究反射对平衡的贡献，并研究人体平衡控制之间常见的肌肉协同作用[7]。相关研究表明，反射控制器能重现正常的人类行走[8]，可以被应用于控制外骨骼[9]、假肢[10]和人形机器人[11]，并且得到了扩展。基于上述分析，本文选择采用肌肉骨骼(MSK)建模结合反射控制的路线，因为其在肌肉层面机制解释、时延效应剖析与连续扰动情境的可重复标定上更具针对性。

我们想研究类似的反射模型能在多大程度上解释在站立扰动实验中观察到的行为。此前已有大量研究将本体感觉、前庭、视觉与其他传感器相融合[12] [13]，但各传感系统的个体贡献仍缺乏系统量化。Suzuki等人[14]的研究在脊髓上位系统框架中引入反射控制器以刻画受扰站立的平衡调控，并能表征出一定的人体平衡能力。Suzuki等人使用的时间延迟高达20毫秒，而由于反馈引起的肌肉激活延迟因肌肉而异，对于下肢肌肉，通常在35 ms~100 ms。此外，所得到的运动并没有直接与具有相同扰动的实验运动进行比较。

本文旨在探索感觉反馈时延在人体平衡控制过程中的影响以及如何影响平衡策略。首先介绍一种用于肌肉反射的改进的三通路反馈控制器方法，继而，我们给出相应的优化建模方法。最后，我们将在不同的扰动条件和3个CoM反馈时延条件下，系统展示正弦平移扰动下的关节角度、力矩以及肌肉激活、CoM与CoP相对运动。同时阐明CoM反馈时延如何耦合扰动幅度共同影响肌肉协同与策略选择。

2. 肌肉骨骼模型建立与肌肉控制系统

本研究采用一个通用二维肌肉骨骼(MSK)模型，该模型基于OpenSim中的gait10dof18musc模型改编而成，其结构如图1所示。该模型体重74.5 kg，身高1.8 m，包含9个自由度和18块Hill型肌肉，所有肌肉均在矢状面内进行驱动。模型中增设了一个质量为10 kg的运动平台。通过控制该平台，可生成前后方向的正弦运动(幅值0 mm~80 mm)。采用1 HZ的扰动频率用于连续扰动方式的平衡研究。足部与平台之间的接触采用Hunt-Crossley接触力模型，在双脚脚趾和脚跟处均设置了接触球体，平台表面则覆盖一个接触半空间。接触力建模所用参数见表1。

Table 1. Contact coefficients for the model’s contact forces

表1. 模型接触力的接触参数

参数	值
刚度	2 × 10^6N/m^2
耗散系数	1
静摩擦系数	0.9
动摩擦系数	0.6
粘性摩擦系数	0.6

2.1. 三通路反馈控制器

Figure 1. (a) Overview of the musculoskeletal model used in the simulations, comprising nine muscles (each labeled); the platform is controller-driven and can impose sinusoidal perturbations of varying amplitudes; (b) Frontal view of the adapted human musculoskeletal model; (c) Lateral view of the adapted human musculoskeletal model

图1. (a) 模拟使用的肌肉骨骼模型整体图，包括9块肌肉(每块肌肉贴有标签)，平台由控制器驱动，可施加不同振幅的正弦扰动；(b) 改编的人体肌肉骨骼模型正面；(c) 改编的人体肌肉骨骼模型侧面

MSK模型采用三向反馈控制方法进行控制。模型肌肉的激活是通过三个不同的反馈控制器计算的：肌肉长度–速度反馈控制器、肌肉力量反馈控制器和CoM反馈控制器。使用多个控制器来平衡模型是为了尽可能真实地模拟人体在维持平衡时依赖的多层次、多来源的神经控制机制。

2.1.1. 肌肉长度–速度反馈控制器

这个控制器基于肌肉纤维长度的变化来激活肌肉。当肌肉被拉伸到一定长度时，这个控制器通过激活相应的肌肉来防止过度拉伸，从而保持平衡。肌纤维长度阈值和激活延迟基于先前的研究[15] [16]。这种控制方式模仿了肌肉纺锤体的反射机制。控制器使用以下方程式计算激活：

$\begin{array}{c}{u}_{LV}=K_L\left(l\left(t-t_{del}\right)-l_0\right)+K_V\dot{l}\left(t-t_{del}\right)\end{array}$ (1)

其中，$u_{LV}$ 表示由肌肉长度–速度反射引发的肌肉兴奋信号，$l$ 为当前归一化肌纤维长度，$l_0$ 表示长度反馈偏移量，$\dot{l}$ 表示肌肉实时收缩速率，$K_L$ 为肌肉长度增益系数，$K_V$ 为肌肉收缩速度增益系数，$t$ 表示当前时间，$t_{del}$ 表示第m个肌肉的反射时间延迟(m为肌肉编号)，其取值范围为35 ms至100 ms。表2汇总了下肢主要肌肉的时间延迟。

Table 2. Muscle-specific delays employed in the controller

表2. 控制器中使用的肌肉延迟

肌肉	时间延迟(ms)
胫骨前肌	35
比目鱼肌	44
腓肠肌	54
股二头肌	60
髂腰肌	40
臀肌群	40
阔肌群	100
股直肌	82

2.1.2. 肌肉力量反馈控制器

肌肉力量反馈控制器基于肌肉的力量反馈来激活肌肉，用于模拟人体中的协同收缩现象，这是一种重要的稳定机制。肌肉的力量变化是平衡控制的关键因素之一，该控制器通过检测某块肌肉的力输出，进而激活相应的拮抗肌，以平衡肢体两侧的力生成。该控制器通过以下公式计算肌肉激活输出：

$\begin{array}{c}{u}_F=K_F{F}_{ant}\left(t-t_{del}\right)\end{array}$ (2)

其中，$F_{ant}$ 是拮抗肌的肌力，$K_F$ 为肌肉力量反馈控制的增益系数。为有效实现肌力反馈机制，本研究已建立多组主动肌–拮抗肌对，并汇总于表3中。

Table 3. Muscle co-contraction pairs

表3. 肌肉共收缩对

主动肌	拮抗肌
臀大肌与腘绳肌群	髂腰肌
腘绳肌群与股二头肌	股直肌与股中间肌
腓肠肌与比目鱼肌	胫骨前肌

2.1.3. CoM反馈控制器

CoM反馈控制器根据模型全身质心(CoM)相对于CoP在矢状面内的位置、速度和加速度信息激活相应肌肉。该控制器的作用在于维持模型质心始终位于支撑基底(BoS)范围内，模拟人体前庭系统在平衡控制中的功能。基于CoM反射控制器所确定的肌肉激活由以下公式定义：

$\begin{array}{c}{u}_{COM}=K_{CP}\left[p_{COP}\left(t-t_D\right)-p_{COM}\left(t-t_D\right)\right]+K_{CV}\left[p_{\dot{C}OP}\left(t-t_D\right)-p_{\dot{C}OM}\left(t-t_D\right)\right]\\ \begin{array}{c}+K_{CA}\left[p_{\dot{C}\dot{O}P}\left(t-t_D\right)-p_{\dot{C}\dot{O}M}\left(t-t_D\right)\right]\end{array}\end{array}$ (3)

其中，$p$ 表示CoM和CoP的位移，$\dot{p}$ 表示速度，$\ddot{p}$ 表示加速度，$K_{CP}$ 、$K_{CV}$ 和$K_{CA}$ 分别表示位移、速度和加速度的反馈控制增益系数，$t_D$ 为CoM反馈延迟时间。本文所采用的延迟时间范围参考McKay等的研究[17]，该研究提出了一个感觉运动反馈模型，成功复现了健康青年、老年人群、帕金森病患者在实验中的肌肉激活模式。与Jones等人[18]研究中采用的公式不同，本文公式(3)新增$K_{CA}$ 加速度项，利用相位超前抵消时延带来的相位损失。

2.2. 目标函数

在每次仿真中，均通过优化算法确定一组最优参数(包括$K_F$ 、$K_L$ 、$K_{CP}$ 、$K_{CV}$ 、$l_0$ )，以使目标函数最小化。目标函数定义如下：

$\begin{array}{c}J=w_{effort}\underset{0}{\overset{t_{term}}{{\displaystyle \int }}}\left[{\displaystyle \sum }_{m=1}^{18}{a}_m{(t)}^2\right]dt+w_{fall}\left(1-\frac{t_{term}}{t_{cnd}}\right)+w_{dof}{J}_{dof}\\ +w_{stab}\underset{0}{\overset{t_{term}}{{\displaystyle \int }}}{\Vert p_{COP}\left(t\right)-p_{COM}\left(t\right)\Vert }^{2}dt\end{array}$ (4)

该目标函数旨在最小化肌肉激活的平方和以降低能耗，同时引入多项稳定性惩罚项：当质心(CoM)高度低于初始高度的0.9倍时仿真终止以模拟跌倒；并对矢状面上膝关节角度 > 30˚、骨盆倾斜 > ±15˚、最小化CoM与CoP的距离的平方以及骨盆前后平移 > ±0.2 m的情况进行惩罚。其中，$w$ 为权重系数($w_{\text{fall}}=100$ ，$w_{\text{effort}}=w_{\text{dof}}=1$ ，$w_{\text{stab}}=1$ )，$t_{cnd}$ 为仿真总时长，设为10 s，$t_{\text{term}}$ 为仿真终止时间，当人体质心(CoM)高度低于初始高度0.9倍时提前终止(如果模型未曾跌落至初始0.9倍高度，则$t_{\text{term}}=t_{cnd}$ )，$J_{dof}$ 表示关节角度超出预设范围时所施加的惩罚项，用于表征失衡趋势：

$\begin{array}{c}{J}_{dof}=\underset{0}{\overset{t_{\text{term}}}{{\displaystyle \int }}}\left(P\left(\text{knee}\left(t\right)\right)+P\left(\text{pelvistilt}\left(t\right)\right)+P\left(\text{pelvistran}\left(t\right)\right)\right)dt\end{array}$ (5)

其中，惩罚函数$P\left(\text{knee}\left(t\right)\right)$ 、$P\left(\text{pelvistilt}\left(t\right)\right)$ 、$P\left(\text{pelvistran}\left(t\right)\right)$ 的取值规则如下：若膝关节屈曲角度 > 30˚、骨盆倾斜角度 > ±15˚或骨盆平移距离 > ±20 cm，则相应函数取值为10，否则为0。

所有优化过程均采用SCONE软件的商用Hyfydy版本以提高计算效率。该软件使用协方差矩阵自适应进化策略(Covariance Matrix Adaptation Evolution Strategy, CMA-ES)对模型参数进行优化。在迭代过程中，每种仿真情景均使用不同随机种子并行运行60次优化，具有更优参数的解将被保留作为下一次优化的起始参数，直至获得最优参数组合。为提高效率，每次扰动幅度优化得到的最优结果将作为下一更高幅度优化初始值。单个情景的优化通常可在一小时内收敛。

2.3. 数据处理

在数据处理过程中，首先将关节角度和肌肉激活数据划分为每秒一个周期，然后计算所有周期数据的均值和标准差。为了消除模拟初期肌肉激活的不稳定性，分析中剔除了前两秒的数据，因为在模拟的开始阶段，肌肉激活变化较大，较为波动。而当模型进入更加稳定的姿态模式后，数据的波动性减小，因此排除初始姿势对结果的影响，可以更准确地反映模型的稳定后期状态。

采用Rudolph [19]等人提出的公式计算拮抗肌对的共收缩指数CCI (Co-Contraction Index)：

$\begin{array}{c}CCI\left(t\right)=\frac{EM{G}_L\left(t\right)}{EM{G}_H\left(t\right)}\left(EM{G}_L\left(t\right)+EM{G}_H\left(t\right)\right)\end{array}$ (6)

其中，$EM{G}_L\left(t\right)$ 表示同一时刻拮抗肌对中较小的肌肉激活，$EM{G}_H\left(t\right)$ 表示同一时刻拮抗肌中较大的肌肉激活。为了方便比较改进后的控制器与原控制器在同幅度同时延下的相同拮抗肌组之间的CCI，我们整合了CCI，并给出了最终标量化后的CCI。相关研究表明[19] [20]，CCI可能与关节稳定性有关，CCI越大，关节稳定性越低。

对均值和峰值肌肉激活数据进行了配对t检验，比较了不同质心反馈延迟组之间的差异，且对每个幅度进行分析。所有统计分析均在MATLAB中进行，显著性水平设定为0.05。

3. 结果与讨论

3.1. 关节角度

Figure 2. Time-normalized joint angles of the pelvis, hip, knee, and ankle over the cycle during an anteroposterior (AP) horizontal sinusoidal perturbation (amplitude ≥ 10 mm), compared across three delay conditions (100 ms, 150 ms, and 200 ms)

图2. 不同时延(100 ms、150 ms、200 ms)条件下，AP方向水平正弦扰动(幅值 ≥ 10 mm)周期内的关节(盆骨、髋、膝、踝)角度

图2展示了在所有前后(AP)方向水平正弦扰动(10 mm至80 mm)和所有质心反馈延迟条件下的关节角度数据。扰动由10 mm增至80 mm时，关节角度响应呈现出清晰的幅值依赖性与时延相关性。总体而言，扰动由10 mm增至80 mm时，各关节的峰摆幅显著扩大，踝与髋对幅值最为敏感；小幅扰动以“踝策略”为主，骨盆与髋部摆动较小；当扰动达到中大幅(≥ 40 mm)或反馈时延加长时，髋、膝参与度明显增强，曲线出现更大的相位滞后与波形畸变(如双峰、拐点)。同时，骨盆角度随幅值与时延增加而摆幅上升，反映近端躯干稳定性需求增大。整体上，这种非线性幅值放大与相位漂移表明：较大扰动与较长时延会放大关节波动、降低周期一致性与稳定裕度。

3.2. 肌肉激活度

多数肌肉的激活幅值随扰动幅度单调增加并呈非线性放大，如图3所示。单关节肌肉如臀大肌和股二头肌由于激活水平较低，所以不包括在内。质心反馈延迟似乎会影响峰值肌肉激活。这个效应在前侧肌肉(股直肌、股中间肌和胫前肌)中尤为明显。此外，腓肠肌在较大的延迟时间下表现出更高的激活，而比目鱼肌则呈现出更为混乱的模式。与此相反，股直肌、腓肠肌和股后肌群等双关节肌肉则表现出较高的激活。对平均肌肉激活的统计分析显示，在所有时间延迟情况下，臀大肌、髂腰肌、股直肌、比目鱼肌和胫前肌这些肌肉的峰值激活存在显著差异(p < 0.05)。

Figure 3. Time-normalized muscle activations over the cycle during an anteroposterior (AP) horizontal sinusoidal perturbation (amplitude ≥ 10 mm), compared across three feedback-delay conditions (100 ms, 150 ms, and 200 ms)

图3. 不同时延(100 ms、150 ms、200 ms)条件下，AP方向水平正弦扰动(幅值 ≥ 10 mm)周期内的归一化肌肉激活

如图4所示，在100 ms时延下，改进的控制器在大多数扰动幅度与拮抗肌对上均表现出更低的CCI (标量及统一化后的值)，尤其在踝关节拮抗肌对与大幅度扰动时下降最为明显，一定程度上可以说明在相同稳定性需求下能以更少的共收缩完成姿态调节，具有更高的控制效率与潜在能耗优势。但同时在较小幅度与股直肌–腘绳肌对中，CCI也有升高。在所有时延与扰动幅度下，改进后的控制器CCI下降的拮抗肌对占全部拮抗肌对比值为61.7%。表明改进的控制器能更好地抑制不必要的共收缩。

Figure 4. Comparison of scalar CCI between the original and improved controllers under a 100-ms delay and perturbation amplitudes > 10 mm

图4. 100 ms时延与大于10 mm幅度扰动条件下，原控制器方法与改进的控制器方法之间的标量化CCI的对比

3.3. CoM和CoP

质心(CoM)与足底压力中心(CoP)的差异可以作为衡量个体稳定性的一种指标[21]。当CoP与CoM之间的差异接近零时，表示较高的稳定性；而差异越大，则表示稳定性较差。如图5所示，足底压力中心在前后方向上的范围大于质心的范围，并且这一范围随着扰动幅度增大而增大。图6展示了在各个反馈延迟条件下，对于较大幅度扰动(大于10 mm)，根据CoP与CoM的差值计算的平均稳定性。总体而言，较大扰动幅度叠加较长时延条件下，CoP和CoM的差值会显著增大，表明姿态稳定性下降且维持姿势稳定的难度提升。然而，在100 ms延迟条件下，不同扰动幅度之间的差异较小。

如图7所示，在绝大多数“时延与幅度”组合下，改进控制器的单周期CoP与CoM偏差均方根低于原控制器。随扰动幅度与反馈时延增大，原控制器的CoP与CoM偏差的RMS值呈上升趋势，而改进控制器的增长更缓，显示出改进控制器在高幅度与长时延情形下的鲁棒性优势。综合来看，改进控制器在更广的时延范围和扰动幅度内维持了更小的CoM与CoP偏差，表明其在扰动跟踪与姿态稳定方面更为有效。

Figure 5. Horizontal displacement between the center of mass (CoM) and the center of pressure (CoP) during anterior-posterior (AP) sinusoidal perturbations (amplitude > 20 mm), under a 100 ms feedback delay

图5. 100 ms时间延迟下，AP方向水平正弦扰动(振幅 ≥ 20 mm)的CoM与CoP的水平位移

Figure 6. COM-COP horizontal displacement difference during anterior-posterior sinusoidal perturbations (amplitudes ≥ 10 mm) under feedback delays of 100, 150, and 200 ms. The x-axis denotes cycle progress (%), and the y-axis shows the COM-COP difference (m)

图6. 不同时延(100、150、200 ms)条件下，AP方向水平正弦扰动(振幅 ≥ 10 mm)周期内的COM-COP水平位移差。横轴为周期进度(%)，纵轴为COM与COP的差值

Figure 7. Single-cycle CoP-CoM difference across feedback delays and perturbation amplitudes. RMS denotes the root-mean-square of the CoP-CoM difference within the cycle; RMS1 corresponds to the original controller, and RMS2 to the improved controller

图7. 不同时延与扰动幅度下单周期CoP与CoM差值对比，RMS表示该周期内CoP与CoM差值的均方根值，RMS1代表原控制器方法，RMS2代表改进的控制器方法

4. 总结

在Jones等人[18]提出的三通路反射的平衡控制框架上，本研究在原有的肌肉长度反馈控制器中引入了速度项，同时将CoM通路由“位置–速度”扩展为“位置–速度–加速度”，并在优化目标中显式加入CoM-CoP位置差的二次惩罚，以稳定性指标直接约束求解。理论上，加速度项为连续正弦类扰动提供额外的相位超前，有助于在存在感觉运动时延时更快速地校正CoM偏移。在不同幅度(0 mm~80 mm)与3个CoM反馈时延(100/150/200 ms)的系统仿真中，我们观察到：1) 关节响应呈现典型的“幅值依赖与相位滞后”特征。较长的时延与较大的幅度会导致髋、膝摆幅增大、参与显著增强，整体周期一致性与稳定裕度下降；2) 多数肌肉的峰值激活随扰动幅度单调且非线性上升，延迟对前侧肌峰值影响更敏感，比目鱼肌呈现更复杂模式。此外，相较于仅含“位置–速度”的原始CoM控制器，加入加速度反馈与稳定性惩罚后，模型在更广的时延范围内维持了较小的CoM-CoP偏差，并在不牺牲跟踪能力的前提下降低了对拮抗肌共收缩的依赖，关节稳定性维持更高效。这些趋势与既往延迟导致肌激活升高、策略占比改变的结论一致，改进后的控制器在延迟环境下具备更好的鲁棒性与生理可解释性。本研究仍保留了肌长反馈与反射通路的并联兼容性，使改进的CoM控制器可与其他控制器组合，以探索年龄与病理对姿势平衡的影响。总而言之，改进后的三向控制器在保持模型简洁性的同时提升了延迟条件下的稳定性与效率，为与更复杂神经控制模型的对比研究提供了稳健的计算基础。