Estimation of Ground Reaction Forces during Sports Movements by Sensor Fusion from Inertial Measurement Units with 3D Forward Dynamics Model

1. Introduction

In sports such as figure skating, gymnastics, and ballet, where artistic expression plays a significant role, rotational jumps involving twists are indispensable for achieving success and obtaining high scores. One crucial element in these jumping maneuvers is the force generated during takeoff. This takeoff force determines factors such as the jump distance [1], jump height [2], and mid-air posture [3]. The analysis of this takeoff force involves using ground reaction force (GRF) and its moment (ground reaction moment, GRM). Particularly in rotational jumps involving twists, conducting a 3D analysis along with joint torques and muscle forces derived from body movements enable the evaluation of control strategies for optimizing sports movements [4] and injury risks in sports [5], among other aspects. The measurement of these kinematic indicators often involves the use of a force plate fixed to the ground. While force plates can accurately measure GRFs during movements, they require specialized measurement facilities. To ensure successful measurements, participants need to perform movements with their entire foot on the force plate, which may introduce unnatural movements. Consequently, there is a possibility that the measurements may not accurately reflect the same movements as those performed during actual athletic competition. Additionally, in facilities where it is not possible to install multiple expensive force plates, continuous tracking of movements becomes challenging. Moreover, in places like surfaces covered with ice, snow, or spring-based floors where it is necessary to install equipment underneath, introducing measurement facilities can be challenging. The cost of measuring GRFs in sports performed in such environments is substantial, making practical measurements difficult. As a solution, various lightweight and compact measurement devices have been developed. However, these may not be sufficient for sports motion analysis due to concerns about durability, constraints from wearing, and other factors. Wearable force plates, while avoiding facility installation costs, add weight and height to the shoes, limiting their applicability for high-difficulty movements. Insole-type force sensors address some of these challenges but may raise concerns about durability and reproducibility, affecting the credibility of their accuracy [6].

To address these challenges, methods have been studied to measure kinematic elements using only motion capture systems or inertial measurement units (IMUs) and to estimate kinetic elements such as GRF and GRM [7]. In this study, we focused on using IMUs to estimate GRFs. In sports motion measurements, it is advisable to use small, lightweight, and self-contained measurement methods to avoid restricting movement [8]. IMUs meet these conditions, allowing measurements with minimal location constraints and in vast areas, such as outdoors or on a skating rink. Previous studies have explored methods such as the use of inertial sensors and machine learning to estimate the GRFs of gait motion in three dimensions [9], as well as research estimating GRFs during sagittal-plane double limb jumping [10]. While this method demonstrates high prediction accuracy, its application to sports with a limited participant pool, making it difficult to obtain sufficient machine learning training data, poses challenges. Another method involves calculating the posture relative to the absolute coordinate system of each body segment from the measurements of IMUs, and estimating GRFs through inverse dynamics calculations using these values as inputs [11]. However, this method requires the use of 17 IMUs, which may potentially hinder movement in sports motion analysis due to the high number of IMUs. The estimation method that can reduce the number of sensors attached to the body may contribute to advancing sports motion analysis, where the effects of sensor interference on motion should be minimized.

This study focuses on a method for estimating GRFs based on forward dynamics calculations using joint torque [12]. This method generates and optimizes body movements through a 3D human model and cost functions defined by motion measurements and internal load. By utilizing a 3D biomechanical model, this approach can target sports movements that are asymmetrical and non-periodic. Moreover, by generating movements based on the cost function, it is possible to calculate biomechanically valid GRFs while following the measured movements even if not all joints are covered by IMUs, thus enabling the estimation of biomechanically valid GRFs with a small number of IMUs. Haraguchi and Hase [12] directly used acceleration and angular velocity data obtained from IMUs in this cost function. However, in the case of sports movements, there is a possibility of significant noise in the IMU data due to the intensity of the movement. Therefore, in this study, the 3D orientation using sensor fusion composed of acceleration and angular velocity data obtained from IMUs and an extended Kalman filter was estimated, aiming to reduce the influence of this noise by using the estimated values in the optimization calculation.

This study aims to estimate GRFs and body loads using a 3D forward dynamics approach incorporating sensor fusion with an extended Kalman filter, based on measurements of rotational jump movements on land using IMUs.

2. Ground Reaction Force Estimation Methods

This study conducted a biomechanical simulation system that combines forward dynamics simulation using a human model with motion measurement using IMUs. From the generated body movements obtained through forward dynamics simulation, the system estimates GRF, GRM, joint angle and joint torque. The simulation consists of a forward dynamics model and optimization calculations, aiming to minimize the cost function composed of errors between the model and IMUs measurements and evaluations of internal biomechanical loads of the model. This allows the model to reproduce biomechanically plausible body movements while following the measurements from the IMUs. Based on the optimized body movements, GRFs are estimated using a contact model between the human and the ground.

The overall flow of the GRF estimation method is as follows:

• The height and weight are input into the system, and a human model is constructed based on the inertia parameters calculated from them.
• Based on the reference joint angles formed by the node points, the forward dynamics model generates body movements.
• The body movement is evaluated by a cost function comprising errors between the generated body movement and measurements from IMUs, as well as internal biomechanical loads.
• Steps 2 and 3 are repeated while adjusting the node points of the reference joint angles to minimize the value of the cost function.
• Estimate GRFs, GRMs, and joint motion using the optimized node points in the forward dynamics model.

2.1. Forward Dynamics Model

2.1.1. Human Model

The human model is illustrated in Figure 1. The model comprises rigid links representing each segment of the body, with rotational degrees of freedom (DOFs) defined at the junctions to represent joints. The segments include the trunk, pelvis, upper arms, forearms, thighs, shanks, and feet, totaling 12 links. The joints have rotational DOFs: three DOFs for hip, one DOF for knee, one DOF for ankle, three DOFs for lumbar, three DOFs for shoulder, and one DOF for elbow joints, totaling 21 DOFs. Additionally, there are six virtual DOFs accounting for rotational and translational motion between the pelvis segment and the global coordinate system, bringing the total to 27 DOFs. The knee and elbow joints exhibit flexion-extension, while the ankle joint demonstrates plantarflexion-dorsiflexion. The hip, lumbar, shoulder joints and pelvis segment are represented by Euler angles in the X (flexion-extension), Y (adduction-abduction), and Z (rotation) sequence. The lengths, moments of inertia, and center of mass are determined using estimation equations based on height and weight in Japanese [13,14]. The neutral position for each joint was set as shown in Figure 1 and determined by the posture estimation system described in Section 2.2.1

In addition, sensor coordinate systems were defined on the pelvis segment of the human model. This allows the output of relative accelerations to the global coordinate system, enabling a comparison between the motion generated by forward dynamics simulation and the motion obtained from IMU.

In this simulation, the motion is generated by calculating the joint motion when torque is applied to the joints defined in Figure 1 using forward dynamics calculations. At each time step of the forward dynamics calculation, the calculation flow is shown in Figure 2. Additionally, the motion equations of the forward dynamics model are shown in Equation (1).

$$\mathit{M}\ddot{\mathit{q}}+\mathit{\Gamma }\left(\mathit{q},\dot{\mathit{q}},{\mathit{F}}_{\mathrm{v}}, {\mathit{F}}_{\mathrm{h}}, {\mathit{\tau }}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}\right)=\mathit{\tau },$$

(1)

where $\mathit{M}$ is the inertia matrix, $\mathit{\Gamma }$ is the vector consisting of Coriolis forces, centrifugal forces, gravity, and external forces, $\mathit{\tau }$ is the joint torque. $\mathit{q}$ is the state variable vector, including joint angles and the six DOFs between the pelvis and the global coordinate system. ${\mathit{F}}_{\mathrm{v}} \mathrm{and} {\mathit{F}}_{\mathrm{h}}$ are the GRF vectors in the vertical and horizontal directions, respectively, and ${\mathit{\tau }}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}$ is an external torque vector preventing the model from falling. By inputting joint torque $\mathit{\tau }$ into Equation (1), the state variable vector is calculated, generating the body movement of the model. The human model was developed using MATLAB/Simulink (MathWorks, Inc., USA) and employed the ode45 variable-step solver, which implements a Runge-Kutta method.

2.1.2. Joint Torque Model

The joint torque $\mathit{\tau }$ is composed of the active torque ${\mathit{\tau }}_{\mathrm{a}}$ generated by muscle forces and the passive torque ${\mathit{\tau }}_{\mathrm{p}}$ acting prominently near the joint's range of motion limits (Equation (2)).

$$\mathit{\tau }={\mathit{\tau }}_{\mathrm{a}}\left(\mathit{q},\dot{\mathit{q}},{\mathit{q}}_{\mathbf{r}},\dot{{\mathit{q}}_{\mathbf{r}}}\right)+{\mathit{\tau }}_{\mathrm{p}}\left(\mathit{q},\dot{\mathit{q}}\right).$$

(2)

The passive torque is represented by a double exponential function to increase the torque near the anatomical joint limits, which are defined by related research [15,16,17,18]. The active torque ${\tau }_{\mathrm{a}}^j$ for the joint $j$, a component of ${\mathit{\tau }}_{\mathrm{a}}$, is calculated using proportional-derivative (PD) control torque based on the joint angle $q^j$ for the joint $j$, a component of $\mathit{q}$, and reference joint angle $q_{\mathrm{r}}^j$ for the joint $j$, as given by Equation (3).

$${\tau }_{\mathrm{a}}^j=K_{\mathrm{P}\mathrm{D}}\left(q_{\mathrm{r}}^j-q^j\right)+D_{\mathrm{PD}}\left({\dot{q}}_r^j-{\dot{q}}^j\right),$$

(3)

where $K_{\mathrm{PD}}$ and $D_{\mathrm{PD}}$ are gains in the PD control, and they were empirically set to $K_{\mathrm{PD}}=1000 \mathrm{N}·\mathrm{m}/\mathrm{r}\mathrm{a}\mathrm{d}$ and $D_{\mathrm{PD}}=10 \mathrm{N}·\mathrm{m}/\left(\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s}\right)$. The reference joint angle $q_{\mathrm{r}}^j$ is determined through the optimization calculation described in Section 2.2.

2.1.3. External Force Model

Regarding the forces ${\mathit{F}}_{\mathrm{v}}$ and ${\mathit{F}}_{\mathit{h}}$ in Equation (1), they are defined as external forces acting on the foot of the human model and modeled by the contact model. In the human model, there are challenges in the dynamic calculations when both feet are in contact and form a closed-loop structure. In gait analysis, statistical models such as transition functions and distribution functions for GRF and GRM during both-foot support, based on past experiments and data, are often used [9,19]. However, in sports movements, there is limited existing data, and statistical models like those used in gait analysis cannot be defined. Therefore, in this study, a contact model was employed to calculate vertical reaction forces by solving the contact problem between the ground and the foot. The horizontal forces were determined by the frictional forces calculated based on the vertical reaction forces. The contact points are set at 22 points for each foot segment, same as the model by Haraguchi and Hase [12], including distances from the ankle to the heel, metatarsophalangeal joint, and the tip of the second proximal phalanx, as well as heel width, metatarsophalangeal joint width, width between the first and second proximal phalanx, and heights from the ground to the ankle and the first proximal phalanx. The vertical GRF $F_{\mathrm{v}}^i$ at contact point $i$, a component of ${\mathit{F}}_{\mathrm{v}}$, and the horizontal GRF $F_{\mathrm{h}}^i$ at contact point $i$, a component of ${\mathit{F}}_{\mathrm{h}}$, were calculated using Equations (4) and (5).

$$F_{\mathrm{v}}^i=-K_{\mathrm{G}}\left(r_z^i-r_{z_0}\right)-D_{\mathrm{G}}{\dot{r}}_z^i,$$

(4)

$$F_{\mathrm{h}}^i=-\mu \left({\dot{r}}_x^i, {\dot{r}}_y^i,{\mu }_{\mathrm{s}}, {\mu }_{\mathrm{d}}\right)F_{\mathrm{v}}^i,$$

(5)

where $r_x^i, r_y^i,$ and $r_z^i$ represent the position of the contact point, and $r_{z_0}$ is set as the height of the ground ($r_{z_0}=0$). Additionally, $K_{\mathrm{G}}$ and $D_{\mathrm{G}}$ are the elastic and viscous coefficients of the ground, respectively. The friction coefficient $\mu$ is defined as a function of the penetration velocities ${\dot{r}}_x^i$ and ${\dot{r}}_y^i$ at each contact point, along with the static friction coefficient ${\mu }_{\mathrm{s}}$ and dynamic friction coefficient ${\mu }_{\mathrm{d}}$. The values for $K_{\mathrm{G}}$, $D_{\mathrm{G}}$, ${\mu }_{\mathrm{s}}$, and ${\mu }_{\mathrm{d}}$ were determined through trial and error and set to $K_{\mathrm{G}}=1.0\times {10}^3 \mathrm{N}/\mathrm{m}, D_{\mathrm{G}}=1.0\times {10}^4 \mathrm{N}/\left(\mathrm{m}/\mathrm{s}\right), {\mu }_{\mathrm{s}}=0.7, {\mu }_{\mathrm{d}}=0.5$.

Furthermore, to prevent the model from tipping over during the simulation, a virtual torque ${\mathit{\tau }}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}$ was defined to constrain the rotation of the model. In forward dynamics analysis, to maintain the balance of the model is essential. When the model attempted to rotate forward, backward, or sideways to the point of tipping over, the virtual torques are applied to a virtual joint. The virtual joint was set between the mass center position of the pelvis segment and the global coordinate system.

$${\mathit{\tau }}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}=-{\mathit{K}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}} {\mathit{q}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}-{\mathit{D}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}} {\dot{\mathit{q}}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}},$$

(6)

where ${\mathit{q}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}={\left[\begin{array}{ccc}{q}_x& q_y& q_z\end{array}\right]}^T$is the virtual joint angle represented by Euler angles of the pelvis relative to the global coordinates. ${\mathit{K}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}$ and ${\mathit{D}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}$ are gains, determined through trial and error and set as ${\mathit{K}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}={\left[\begin{array}{ccc}400& 400& 0\end{array}\right]}^{\mathrm{T}} \mathrm{N}·\mathrm{m}/\mathrm{r}\mathrm{a}\mathrm{d}$ and ${\mathit{D}}_{\mathrm{v}\mathrm{i}\mathrm{r}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{l}}={\left[\begin{array}{ccc}40& 40& 0\end{array}\right]}^{\mathrm{T}} \mathrm{N}·\mathrm{m}/\left(\mathrm{r}\mathrm{a}\mathrm{d}/\mathrm{s}\right)$.

2.2. Optimization Calculation

The body movement of the forward dynamics model is generated based on the reference joint angle $q_{\mathrm{r}}^j$ used in the joint torque model. Typically, optimization calculations for the reference joint angle are performed at each time step of the simulation. By performing optimization calculations at high frequencies over time steps, the generated motion becomes closer to the measured motion. However, this approach has drawbacks, such as requiring a significant computation time or being susceptible to measurement errors due to large noise of IMUs during certain time steps, like impacts. In this study, to generate motion closer to the measured motion while reducing the frequency of optimization calculations, a method was employed where the reference joint angles were optimized not at each time step but over the entire measurement period. Multiple nodes were used as inputs to the simulation, and their values were determined through optimization calculations to create waveforms of reference joint angles by cubic spline interpolation. The reference joint angles $q_{\mathrm{r}}^j$ were determined as follow:

$$q_{\mathrm{r}}^j\left(t\right)=spline(N_1^j,\dots ,N_{\mathrm{m}}^j,t),$$

(7)

where, $spline$ is a spline function defined by multiple nodes $N_1^j,\dots ,N_{\mathrm{m}}^j$ and the sampling time $t$. The number of nodes $\mathrm{m}$ was determined through trial and error: each upper body joint (shoulder, elbow, lumber) had two nodes at the start and end of the motion, and each joint of the leg involved in the jump (hip, knee, ankle) had five nodes evenly spaced from the start to the end of the motion. Other joints had three nodes set from the start to the end of the motion. The gait motion model proposed by Haraguchi and Hase [12] set node points periodically and symmetrically. However, for the generation of asymmetric and non-periodic movements, node points were set as variable number for all joints with DOFs.

Cost function was established to assess the body movements generated by the forward dynamics simulation, and optimization calculations were performed to minimize this value. Genetic algorithm optimization methods were employed in this study. Cost functions defined by Equation (8).

$$I_{\mathrm{all}}={\xi }_1{I}_{\mathrm{Q}}+{\xi }_2{I}_{\mathrm{a}, \mathrm{pelvis}}+{\xi }_3{I}_{\mathrm{muscle}}.$$

(8)

$I_{\mathrm{Q}}$ and $I_{\mathrm{a}, \mathrm{pelvis}}$ represent the cost functions for evaluating errors in the 3D orientations of each segment and acceleration of the pelvis segment, respectively, between the simulation model and the experimental record by IMUs, as follow:

$$I_{\mathrm{Q}}=\sum _{s=1}^{12}\int \left\{{\left({\varphi }_{\mathrm{model}}^s-{\varphi }_{\exp }^s\right)}^2+{\left({\theta }_{\mathrm{model}}^s-{\theta }_{\exp }^s\right)}^2+{\left({\psi }_{\mathrm{model}}^s-{\psi }_{\exp }^s\right)}^2\right\}\mathrm{d}t,$$

(9)

$$I_{\mathrm{a},\mathrm{pelvis}}=\int \left\{{\left(a_{\mathrm{model},x}^{\mathrm{pelvis}}-a_{\exp ,x}^{\mathrm{pelvis}}\right)}^2+{\left(a_{\mathrm{model},y}^{\mathrm{pelvis}}-a_{\exp ,y}^{\mathrm{pelvis}}\right)}^2+{\left(a_{\mathrm{model},z}^{\mathrm{pelvis}}-a_{\exp ,z}^{\mathrm{pelvis}}\right)}^2\right\}\mathrm{d}t,$$

(10)

where ${\varphi }_{\mathrm{model}}^s$, ${\theta }_{\mathrm{model}}^s$, and ${\psi }_{\mathrm{model}}^s$ are the simulated orientations of sth segment with respect to the global coordinates, represented by Euler angles in the X, Y, and Z sequence, determined by the ${\mathit{Q}}_{model}^s$. ${\varphi }_{\exp }^s$, ${\theta }_{exp}^s$, and ${\psi }_{exp}^s$ are the measured orientations of sth segment with respect to the global coordinates, represented by Euler angles in the X, Y, and Z sequence. These Euler angles are determined by the measured 3D orientations ${\mathit{Q}}_{\exp }^s$ expressed by quaternions in the global coordinate system recorded by IMUs, which is calculated by the extended Kalman filter, as described in Section 2.2.1. $a_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l},x}^{\mathrm{p}\mathrm{e}\mathrm{l}\mathrm{v}\mathrm{i}\mathrm{s}}$, $a_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l},y}^{\mathrm{p}\mathrm{e}\mathrm{l}\mathrm{v}\mathrm{i}\mathrm{s}}$, and $a_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l},z}^{\mathrm{p}\mathrm{e}\mathrm{l}\mathrm{v}\mathrm{i}\mathrm{s}}$are the simulated acceleration of the pelvis segment. $a_{\mathrm{e}\mathrm{x}\mathrm{p},x}^{\mathrm{p}\mathrm{e}\mathrm{l}\mathrm{v}\mathrm{i}\mathrm{s}}$, $a_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l},y}^{\mathrm{p}\mathrm{e}\mathrm{l}\mathrm{v}\mathrm{i}\mathrm{s}}$, and $a_{\mathrm{e}\mathrm{x}\mathrm{p},z}^{\mathrm{p}\mathrm{e}\mathrm{l}\mathrm{v}\mathrm{i}\mathrm{s}}$ are the measured acceleration of the pelvis by the IMUs. ${\mathit{Q}}_{\mathrm{model}}^s$, $a_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l},x}^{\mathrm{p}\mathrm{e}\mathrm{l}\mathrm{v}\mathrm{i}\mathrm{s}}$, $a_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l},y}^{\mathrm{p}\mathrm{e}\mathrm{l}\mathrm{v}\mathrm{i}\mathrm{s}}$, and $a_{\mathrm{m}\mathrm{o}\mathrm{d}\mathrm{e}\mathrm{l},z}^{\mathrm{p}\mathrm{e}\mathrm{l}\mathrm{v}\mathrm{i}\mathrm{s}}$ are calculated after computing the $\mathit{q}$, $\dot{\mathit{q}}$, and $\ddot{\mathit{q}}$ based on the forward dynamics calculation. $I_{\mathrm{Q}}$ aims for alignment with the measurements from IMUs regarding joint motion, $I_{\mathrm{a}, \mathrm{pelvis}}$ aims to match the acceleration from the impact at the moment the foot makes contact with the ground and translational motion of the body. $I_{\mathrm{muscle}}$ is the cost function for the overall muscle load, as follows:

$$I_{\mathrm{muscle}}=\sum _{j=1}^{21}\int \left({\left|{\sigma }_{+}^j\right|}^3+{\left|{\sigma }_{-}^j\right|}^3\right)\mathrm{d}t,$$

(11)

$${\sigma }_{+}^j=\left\{\begin{array}{c}0 :{\tau }_{\mathrm{a}}^j<0\\ \frac{{\tau }_{\mathrm{a}}^j}{{\tau }_{+,\max }^j} :0\le {\tau }_{\mathrm{a}}^j\end{array}\right.,$$

(12)

$${\sigma }_{-}^j=\left\{\begin{array}{c}\frac{{\tau }_{\mathrm{a}}^j}{{\tau }_{-,\max }^j} :{\tau }_{\mathrm{a}}^j<0\\ 0 :0\le {\tau }_{\mathrm{a}}^j\end{array}\right.,$$

(13)

where ${\sigma }_{+}^j$ and ${\sigma }_{-}^j$ are the active torques normalized by the maximum positive joint torque ${\tau }_{+,\mathrm{m}\mathrm{a}\mathrm{x}}^j$ and maximum negative joint torque ${\tau }_{-,\mathrm{m}\mathrm{a}\mathrm{x}}^j$, respectively, for joint $j$. For instance, the maximum flexion torque at the hip joint is ${\tau }_{+,\mathrm{m}\mathrm{a}\mathrm{x}}$ and the maximum extension torque at the hip joint is ${\tau }_{-,\mathrm{m}\mathrm{a}\mathrm{x}}$, which were determined from previous studies across whole joints at neck [20], lumbar [21,22,23], shoulder [24,25,26], elbow [27], hip [28,29,30], knee [24], and ankle [31]. $I_{\mathrm{muscle}}$ contributes to generating biomechanically valid movements while reducing overall physical load. Additionally, $\mathsf{\xi }$ represents the weight factors for each term of the cost function, and they were determined through trial and error as follows: ${\mathsf{\xi }}_1={10}^3 , {\mathsf{\xi }}_2={10}^3 \left(\mathrm{r}\mathrm{a}{\mathrm{d}}^2\right), {\mathsf{\xi }}_3=10 {(\left(\mathrm{m}/{\mathrm{s}}^2\right)}^2)$.

2.2.1. Orientation Estimation

The algorithm for orientation estimation during motion is shown in Figure 3. IMU orientation estimation algorithms are developed, such as utilizing data from 3-axis accelerometers, 3-axis gyroscopes, and 3-axis magnetometers to estimate orientation [32]. However, these methods often struggle to achieve sufficient accuracy when dealing with large dynamic accelerations in sports movements due to compensating for dynamic acceleration as white noise. Additionally, they rely on the magnetic field information obtained from magnetometers, making them challenging to use in indoor sports settings prone to magnetic disturbances. Therefore, in this proposed method, 3D orientations were estimated using sensor fusion of the 3-axis accelerometer and 3-axis gyroscope data from the IMUs [33,34]. It thereby enables estimation of 3D orientation without the effects of drift errors due to integration of angular velocities, dynamic accelerations, and geomagnetic disturbances. This sensor fusion approach is constructed using an extended Kalman filter to address the nonlinearity inherent in the relationship between 3D orientations and acceleration.

The output of acceleration in an IMU ${\mathit{a}}_{\exp }^s$ attached to the segment is represented as the sum of translational, centripetal, tangential, and gravitational, and Coriolis accelerations. In this study, the Coriolis acceleration can be neglected due to the constant distance from the joint to the IMU. Additionally, the gravitational and translational accelerations can be obtained by subtracting the centrifugal and tangential accelerations from the sensor output. Considering the entire system, IMU sensors provide the same gravitational and translational accelerations. Therefore, expressing centrifugal and tangential accelerations vector of sth segment at the sampling time $t$ as ${\mathit{a}}_{\exp ,\mathrm{ct}}^s$, this relationship was expressed as in the observation equation of the extended Kalman filter (Equation (15)). Furthermore, the state values ${\mathit{x}}^s={\left[\begin{array}{cc}{\mathit{Q}}_{\mathit{e}\mathit{x}\mathit{p}}^s& {\mathit{Q}}_{\mathit{e}\mathit{x}\mathit{p}}^{s+1}\end{array}\right]}^{\mathrm{T}}$ were calculated by integrating angular velocity at the sampling time $t$ to calculate the 3D orientations at the next sampling time [35] (Equation (14)). In this study, the 3D orientations ware described using quaternions that are not affected by gimbal lock. The quaternion represents the real part as the first element and the vector of imaginary parts as the subsequent elements.

$${\mathit{x}}^s(t+1)=\left(I+\frac{dt}{2}\left[\begin{array}{cc}{\mathit{\Omega }}^s& \mathit{O}\\ \mathit{O}& {\mathit{\Omega }}^{s+1}\end{array}\right]\right){\mathit{x}}^s\left(t\right)+{\mathit{w}}^{\mathit{s}},$$

(14)

$$\left[\begin{array}{c}{\mathit{a}}_{\exp }^s-{\mathit{a}}_{\exp ,\mathrm{ct}}^s\\ {\mathit{a}}_{\exp }^{s+1}-{\mathit{a}}_{\exp ,\mathrm{ct}}^{s+1}\end{array}\right]=\left[\begin{array}{c}{}_{s+1}{}^s\mathit{R}\left({\mathit{x}}^s\left(t\right)\right)\left\{{\mathit{a}}_{\exp }^{s+1}-{\mathit{a}}_{\exp ,\mathrm{ct}}^{s+1}\right\}\\ {}_s{}^{s+1}\mathit{R}\left({\mathit{x}}^s\left(t\right)\right)\left\{{\mathit{a}}_{\exp }^s-{\mathit{a}}_{\exp ,\mathrm{ct}}^s\right\}\end{array}\right]+{\mathit{v}}^{\mathit{s}},$$

(15)

$${\mathit{a}}_{\exp ,\mathrm{ct}}^s={\mathit{\omega }}_{\exp }^s\times \left({\mathit{\omega }}_{\exp }^s\times {\mathit{l}}^s\right)+{\dot{\mathit{\omega }}}_{\exp }^s\times {\mathit{l}}^s,$$

(16)

where ${}_{s+1}{}^s\mathit{R}$ is the rotation matrix from segment $s$ to $s+1$. ${\mathit{\omega }}_{\exp }^s$ is the angular velocity vector of sth segment at the sampling time $t$, expressed as ${\mathit{\omega }}_{\exp }^s={\left[\begin{array}{ccc}{\omega }_{\exp ,x}^s& {\omega }_{\exp ,y}^s& {\omega }_{\exp ,z}^s\end{array}\right]}^T$, and ${\mathit{l}}^s$ represents the position vector indicating the location of the IMU attached sth segment from the proximal joint $j$. ${\mathit{\Omega }}^s$ is expressed in terms of angular velocities, as follows:

$${\mathit{\Omega }}^s=\left[\begin{array}{cccc}0& -{\omega }_{\exp ,x}^s& -{\omega }_{\exp ,y}^s& {-\omega }_{\exp ,z}^s\\ {\omega }_{\exp ,x}^s& 0& {\omega }_{\exp ,z}^s& {-\omega }_{\exp ,y}^s\\ {\omega }_{\exp ,y}^s& -{\omega }_{\exp ,z}^s& 0& {\omega }_{\exp ,x}^s\\ {\omega }_{\exp ,z}^s& {\omega }_{\exp ,y}^s& -{\omega }_{\exp ,x}^s& 0\end{array}\right].$$

(17)

In this study, these orientation estimation systems were adapted to the adjacent IMUs across each joint. Additionally, the process noise in the extended Kalman filter was set to ${\mathit{w}}^s={\left[{10}^{-3} {10}^{-3} {10}^{-3} {10}^{-3} {10}^{-3} {10}^{-3} {10}^{-3} {10}^{-3}\right]}^{\mathrm{T}}$, and the observation noise to ${\mathit{v}}^{\mathit{s}}={\left[{10}^{-1} {10}^{-1} {10}^{-1} {10}^{-1} {10}^{-1} {10}^{-1}\right]}^{\mathit{T}}\left(\mathrm{m}/{\mathrm{s}}^2\right)$. The initial value of the posterior error covariance matrix ${\mathit{P}}_t$ was set to ${\mathit{P}}_{t_0}=\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\{{10}^2, {10}^2, {10}^2, {10}^2, {10}^2, {10}^2, {10}^2, {10}^2\}$. These values were determined through trial and error, taking into account the magnitude of noise from the accelerometers and gyroscope sensors.

The initial state values of Equation (14) and the neutral position of the model were determined using IMUs data from the calibration posture and this posture estimation system.

3. Experiment

3.1. Participants

The participants for measurement included seven healthy adult males and three healthy adult females (average body height: $1.69\pm 0.10 \mathrm{m}$, average body weight: $63.3\pm 12.7 \mathrm{k}\mathrm{g}$, average age: $23.5\pm 2.5$ years). This study has been approved by the Ethics Committee of Tokyo Metropolitan University. Participants were provided with written and verbal explanations of the measurement details before the measurements, and their consent was obtained through signed agreement documents before conducting the measurements.

3.2. Conditions

The jumping movement involved a run-up followed by a jump with a takeoff using the left leg. The run-up motion began in the direction of the jump, and on the second step of the run-up, a takeoff motion was performed, followed by a counterclockwise half-turn jump. The distance of the run-up was adjusted so that the participants could jump naturally. The jumping distance was determined with consideration for safety and the participant's athletic ability, with an average ± standard division (SD) value of $0.52\pm 0.19 \mathrm{m}$. Participants familiarized themselves with the movement before the experiment.

3.3. Measurements

The IMUs (TSND151, ATR-Promotions Inc., Japan) was used to measure 3-axis acceleration and 3-axis angular velocity (sampling frequency: 1,000 Hz, accelerator measurement range: ±16 G (resolution: 0.48 mG), gyroscope measurement range: ±2,000 dps (resolution: 0.061 dps)). 12 IMUs were attached to the participant’s trunk, pelvis, upper arms, forearms, thighs, shanks, and feet. The attachment points were as follows: sensors on the trunk were positioned on the upper thoracic vertebrae, those on the pelvis were at the midpoint of the bilateral posterior superior iliac spines, those on the upper arms were at a quarter of the distance from the elbow joint along the line connecting the shoulder and elbow joints, those on the forearms were at the midpoint along the line connecting the elbow joint and the styloid process of the ulna, those on the thighs were at a quarter of the distance from the lateral condyle of the femur along the line connecting the greater trochanter and the lateral condyle of the femur, those on the shanks were at the midpoint along the line connecting the lateral malleolus and the medial malleolus of the tibia, and those on the feet were at the midpoint along the line connecting the calcaneus and the base of the first metatarsal bone. Two cylindrical rods were attached to the IMUs on the forearms, and they were worn to sandwich the radius and ulna. All IMUs were secured to each segment using elastic belts.

To validate the simulation, GRFs were measured using a force plate (TF-4060-D, Tech Giken Co., Ltd., Japan) at a sampling frequency of 100 Hz. Additionally, the 3D coordinates of markers attached to the entire body were measured at a sampling frequency of 100 Hz using an optical motion capture system (OptiTrack Flex 3, Natural Point Inc., USA). Marker positions and calibration poses were aligned with a full-body musculoskeletal model [36] used for calculating joint angles and joint torques. All measurements underwent processing with a low-pass filter (Butterworth fourth-order type, -3 dB at 18 Hz).

3.4. Data Analysis

The measurements from the IMUs were input into the forward dynamics model. To validate the model, GRF, GRM, joint angle and joint torque were calculated using force plates and an optical motion capture system as the reference values. GRF and GRM were computed from the force plate measurements. The point projected onto the ground as the origin was the ankle joint center, and the calculations were performed in a coordinate system aligned with the foot's horizontal direction. The ankle joint center was determined as the midpoint between markers on the medial and lateral sides of the heel, obtained from the motion capture system. Additionally, these values were normalized by body weight (BW) and the product of BW and body height (BH). Furthermore, using the motion capture system and the open-source musculoskeletal modeling software OpenSim 4.4 [37,38], along with the full-body musculoskeletal model [36], inverse kinematics and inverse dynamics calculations were performed to compute joint angles and joint torques, including virtual joints. The joint torques were normalized by the product of BW and BH. These results are normalized from the moment the stepping foot makes contact with the ground until it leaves the ground.

The concordance between the estimated and true values is assessed using the Pearson correlation coefficient ($\rho$), categorized as weak correlation for $\rho \le 0.35$, moderate correlation for $0.35<\rho \le 0.67$, strong correlation for $0.67<\rho \le 0.9$, and excellent correlation for $0.9<\rho$ [39]. Additionally, the root mean square error (RMSE) and relative RMSE (rRMSE) were calculated. rRMSE normalizes RMSE by the average values of the two data ranges.

4. Results

Figure 4 shows the time-series data of GRFs ((a) the medial direction, (b) the anterior direction, and (c) the vertical direction) and GRM ((d) the sagittal plane, (e) the frontal plane, and (f) the transverse plane). Strong correlations were observed in the vertical GRF ($\rho =0.749, \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=63.5 \%\mathrm{B}\mathrm{W},r\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=27.4 \%$), moderate correlations in the sagittal GRM ($\rho =0.520, \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=5.67 \%\mathrm{B}\mathrm{W}\mathrm{･}\mathrm{B}\mathrm{H},r\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=61.6 \%$), weak correlations in the medial GRF ($\rho =-0.0326,\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=32.1\%\mathrm{B}\mathrm{W},r\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=72.9 \%$), anterior GRF ($\rho =-0.206, \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=44.5 \%\mathrm{B}\mathrm{W},r\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=158 \%$), frontal GRM ($\rho =-0.481, \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=9.38 \%\mathrm{B}\mathrm{W}\mathrm{･}\mathrm{B}\mathrm{H},r\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=84.0 \%$), and transverse GRM ($\rho =0.157, \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=168 \%\mathrm{B}\mathrm{W}\mathrm{･}\mathrm{B}\mathrm{H},r\mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=73.7 \%$).

Table 1 shows the estimation accuracy for joint angles and joint torque. Strong correlations were observed in the flexion angles of the hip, knee, and ankle joints. Weak correlations were observed in the abduction and external rotation angles of the hip, and all joint torques.

5. Discussion

5.1. Advantages of the Proposed Approach

In this estimation method, accurate estimation was achieved for joint angles (hip, knee, and ankle flexion angles) with relatively large displacement. This was accomplished by optimizing biomechanically plausible flexion angles using forward dynamics calculations and the cost function defined by motion measurements and internal biomechanical loads. The correct generation of vertical foot movement resulted in high precision in vertical GRF estimation and slightly improved accuracy in sagittal GRM estimation.

In prior research conducted on estimating GRF in jump movements using IMUs and biomechanical models, 17 IMUs were utilized [40]. Additionally, methodologies employing machine learning require a substantial amount of experimental data for training purposes [10]. However, employing a hybrid cost function including internal biomechanical loads and less noise data estimated from sensor fusion with extended Kalman filter, our system requires fewer IMUs and eliminates the need for extensive training data. It enabled cost-effective estimation of GRF, GRM, and joint movements in sports movement. By reducing the hindrance caused by sensor attachment, our proposed estimation approach contributes to sports movement analysis. Moreover, employing a hybrid cost function including internal biomechanical loads necessitates accuracy validation, but also offers the potential to further reduce the number of required IMUs.

5.2. Accuracy

Although this study achieved high accuracy in estimating GRFs and joint movements during rotational jump actions, the precision attained was not sufficient for practical application in sports competitions.

This study employed IMUs to measure motion. IMUs are susceptible to noise in measurement and require robust attachment and careful calibration to avoid inhibiting movement. In this study, 12 IMUs were manually attached according to predefined positions, which may have introduced errors in the estimated GRFs due to positional inaccuracies. Furthermore, movement of segments during jumping and noise from soft tissues could have adversely affected the estimation of GRFs. IMUs measure acceleration and angular velocity of soft tissues relative to bones [41]. During measurement, attachment positions were chosen to minimize the influence of soft tissues and were secured with elastic straps. However, in vigorous activities such as sports movements, IMUs within soft tissues or within the belt might experience significant vibrations, potentially leading to considerable errors in measurements.

When compared to the gait analysis by Haraguchi and Hase [12], the accuracy of anterior and medial GRFs and frontal and transverse GRMs decreased. This is likely due to the lack of consideration for the attenuation of impact acceleration by joints and the virtual torque for fall prevention. The 3D orientation of measured body movement used in the cost function was estimated by sensor fusion of acceleration and angular velocity information obtained from IMUs. The extended Kalman filter for sensor fusion was constructed using the relationship between the two IMUs attached to the two-link system, enabling it to track high-speed movements such as jumping. However, the lack of consideration for the attenuation of impact acceleration by joints may lead to errors in the sum of translational acceleration and gravitational acceleration for each link, particularly during instances of high impact such as ground contact or takeoff in sports movements, resulting in potential errors in the estimated 3D orientation and adversely affecting the optimization of body movements in the model. In joints with a small range of motion, errors due to impact acceleration can have a significant effect on the displacement of the 3D orientation. Consequently, errors may occur in the estimation of non-flexion joint angles, leading to inaccurate foot contact and potentially reducing the accuracy of medial GRF and transverse GRM. Furthermore, while Haraguchi and Hase added virtual forces for fall prevention to the lateral direction based on the relationship between the overall center of gravity position and the support base, this study introduced virtual torques corresponding to the rotation angles of virtual joints to enable the generation of unstable movements such as jumping with a single foot while preventing falls. The addition of virtual torques in the sagittal and frontal plane rotations of the body might have limited movements. The decrease in accuracy of anterior GRF and frontal GRM may be attributed to the application of forces akin to brakes in those directions through virtual torques. The system cannot eliminate the virtual torques for fall prevention of the model. Therefore, reconsidering the reduction of the impact of virtual torques is necessary to improve the accuracy of body movements and joint torques.

5.3. Limitations

This study utilized healthy non-athlete participants. When adapting the system to athletes, there is a potential for errors in the estimated GRFs. Currently, the model is based on inertial parameters determined using estimation formulas based on height and weight and internal values such as joint torques and muscle forces. They were obtained from a population that does not represent athletes. However, when targeting sports movements with athletes as participant, it may be necessary to use parameters specific to the participants or the average values of athlete populations in that sport [42].

In this study, focusing on rotational jump movements, a total of 55 node points were set as parameters for the optimization calculation. These node points are used to calculate the reference joint angle $q_{\mathrm{r}}^j$, which serve as inputs to the simulation. Increasing the number of node points may enable the generation of more complex body movements and enhance the accuracy of ground reaction force estimation. However, increasing the variables, represented by node points, in the optimization calculation may negatively impact convergence speed and precision. To reduce measurement and computation costs while improving measurement accuracy and optimization calculation precision, kinematic synergy can be considered. Humans smoothly control the redundant multi-degree-of-freedom musculoskeletal system through a mechanism known as kinematic synergy, reducing the controlled DOFs [43]. This control mechanism flexibly manages joint motion through groups with synergies, allowing for the generation of diverse body movements. By introducing kinematic synergy into this system and controlling reference joint angle $q_{\mathrm{r}}^j$ for groups with synergies, it may be possible to achieve the generation of body movements and estimate GRFs even with fewer IMUs.

While the accuracy of vertical GRF was high in this system, the precision of anterior GRF and frontal GRM was lower. Although it may not predict the dynamics of all sports, it could use in sports involving sliding, such as those on ice or snow. For example on ice, anterior GRF tends to be small, and the frontal GRM is also diminished due to the narrow width of the blade [44]. Moreover, considering the possibility that virtual torque may have contributed to these low estimation accuracies, this method could be applicable to figure skating rotation jump movements, where these values are small. If this estimation method can be applied to figure skating, it would enable the estimation of forces exerted on the ice surface and physical loads with simple measurements, facilitating their application in movement analysis.

This system requires approximately 24 hours for a single trial (CPU: Intel® Core i9-13900KF, with 24 cores, 32 logical processors, and an average speed of 3.0GHz; memory: 128GB; software: MATLAB R2023a (9.14); operating system: Windows 11 Pro). Therefore, real-time analysis of sports movements is not feasible.

6. Conclusions

We estimated GRFs and body loads using a 3D forward dynamics approach incorporating sensor fusion with an extended Kalman filter, based on measurements of rotational jump movements on land using IMUs. By generating movements based on the cost function, it was possible to calculate biomechanically valid GRFs while following the measured movements even if not all joints are covered by IMUs, thus enabling the estimation of biomechanically valid GRFs with a small number of IMUs. To reduce the influence of the significant noise of IMUs due to sports movements in the optimization calculation, the 3D orientation using sensor fusion composed of acceleration and angular velocity data obtained from IMUs and an extended Kalman filter was estimated. This estimation method enables 3D motion analysis independent of measurement environment and machine learning training data.