Dynamic assessment for low back-support exoskeletons during manual handling tasks

1 Introduction

Back-support exoskeletons help farmers, nurses, and industrial workers reduce lumbar burden and improve working efficiency (Kobayashi et al., 2009; Hasegawa and Muramatsu, 2013; Upasani et al., 2019). As one of the specific human-robot collaboration solutions for manual handling tasks, back-support exoskeleton should satisfy the safety requirement for both robot and humans (De Looze et al., 2016; Ajodani et al., 2018). Thus, developing safe and effective exoskeletons will enable a broad range of applications that could benefit users.

ISO 13482 states that the purpose of exoskeletons is to reduce physical workload (ISO, 2014). The safe limit for human lumbar workload averages 3.4 kN (Water et al., 1993). However, accurately evaluating exoskeleton safety is difficult owing to the complexity of human–robot interactions and the unpredictability in user movements. Moreover, exoskeletons can constrain human movement, leading to discomfort (Baltrusch et al., 2018). Performance variations at actual rehabilitation, industrial, and agriculture work sites highlight the need for a standardized and dynamic assessment method for exoskeletons (De looze et al., 2016; Omoniyi et al., 2020; Zheng et al., 2022).

Exoskeleton assessment requires obtaining lumbar load and human movement data. Human movement can be measured by optical or inertial measurement unit (IMU) motion capture systems. Several methods have been used to measure exoskeletons’ lumbar load, including directly recording the assistive force by inserting additional load cells into the exoskeleton (Abdoli-Eramaki et al., 2007; Abdoli-Eramaki et al., 2008); assuming a relationship between the magnitude of electromyographic signals and assistive forces, and then estimating the assistive forces by recording trunk muscles’ activities (Marras et al., 2000; Lamers et al., 2018; Weston et al., 2018; Koopman et al., 2019b); and establishing an exoskeleton by testing its characteristic performance beforehand (Koopman et al., 2019a). Nabeshima et al. (2018) developed a non-human testing framework to obtain lumbar torque. These methods, which estimate assistive forces, can be used with an inverse dynamic human model to calculate the actual lumbar load when using exoskeletons.

Further, several statistical methods have been employed to examine lumbar load and other key variables. These methods can be flexibly employed to scrutinize data across various temporal frames—individually, collectively, or even utilizing time-weighted averages—based on the specific demands of their research objectives. For instance, the effect of using exoskeleton on a biomechanical variable can be easily examined using the t-test (Abdoli-Eramaki et al., 2007; Whitfield et al., 2014; Lamer et al., 2018). Analysis of variance (ANOVA) would be suitable for evaluating the mixed effect between the testing conditions and exoskeleton modes on users (Marras et al., 2000; Abdoli-Eramaki et al., 2008; Ulrey and Fathallah, 2013; Weston et al., 2018; Koopman et al., 2019a; Koopman et al., 2019b; Poliero et al., 2020). Principle component analysis was used to select the important features and identify the differences of using exoskeletons (Sadler et al., 2011).

Except peak burden, when it is effective to use exoskeletons, and when motion is restricted are also concerned and can contribute to the risk of lumbar injury (Upasani et al., 2019; Omoniyi et al., 2020). Thus it is necessary to consider the effect of the exoskeletons not only on the timing when peak lumbar burden occurs but on the lumbar burden across the whole task. In addition it is also found the phase shift.

To the author’s knowledge, the effect of exoskeleton at the time other than the peak value in manual handling tasks has not been well-investigated. Therefore, we employed a time-series analysis method Functional data analysis (FDA), specifically functional ANOVA (FANOVA), which is designed to handle functional data such as biomechanical data, accounting for their continuous nature and temporal dependencies (Ramsey and Silverman, 2005). The FANOVA could be separated into a few steps. First, time-series biomechanical data such as the lumbar load and flexion angle were collected (Section 2.3 and Section 2.4). Second, functions were used to present the waveforms of the biomechanical variables, with B-spline the most commonly used (Ramsey and Siverman, 2005) (Section 2.5.1). Third, the obtained functions were aligned at identical timing points to obtain a representative comparison using a data registration (or data alignment) technique (Godwin et al., 2010) (Section 2.5.1). Finally, the FANOVA model was used to calculate the F-ratio between using and not using the exoskeleton (Section 2.5.2). In agriculture, FANOVA has been used to accurately estimate continuous growth trends (Xu et al., 2018) and demonstrate significant differences in various biomechanical contexts, such as lip kinematics and fatigue-induced kinematics changes (Ramsay et al., 1996; Godwin et al., 2010).

Compared to the traditional t-test, ANOVA, and PCA methods, FDA is better at dealing with the time-series dataset. The traditional methods usually identify vital features related to the performance of motions and injury mechanisms from the waveforms of the biomechanical data by referring to some individual points and reducing the dimensionality of the waveforms (Moudy et al., 2018). These methods are limited in that the important features are identified before the data have been analyzed. In contrast, FDA can be applied to multidimensional signals and eliminates the need for the prior identification of the relevant features (Donoghue et al., 2008). In addition, the traditional methods have difficulty finding the differences between groups, and individual differences may produce conditions that will cause timing/phase variability in the waveform (Moudy et al., 2018). For example, in different trials, the subject may not be able to reach the maximal flexions at the same time, while the FDA can help us to minimize the time difference (data alignment) among the maximal flexions, and to maintain the shape and amplitude of each curve. Using this technique, FDA can highlight the features of waveforms to reduce the timing/phase variability so that we can analyze the effect of exoskeletons on the magnitude at all timings, and it is not necessary to identify the peak value or mean value (Moudy et al., 2018). Thus, using FDA to assess the performance of an exoskeleton makes it possible to determine when the exoskeleton can significantly reduce the lumbar load or restrict human motion during a task.

Industrial exoskeleton usage requires safety considerations in high-risk scenarios, such as dynamic symmetrical and asymmetrical lifting tasks (De Looze et al., 2016; Huysamen et al., 2018). FANVOA can be a suitable method for revealing exoskeletons’ effect on key variables in each phase of these tasks, thus promoting their standardization.

Although important discrete peak and mean values of the lumbar load and kinematics factors have been studied in exoskeleton evaluations, no method is available to evaluate the effectiveness of various postures during lifting-flexion movement of the exoskeleton. Thus, this study was novel because it not only evaluated the performance of the exoskeleton when variables reached their peak values but also at all other times. This made it possible to identify the variability when using the exoskeleton throughout the entire lifting-flexion motion and assess its applicability to the entire manual handling movement.

In this study, we mainly focus on developing a FANOVA-based method to evaluate the effects of an exoskeleton on the biomechanical variables at every time point during the manual handling tasks. This analysis allows us to examine when exoskeletons can significantly reduce lumbar load and restrict motions during the task. In addition, we also proposed an exoskeleton–human model to estimate the dynamic lumbar load with the exoskeleton’s assistance. It was hypothesized that the FANOVA would demonstrate the exoskeleton, and significantly affect not only the peak biomechanical variables but also that at the other period.

2 Materials and methods

In this study, we used the FANOVA to investigate the effect of the exoskeleton. In this section, 2.1 and 2.2 show the tasks (Figure 1), and instrumentation of this experiment. In the assessment of the exoskeleton, the details are developed as two steps: estimation of biomechanical variables (Section 2.3 and Section 2.4) and FANOVA assessment (Section 2.5 and Section 2.6). As shown in Figure 2, in the first step, using a 3D human model and exoskeleton model to estimate the biomechanical variables, which will be analyzed in the second step using FANOVA after data smoothing and registration, and finally, the effect of the exoskeleton on each variable can be assessed using F-ratio of time points.

FIGURE 1. Manual handling task of a 10 kg box with an exoskeleton under (A) symmetrical and (B) asymmetrical conditions.

FIGURE 2. The procedure for assessing the effectiveness of the exoskeleton using FANOVA; the movement data are imported into an exoskeleton-human model to estimate the biomechanical variables, which are defined in Section 2.3; then smoothing and registration ensures the data can be presented as functions and aligned with each other; finally, conducting FANOVA and estimate F-ratio to assess the exoskeleton’s effect at the time points.

2.1 Participants and tasks

Ten male participants (height: 1.72 ± 0.08 m, body mass: 68.1 ± 8.8 kg, age: 30.9 ± 7.7 years) were recruited and all of them consented to join this experiment. This study was approved by the institutional review board of the Institute of Agricultural Machinery, National Agriculture and Food Research Organization (approval no. Kakushin-ken_Rinri_R03-02). Participants aged 20–40 years were selected because manual handling tasks pose a high risk of low back pain, and young individuals can tolerate relatively high lumbar loads (Kudo et al., 2019).

The manual handling tasks comprised three phases: free-flexion, box-lifting, and box-laying. Figures 1A, B show the beginning of symmetrical and asymmetrical lifting, respectively. The participants performed 2 × 2 non-repetitive tasks (asymmetrical/symmetrical condition with/without the exoskeleton) five times each, the time interval between two trials was around 30 s to reduce the effect of muscle fatigue. The participants’ feet were always pointing forward. Prior to the experiment, participants were instructed to perform manual handling tasks at their preferred speed to test their strength.

2.2 Instrumentation

A motion capture system (Xsens MVN Analyze, Xsens, Inc., Enschede, Netherlands) was used to reconstruct motion using IMUs. As shown in Figure 1, 15 IMUs were attached to the participants’ heads, shoulders, L5/S1, upper arms, forearms, thighs, shanks, and feet. Two identical force-measuring devices recorded external forces on the hands, each consisting of two three-axis force sensors (USL08-H6, Tec Gihan Co., Ltd., Kyoto, Japan). The box size was 57 cm × 28 cm × 10 cm, with a total mass of 10 kg, which is the limit mass for a normal adult in a one-time lift (ISO 11228-1, 2021). Data recorded at 60 Hz were filtered using a low-pass filter with a 4 Hz cut-off frequency.

2.3 Exoskeleton’s effect on biomechanical variables

Lumbar load reduction is the biomechanical criterion for relieving the lumbar burden (Waters et al., 1993). The kinematics variables such as trunk angle, trunk angular velocity, trunk angular acceleration, and the horizontal displacement between the wrist and the lumbar represent the effect of the exoskeleton on the motion restrictions (Potvin, 1997; Marras et al., 2000).

The five representative biomechanical variables are presented as lumbar load ($F_c$), trunk angle ($\theta$), trunk angular velocity ($\omega$), trunk angular acceleration ($\alpha$), and horizontal displacement between the wrist and lumbar ($H$). The trunk angle $\theta$ = ∠C7·L5·C is shown in Figure 1A, and the trunk angular velocity and acceleration are the derivatives of the trunk angle and angular velocity, respectively. In Figure 1A, the C7 and L5 positions are determined by virtual markers generated by the motion reconstructed using Xsens MVN Analyze. The wrist (D) and knee (C) positions were determined on the basis of the center of the virtual markers of the left and right wrist (D1, D2), as well as the left and right knee (C1, C2), respectively. During the asymmetrical task, the angle (γ) between the box and the table (65 cm height) in the horizontal plane was 90°, as shown in Figure 1B.

2.4 Lumbar load estimation

2.4.1 Exoskeleton model

The “Muscle Suit Every” (Innophys, Inc., Tokyo, Japan) exoskeleton was used in this study. This exoskeleton uses artificial muscle to provide assistive torque at different trunk angles. In addition, angular velocity influences assistive torque by altering artificial muscle extension speed, which affects the friction force of the artificial muscle and assistive torque generation (Sugimoto et al., 2011; Tondu, 2012). Thus, we used the thin-plate spline (TPS) model for estimating assistive torque based on trunk angle and angular velocity. TPS can provide a robust estimation for spatial data interpolation and surface fitting (Bookstein, 1989; Donato Belongie, 2002).

To establish the TPS model, we used a testing machine to record the assistive torque of the exoskeleton when the trunk angle rotated from 0° to 90° at five speeds (10, 30, 45, 60, and 90°/s) at a sampling rate of 1 kHz (Tanaka et al., 2020). In every test, air pressure of the exoskeleton was set to 0.1 MPa. The tests were repeated at each speed 10 times. Subsequently, the recorded data were processed with a low-pass filter of 10 Hz. The results, shown in Figure 3A, are the relationships between the extension angle and average assistive torque in the 10 trials under each speed condition. Finally, utilizing the TPS method, which was established using MATLAB (version 2022a), the estimated assistive torque could be presented in terms of trunk angle and angular velocity, as shown in Figure 3B. The TPS model’s R2 was 0.955 with a root mean square error of 3.70 for all testing data.

FIGURE 3. Procedure to estimate exoskeleton’s assistive torque based on the TPS method. (A) The measured assistive torque and trunk angle at different speeds. (B) The TPS model estimation of assistive torque using trunk angle and angular velocity. Abbreviations: TPS, thin-plate spline.

2.4.2 Actual lumbar load estimation

As shown in Figure 4, A link-chain human three-dimensional (3D) model was established. The coordinate, movement and center of mass for each segment, and motion reconstruction were recorded using the IMU motion capture system attached to the subject body (data collection). The dynamic link-chain model had nine segments, involving the forearms, upper arms, head, shoulder, thoracic spine, lumbar spine, and pelvis (the lower limbs are not included in the dynamic calculation). The mass of each segment was a proportion of the total body mass, as estimated by a previous study (Winter et al., 2009; Hof, 1992). In addition, the mass of the exoskeleton was added to the center mass at the lumbar segment.

FIGURE 4. Schematic of estimating biomechanical variables using the human 3D model. After obtaining the exoskeleton model, the motion and the external force data are imported into the human 3D model, and then the biomechanical variables are estimated.

To calculate the lumbar torque, we initially estimated the non-assisted lumbar torque using the top–down inverse dynamic method based on the link-chain human model. Given that the assistive torque only acts in the flexion plane (plane C7L5C in Figure 1A), we subtracted the assistive torque in this plane from the non-assisted torque to obtain the actual torque at the lumbar joint. The geometric model of the trunk muscles was determined using previously reported data (Schultz et al., 1982; Granata and Marras, 1993; Gagnon et al., 2001), and the muscular forces were estimated by minimizing the sum of the trunk muscles’ stress square (Anderson and Pandy, 2001). Finally, the lumbar load was obtained by the force resulting from the muscular forces and the upper body load in the direction perpendicular to the lumbar vertebra. The inverse dynamic computation was completed in MATLAB (version 2022a), and the optimization procedure was completed by the quadratic programming algorithm (Stellato et al., 2020).

In order to investigate the accuracy of the 3D human model, we compared the estimated lumbar load from this 3D human model with the in vivo data reported by Wilke (Wilke et al., 1999; Wilke et al., 2001) under several body conditions. The estimated lumbar load for each posture was obtained from the average estimation among eight volunteers with a mean body height of 1.72 m and mean body mass of 69 kg, similar to the participant (1.68 m, 70 kg) in Wilke’s study.

The accuracy estimation procedure was presented in Figure 5, where the lumbar loads estimated by the human model at 12 postures including standing, 36° flexion, maximal flexion, extension, lateral flexion, twisting, stoop lifting, squat lifting, one-hand carrying, and close-to-chest handling were compared with these obtained from the reported in vivo experiments. The in vivo lumbar intradiscal pressure (MPa) estimated by Wilke was converted to lumbar load (N) using the correction factor proposed by Dreischarf et al. (2013). Finally, calculating the R of Intraclass Correlation Coefficient (ICC) between the estimated lumbar loads of the human model and the invasive data for the postures [type of ICC (1,1), and on the 95% confidence interval]. The estimated R is 0.93, which can be explained as excellent reliability (>0.9) (Koo and Li, 2016). The result indicates that the lumbar load estimated by the human model has a high consistency with the in vivo data.

FIGURE 5. Comparison between reported the in vivo data and the estimated lumbar load using the 3D human model. The mass of the box is simulated as 20 kg.

2.5 Functional data analysis

2.5.1 Smoothing and processing

This study’s biomechanical time-series discrete data, including tracks, angles, and lumbar moment, were converted to functional data in the FDA. In this step, we use a series of the basis functions to fit the recorded biomechanical data and one advantage of this step is to reduce the influence of noise. Since non-repetitive tasks were conducted in this experiment, B-spline is used as the basis function system. Each obtained function can be presented as follows (Ramsay and Silverman, 2005):

where $x_i\left(t\right)$ represents the function converted from the i-th observed data series; t represents the number of time points; $c_{ik}$ represents the coefficients; and ${\phi }_k\left(t\right)$ are B-spline basis functions with the number, K.

The B-spline fitting functions for the recorded data can be determined using the least square method. The residual sum of squares and a penalty term based on the second derivative of the fitted curve were minimized. The following equation expresses the minimized penalized least squares problem (Ramsay and Siverman, 2005):

where $X\left(t_i\right)$ represents the observed data points at the i-th time point (i = 1, 2,…, N); $x\left(t_i\right)$ represents the estimated data at the i-th time point using the function obtained with Eq. 1; $\lambda$ is the smoothing variable, which is a non-negative constant that controls the smoothness of the function, with a larger value leading to a smoother function; and $D^{2}x\left(t\right)$ represents the second derivative of the fitted function in time-series. The first part measures goodness of fit between the data points and the fitted function, and its objective is to minimize the residuals. The second part is the penalty term that encourages smoothness in the fitted function by penalizing abrupt changes in curvature; the penalty term is proportional to the integral of the squared second derivative of the function $x\left(t\right)$ over the domain (a, b). Thus, this objective function allows us to control the trade-off between the goodness of fit and the smoothness of the function.

Data registration is a technique that aligns generated functional data that might be misaligned. This can improve the ANOVA before estimating the main effects and interactions. We applied data registration to all observed data for each condition and variable using continuous registration (Ramsay and Silverman, 2005).

2.5.2 FANOVA

After obtaining the functional data, we performed FANOVA to investigate the effect of exoskeletons on each variable under symmetrical and asymmetrical conditions. The model for each variable, $f\left(t\right)$, with time history can be presented as follows:

where $\mu \left(t\right)$ is the mean function indicating the average value of all trials under the symmetrical or asymmetrical condition; ${\alpha }_i\left(t\right)$ represents the effect of using (i = 1) and not (i = 2) exoskeletons on the variable; and $\epsilon \left(t\right)$ is the unexplained variation. We identified the specific effects of using exoskeletons; the constraint added for all t, as ${\alpha }_1\left(t\right)+{\alpha }_2\left(t\right)=0$.

Then, the model for each variable for each t can be rewritten as a matrix form $\mathit{f}\left(t\right)$ as follows:

where $\mathit{f}\left(t\right)$ is the 20 × 1 function vector and $\mathit{Z}$ is the 20 × 3 design matrix, with the 20 rows corresponding to the 10 participants, each contributing with two curves: one when using exoskeletons, and the other when not using them. The first column has ones; the second column has zeros in the first 10 rows, followed by ones; and the final column has ones in the first 10 rows, followed by zeros. $\mathit{\beta }\left(t\right)$ is the 3 × 1 vector of parameter functions, with $\mathit{\beta }\left(t\right)={\left[\mu \left(t\right),{\alpha }_1\left(t\right),{\alpha }_2\left(t\right)\right]}^{\prime }$; $\mathit{\epsilon }\left(t\right)$ is the 20 × 1 vector of residual functions.

The vector $\mathit{\beta }\left(t\right)$ can be estimated by minimizing the linear minimum mean square error (LMSSE):

$\text{LMSSE}\left(\mathit{\beta }\right)$ should be minimized under the condition ${\alpha }_1\left(t\right)+{\alpha }_2\left(t\right)=0$.

As with traditional ANOVA, the error sum of squares for the residual function and the mean curve were evaluated as a function for each time point in (t) with the following. SSE (sum of squared errors) and SSY (sum of squared in a total) of the model at time t can be calculated as:

where $\widehat{\mathit{\beta }}\left(t\right)$ is the 20 × 1 vector of the predicted parameters function; $\widehat{\mathit{\mu }}\left(t\right)$ is the 20 × 1 vector of the predicted mean function of all trials.

The F-ratio determines whether the variance between two data sets is equal, and the FANOVA can evaluate the F-ratio across time t. The F-ratio can be presented as follows:

where df(error) is the degree of freedom for error [df(error) = 1], and the df(regression) is the difference in degrees of freedom [df(regression) = 18]. Based on the F-criterion is 4.41. Any time the functional curve F-ratio > F-criterion, the effect of time had reached significance at the chosen α level of 0.05.

2.6 Statistical analysis

For all tasks, the paired t-test (at an α level of 0.05) was used to compare the peak values of the lumbar load, trunk angle, angular velocity, angular acceleration, and horizontal distance between the wrist and lumbar spine under symmetrical and asymmetrical conditions. Subsequently, the result of t-test was compared to that of FANOVA, which was utilized to examine the continuous effect of the exoskeletons on these variables throughout the normalized manual handling tasks. The F-ratio, obtained as time history, was compared to the F-criterion at an α level of 0.05.

The functional analysis for time-series data was performed using the package developed by Ramsey and Silverman, available at https://www.psych.mcgill.ca/misc/fda/downloads/FDAfuns/. All the statistical analyses were conducted using MATLAB (version 2022a).

3 Results

3.1 T-test analysis: effect on peak values of each variable

The estimated peak variables of lumbar load, trunk angle, angular velocity, angular acceleration, and horizontal displacement between the wrist and lumbar are shown in Figure 6, which shows the average peak variables and standard deviations for all participants. Peak lumbar load and angular velocity were significantly lower (p < 0.05) when using the exoskeleton during both tasks; however, peak horizontal displacement was significantly lower only during asymmetrical tasks (p < 0.01). Compared with not using exoskeletons, using them reduced peak lumbar load by 388 N (14%) and 427 N (17%) during both tasks, respectively. Similarly, the peak trunk angular velocity reduction was 23°/s (24%) and 25°/s (24%) under both conditions. During asymmetrical tasks, the peak horizontal displacement was reduced by 0.02 m (5%). No significant difference was observed in peak trunk angles and angular acceleration variables for both tasks.

FIGURE 6. The t-test result: exoskeletons’ influence on the peak value of representative variables during manual handling tasks. Abbreviations: acc., acceleration; disp., displacement; vel., velocity; w, with; w/o, without (*p < 0.05; **p < 0.01; ***p < 0.005).

3.2 FANOVA: effect on time-series values of each variable

The results of the functional analysis, as presented in Figures 7–9, demonstrate the effect of using exoskeletons on different variables during asymmetrical and symmetrical tasks.

FIGURE 7. Effect of using exoskeletons on the variables during asymmetrical tasks. (A) Absolute values of the variables when using (solid line) and not using (dash line) exoskeletons, where positive (+) and negative (−) for angular velocity and acceleration represent flexion and lifting, respectively. (B) The corresponding F-ratio generated by FANOVA (solid line) and Fcri = 4.41 with statistical significance at an α level of 0.05 (red dash line). The whole normalized time is separated into flexion (0%–33%), lifting (33%–67%), and laying (67%–100%) phases. Abbreviations: FANOVA, functional analysis of variance; Fcri, F-criterion; acc., acceleration; disp., displacement; vel., velocity; w, with; w/o, without.

During asymmetrical tasks (Figure 7), the exoskeleton reduced lumbar load by 412, 393, and 383 N with the greatest significant differences during the flexion, lifting, and laying phases, respectively. In addition, the most significant lumbar load reduction (412 N, F-ratio = 33.8) for all phases occurred during the flexion phase (31.7%), as shown in Figure 9. Reductions in the kinematic variables of trunk angle, angular velocity, and angular acceleration were observed in all phases (Figure 7), indicating the exoskeletons’ impact on posture and movement restrictions during asymmetrical tasks. As shown in Figure 9, the most significant reduction (5.5°, F-ratio = 18.6) in trunk angle occurred during the laying phase (86.1%); for trunk angular velocity, the reduction (20.0°/s, F-ratio = 18.1) occurred during the flexion phase (14.1%); and for trunk angular acceleration, the reduction (61.5°/s², F-ratio = 18.6) occurred during the flexion phase (20.0%). Horizontal displacement was only significantly reduced during the lifting phase (61.9%), by 0.02 m (F-ratio = 11.1).

During symmetrical tasks (Figure 8), the exoskeleton reduced lumbar load by 300, 672, and 280 N with the greatest significant differences during the flexion, lifting, and laying phases, respectively. The most significant lumbar load reduction (300 N, F-ratio = 35.9) for all phases occurred during the flexion phase (26.7%) (Figure 9). Similarly, in asymmetrical tasks, reductions in the kinematic variables such as trunk angle, angular velocity, and angular acceleration were observed in all phases (Figure 8). As shown in Figure 9, the most significant reduction in trunk angle (5.9°, F-ratio = 20.6) occurred during the laying phase (83.2%) and that in trunk angular velocity (22.5°/s, F-ratio = 24.1) during the flexion phase (13.4%). For trunk angular acceleration, the reduction (79.4°/s², F-ratio = 16.8) occurred during the flexion phase (10.6%). The patterns of the above kinematic variables were also similar to those in asymmetrical tasks. Horizontal displacement was significantly reduced during the lifting (15.8%) and flexion (42.6%) phases by 0.02 m with similar F-ratios (5.3 vs. 4.8).

FIGURE 8. Effect of using exoskeletons on the variables during symmetrical tasks. (A) Absolute values of the variables when using (solid line) and not using (dash line) exoskeletons, where positive (+) and negative (−) for angular velocity and acceleration represent flexion and lifting, respectively. (B) The corresponding F-ratio obtained by FANOVA (solid line) and Fcri = 4.41 (red dash line). The whole normalized time is separated into flexion (0%–30%), lifting (30%–60%), and laying (60%–100%) phases. Abbreviations: FANOVA, functional analysis of variance; Fcri, F-criterion; acc., acceleration; disp., displacement; vel., velocity; w, with; w/o, without.

FIGURE 9. FANOVA result: The reduced variables at the peak F-ratio timing for each phase during the manual handling task between w, and w/o assistance; The bar values, representing the reduction of each variable (blue: asymmetrical task; red: symmetrical task), and are obtained at the peak F-ratio of flexion, lifting, and laying phase, respectively. As F-ratio is smaller than 4.41, it is considered no significant difference is found between w, and w/o assistance. Abbreviations: acc., acceleration; disp., displacement; vel., velocity.

Except for providing assistive torque, exoskeletons can reduce the lumbar load by restricting movements. As shown in Figures 7, 8, some significant overlaps were observed between the F-ratio of the lumbar load and the kinematic variables in each phase. During the flexion phase, these overlaps of the significant F-ratios were as follows: lumbar load: 13%–31%; trunk angle: 13%–24%; angular velocity: 13%–18%, 25%–31% (in the flexion phase, both periods were where F-ratios of angular velocity were significant); angular acceleration: 15%–30%; and horizontal displacement: 14%–17%. During the lifting phase, lumbar load: 35%–45%; trunk angle: 39%–49%; angular velocity: 41%–45%; and horizontal displacement: 41%–43%. During the laying phase, lumbar load: 56%–100%; trunk angle: 73%–100%; angular velocity: 56%–58%; and angular acceleration: 79%–82%. These overlaps show the consistency between the reduction of lumbar load and that of kinematic variables, which indicate that the lumbar load reduction is not only affected by lumbar moment reduction but also by restricting the kinematic variables when using the exoskeletons.

4 Discussion

4.1 T-test analysis: effect on peak values of each variable

The findings presented in Figure 6 show that exoskeletons significantly affect the representative biomechanical variables during symmetrical and asymmetrical tasks. Exoskeleton use substantially reduced peak lumbar load and motion speed, consistent with previous findings (Koopman et al., 2019b). Increasing trunk angular velocity requires greater trunk muscle activation (Dolan and Adams, 1993). This result aligns with those of previous studies, suggesting that exoskeletons help reduce lumbar load and muscle activation (Huysamen et al., 2018; Lamer et al., 2018), thus, lowering the risk of musculoskeletal disorders and enhancing worker comfort.

No significant differences were detected in peak trunk angles and angular accelerations, in contrast to angular velocity, for both tasks. Thus, it results difficult to infer from peak values whether exoskeletons limit the range of motion or hinder human movement while providing support and reducing lumbar load. Previous studies have yielded mixed results on the exoskeletons’ influence on trunk angle: reduction was observed on the nylon elastic support, while no reduction was shown for another passive exoskeleton (Laevo V2.4 Delft, Netherlands) (Marras et al., 2000; Koopman et al., 2019a). The differing results from our experiment and previous studies regarding the influence of exoskeleton on trunk angle can be attributed to variations in design.

The lumbar reduction of the current exoskeleton can meet the 3.4 kN average lumbar load criterion recommended by the National Institute of Occupational Safety and Health (Waters et al., 1993), as shown in Figure 6. However, the lumbar load limit will vary across age groups and sexes (Genaidy et al., 1993). Aging societies have a growing population of older workers, whose lumbar load limit is 1.69 kN lower than that of younger workers (Kudo et al., 2019). Women have lower lumbar load limits than men (Genaidy et al., 1993; Kudo et al., 2019). Thus, deterministic assistive forces should consider lumbar load limits for different age groups and sexes.

4.2 FANOVA: effect on time-series values of each variable

FANOVA enables a more comprehensive examination of the effects of exoskeletons on variables throughout an entire task rather than focusing solely on peak values. The FANOVA results suggested that exoskeletons can alleviate lumbar burden at peak load timings and throughout all task phases (Figure 9). The lumbar load reduction in Figures 7, 8 implies that exoskeletons are probably effective in reducing the lumbar load in all phases during asymmetrical tasks while being much more effective in the lifting phase during symmetrical tasks. Asymmetrical lifting reportedly results in a higher lumbar load than symmetrical lifting, implying that workers are more easily prone to getting lumbar injuries during asymmetrical lifting (Kim and Zhang, 2017). Therefore, improving exoskeletons’ performance during asymmetrical lifting can reduce users’ lumbar injury risks.

The reduced trunk angle, angular velocity, and acceleration during both manual handling tasks suggest a similar tendency when exoskeletons assist participants. The results indicated that using the exoskeleton imposes significant restrictions on body movement at larger flexion postures (>25°), in the vicinity of peak lumbar load occurrence, and throughout most of the box-laying movements in both tasks. The most significant reduction in angular velocity and acceleration occurred in the middle of flexion, implying that the human body may experience greater kinematic restrictions from exoskeleton use when not under loading conditions. This could be a reason for exoskeletons not supporting the lumbar load in a small trunk angle, which will be discussed subsequently. However, the restrictions on horizontal displacement differ between symmetrical and asymmetrical conditions. Under asymmetrical conditions, a significant reduction occurred at the end of lifting. In contrast, in symmetrical tasks, restrictions occurred in the middle of flexion and the vicinity of peak lumbar load occurrence. Determining the relationship between these changes using other variables in this study was challenging. This discrepancy in horizontal displacement restrictions may originate from the interaction between the trunk and upper limbs.

Comparing lumbar load and trunk angle in Figures 7, 8, significant reductions when using the exoskeletons were difficult to observe when the trunk angle was less than 25°. The only exception was during the asymmetrical lifting phase, when a significant lumbar load reduction occurred even with a trunk angle smaller than 25°. This lumbar load reduction was likely the result of reduced horizontal displacement rather than trunk angle changes. These results suggest that exoskeletons at smaller trunk angles (<25°) neither significantly reduce lumbar load nor restrict trunk movement. Since the greatest lumbar load typically occurs at larger trunk angles, this design approach may prioritize allowing users a greater range of motion and, consequently, sacrifice effectiveness at smaller trunk angles. This trade-off between range of motion and lumbar load reduction at smaller trunk angles may be a deliberate design choice to improve user comfort and task efficiency while offering lumbar support when needed.

Although restricting body movement may lead to inconvenience, adopting a low-speed or small-inclination posture will lead to lower lumbar injury risk (Marras et al., 1993; Waters et al., 1993). In addition, the lumbar load can also be reduced by decreasing the kinematic variables. The overlap of F-ratios with significance (F-ratio > 4.41) suggests that the reduction in lumbar load is affected by the kinematic variables in all phases. Significant trunk angle and horizontal displacement reductions were observed near peak lumbar load during symmetrical tasks. During asymmetrical tasks, only trunk angle reduction was noted. Quantitatively evaluating the impact of each kinematic factor on lumbar load reduction using FANOVA in conjunction with established ergonomic equations is possible (Potvin, 1997; Merryweather et al., 2009; Arjmand et al., 2011; Arjmand et al., 2012). Although quantitative assessment methods for the inconvenience of exoskeletons are lacking, the discomfort in kinematic restrictions probably results from contact pressure and friction (Baltrusch et al., 2018; Huysamen et al., 2018). Further investigation into the complex interplay between the various biomechanical variables and the different conditions could improve the understanding of the effects of exoskeleton use on human kinematics.

4.3 Exoskeleton simulation method

Direct measurement methods, such as electromyography or extension force, may introduce errors due to individual differences (Marras et al., 2000; Abdoli-Eramaki et al., 2007; Lamers et al., 2018). To address this, we employed a 3D interpolation method, the TPS, to represent the relationship between extension angle, angular velocity, and assistive torque. Comparing to torque–angle relationship (Koopman et al., 2019a), our approach accounts for certain dynamic factors affecting assistive torque. The advantages of TPS include efficient overfitting reduction through regularization, reduced computational complexity once the analytical form of the TPS model can be obtained, and adaptability to multidimensional datasets (Bookstein, 1989; Donato and Belongie, 2002), making it suitable for modeling various exoskeletons whose assistive torques may be influenced by different factors. Consequently, TPS is expected to contribute to exoskeleton standardization. By incorporating TPS in exoskeleton testing and refining FANOVA assessment methods, researchers can develop more effective, comfortable, and efficient exoskeletons that cater to workers’ diverse needs in various industries.

4.4 Comparison between FANOVA and other methods

Compared to the traditional methods such as a t-test or ANOVA, FANOVA includes a smoothing and data registration procedure that reduces the noise and variability of the timing/phases. Therefore, the biomechanical variables can be analyzed at each time point.

In this study, we conducted a t-test and determined how the peak value was affected by the dynamic movement of an exoskeleton. However, it was difficult to investigate the significant differences in the variables at the other time points because of the variability of the timing/phases of subjects. To reduce this variability, previous studies usually had subjects maintain specific postures before applying a t-test or ANOVA (Lamers et al., 2018; Koopman et al., 2019a). For example, Lamers et al. (2018) found that using an exoskeleton could significantly reduce a user’s lumbar load under flexion angles of 30°, 60°, and 90°. However, the static assumption may cause the lumbar load to be underestimated compared with those of actual tasks, which are usually dynamic (van Dieën et al., 2010). Moreover, it is a significant task to evaluate the effectiveness of an exoskeleton at all flexion angles, users are normally concerned about how to best use it (Upasani et al., 2019). In contrast, FANOVA makes it possible to evaluate the reduction in the dynamic lumbar load resulting from the use of an exoskeleton, and we found that the lumbar load was reduced at the majority of the flexion angles, with the exception of small angles (<25°) or during a flexion-lifting shift. Thus, compared to the traditional t-test and ANOVA, FANOVA can find the most effective conditions for using an exoskeleton.

FANOVA can also find the difference in the flexion angle at times other than under peak or specific conditions. Using a t-test, no significant difference can be found in the peak flexion angle under either a dynamic movement (Figure 6) or in static postures (Lamers et al., 2018). However, FANOVA shows that a significant restriction of the trunk angle can be found in all phases when an exoskeleton is used (Figures 7, 8). This could be because the exoskeleton mainly affects the flexion angle not at the beginning or end of lifting or flexion, but during the middle of the task.

The results of this study are expected to contribute to safety standards for exoskeletons. The safety requirements set forth in the international standard for wearable robots, ISO 13482, are limited only to conceptual design guidelines. This study investigated a method that is expected to assist manufacturers in quantitatively evaluating their products throughout the entire movement process and guide users in the appropriate use of an exoskeleton in lifting-flexion tasks.

4.5 Limitations

This study had a few limitations. First, this study relied on a single type of exoskeleton, although design can significantly influence an exoskeleton’s effectiveness in reducing lumbar load and modifying kinematic variables (Baltrusch et al., 2018; Kozinc et al., 2020; Luger et al., 2021). Second, the complex and time-intensive computations involved in FDA, such as the data registration procedure, may present challenges as the curves involve a considerable amount of data from lengthy experiments.

Furthermore, multiple comparisons may pose limitations. As the number of exoskeleton modes increases and significant differences are assessed at each time point, methods such as the Bonferroni adjustment, which may reduce the statistical power of the analysis, can be used (Khalaf et al., 1999). An alternative approach is establishing a critical number of simultaneous F-ratios that must exceed the F-crit to be considered significant for that period. However, the optimal number of F-ratios for this method remains unestablished (Godwin et al., 2010).

Other than the criteria mentioned, cumulative load can also contribute to low back pain. However, the lack of a safety criterion for cumulative load makes evaluating the safety of exoskeletons’ assistance on this variable difficult, and the value of the safety criterion will affect the evaluation of the exoskeletons’ assistance.

This experiment took an interval of around 30 s between two tasks. However, considering individual differences, taking a maximal voluntary contraction (MVC) test after resting would be a better way to confirm whether the rest time was sufficient to reduce muscle fatigue, which was a limitation of our work.

Despite these limitations, the study demonstrates the potential of using FANOVA to assess the effectiveness of exoskeletons in various manual handling tasks. Future research could focus on refining measurement techniques, developing more accessible tools for FANOVA, and exploring alternative statistical approaches for handling multiple comparisons.

5 Conclusion

A dynamic assessment method based on FANOVA was used to investigate the effect of utilizing exoskeletons on five representative biomechanical variables. The result implied that exoskeletons could reduce the lumbar load during manual handling tasks, particularly under symmetrical lifting conditions. The significant reductions in lumbar load and kinematic variables indicate that exoskeletons are crucial in protecting users’ lumbar spine and reducing the risk of low back injury. Furthermore, the exoskeletons achieve this reduction by restricting movement, which helps to maintain proper posture during handling tasks. The results also showed how exoskeletons indirectly affect the lumbar load, influencing other kinematic variables in time history. These findings contribute to developing safer and more effective exoskeleton designs, ultimately enhancing the practical adoption of exoskeletons in various scenarios involving manual handling tasks such as in agriculture, industry, and physical rehabilitation.

Data availability statement

The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

Ethics statement

The studies involving humans were approved by the institutional review board of the Institute of Agricultural Machinery, National Agriculture and Food Research Organization. The studies were conducted in accordance with the local legislation and institutional requirements. The participants provided their written informed consent to participate in this study.

Author contributions

XX: Conceptualization, Data curation, Formal Analysis, Methodology, Validation, Writing–original draft. MT: Conceptualization, Data curation, Funding acquisition, Methodology, Writing–original draft. SU: Methodology, Resources, Writing–review and editing. YuK: Conceptualization, Funding acquisition, Writing–review and editing. YoK: Resources, Writing–original draft.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. This research was supported by the Japan Society for the Promotion of Science (JSPS) KAKENHI (Grant Number JP21K05865).

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References