A total of 25 experienced male badminton athletes (age: 24.3 ± 4.5 years; height: 175 ± 3.6 cm; weight: 71 ± 4.2 kg; years of badminton playing: 7.1 ± 3.2 years; all were right-hand dominant) participated in the lab test, and the synchronized collection of motion data and IMU data (fist session) was performed. This study was approved by the Ethical Committee of the Research Institute in Ningbo University (RAGH20190901). All athletes were informed of the requirements, objectives, and procedures of the lab and on-court tests, and written consent was obtained.

The first session of the lab test synchronized the 3D motion capture and wearable sensors. The test involved a 12-camera Vicon system and Vicon IMU wearable sensors (Vicon Metrics Ltd., Oxford, United Kingdom) and an AMTI 3D force plate (AMTI, Watertown, MA, United States) (Figure 1). The collection frequency of the Vicon camera system was set at 200 Hz, and the IMU and force plate were set at 1000 Hz.

In order to mimic the real movements, we employed an established full-body marker-set model during motion capture (Rajagopal et al., 2016). The IMU sensors for the ankle and knee joints are illustrated in Figure 1B. In particular, the knee IMU was placed 2 cm below the medial condyle of the proximal tibia, and the ankle IMU was placed 3 cm above the medial malleolus of the distal tibia (Sheerin et al., 2020). The vertical axis of the IMU sensor was axially aligned with the midpoint of the ankle and knee joints in the tibia, which was validated in a clinical protocol of our previous study (Yeung et al., 2022). To follow up on our previous studies (Mei et al., 2017; Yu et al., 2021a), the four directional sub-maximal and maximal lunges were performed with a synchronous collection of motion, ground reaction force, and IMU data for further processing and analysis. Particularly, the sub- and maximal-right-forward forehand lunges ( Sub-RF and Max-RF ), sub- and maximal-left-forward backhand lunges ( Sub-LF and Max-LF ), right-backward backhand sub-maximal and maximal lunges ( Sub-RB and Max-RB ), and left-backward backhand sub-maximal and maximal lunges ( Sub-LB and Max-RB ) were performed with 80% (for sub-maximal) and 100% (for maximal) efforts, following the previous protocol (Lam et al., 2017; Lam et al., 2018; Yu et al., 2021a).

The joint kinematics, kinetics, and contact forces were calculated, following the previously established protocols of musculoskeletal OpenSim modelling. Machine learning models were also developed and tested using motion capture data against the acceleration and angular velocity data from wearable sensors.

First, the static marker positions and body mass were used to “scale” the generic model to match subject-specific musculoskeletal models (Figure 2), as per the standardized workflow (Delp et al., 2007), which was validated in our recent studies (Mei et al., 2019; Yu et al., 2021a). The “inverse kinematics” ( IK ) algorithm, which minimized errors between virtual markers in the model and experimental marker trajectories, was applied to compute joint angles. Then, the “inverse dynamics” ( ID ) algorithm was performed to calculate joint moments, and the “static optimization” ( SO ) algorithm was applied to compute muscle activation and forces. The estimated muscle activation was compared with the measured surface EMG signals to validate the model (Yu et al., 2021a). Lastly, the contact forces towards the knee and ankle joints in the anterior/posterior (x), superior/inferior (y), and medial/lateral (z) directions were computed using “joint reaction” (JR) analysis.

In addition to the biomechanical variables, we also calculated the parameters of contact time, approaching velocity, and peak IMU accelerations (G) in the knee and ankle joints, as shown in Figure 3. The loading rate was calculated following the previous established protocol (Mei et al., 2019; Yu et al., 2021b). The axial (y-axis) acceleration of particular interest was taken for analysis to quantify the accumulation of impact in the lower extremity (tibia) (Rice et al., 2019; Tenforde et al., 2020), which was normalized by gravitational acceleration (G = 9.8 m/s²).

The processed time-varying moment and contact force parameters during the lunging stance (Figure 4) were interpolated (normalized) into a 101 datapoint for statistical modelling (Yu et al., 2021b; Mei et al., 2021). Particularly, the biomechanical parameters included the knee flexion/extension moment, knee varus/valgus moment, knee int/ext rotation moment, ankle dorsi/plantar flexion moment, ankle inversion/eversion moment, knee ant-post/med-lat/vertical contact forces, and ankle ant-post/med-lat/vertical contact forces. The joint moment was normalized to body mass in kg (unit: Nm/kg), and the joint force was normalized to body weight (unit: times Newton in BW).

The discrete values of approaching velocity, stance time, knee and ankle acceleration, peak knee and ankle moments, and peak knee and ankle joint contact forces were first checked for the normality distribution and were compared using the paired sample t-test with a significance level at 0.05. The time-varying joint moment and force over stance were then modelled with the multivariate statistical models. In the current study, the statistical models were developed and validated, as described in our previous studies, using MATLAB software (R2019a, MathWorks Inc., MA, United States of America), specifically PCA (Yu et al., 2021b) and PLSR modelling (Mei et al., 2020).

In this study, the PCA multivariate technique (Wold et al., 1987; Lever et al., 2017; Yu et al., 2021b) was used to reduce the high-dimensional data matrices into orthogonal principal components (PCs), which explained major variations within the dataset (Deluzio et al., 1997; Lever et al., 2017). Each variation reported in the PCA modelling was a feature extraction applied in the machine learning (PLSR) technique (Phinyomark et al., 2018).

As presented in Eq. 1, the original matrices (X = x ¹ , x ² , x ³ , …, x ⁹⁹ , x ¹⁰⁰ , x ¹⁰¹ ) *m were orthogonally transformed into uncorrelated principal components (Z = z ¹ , z ² , z ³ , … , z ^p ) (p < 101), corresponding loading vectors (T ² = T ₁ , T ₂ , T ₃ , …, T _m ), and residuals (Q), which was defined as Z = X*T ² (Deluzio et al., 1997). x 1 1 x 1 2 ⋯ x 1 100   x 1 101 ⋮ ⋮ ⋱ ⋮ x m 1 x m 2 ⋯ x m 100   x m 101 = z 1 1 z 1 2 z 1 3 ⋮ ⋮ ⋮ z m 1 z m 2 z m 3 T 1 2 ⋮ T m 2 Q 1 ⋮ Q m . (1)

The m equals 200 (4*2*101 matrices) for PCA modelling of the four lunges (RF, LF, RB, and LB) with sub-maximal and maximal performance. This study mainly considered the main variations in the first three PCs ( z ¹ , z ² , and z ³ ), which accounted for approximately 85–90% of the variation (Yu et al., 2021b). The variations in the vertical ground reaction force, knee and ankle moments, and contact forces of the first three PCs were then plotted against the mean for the visualization of the key features of variances, with “ + ” and “▽” representing the upper and lower limits, respectively.

In terms of PLSR statistical modelling (Wold et al., 1984; Mei et al., 2020), the two fundamental equations are the predictor matrix ( X _NM ) and the response matrix ( Y _NP ), which are expressed as follows: X N M = T N L   P M L T + E N M , (2) Y N P = U N L   Q P L T + F N P . (3)

The subscript N represents the number of datasets (25*4 training samples in this study). The subscript M represents the number of predictor variables (four metrics, namely, contact times, velocity, peak knee G, and peak ankle G). The subscript P represents the number of response variables (12 loading variables, such as loading rate, knee flex-extension/varus-valgus/int-ext rotation moments, ankle dorsi-plantar flexion/inversion-eversion moments, knee ant-post/med-lat/vertical contact forces, and ankle ant-post/med-lat/vertical contact forces), and the subscript L represents the number of components.

T and U are the projection matrices (also called the scores); P and Q are the transposed orthogonal loading matrices (where the rows are created from eigenvectors or principal components); and E and F are the error or residual terms. The score vectors are related using a linear function, U = f(T) + H, where H is the vector of residuals.

The vertical ground reaction force (Figure 5A) was classified into four key phases: initial impact peak (I), secondary impact peak (II), weight acceptance (III), and drive-off (IV) phases (Figure 5B). Following the PCA modelling of the vertical GRF, the first mode (PC1, Figure 5C) showed the main variations (percentage of variation explained: 31.67%) from landing to the initial impact peak (Phase-I), where the loading rate (LR) was calculated as a key impact parameter. The second mode (PC2, Figure 5D) occurred in the secondary impact peak (Phase-II) (27.58%), and the third mode (PC3, Figure 5E) was the combination of variations (14.1%) in both the initial and secondary impact peaks (Phase-I and Phase-II), which was the impact transient.

Consistent with the variations in vertical GRF, the knee flexion/extension, varus/valgus and int/ext rotation moments showed great variances mainly during the landing (impact absorption) phase. As a particular interest to illustrate key variations, the PC1 of knee flexion/extension (31.45%, impact phase), varus/valgus (32.25%, landing and drive-off phases), and int/ext rotation (52.91%, over the stance) moments is presented in Figure 6, with highlighted regions.

As shown in Figure 7, the key variations in the knee contact forces were observed over the stance with PC1 of ant-post (48.82%) and med-lat (48.36%) forces, while the axial contact force was mainly observed during the mid-stance (39.91%), especially the weight acceptance phase showing the difference between sub-maximal and maximal lunging steps.

The ankle dorsi/plantar flexion and inversion–eversion moments had a principal variance over the stance, especially PC1 (59.86% and 35.03%, respectively) during impact absorption and drive-off phases (Figure 8). It is worth noting the secondary variation in the dorsi/plantar flexion moment during the drive-off phase (PC2: 17.38%), which may explain the difference in ankle contributions during the push-off phase to the return phase. The inversion–eversion moments mainly varied during mid-stance (PC: 24.02%), which was another difference between sub-maximal and maximal lunges.

Similarly for the ankle contact forces, great variations (PC1) in the impact peaks (initial and secondary) and drive-off phases were observed in the ant-post (38.37%), axial (49.45%), and med-lat (49.64%) vectors (Figure 9). Secondary (PC2) variations in the impact peaks were found in the ant-post (29.13%) and axial (20.86%) forces, while the med-lat forces varied across the stance (24.86%).

Following PCA modelling, statistical comparisons of contact times, approaching velocities, peak acceleration (knee), and peak acceleration (ankle) were conducted, as presented in Table 1, which were then input as predictors in the PLSR machine learning model. Sub-maximal lunges showed shorter contact times, smaller approaching velocity, and smaller ankle impact than those in maximal lunges in four directions.

LR, peak knee flexion/extension, varus/valgus, and int/ext rotation moments, peak ankle dorsi/plantar flexion and inver/eversion moments, peak knee ant-post/med-lat/axial (vertical) contact forces, and peak ankle ant-post/med-lat/axial (vertical) contact forces were analyzed, as shown in Table 2. The maximal lunges had greater LR, joint moments, and contact forces overall than those in submaximal lunges. These discrete parameters were then used as response metrics to train the PLSR machine learning model.

Together with the four predictors, a prediction accuracy of 94.52% was observed for the moments and contact forces in the knee and ankle joints (Figure 10). To test the sensitivity of the knee and ankle peak acceleration, we performed a “leave-one-out” cross validation and found that both knee and ankle peak acceleration could predict 93.72% of the loadings. In particular, the knee peak acceleration had an 88.76% prediction accuracy and the ankle peak acceleration had a 93% prediction accuracy.