1. Introduction
In the fields of rehabilitation and sport, analyzing ground reaction forces (GRFs) and ground reaction moments (GRMs) offers advantages in assessing the risk of musculoskeletal disorders and evaluating physical performance [
1,
2]. To address the challenges associated with the complexity of recording GRFs and GRMs, previous studies have developed simplified measurement techniques that do not require a force plate. One such technique involves wearable systems for recording GRFs and GRMs, such as pressure insoles or instrumented force shoes [
3]. These systems offer advantages in recording GRFs and GRMs with fewer limitations on measurement location. However, concerns remain regarding their low durability and the impact of the weight and height of the wearable instruments on the motion of participants.
Another technique for measuring GRFs involves a prediction approach using kinematic data recorded by optical motion capture (OMC) or inertial measurement units (IMUs) [
3,
4]. The OMC-based method contributes to high accuracy in GRF estimation because it predicts GRFs based on the body-segment positions obtained with high precision by OMC which correctly records three-dimensional positions. While the IMU-based method faces challenges in accurately estimating GRFs from kinematics data, as IMUs cannot directly record body-segment positions, they have the advantage of fewer limitations on measurement location. Researchers choose measurement systems based on the accuracy and applicability required for their estimation system, as well as the computational approach. Many GRF estimation techniques include two computational approaches: the statistical-model-based method or the biomechanical-model-based method. The statistical-model-based approach enables the prediction of GRFs and GRMs with high accuracy, relying on statistical algorithms such as machine learning [
5,
6,
7,
8,
9,
10]. Recently, a GRF estimation method combining electromyography measurements and machine learning has also been proposed [
11], and the statistical-model-based approach is rapidly developing as an estimation technique. However, this method has a disadvantage in its application to situations where obtaining sufficient training data is challenging, thus limiting the range of movements that can be analyzed.
The biomechanical-model-based approach enables the estimation of GRFs and GRMs using multibody dynamics calculations [
3,
4]. Although this method has the potential to be applied in various movements, it faces limitations in making accurate predictions during the double-support phase due to the closed-loop structure between humans and the ground. To address this challenge, several studies have proposed computational techniques within the biomechanical-model-based approach. While artificial neural networks and smooth transition assumptions have been employed as accurate estimation methods [
12,
13], they lack applicability due to the need for training data and gait assumptions. In contrast, some studies have proposed computational techniques that incorporate a contact model in which the contact state between the foot and the ground is represented by a mathematical model, such as a spring-damper model [
14,
15]. Since this method calculates GRFs and GRMs from the measured contact state between the foot and the ground, it does not require any training data or gait assumptions, making it applicable to activities other than walking [
16,
17].
Previously, we developed the contact-model-based approach for predicting GRFs, GRMs, and joint torques [
18]. This method employs multibody dynamics calculations based on kinematic data for GRF prediction without requiring training data or making gait assumptions, thus providing advantages in terms of versatility and applicability. However, our study faced a limitation of significant computational cost due to the optimization technique to tune the parameters necessary for estimating GRFs and GRMs. Similarly, other studies using the contact-model-based approach also rely on optimization processes to tune parameters such as the coefficients of the spring-damper model [
14,
15,
16,
17]. While alternative methods exist for estimating GRFs and GRMs using different computational techniques [
19,
20], these methods also incorporate optimization techniques for GRF prediction. Therefore, a prediction method focused on reducing computational time could significantly contribute to the clinical applications of GRF estimation, such as improving the efficiency of the biomechanical analysis and rapidly presenting analysis results to medical staff and patients.
This study developed a new prediction method for GRFs, GRMs, and joint torques with low computational cost by focusing on deformations of the foot alignment. Previous studies have reported that the foot arch and soft tissues contribute to foot deformations under load, reaching magnitudes of millimeters or centimeters under the load of the human body mass [
21,
22]. Based on this fact, we hypothesized that foot deformations can be recorded with OMC, which correctly obtains three-dimensional positions, and that the load condition of the foot can be predicted from the recorded foot deformations. Furthermore, if the load condition can be estimated from foot deformations, GRFs and GRMs during the double support phase can be predicted without optimization techniques by distributing external forces calculated by the biomechanical model according to the load condition. Therefore, this study developed a new OMC-based prediction method for GRFs, GRMs, and joint torques using a hybrid approach of a multibody-dynamics model and foot deformations and evaluated the estimation accuracy and its advantages for biomechanical analysis.
2. Materials and Methods
2.1. Participants
The study participants were 18 volunteers (10 males, 8 females; mean age: 23.3 ± 2.5 years; mean height: 1.67 ± 0.10 m; mean weight: 55.2 ± 7.46 kg) with no history of musculoskeletal disorders. Approval for the study was obtained from the Ethics Committee of Tokyo Metropolitan University. Prior to the experiment, all participants received both verbal and written explanations regarding the study’s content, and they provided written informed consent.
2.2. Conditions
As the initial step in developing the prediction method, this study evaluated the accuracy of estimation during level walking with bare feet, a fundamental human motion. The experiment included three walking speeds: normal, fast, and slow. The normal speed was adjusted to each participant’s preferred walking pace, while the fast and slow speeds were approximately +20% and −20% of the normal walking speed, respectively. Before the experiment, participants were familiarized with each walking speed.
2.3. Measurements
The three-dimensional coordinates, denoted as
x for the anterior−posterior axis,
y for the vertical axis, and
z for the medio−lateral axis, of markers attached to the whole body were measured using an OMC system (OptiTrack Flex3; Natural Point Inc., Corvallis, OR, USA). A total of 49 reflective markers were placed on various body locations based on the marker set of the Gait 2392 musculoskeletal model in OpenSim [
23,
24,
25,
26]. The marker data were digitally filtered (low-pass filter, Butterworth fourth-order type, −3 dB at 6 Hz) and were sampled at 100 Hz.
To validate the OMC results, the GRF and GRM were measured using a force plate (TF-4060-D; Tec Gihan Co. Ltd., Kyoto, Japan). The force data were digitally filtered (low-pass filter, Butterworth fourth-order type, −3 dB at 18 Hz) and were sampled at 1000 Hz.
2.4. Data Analysis
2.4.1. Kinematics
At the beginning of the prediction procedure, the kinematics of each joint, segment, and contact point were computed based on the three-dimensional coordinates of markers recorded by the OMC system. The joint and segment kinematics were calculated using the Gait 2392 musculoskeletal model in OpenSim [
23,
24,
25,
26]. This model consists of 12 segments, including the torso, pelvis, femurs, tibias, taluses, calcanei, and toes, and had 23 degrees of freedom (DOFs). The segment reference frame and joint coordinate system were defined based on previous studies in the lower extremity [
23], plantar knee [
24], and low back model [
25]. The hip and back joints had three DOFs expressed in Euler angles with the sequence
Z (flexion−extension)
X (abduction−adduction)
Y (rotation). The knee, ankle, subtalar, and metatarsophalangeal joints each had one DOF in flexion−extension or dorsiflexion−plantarflexion. The pelvis had six DOFs for global coordinates, which included three translational DOFs and three rotational DOFs expressed in Euler angles with the sequence
Z (flexion−extension)
X (abduction−adduction)
Y (rotation). The inverse kinematics using the model in OpenSim computed the joint kinematics vector
, which includes the joint angles of all joints, and the segment kinematics vectors
,
, and
, which include the translational position of the center of mass (COM) of all segments [
27].
To compute GRFs and GRMs based on foot deformations, a total of 20 contact points (10 per foot) were positioned on the calcaneus and toe segments, as shown in
Figure 1. The segmentations of the calcaneus and toe was defined based on the Gait 2392 musculoskeletal model in OpenSim [
23]. The horizontal locations of the contact points were derived from the marker positions based on the marker set of the Gait 2392 musculoskeletal model in OpenSim. The vertical locations of the contact points were set to have an offset of 30 mm toward the ground during static standing, using the OMC data of the static standing calibration with weight on the heel side. Subsequently, the contact point computation yielded translational positions of the contact points on the foot, represented by
,
, and
.
2.4.2. Ground Reaction Forces and Moments
Before computing GRFs and GRMs, the total external forces acting on the human body were calculated using the translational equations of motion. The external forces
,
, and
were expressed as follows:
where
,
, and
represent the anterior, vertical, and medial position of the COM of the
sth segment, which are components of
,
, and
, respectively;
is the mass of the
sth segment as defined in the previous study [
25]; and
is gravitational acceleration.
Because the external forces, calculated from the equation of motion, represent the total forces applied to the human body, it is necessary to allocate these forces to each foot to calculate the GRFs and GRMs. The present study determined the amount of sinkage to the ground at the
ith contact-point position
based on the contact-point positions, and the total external force was distributed to each contact point based on the ratio of
to the overall sinkage. The vertical GRF
, frontal GRM
, and sagittal GRM
applied to one foot were computed as follows:
where
,
, and
represent the anterior, vertical, and medial position of the
ith contact point, which are components of
,
, and
, respectively, and
and
are the critical position and velocity of the contact point, respectively, which were set to
= 0 and
= 0.05, as determined empirically. In this study, the ground height was set to
= 0. The origin of the horizontal directions
and
was positioned at the ankle joint, as measured by OMC (the midpoint of the markers of the lateral and medial malleolus).
During the single-support phase, the anterior GRF
and medial GRF
applied to one foot were determined using the horizontal external forces calculated from the equation of motion, as follows:
During the double support phase, the anterior and medial GRFs were calculated using a different method compared to the single-support phase. Previous studies have reported that the GRF vector intersects a point located above the COM of the system, referred to as the virtual pivot point (VPP), in a dynamic system such as a human in motion [
28]. Accordingly, the present study computed the horizontal GRF based on the assumption that the GRF vector intersects the VPP of the human body. First, three-dimensional positions of the center of pressure (COP) of the GRF, which is the starting point of the GRF vector, are calculated by the vertical GRF and the GRM around the horizontal axis, as follows:
where the vertical position of the COP was set to the ground height, that is,
= 0. Then, if the GRF vector is assumed to intersect the VPP, the GRF can be expressed as a vector connecting the COP and the VPP multiplied by a real number. Consequently, the anterior GRF
and medial GRF
during the double support phase are calculated using the real number
, as follows:
where
,
, and
represent the anterior, vertical, and medial positions of the VPP of the entire body, respectively. The VPP position was defined as 37.5 mm upward along the axis direction represented in the trunk segment coordinate system from the COM of the whole body, after a previous experimental study [
28]. The whole-body COM was determined from the segment kinematics vectors
,
, and
. Subsequently, the transverse GRM
was computed using the horizontal GRFs and the COP of the GRF, as follows:
The computational flow of GRFs and GRMs is summarized as follows. First, the positions of the COM of each segment and the contact points on the foot segment are derived from the OMC data. Based on these COM kinematics, the equation of motion determines the external forces acting on the human body. Through the foot deformation approach, these external forces are distributed by the positional relationship between the contact point and the ground, obtaining the vertical GRF, sagittal GRM, and frontal GRM. Following this computation, the anterior GRF and medial GRF, and COPs are calculated using different approaches for the single-support phase and double support phase. Finally, the transverse GRM is derived from the computed GRFs and COPs.
2.4.3. Joint Torques
Joint torques were determined through inverse dynamics analysis using the estimated GRFs and GRMs. The joint torque vector
was computed as follows:
where
represents the inertia matrix as defined in the previous study [
25], and
is a vector consisting of Coriolis, centrifugal, gravitational, and external forces. The inverse dynamics analysis was conducted with the Gait 2392 musculoskeletal model in OpenSim [
27], which is the same model used in the inverse kinematics analysis.
2.5. Accuracy and Sensitivity Analysis
To evaluate the prediction accuracy, the predicted GRFs and GRMs were compared with the force-plate data. The predicted joint torques were compared with an inverse dynamics solution using the Gait 2392 musculoskeletal model in OpenSim with the input measurement data [
27]. The agreement between the prediction and measurement was evaluated using Pearson’s correlation coefficient (
), which was classified as follows:
for weak,
for moderate,
for strong, and
for excellent correlation [
29]. In addition, we computed the root-mean-square error (
) and relative
(
), which normalized the
by the average peak-to-peak amplitude for the two solutions [
13]. All statistical analyses were performed using MATLAB 9.9 (MathWorks, Inc., Natick, MA, USA).
In the proposed method, which relies on foot deformation, the measurement error of the contact point may have a significant impact on the estimation accuracy. Therefore, a sensitivity analysis was conducted to investigate the influence of the measurement error of the contact point on prediction accuracy. To simulate data with measurement errors, white noise was added to all contact-point positions, , , and using the rand function in MATLAB. This study compared the RMSE and rRMSE of the GRFs and GRMs under three conditions with maximum white noise amplitudes of 1, 10, and 100 mm.
3. Results
The GRF and GRM curves throughout a gait cycle during normal-speed walking are shown in
Figure 2. At all three experimental walking speeds, excellent or strong correlations were observed for all GRFs and GRMs. A comparison of the estimation accuracies in GRFs and GRMs with those of previous studies is presented in
Table 1.
The joint torque curves throughout a gait cycle during normal-speed walking are shown in
Figure 3. At all three experimental walking speeds, excellent or strong correlations were observed for all joint torques. A comparison of estimation accuracies in joint torques with those of previous studies is presented in
Table 2.
The RMSEs and rRMSEs of the GRFs and GRMs with white noise in contact-point position data during normal-speed walking are presented in
Table 3. The sensitivity analysis confirms that the estimation accuracies of GRFs and GRMs tended to deteriorate as the amplitude of the white noise increased.
The computational time from OMC data input to GRFs and GRMs output was approximately 4 s for one trial (CPU: Intel Core i7-10700 with 4.0 GHz of average speed; memory: 32 GB; software: MATLAB 9.9 and OpenSim 4.1; OS: Windows 11 Home).
The results of the statistical analysis for each speed are summarized in
Table S1. The time-series data for each speed are included in
Dataset S1.
4. Discussion
The proposed method demonstrated estimation accuracy comparable to that of previous studies, as shown in
Table 1 and
Table 2 [
12,
13,
16]. Thus, the proposed approach, which does not rely on optimization techniques, succeeded in achieving these estimation accuracies for GRFs, GRMs, and joint torques at a low computational cost of a few seconds. In previous studies that employed optimization techniques for prediction, the computational time was relatively long because the parameters for computing the estimated values had to be determined through exploratory calculations [
18]. Although the computational time depends on various factors, it is assumed that optimization approaches requiring exploratory calculations involve a computational cost at least in the order of several minutes. Therefore, the proposed method offers advantages over conventional methods, such as rapid feedback of the analysis results, thereby contributing to the clinical application of these kinetic factor estimation methods.
Another feature of the proposed method is its potential applicability as a GRF estimation technique for various activities. In previous studies, the smooth transition assumption approach has been applied to GRF estimation only during walking [
13,
30]. Although the artificial neural network approach can be extended to GRF estimation for activities other than walking [
31], its applicability might be constrained in situations where sufficient training data cannot be obtained. In contrast, because the present method operates independently of statistical models or walking assumptions, it has the potential to be applied to GRF estimation across a wide range of activities, such as sports, daily activities, and orthopedic or prosthetic conditions.
Here, we provide a detailed description of the prediction accuracy of the proposed method. Although previous studies estimated GRFs during the single-support phase using an approach relying on the equation of motion as used in our study, GRFs during the double support phase were calculated using complex computational techniques, such as a statistical model, walking assumptions, and optimization [
12,
13,
14,
15,
16,
17,
18,
19,
20]. In contrast, our proposed estimation method calculated GRFs using only the external force derived from the equation of motion and the amount of sinkage at the contact point. Consequently, the prediction results of this study were more susceptible to the influence of measurement errors compared to conventional approaches. Hence, the simplicity of the computational method employed in this study had a detrimental impact by producing a lower prediction accuracy for the GRFs.
Previous studies estimated the GRMs using the equation of motion and some computational techniques, similar to their GRF prediction [
12,
13,
14,
15,
16,
17,
18,
19,
20]. However, the estimated GRMs based on the rotational equation of motion could include modeling errors related to the mass, moment of inertia, and COM of each body segment. These errors potentially degrade the estimation accuracy compared to the GRF prediction based on the translational equation of motion, which only incorporates modeling errors of the mass of each body segment. On the contrary, the proposed estimation method computed the GRMs using the estimated GRFs and the position of the contact point, thereby eliminating the modeling errors associated with the moment of inertia and COM of each body segment from the GRM prediction. Consequently, the accuracy of GRM estimation for the entire gait cycle was improved compared to previous studies.
The proposed method provided estimation accuracies in joint torques comparable to those of the estimation method based on the smooth transition assumption [
13]. In addition, a previous study has reported that the optimization approach provides joint torque prediction accuracy levels similar to those of the smooth transition assumption [
16]. Therefore, while the proposed method does not achieve an estimation accuracy comparable to the artificial neural network approach, which provides superior estimation accuracy compared to other methods, the present study demonstrates a practical accuracy of the proposed approach comparable to that of the smooth transition assumption and the optimization approach.
While the GRF prediction was not significantly affected by white noise up to 10 mm, the GRM prediction showed a deterioration in accuracy at 10 mm of white noise. In the proposed prediction method, noise in the contact-point positions influenced the distribution of the external force computed by the equation of motion, but not for the magnitude of the external force. Therefore, the noise in the contact-point positions had less impact on the GRFs, which are calculated as the sum of the forces acting on the contact points of each foot. In contrast, the distribution of the external force at each contact point had a significant effect on the COP, leading to a deterioration in estimation accuracy for the GRMs. Based on the results of the sensitivity analysis, it is expected that the estimation accuracy of GRMs will decline if the proposed method is applied to IMU-based prediction, which has fewer limitations for measurement locations and lower measurement accuracy than that of the OMC-based approach. Nevertheless, the present method has the potential for application in IMU-based methods with practical accuracy because the estimation accuracy with white noise at 10 mm was not greatly different from the estimation accuracy in some previous studies, as shown in
Table 3 [
13,
16].
Several limitations of this study are noted. Although the present study developed the GRF estimation method to address the location limitations of force plates, the proposed method can only be applied in situations where an OMC system can be installed. While the location limitations of the present method can potentially be relieved by employing IMUs instead of an OMC system, a critical challenge remains in determining the position of the foot contact points when using IMUs, as IMUs cannot directly measure three-dimensional coordinates.
In the present study, the experiment was conducted with participants in bare feet to accurately capture foot deformations. Consequently, this study cannot ensure prediction accuracy under conditions other than bare feet, such as when shoes are worn. To extend the applicability of the method to various motions in the future, the accuracy of this method should be confirmed under conditions other than bare feet.
This study estimated GRFs and GRMs for movement on level ground by assuming the critical position shown in Equation (4) to be constant. While the prediction accuracy needs confirmation, the proposed method may be applicable beyond level ground to situations where the critical position can be formulated, such as when the ground has a constant slope. However, it is challenging to apply this method to a rough walking surface with irregularities in the ground.
The estimation accuracy in this study was confirmed only for healthy young adults walking at three different speeds. Although the proposed method is expected to be applicable to various activities other than walking, further accuracy verification is necessary to apply the proposed method to those activities.