Combining forces and kinematics for calculating consistent centre of mass trajectories

For biomechanical studies on animal and human locomotion as well as for clinical gait assessment, a precise estimation of the body centre of mass (CoM) trajectory is crucial. In general, the CoM motion is calculated either from kinematic data (kinematic method), often by incorporating an anthropometric model, or by double integration of the acceleration obtained from ground reaction forces (GRF; dynamic method).

The kinematic methods can be subdivided into pure marker methods and segmental analysis methods. Pure marker methods use a single marker (Saini et al., 1998; Thirunarayan et al., 1996) or a minimalistic marker set (Forsell and Halvorsen, 2009; Halvorsen et al., 2009) to estimate the CoM coordinates. However, some uncertainty arises in these minimalistic estimates because limb motion is not accounted for (Whittle, 1997).

The segmental analysis methods rely on full-body marker sets, e.g. Helen–Hayes (cf. Castagno et al., 1995; Kadaba et al., 1990; Sutherland, 2002), and calculates the CoM trajectory by assuming the masses and CoM locations of each segment (Eames et al., 1999). For humans, segmental data are usually obtained from anthropometric literature (Dempster, 1955) or individually determined from a reference measurement (Forsell and Halvorsen, 2009). However, for animals biometric data are much less available.

The accuracy of the segmental analyses relies on the correct marker placement and the correct estimation of segment properties (cf. Shan and Bohn, 2003). In addition, potential errors arise from unrecorded dynamics of segment masses, i.e. wobbling masses (Gruber et al., 1998; Günther et al., 2003; Schmitt and Günther, 2010) and motion of the viscera (Minetti and Belli, 1994).

In the dynamic method these problems do not occur. According to Newton's second law, the CoM motion is fully determined by the body's mass, external forces, the initial velocity and the initial position. Therefore, the accuracy of the dynamic method is limited by the precision of the measured forces and the precision of the integration constants (Cavagna, 1975). Although the initial position can be estimated quite well, the determination of the initial velocity is a common problem in practice. Gutierrez-Farewik et al. (Gutierrez-Farewik et al., 2006) and Günther and Blickhan (Günther and Blickhan, 2002) showed that the calculated trajectory strongly depends on the body's mass measurement and the estimated initial velocity.

It is possible to optimize the guess for the initial velocity with the path matching method, i.e. the CoM initial conditions are estimated from a reference marker (Daley et al., 2006; McGowan et al., 2005). Yet, as this maker-based criterion itself is only a guess, the optimisation is prone to error. In addition, because of systematic errors in the force signals (Mack, 2007), high-pass filters have to be applied to the force signals for long-term integration of acceleration. Another method to correct the integration constants on a per-stride basis has been proposed by Saibene and Minetti (Saibene and Minetti, 2003). Here, the mean velocity of each stride is replaced by a kinematic estimate, thus attenuating the long-term drift. However, it remains unclear to what extend typical systematic errors in both kinematic and dynamic measurements bias the results.

The sources of errors in kinematic and dynamic methods are of a different kind. Although the dynamic method has drawbacks on long-term scales, the kinematic methods do not fully capture some effects on short time scales, e.g. wobbling masses. Thus, typical inherent errors of both methods can be regarded as opposed, resulting in systematic mismatches of kinematic and dynamic CoM estimates (Gard et al., 2004). Gard et al. (Gard et al., 2004) further proposed that a combination of kinematic and dynamic methods could potentially lead to more accurate results.

Here, we present a simple yet powerful approach to calculate CoM trajectories from measured kinematic and dynamic data. It is based on the combination of the low-frequency content of kinematic data and the high-frequency content of GRF to calculate a more accurate CoM velocity. From this velocity, the CoM trajectory and a corrected GRF are calculated by integration and derivation, respectively. Figuratively speaking, this method combines the strengths of both approaches by taking the coarse motion from the kinematic data and adding the fine structure from the GRF. The proposed method is tested on simulation data and subsequently applied to real measurements. The resulting CoM trajectory and GRF are physically consistent and closely resemble the original measurements.

Here, we present the new CoM algorithm and the simulation model used for verification. Further, we describe real measurements to which the algorithm is applied.

The main idea of the algorithm is to compute the CoM velocity based on kinematic and dynamic velocity estimates, taking only the frequencies of each estimate into account that we consider to be reliable. From this combined velocity, the GRF and CoM trajectory are calculated. The algorithm is summarized in Fig. 1. The individual steps are as follows and hold for each coordinate separately:

• We computed two estimates of the CoM velocity, a kinematic estimate v_k, by differentiating the kinematic CoM estimate, and a dynamic estimate v_d, by integrating the GRF (taking gravity and body mass into account). Here, it is important that the mean force is accurate. For a subject standing at the beginning and the end of a trial, we set the mean force to exactly zero, because here the mean GRF exactly compensates gravity.

• We computed the Fourier transform v_k and v_d of the kinematic and dynamic velocity estimates, respectively.

• We selected a weighting factor w between 0 and 1 as a function of the frequency, i.e. w(f) expressing the confidence in each method with respect to the frequency. The choice of w can be dependent on the experiment and the equipment. Here, we used a sigmoid function that is close to 0 for low frequencies and close to 1 for high frequencies (Fig. 2). Because the Fourier spectrum is symmetric, this function must also be symmetric. For the appropriate selection of the threshold separating high and low frequencies, see Selection of the weighting function and Fig. 3.

• We used the weighting factor w to create a combined spectrum , which is a weighted sum of the low frequencies of the kinematic estimate and the high frequencies of the dynamic estimate.

• We computed the combined velocity v_c as the inverse Fourier transformation of the combined spectrum v_c.

• We computed the combined CoM position and GRF from the combined velocity by integrating or differentiating, respectively.

Some parts of the segment dynamics, primarily the composition of soft and rigid body structures, are not well recorded in the kinematics. As this segment-internal motion, i.e. motion of the soft tissue with respect to the bone, mainly affects harmonics of the stepping frequency and higher frequencies (Günther et al., 2003), we expect the part of motion concerning frequencies below a certain threshold to be well captured, and that this threshold is not substantially lower than the stepping frequency.

In contrast, errors in force measurement usually are of a type that mainly affects low frequencies, e.g. drift and slowly varying offset (Mack, 2007; Nigg and Herzog, 1999). This is why integration of acceleration obtained from measured GRF without filtering gives reasonable results only if trials are short, i.e. if very low frequencies are not present. However, we have no reason to assume systematic errors in the force measurement at high frequencies. An exception may be the eigenfrequency of the measurement system. In our system, the eigenfrequency is 120 Hz, which corresponds roughly to the 40th harmonic of the stepping frequency and thus is far out of the region of interest (Racic et al., 2010).

These arguments lead us to use a weighting function that takes low frequencies from the kinematics and higher frequencies from the dynamics. As we expect that the wobbling masses move approximately with the frequency of the driving force (i.e. the stepping frequency) and their harmonics, we assume that a threshold frequency slightly below the dominant frequency of the motion would be optimal. The validity of this argument is strongly supported by our simulation results, as shown in Fig. 3.

In order to verify the accuracy of the CoM estimate, we used a simulation model of a bipedal walker to create an artificial dataset with precisely known dynamics. Subsequently, we added virtual measurement errors that resemble errors we expect from real-world measurements, including errors concerning the wobbling masses.

We used the model of a human walker (mass=80 kg) from Geyer and Herr (Geyer and Herr, 2010), which is able to predict typical human-like GRF, CoM motion and even muscle activation patterns. We modified this model by setting ¼ of each segment's mass to a wobbling mass, which is connected to the segment by a spring-damper element [roughly according to Minetti and Belli (Minetti and Belli, 1994)]. The spring and damping constants are chosen such that the eigenfrequency is ∼10 Hz for the limbs and ∼3 Hz for the head–arms–trunk segment, with a decay time of 0.4 s. Our modified model walked in an aperiodic manner with a dominant motion frequency of 1.64 Hz. Data were sampled every 1 ms simulation time.

The kinematic CoM estimate of the model was computed using standard segmental analysis (Winter, 2009), including only measurements of the position of the rigid segments, and neglecting wobbling masses. Additionally, white noise with a r.m.s. of 0.5 mm was added.

The force data were modified by: (1) adding a non-stationary noise signal with amplitude of 5 N; (2) adding white noise with a r.m.s. of 5 N; and (3) applying a small nonlinear scaling. In detail, these corruptions were calculated to: (1) , (2) R_i =Aη′_i and (3) F_i′=0.98(F⁰_i)^1.02, where R¹_i and R²_i are random non-stationary and stationary noise, respectively; A is amplitude (5 N); η_i and η′_i are gaussian random numbers with zero mean and unit variance; F_i′ is the scaled GRF; and F_i⁰ is the real GRF in units of body weight. The index i denotes the particular sampling frame. Then, the virtually measured force is F_i^m=F_i′ + R_i¹ + R_i².

The dynamic estimate of the CoM was computed by a double integration of the acceleration obtained from GRF, thereby applying a 0.35 Hz high-pass filter for the force and obtained velocity. The cut-off at 0.35 Hz was chosen because in this model it resulted in the most accurate CoM reconstruction. Lower cut-offs led to increased long-term oscillation whereas higher cut-offs led to an underestimation of the oscillation amplitude.

For the CoM estimation using the proposed algorithm, the threshold frequency was set to 1.5 Hz, which is below the dominant frequency of the motion at 1.64 Hz.

In order to demonstrate the new method, we calculated the CoM trajectories of human running and walking and compared these with kinematic and dynamic CoM estimates. We further calculated the ‘kinematic GRF’ as second derivative of the kinematic CoM estimate to compare kinematic GRF, measured GRF and calculated GRF.

All measurements were conducted on an instrumented custom-built treadmill (ADAL, HEF Medical Development, Andrezieux-Boutheon, France). We further used a marker-based kinematic system (Qualisys, Gothenburg, Sweden) to capture the subject's motion. Force data were sampled at 1000 Hz; kinematics were sampled at 240 Hz.

The subject walked at 1.8 m s^–1 and ran at 2.7 m s^–1 on the treadmill. Each trial started and ended with 5 s of quiet standing. The total trial duration was 50 and 55 s, respectively.

A linear drift in the total vertical and total horizontal force was removed from the raw data. Kinematic data were linearly interpolated from 240 to 1000 Hz to match the treadmill sampling frequency. We do not expect relevant numerical errors from the interpolation because we used only kinematic frequencies up to ∼1.5 Hz, which are well over-sampled and thus hardly affected by this interpolation.

The kinematic CoM estimate was obtained using a standard segmental analysis (Winter, 2009). Corresponding GRFs were calculated by computing the second derivative of the CoM estimate after applying a second-order Butterworth filter with a cut-off frequency of 15 Hz.

The dynamic CoM estimate was obtained by twice integrating the GRF. To avoid spurious drifts in the resulting trajectory, components below 0.35 Hz were removed before integration by applying a first-order Butterworth filter to GRFs and computed velocity.

Finally, we applied the proposed algorithm to calculate the combined CoM and GRF. We used the measured GRF and the kinematic CoM estimate as inputs for the algorithm. The threshold frequency was set to 1.5 Hz.

The analysis of the simulated model shows that the proposed method reconstructs the CoM motion with greater accuracy than the pure kinematic or dynamic methods (Fig. 4). For better visual comparison, only the displacement relative to the first apex position is shown in Fig. 4A. Also, we applied the new method to the forward motion. The results show a proper tracking of the motion (Fig. 4C). The differences of the kinematic estimate are apparent in every step, whereas, because of the tracking problem, the dynamic method accumulates the error in the long term (Fig. 4B). These slow drifts are not present in the combined estimate.

The results of applying the algorithm to walking and running data are shown in Fig. 5. Here, we focused on the vertical CoM component, but the method is applicable to any CoM component. The calculated data resemble both the kinematic CoM estimate and the measured GRF. Because of this and its inherent physical consistency (i.e. the combined GRF equals the second derivative of the combined CoM minus gravity), this method provides a suitable enhancement to common kinematic CoM estimations. These consistent data can then be used for further analysis, e.g. of external work and external power, with greater confidence.

Our results show systematic deviations of the GRF obtained directly from kinematics compared with the measured GRF. This comprises mainly an underestimation of the GRF after lift-off in running, which is accompanied by an overestimation of the GRF during stance (Fig. 5). Similar results were shown by Racic et al. (Racic et al., 2010). As negative vertical GRF cannot occur in typical running and hopping experiments, this indicates a systematic shortcoming of kinematic GRF estimates, which also extends to the corresponding CoM estimate according to Newton's second law. The combined CoM/GRF estimates do not show this systematic deviation but closely resemble the measured forces. In the Appendix, we show that this property also holds for very simple kinematic CoM estimations, such as a single marker, both in human walking and dog trotting. Thus, accurate CoM trajectories can be obtained without knowing segment properties using this algorithm. Further, this also allows highly reduced experimental effort for some gait analyses.

The proposed method also accounts for wobbling masses (Gruber et al., 1998; Günther et al., 2003; Schmitt and Günther, 2010) because their motion is included in the GRF. Wobbling masses are not necessarily captured by the kinematic measurement. This can lead to differences especially when the wobbling mass motion is substantial or out of phase with respect to the skeletal motion (Minetti and Belli, 1994). Conversely, it appears plausible that this soft tissue motion could be estimated by calculating the difference of the calculated CoM and the kinematic CoM estimate. Thus, this method could also provide a basis for analyzing soft tissue motion.

A standard method in engineering to combine both kinematic and dynamic input to obtain a more reliable CoM estimate is the Kalman filter (Kalman, 1960). Roughly speaking, the main idea is to propagate the trajectory based on a system model using the GRF as input [so-called estimate x_i^–=f(x_i–1⁺, GRF)], and updating this estimate with each measurement y_i of the CoM trajectory, x_i⁺=x^–+K(y_i–x_i^–). The relative weight K of the update is based on the relative confidence in the estimate and the measurement. Under certain restrictions of the measurement errors, the Kalman filter or its modifications give an optimal estimation. However, this does not apply here as both the kinematic and GRF measurement errors have a systematic structure that renders them very different from these restrictions. In order to use Kalman filtering here, a model of the measurement errors would have to be included. In contrast, the proposed method offers a convenient and intuitive way, namely the selection of a threshold frequency, to account for the typical structure in the measurement errors.

Physical consistency in long trials can also be obtained by the dynamic method, when low frequencies are discarded. When low frequencies are not discarded, spurious slow oscillations with high amplitude occur. However, when low-frequency components are discarded, slow changes such as the descent of CoM in walking compared with standing cannot be captured (Fig. 6). Our approach solves this problem by adequately replacing the low-frequency components with those obtained from the kinematics. This is similar to cutting the long trial into short trials and taking the initial conditions for each short trial from kinematics [see also the method of Saibene and Minetti (Saibene and Minetti, 2003)]. In our case, the length of each short trial would then be 1 divided by the threshold frequency. However, in our method the kinematic estimate not only influences the initial value of each short trial, but also reshapes the result equally over time. This demonstrates the wider applicability of this method to accelerated (Segers et al., 2007; Williams et al., 2009) and unsteady (McGowan et al., 2005) conditions. Further, the need for repeated estimates of the CoM velocity from kinematics is omitted. Thus, the proposed method can be regarded as an extension of the dynamic method that enables its application especially for long trials, when low-frequency components become important, or non-periodic movements, where estimation of integration constants is difficult.

When multi-segment kinematics are not available or the segment properties are unknown, the proposed algorithm still can be applied to improve the resulting CoM estimate. Here, we demonstrate the result when we take the sacral marker as the kinematic CoM estimate and compare it with the CoM and GRF obtained when using the full maker set. Further, we show the results of computing the CoM and GRF using a single marker on the back of a trotting dog (Fig. A1). For human data (Fig. A1, left), a substantial overestimation of the CoM amplitude by the sacral marker estimation is visible, which is also reflected in substantial differences in the GRF. For animal data (Fig. A1, right), differences between the single-marker estimate and the combined estimate are apparent. However, in both situations, taking the single marker as input for the algorithm results in reliable CoM and GRF estimates, which is reflected in both cases by close resemblance to the measured GRF.

However, it is clear that especially slow or permanent, unrecorded motions (e.g. a change of posture when using a single marker) induce an offset in the simplified kinematic CoM estimate that transfers to a corresponding offset of the combined estimate. This has to be taken into account when designing the experiment.

 
• CoM

  centre of mass

 
• GRF

  ground reaction force

 
• r.m.s.

  root mean square

 
• v_c

  CoM velocity calculated from kinematics and GRF

 
• v_d

  CoM velocity, estimated from GRF

 
• v_k

  CoM velocity, estimated from kinematic measurement

 
• 

  Fourier spectrum of v_k, v_d and v_c, respectively

 
• w(f)

  weighting parameter as a function of frequency f, with values in [0, 1]

The authors thank Margrit Schaarschmidt for assisting in the data acquisition and Hartmut Geyer for providing his walking model (which can be downloaded from http://www.cs.cmu.edu/~hgeyer/Research_MotorControl.html). The authors further thank the anonymous reviewers whose comments led to great improvements in this paper.