复杂装备具有结构组成复杂、尺寸大、耦合紧密等特点，在其维修过程中常须作业人员长时间手持工具过头作业，易引起上肢及肩部肌肉酸痛、僵硬，不仅导致作业时间延长与维修效率降低，长期还可能引发与工作相关的肌肉骨骼损伤（work-related musculoskeletal disorders，WMSDs）^[1]，这已成为我国职业病防控的重点问题之一. 因此，明确过头作业环境下的上肢肌肉疲劳特性，对于提高装备维修效率、降低作业人员的WMSDs风险及保护作业人员健康具有重要现实意义^[2].

人体疲劳包括生理上的疲劳和心理上的疲劳，针对上肢作业一般结合主、客观2种方法进行肌肉性疲劳、心理性疲劳、脑力性疲劳等方面的研究. 主观方法常用Borg量表，如Frasie等^[3]证实了Borg CR10量表在提升作业中评估肩部疲劳的有效性；Chin等^[4]将Borg RPE量表与NOISH手工提举分析、北欧身体地图问卷相结合对钣金工业中的提举作业疲劳进行评估. 客观测量方法包括表面肌电(surface electromyography，sEMG)、心率、呼吸频率和血压等生理数据的采集分析，以量化分析工作负荷和疲劳程度. 其中，sEMG技术以其高信噪比和对环境因素的低敏感性而被广泛应用，如Zhou等^[5]基于sEMG研究了战斗机飞行员任务执行中的上肢肌肉疲劳状况及与作业效能的关系；徐兆等^[6]结合sEMG与动作捕捉监测上肢运动疲劳；La Delfa等^[7]利用sEMG评估不同手部垂直施力高度对肩部肌肉疲劳反应的影响，发现过头作业使肩部力量最大降低18%；Yang等^[8]利用sEMG和上肢运动学数据评估使用坐、立两用座椅作业时的上肢疲劳改善；Musso等^[9]研究了使用上肢外骨骼进行过头装配作业、砌砖和移箱作业时的上肢sEMG变化，发现肩屈肌激活减少，而非过头作业时伸肌激活增加.

以上研究探讨了不同作业场景下的上肢疲劳特性，发现其中过头作业更易增加上肢WMSDs风险. 然而，针对过头作业所致上肢肌肉疲劳的有效分类及其非线性动态演变规律的解析，仍须进一步探索. 当前，有监督学习的机器学习方法，如支持向量机（support vector machine, SVM）、k最近邻(k-nearest neighbor, KNN)、随机森林（random forest，RF）等，仍须人为对肌肉疲劳分类^[10-11]. 以K-means聚类、层次聚类为代表的无监督学习虽不依赖标注数据，但前者随机的初始划分影响聚类结果且要求指定聚类数量^[12]，后者则对异常值和数据点密度波动较为敏感^[13]. 另外，仅通过分类难以捕捉肌肉疲劳的非线性动态特征，而Hammerstein模型则能够有效描述非线性动态过程，因而被广泛应用^[14-15]. Li等^[15]提出基于Hammerstein结构的动态神经网络模型，应用于上肢神经康复的功能性电刺激中，实现肘关节运动的精确模拟和跟踪；Zhang等^[16]针对帕金森综合征患者手腕震颤抑制性能受限问题，建立了腕部肌骨系统的Hammerstein模型，提出了自动、快速的参数辨识方法并设计了震颤抑制控制方案；杨坤等^[17]基于Hammerstein系统提出航空发动机部件级辨识模型，以预测关键参数并描述航空发动机的动态特性.

针对上述问题，本研究融合Hammerstein模型和sEMG提出新的过头作业上肢肌肉疲劳评估方法. 以sEMG的中位频率(media frequency，MF)和均方根值(root mean square，RMS)作为Hammerstein模型的输入和输出，辨识动态参数，借助K-means++聚类算法对肌肉疲劳分类评估，同时利用肌电疲劳阈值（electromyographic fatigue threshold，EMGFT）确定具体肌肉疲劳时间. 设计过头作业实验并采集上肢sEMG数据，通过关键参数辨识与疲劳分类深入分析上肢肌肉疲劳的动态变化规律.

为了从复杂的上肢肌肉sEMG信号中获得关键信息，须进行sEMG特征提取，常用特征包括时域特征和频域特征.

时域特征常用的有积分肌电值(integrated electromyography，IEMG)与RMS. IEMG表达式如下：

(1)$ {\rho _{{\text{IEMG}}}} = \sum\limits_{i = 1}^N {\left| {{x_i}} \right|} . $

式中：N为sEMG总长度，x_i为第i个窗口的sEMG幅值.

RMS表达式如下：

(2)$ {\rho _{{\text{RMS}}}} = \left[\frac{1}{N}\sum\limits_{i = 1}^N {x_i^2} \right]^{1/2}. $

通常随着肌肉疲劳加深，IEMG幅值会增加，而RMS对信号波动和瞬时变化响应较平滑，适用于分析疲劳过程中肌肉活动的变化趋势^[18].

频域特征主要涉及平均功率频率(mean power frequency，MPF)与MF. MPF表达式如下：

(3)$ {\gamma _{{\text{MPF}}}} = \sum\limits_{j = 1}^M {{f_j}{P_j}} \bigg/\sum\limits_{j = 1}^M {{P_j}} . $

式中：P_j为第j个频率点对应的功率谱密度，M为功率谱中的总频率点数，f_j为第j个频率点对应的频率.

MF的取值$\gamma _{{\text{MF}}} $表达式如下：

(4)$ \sum\limits_{j = 1}^{{\gamma _{{\text{MF}}}}} {{P_j}} = \sum\limits_{j = {\gamma _{{\text{MF}}}}}^M {{P_j}} = \dfrac{1}{2}\sum\limits_{j = 1}^M {{P_j}} . $

在肌肉疲劳过程中，MPF和MF都表现为下降趋势^[19]；但MF相较于MPF，在监测疲劳期间受噪声影响较小，当信噪比小于10 dB时MPF误差过高，而MF的误差稳定在2~10 Hz，更适于跟踪肌肉疲劳^[20].

肌肉疲劳主要原因之一是力竭性肌肉活动导致肌肉收缩动力学变化. 为了进一步分析肌肉疲劳的动力学过程，引入Hammerstein模型，融合sEMG特征深入分析过头作业上肢肌肉疲劳的非线性规律. Hammerstein模型由2部分组成（见图1）：无记忆的静态非线性模块和动态线性模块. 输入经无记忆线性函数加和转化为中间变量，结合动态线性模块及外界干扰因素输出特定结果^[21].

Hammerstein模型的输入转化如下：

(5)$ \begin{split} v(t) =& F(u(t)) = \sum\limits_{l = 1}^p {{c_l}{f_l}[u(t)]}= \\ & {c_1}u(t)+{c_2}{u^2}(t)+\cdots +{c_p}{u^p}(t)\end{split} . $

式中：u(t)为t时刻系统输入，v(t)为非线性部分输出，p为基函数阶数，c_l为第l个基函数的系数，$f_l(\cdot) $表示第l阶基函数.

图1的动态线性模块表达式如下：

(6)$ A(z)Y(t) = B(z)v(t)+e(t) . $

式中：Y(t)为系统输出，e(t)为均值为0的高斯白噪声.

(7)$ \left.{\begin{array}{*{20}{l}} {A(z) = 1+\displaystyle \sum\limits_{j = 1}^m {{a_j}{z^{ - j}}} }, \\ {B(z) = \displaystyle \sum\limits_{k = 1}^n {{b_k}{z^{ - k}}} } ,\\ {G(z) = B(z)/A(z)}. \end{array}} \right\} $

式中：m、n分别为模型的输出、输入滞后项数目； a_j为第t−j个滞后值Y(t−j)的系数；b_k为第t−k个滞后值v(t−k)的系数；z^−j为滞后算子，并使得z^−jY(t)=Y(t−j).

联合式(5)~(7)，可得

(8)$ Y(t) = \sum\limits_{j = 1}^m {{a_j}Y(t - i)+\sum\limits_{k = 1}^n {\sum\limits_{l = 1}^p {{b_k}{c_l}{u^l}(t - i)+\xi (t)} } } . $

式中：ξ(t)为外界干扰. a、b、c为参数，a用以描述不同阶数的历史输出Y(t−i)对Y(t)的影响；b为非线性静态部分，将u(t)转换为v(t)，并进一步影响Y(t)；c反映不同阶次u(t)对v(t)的影响.

进而将式(8)简化为

(9)$ Y(t)={{\boldsymbol{\theta}}}^{\mathrm{T}}(t) {{\bf\textit{φ}}}+\xi(t) . $

其中，

(10)$ \begin{split} {\bf\textit{φ}}=&\left[a_1, a_2, \cdots, a_m,\left(b_1 c_1\right)_{m+1},\left(b_1 c_2\right)_{m+2},\cdots,\right. \\& \left.\left(b_1 c_p\right)_{m+p}, \cdots,\left(b_n c_p\right)_{m+p n}\right]^{\mathrm{T}},\end{split}$

(11)$ \begin{split} & {\boldsymbol{\theta }}(t) = [Y(t - 1),Y(t - 2),\cdots ,Y(t - m), \\ &\qquad {u(t - 1),{u^2}(t - 1),\cdots ,{u^p}(t - 1),\cdots ,{u^p}(t - n)]^{\mathrm{T}}}.\end{split}$

式中：${\boldsymbol{\theta}}$(t)为历史输出和输入组成的信息矩阵；${\bf\textit{φ}} $为a、b、c组成的参数矩阵；(b₁c₁)_m+1, (b₁c₂)_m+2, $\cdots $, (b₁c_p)_m+p表示b₁与c₁,c₂,$\cdots $,c_p的乘积组合，以此类推，(b_nc₁)_m+p(n−1)+1, (b_nc₂)_m+p(n−1)+2, $\cdots $, (b_nc_p)_m+pn表示b_n与c₁, c₂, $\cdots $, c_p的乘积组合.

采用带遗忘因子的递推最小二乘法（forgetting factor recursive least squares，FFRLS），通过降低历史数值的信息贡献，增加新数据的信息，以解决数据过载问题^[22]. 表达式如下：

(12)$\left.\begin{array}{l}\boldsymbol{K}(t)=\dfrac{\boldsymbol{P}(t-1) \boldsymbol{\theta}(t)}{\lambda+\boldsymbol{\theta}^{\mathrm{T}}(t) \boldsymbol{P}(t-1) \boldsymbol{\theta}(t)}, \\\boldsymbol{P}(t)=\dfrac{1}{\lambda}\left[\boldsymbol{I}-\boldsymbol{K}(t) \boldsymbol{\theta}^{\mathrm{T}}(t)\right] \boldsymbol{P}(t-1), \\{\bf\textit{φ}}(t)={\bf\textit{φ}}(t-1)+\boldsymbol{K}(t)\left[Y(t)-\boldsymbol{\theta}^{\mathrm{T}}(t) {\bf\textit{φ}}(t-1)\right].\end{array} \right\} $

式中：K(t)为卡尔曼增益矩阵，P(t)为协方差矩阵，$ {{\bf\textit{φ}}}(t) $为参数估计矩阵，λ为遗忘因子.

利用FFRLS辨识出参数a，参数b、c则用奇异值分解(singular value decomposition，SVD)方法辨识^[23-24]. 建立增广矩阵A，SVD分解如下：

(13)$ {\boldsymbol{A}} = {\boldsymbol{UQ}}{{\boldsymbol{V}}^{\text{T}}} . $

式中：U和V为不同阶的酉矩阵，Q为对角矩阵.

(14)$ {\boldsymbol{C}} = \frac{1}{{{{V}}(1,1)}} \cdot {\boldsymbol{V}}(:,1) \text{，} $

(15)$ {\boldsymbol{B}} = {{V}}(1,1) \cdot {{Q}}(1,1) \cdot {\boldsymbol{U}}(:,1) . $

式中：C、B分别为参数c、b的向量矩阵.

为了进一步评估肌肉疲劳状态，利用Hammerstein模型识别的参数分类检测肌肉疲劳程度，但当前缺乏对肌肉疲劳的多等级分类研究^[25]. 为此，本研究利用K-means++聚类算法进行过头作业上肢肌肉疲劳分类，通过改进的初始质心选择，避免K-means算法对初始值敏感的现象，从而提高聚类结果质量. 具体步骤如下：随机选择初始质心，计算每个数据点到已选质心的最小欧氏距离，选中下一个质心的概率与该距离的平方成正比. 重复此过程直至选出K个初始质心. 之后利用标准K-means算法，将数据点分配至最近质心并更新质心位置，直至算法收敛从而完成聚类^[26]. 最佳聚类数由肘部法则确定，在聚类过程中，随着K增大，误差平方和(sum of squares due to error，SSE)逐渐降低，当出现明显拐点则可确定为最佳聚类数.

为了量化a、b参数的变化规律，建立其与疲劳分类的关系，引入参数特征矩阵${\boldsymbol{\varPsi }} = [{\boldsymbol{\varGamma }},{{\boldsymbol{K}}_{\boldsymbol{\varGamma }}},\Delta {\boldsymbol{\varGamma }}]$，其中，

(16)$ \left. {\begin{array}{*{20}{l}} {{\boldsymbol{\varGamma }}(i,:) = [{a_1},{a_2}, \cdots ,{a_m},{b_1},{b_2},\cdots ,{b_n}]} ,\\ {{{\boldsymbol{K}}_{\boldsymbol{\varGamma }}}(i,:) = \dfrac{{{\boldsymbol{\varGamma }}(i,:) - {\boldsymbol{\varGamma }}(i - 1,:)}}{{{\boldsymbol{\varGamma }}(i - 1,:)}}}, \\ {\Delta {\boldsymbol{\varGamma }}(i,:) = {\boldsymbol{\varGamma }}(i,:) - \displaystyle \sum\limits_{t = \max \;(m,n)}^i {{\boldsymbol{\varGamma }}(t,:)} /i} .\end{array}} \right\} $

式中：${\boldsymbol{\varGamma }}(i,:)$、${{\boldsymbol{K}}_{\boldsymbol{\varGamma }}}(i,:)$、$\Delta {\boldsymbol{\varGamma }}(i,:)$分别为i时刻a、b参数的增广矩阵、变化率矩阵以及参数值与其瞬时平均值的偏差矩阵.

为了更准确地确定各疲劳阶段的具体疲劳情况，采用改进的肌电疲劳阈值（electromyographic fatigue threshold，EMGFT）方法. 与以往依靠RMS单一特征不同，本研究利用Hammerstein模型参数的量化矩阵，将每个疲劳阶段的动态参数矩阵${\boldsymbol{\varPsi }}$分段拟合，计算各段拟合直线交点，如图2所示. 其具体步骤如下：1）将数据分成2部分，前15个数据点为S₁，剩余为T₁；从第16个数据点起，依次形成S₂、T₂等组，直至T_n与S₁数据点数相同. 2）对每组数据进行一阶线性回归，生成2条拟合直线并计算斜率乘积kS_n·kT_n. 3）选择斜率乘积最大值及其对应的2条直线交点，定义为EMGFT^[27].

选择10名身体健康、年龄为（23±2）岁的志愿者，身高为（174±2.5）cm，体重为（65±5）kg. 所有志愿者均无上肢肌肉骨骼疾病史，均为右利手. 实验前3个月内无肌肉骨骼损伤，且实验前72 h内无剧烈运动. 被试者在实验前都阅读实验指导书并签署同意书. 本实验得到长安大学学术委员会批准.

利用Anybody进行过头作业仿真，发现上肢平均肌肉激活程度大于0.1的肌肉为右臂的肱桡肌（m₁）、肱二头肌（m₂）、三角肌（m₃）和斜方肌（m₄），及左臂的斜方肌（m₅）、肱二头肌（m₆）、三角肌（m₇）和肱桡肌（m₈）. 因此，选取这8块肌肉进行过头作业实验采集sEMG.

主要设备如下：电动螺丝刀、采样频率为1 000 Hz的意大利BTS Bioengineering公司的BTS FreeEMG 1000无线表面肌电系统、杭州迅达无线电器材有限公司生产的一次性银/氯化银(Ag/AgCl)电极片. 过头作业平台如图3所示.

被试者在测量最大自主收缩(maximum voluntary contraction，MVC)前进行肌肉放松，使用磨砂膏和酒精棉擦拭清洁m₁~m₈肌腹表面皮肤，沿肌纤维方向黏贴电极片，并按表1要求完成MVC测量. 每块肌肉测试3次，依次测量左右m₁~m₈肌肉. 每2次测试间休息20 min，以避免肌肉疲劳影响结果，最终取3次测量的平均值作为该肌肉的MVC.

首先，让被试者熟悉Borg CR10主观疲劳等级表(见表2)直至能完全记住各等级含义. 其次，将室温调至26 ℃，被试者进行1 min上肢热身. 然后，按照MVC实验时的相同方法清洁上肢皮肤和黏贴电极片. 最后，被试者站立在工作凳上，右手持内六角电动螺丝刀，按顺时针拆卸图4（b）区域螺栓. 实验中，持续采集8通道sEMG数据并观察采集状况，每隔30 s询问被试者并记录一次疲劳等级. 观察被试者工作状态，如有电极片脱落及时备注、休息并重新实验. 实验过程如图4所示.

将采集的sEMG经50 Hz陷波滤波器去除工频信号干扰，再用20~500 Hz带通滤波器去除噪声保留主要sEMG，进而使用小波阈值去除通道间串扰及高斯噪声. 如图5所示为m₄降噪前、后对比. 其中，Amp表示幅值.

将被试者m₁~m₈的sEMG幅值除以其对应的MVC，对sEMG标准化，以消除个体差异和肌肉激活强度影响，确保sEMG的可比性.

肌肉贡献率是在肌肉群中选择合理肌肉的重要指标，反映肌肉在特定动作中的参与程度. 较高的贡献率意味着该肌肉承担更多负荷，其计算方式为所采集某块肌肉的IEMG占所有采集肌肉IEMG总和的百分比^[28]. 设置时间窗与滑动窗为5 000 ms，并对处理后的sEMG利用滑动时间窗获得每个通道的IEMG，进而得到8个通道肌肉贡献率，如表3所示. 其中P_m为肌肉贡献率. 可以看出，m₄的平均肌肉贡献率最大，为21.28%，因此选择m₄作为目标肌肉进行分析.

以IEMG、RMS、MF、MPF两两组合为自变量，以Borg CR10评分结果为因变量，建立多元回归模型，不同特征组合的判定系数R²如表4所示. 可知，所有被试者的RMS和MF组合的R²最大. 同时，考虑RMS可视为肌肉激活强度和能量变化指标^[29]，而MF相较于MPF更适于追踪疲劳过程^[20]，因此选择RMS和MF作为后续研究的特征组合.

由于Hammerstein模型的输出依赖历史数据且具有记忆性，而输入具有无记忆性，选用MF和RMS分别作为Hammerstein模型的输入与输出进行建模，以描述疲劳过程中频域变化对肌肉激活水平的影响. 为了确保模型的稳定性和有效捕捉动态变化，对Hammerstein模型阶数利用最小描述长度（minimum description length，MDL）准则计算^[30]，当MDL取最小值时，m=3、n=4、p=3为最佳阶数组合. 根据文献[31]和[32]，有λ∈[0.90,1.00]，当模型残差平方和最小时，取λ=0.996. 通过式(12)~(15)辨识a、b、c参数，其中一个被试者的辨识结果如图6所示. 其中，No.为时间窗口序号. a参数随时间波动反映历史RMS对当前RMS的影响程度与趋势. 在初始阶段，系统会经历较大程度的动态调整，表现为a参数的剧烈变化，反映系统初始化过程. 随着时间推移，a参数趋于平稳，表明动态线性部分调整减少. b参数在系统初始时波动较小，表明模型非线性部分在该阶段相对平稳. 随着疲劳加剧，b参数变化变得显著，模型中的非线性效应增强. 这说明在疲劳过程中，MF和RMS的关系变得更加复杂. c参数在初期系统未稳定时表现出波动，而在后期基本保持恒定，其中微小波动主要由即时输入变化引起，符合Hammerstein模型无记忆特性^[23]. 相较于简单的时、频域特征分析，Hammerstein模型参数分别描述了系统的动态与非线性变化，能够更深入地呈现肌肉疲劳过程中的复杂现象.

利用K-means++聚类，SSE变化如图7所示. 可知，最佳聚类数为4，且平均轮廓系数为0.81，表明聚类效果理想. 进而对动态参数矩阵${\boldsymbol{\varPsi }}$进行聚类分析，得到预备疲劳、轻度疲劳、中度疲劳和重度疲劳4类，如图8所示. 其中，L为肌肉疲劳等级.

因每位被试者在完成过头作业实验任务时的Borg CR10评分均未达到9，将10名被试者的Borg CR10评分（0~8）分为4类，共计84种组合，并将每种组合均与被试者的疲劳分类进行皮尔逊相关性分析. 结果显示，6名被试者的最高皮尔逊相关系数均出现在第65组，4类疲劳等级的划分方式为{0~2,3,4,5~8}，且p<0.01. 其他4名被试者，在第65组的皮尔逊相关系数大于0.86，整体平均皮尔逊相关系数为0.904（见表5），表明{0~2,3,4,5~8}的疲劳等级划分具有合理性.

进一步对比{0~2,3,4,5~8}的疲劳等级划分和聚类结果的关系，如图9所示. 结果显示，聚类结果在肌肉疲劳阶段的划分上比主观评分更敏感，能更早检测到疲劳进展. 主观评分可能因感知滞后或判断偏差而延迟，并出现波动，如被试者9，这反映了不同被试者个体对疲劳感知的差异性.

计算其中一名被试者的MF、RMS、a、b、c参数及对应聚类的疲劳等级分段，如图10所示. 参数a、b、c在不同疲劳阶段的波动如图11所示. 其中①表示预备疲劳，②表示轻度疲劳，③表示中度疲劳，④表示重度疲劳.

对于参数a，在分段①中，当肌肉刚开始收缩时，RMS较低，当前神经肌肉招募主要是慢肌纤维. a参数不确定性较高（均值：−0.097，标准差：0.438），系统初始动态响应不稳定. 在分段②中，随着肌肉工作时间延长，快肌纤维开始被动招募以补充所需肌力，a的均值从负转为正，且波动减小（标准差：0.148），方差和极差分别为0.022和0.740，说明输出对历史RMS的依赖性有所增强，动态反馈处于调整状态. 在分段③中，慢肌纤维在长时间工作后疲劳累积，贡献逐渐减少，而快肌纤维招募显著增加，以维持肌肉工作状态. 由于快肌纤维能量消耗较快且抗疲劳能力较低，导致RMS在此阶段出现波动，参数a趋于稳定（均值：0.120），但动态影响有所减弱，表现为标准差（0.055）和方差（0.003）的显著下降. 在分段④，肌肉进入重度疲劳后，慢肌纤维的工作能力明显下降，仅能维持少量基础活动，快肌纤维成为主导. 快肌纤维由于快速抽搐特性，能够维持较高的力量输出，但过多的能量消耗导致疲劳进一步加重. 此时，a的均值进一步升高（0.151），波动幅度显著降低（标准差：0.027），反映系统动态调整能力进一步减弱.

对于参数b，在分段①中MF较高，慢肌纤维为主导，快肌纤维参与较少，b数值整体稳定（均值：0.121，标准差：0.478，方差：1.710，极差：0.252），表明此阶段MF与RMS间的关系较为线性. 在分段②中，随着疲劳累积，MF下降，b的波动性增大（标准差：52.996，方差：240.540，极差：3107.640）. 此时慢肌纤维的招募增加，快肌纤维开始被动招募，非线性增强，RMS与MF的关系表现为较大的波动性. 在分段③中，疲劳进一步加重，快肌纤维成为主导，MF下降明显，非线性特性进一步放大（均值：−123.209，标准差：790.320，方差：248.950，极差：7494.690）. 在分段④中，肌肉逐渐达到重度疲劳，b的波动性相比分段③有所减小（标准差：49.312，方差：167.483，极差：2744.960）. 此时慢肌纤维逐渐减少，快肌纤维继续维持肌肉活动.

对于参数c，在分段①的初始阶段波动较大（均值：−0.276，标准差：11.488，方差：43.010，极差：209.650），说明此时系统尚未稳定. 随着时间推移，模型的动态行为趋于平稳，因而c趋于稳定，如分段②~④的均值接近于0，标准差显著减小，方差和极差也随之降低，该稳定性反映参数c对即时输入的响应特性，与无记忆输入特征对应.

为了进一步探究上肢肌肉疲劳趋势和模型稳定性，利用式(7)的G(z)将参数对应极点绘制在极平面上，再由不同疲劳阶段划分极点. 如图12所示为所有被试者的极平面图. 在整个模型参数跟踪过程中，极点均位于单位圆内，说明当前输入输出下的Hammerstein模型具有良好稳定性. 且不同疲劳阶段的极点都随疲劳等级增加，由单位圆边缘向圆心收敛，说明随着过头作业上肢肌肉疲劳加重，肌肉动态响应能力下降，振荡逐渐减弱，系统响应速度变慢.

过头作业上肢肌肉疲劳状态的4个阶段反映了上肢肌肉疲劳的全局变化趋势，为了进一步确定具体疲劳时间，根据式(16)得到的矩阵${\boldsymbol{\varPsi }}$计算各疲劳阶段的具体EMGFT，结果如表6所示.

如表6所示，每位被试者在完成相同作业任务时的具体操作存在差异，使得疲劳阶段与完成作业时间不同. 利用双因素方差检验被试者的主体间效应，以评估个体差异对EMGFT的影响，结果如表7所示. 其中，SS为平方和，MS为均方. 对于因素“被试者”，显著性水平p>0.05，表明不同被试者之间的疲劳阈值差异不具显著性. 被试者在各疲劳阶段表现趋于一致，且对单个被试者进行相关性分析具有可靠性.

通过Hammerstein模型与sEMG特征，提出新的过头作业肌肉疲劳评估方法，将sEMG的MF、RMS分别作为Hammerstein模型输入与输出，辨识模型参数，解释过头作业上肢肌肉疲劳过程并揭示其动态变化规律. 同时，将动态参数与EMGFT算法融合，识别由聚类获得的不同疲劳分段的具体疲劳阈值. 本研究不仅能够有效分离肌肉疲劳的各个阶段，且与Borg CR10主观分类具有较高一致性，有助于辅助过头作业中的上肢肌肉疲劳监测、预测与评估，为避免WMSDs、过头作业优化及助力装置设计提供理论支持. 然而，当前研究基于健康被试者在受控环境下开展，尚未充分验证所提方法在真实工业场景中长期动态作业下的鲁棒性与适用性. 下一步将拓展至复杂工况下的多任务作业人群，结合运动学数据与长期跟踪分析，进一步提升模型在实际应用中的泛化能力与实用性.