Decomposition of three-dimensional ground-reaction forces under both feet during gait.

Three-dimensional ground reaction forces (3D-GRF) are essential for functional evaluation for rehabilitation. A platform path is required to obtain the 3D-GRF. The main shortcoming of these platform paths is that during double stance phases of gait, both feet can be placed on the same force platform causing the need for decomposing the 3D-GRF under each foot. Despite the high number of studies on force decomposition, there is still no method on the decomposition of 3D-GRF based on data from platforms.
OBJECTIVE: This study aims to present an automatic method using parametric curve fitting modeling to increase the accuracy of decomposition of 3D-GRF during double stances under each foot.
METHODS: The decomposition method was applied to the global 3D-GRF using 3rd order polynomial, sine, and sine-sigmoid functions. The computed 3D-GRF was compared to the 3D-GRF independently recorded by force platforms for each subject.
RESULTS: The relative average error between the computed 3D-GRF and the recorded 3D-GRF were equal to 3.3±1.6%. In details for the vertical, antero-posterior, and medio-lateral GRF, these errors were 2.9±1.6%, 6.3±4.3%, and, 9.5±3.6%, respectively, for 30 subjects.
CONCLUSION: The global error on the GRF is the best one in the literature. This method can be validated on various populations with musculoskeletal disorders.

Introduction

Evaluating Ground Reaction Forces (GRF) under each foot is crucial in different domains ranging from biomechanics, to compute centers of pressure and spatio-temporal parameters[1-3], health monitoring[4-6] and structural dynamics, to model the interaction between the foot and the structure[1,7].

Three-Dimensional Ground Reaction Forces (3D-GRF) are traditionally recorded using separate platforms under each foot. This generally leads to a large number of trials to get cycles where both steps are well inside the force platforms separately. These repeated trials can be both time consuming and difficult to achieve[8]. Moreover, it may also lead to undesired unconscious targeted steps and muscular fatigue[9-13].

Consequently using unique large force platforms emerges as a practical method to facilitate and accelerate the procedure of recording GRF data for a varied population[1,14].However, a main challenge remains in how to decompose the global 3D-GRF under each foot along X, Y and Z axes corresponding to medio-lateral (ML), antero-posterior (AP) and vertical (V) GRF, respectively, since the unique force platform only provides the resultant of GRF for the two feet combined.

Accurately decomposing the GRF has been a challenge for more than two decades, with a dozen proposed methods described the significant ones[9,10,15-20] here below. First Davis and Cavanagh (1993), proposed a method to detect the transition between single stance phase (SS) and double stance phase (DS) and to decompose V GRF profiles[9]. To detect the transition between phases, their proposed method was based on the analysis of the side-to-side oscillations of the global center of pressure (COP). In numerous cases, the path of the COP of ML GRF does not systematically show a clear side-to-side inflection point[10]. In spite of that, this method is the most referred one for the decomposition of V GRF[9,10,14,17,20].

Ballaz et al. (2013) were the first authors to propose an automatic method to detect SS and DS transition based on the difference of COP in the horizontal plane between two given time instants. They also proposed an automatic method to decompose V GRF profiles modeling the first foot leaving the ground by using a spline interpolation-based method. Their study achieved a rather low average error of 3.8% and was tested on 6 typically developed children[10,21]. To our knowledge, this innovative method, is the best and the more accurate method in the literature. It is therefore natural to be used as a measure for assessing the performance of new methods.

It is however, worth noting that their work still presents some limitations. First off, spline interpolation does not consider that the typical GRF patterns are not always smooth. This omission can potentially lead to inaccuracies while decomposing AP and ML GRF[17]. The V GRF under each foot are vertically oriented (opposed to gravity), whereas GRF along X and Y axes change from negative to positive directions during each cycle of gait. Consequently, GRF along these two axes show more complex patterns in comparison to the V GRF. Therefore, a simple spline interpolation cannot accurately model the remaining two components.

Recently, Oh et al. (2013), Karatisdis et al. (2016) and Villeger et al. (2014) decomposed 3D-GRF by applying a genetic algorithm-general regression neural network[16], smooth transition assumption[22] and multiple regression parameter models[15], respectively. The major practical limitation to compute 3D-GRF using the two first methods is the need for motion capture systems to build the learning database. Additionally, Hijazi and Makssoud (2015) decomposed only V GRF using a hyperbolic tangent function[19], but the robustness of the method is not known since they did not validate it by comparing the computed V-GRF to the recorded one, besides that their model was very sensitive to the choice of initial values of the parameters.

In fact, 3D-GRF decomposition has been a challenge for 20 years with different methods. However, most of these methods have been only applied to vertical forces. In other words, it has never been evaluated whether these methods on Z are transferable to the X and Y directions. The hypothesis at the beginning of this study is that these methods are not a priori easily transferable to the X and Y components since their behaviors and profiles are very different from V GRF. Particularly, while the V GRF under each foot are always oriented upwards, the forces in X and Y go from positive to negative during each walking cycle (Figure 1). Moreover, their profile is more variable and of lower magnitude than V-GRF[23]. This makes the problem more complex. Consequently, as of today, no method has been proposed for the decomposition of 3D-GRF solely using force platform measurements. In fact, for functional evaluation in rehabilitation, all three GRF are absolutely essential. For instance, defining ML and AP GRF is required to study the balance of stance phases during gait[24] and distinguish functional compensation from physiological restitution, respectively[25].

Figure 1

Medio-lateral (ML), Antero-posterior (AP) and Vertical (V) GRF during DS of an individual. Medio-lateral (ML) GRF (black), Antero-posterior (AP) GRF (red), Vertical (V) GRF (blue). SS: The single stance phase, DS: Double stance phase.

Furthermore, no method has ever been tested on a large number of subjects to enable large group statistics analysis (30 subjects in this study). Therefore, the objective of this study is to develop an automatic method that decomposes global 3D-GRF under each foot during DS. This method should be usable on any platform type, whether it is long force platforms or juxtaposed platforms, and without the need of a motion capture system.

Methods

Selecting participants

Thirty healthy adults volunteered for this study with an average age of 24.8 (standard deviation (STD): 3.09 year), range: 20-34; 17 males and 13 females.

Set-up and procedure for data collection

The process of collecting participants data was carried out in a motion laboratory, equipped with two separate force platforms 50×50 cm2) (AMTI, USA) embedded in the floor to record GRF at a frequency of 1000 Hz (Figure 2a & b).

Figure 2

a. Focus on the human body walking on the force platforms (VICON Nexus). V: Ground reaction forces (GRF) along with the vertical axis, ML: GRF along with the X axis, AP: GRF along with the Y axis, 3D-GRF: Global GRF, which is the norm of the V, AP and ML components of GRF. b. An individual during the evaluation test in the motion laboratory.

The protocol of the study was approved by the ethic committee of Sainte-Justine University Hospital Center (UHC). The experimental procedure was explained for each participant and all participants provided a written informed consent prior to the test.

To record left and right GRF, they were independently measured on two platforms during untargeted gait. The participants were asked to walk naturally, barefoot, at their own preferred speed, looking straight ahead without targeting their steps on the platforms. The trials with both steps inside the force platforms were recorded and used for this study.

3D-GRF decomposition: identification of functions and parameters

To decompose 3D-GRF, first, it is essential to identify the transition instants between SS and DS, since decomposition is only needed during DS when both feet are in contact with the force platforms. Therefore, we used the automatic method, developed in our research group and proposed by Ballaz et al.[10] (Step 1). At the next step (Step 2), the 3D-GRF of the foot that leaves the ground was modeled, using a parametric curve modeling. Then the 3D-GRF of the second foot was obtained by subtracting the modeled first foot from the total recorded 3D-GRF (Step 3).

Due to the vibration of force platforms, during the evaluation, the recorded 3D-GRF contained noise signals. First, a 4th order low-pass Butterworth filter with a 4 Hz cut-off frequency was applied to the recorded 3D-GRF signals to remove the noise, then Ballaz et al. phase detection method and our decomposition proposed method were applied on the 3D-GRF data.

Step 1: Identification of transition instants between SS and DS

This step aimed to detect the transition between the SS and DS phases, using Ballaz et al. method[10]. This method has the advantage of providing the exact instant of transition between SS and DS. Indeed these instants “are determined by the peaks of a variable called ΔCOP[2] (the squared Euclidian norm in the horizontal plane between two given time instants) that is very sensitive to subject COP displacement during gait”[10]. The ΔCOP[2] is compared with the constant thresholds to detect the phase transition.

Step 2: Computation of GRF under the taking-off foot during DS

This step aimed to compute the GRF under the first taking-off foot during the DS, by using parametric curve fitting model. First the general geometric pattern of the recorded 3D-GRF by force platforms was observed, as shown in Figure 1 with the black line: generally, the 3D-GRF curves present S-type curve patterns[24]. Then we decided to find the functions that can best fit the 3D-GRF curves. Therefore, the two main criteria used to select the best functions were their ability to form an ‘S’ type curve, with a manageably low number of unknown parameters. This number was here limited to 4 since it is related to the 4 main parameters, which are amplitude, codomain, vertical and horizontal shift of 3D-GRF curves during DS. Having too many parameters in functions with a limited amount of information would lead to inaccuracy while calculating them. A preliminary study was performed to examine and compare different options of possible mathematical functions to model GRF along X, Y and Z axes during DS. In this preliminary study, parameters of functions were obtained by curve fitting method during DS. The functions with less error between the fitted curve and measured 3D-GRF by force platforms were selected for this study. Table 1 and Figure A1 in the Appendix shows that the models utilizing sine, polynomial function of degree 3, and arc-tangent achieved a better fit when compared to the functions of Hill[26] and Strejc[27]. The arctangent function was replaced by the sigmoid function since it has a similar S curve pattern; however solving sigmoid functions is less complex. Therefore, the sine, sigmoid, and 3rd degree polynomial functions were chosen to model 3D GRF of the first landed foot during DS. In the following subsections, we present our method for identifying the parameters of each function, using the kinetics data during SS, measured by force platforms. The corresponding results are presented in the Results section.

Table 1

Curve fitting results (comparing functions).

Figure A1

a. Antero-posterior (AP) GRF during DS. b. Medio-lateral (ML) GRF during DS. Recorded trial by force platforms (black), curve fitting using 3rd order polynomial function (red), curve fitting using sine function (green), curve fitting using arctangent function (blue), curve fitting using Strejc function (magenta), curve fitting using Hill function (cyan).

Modeling 3D-GRF of the first landed foot using 3rd order polynomial function

To compute the 3D-GRF during DS using 3rd order polynomial function, the following equation was used.

- A (1×1) is the direction of the polynomial function, defines if it goes from negative to positive or opposite (Obtained from solving Equation 4 & 5 with the values of C and D).

- B (1×1) is the parabola shape of the polynomial degree 3 function (Obtained from solving Equation 4 & 5 with the values of C and D);

- C (1×1) is the slope of the polynomial degree 3 function (Obtained from Equation 3);

- D (1×1) is the vertical shift of the polynomial degree 3 function (Obtained from Equation 2).

The equations 2-5 are used to compute the parameters of the 3rd order polynomial function, as explained above:

where:

- fpolynomial is the function of 3rd order polynomial;

- f’polynomial is the derivative function of 3rd order polynomial.

All variables are defined in Table 2 and shown in Figure 3-5.

Table 2

Variables used in the functions.

Figure 5

Sigmoid and sine-sigmoid functions parameters to compute 3D-GRF.

Sine function parameters to compute the 3D-GRF.

Modeling 3D-GRF of the first landed foot using sine function

To compute the 3D-GRF during DS using sine function, the following equation was used:

And the required parameters were calculated as bellow (see Figure 3):

Figure 3

Polynomial degree 3 function parameters to compute 3D-GRF.

Where,

- A (1×1) is the amplitude of the sine function, which is proportional (a) to the value of LFSS, globally representing the amplitude of the actual DS curve; this was obtained from the experimental data (see Table A1 in Appendix);

Table A1

Variables used in the functions.

- 1/B (1×1) is the period of the sine function, which is proportional (b) to the value of t, globally representing the period of the actual DS curve. This was obtained from the experimental data (see Table A1 in Appendix);

- C (1×1) is the horizontal phase shift of the sine function, defining the starting point of the sine function, i.e. globally representing the starting point of DS for each 3D-GRF in horizontal axis. This was obtained from an equation based on the experimental data (see equation A1 and Table A1 in Appendix);

- D (1×1) is the vertical phase shift of sine function, defines where the computed curve should start vertically in DS. D is proportional (d) to the value of LFSS; This was obtained from the experimental data (see Table A1 in Appendix).

Modeling 3D-GRF of the first landed foot using sigmoid function

To compute the 3D-GRF during DS using a sigmoid function, the following equation was used (see Figure 5, Table 2 and Table A1 in Appendix):

where,

- A (N)(1×1) is the codomain of the sigmoid function, i.e. the vertical range of the sigmoid curves, which corresponds to the vertical range of 3D-GRF curves during DS;

- B (1×1) is a constant value obtained practically on the basis of the geometric curve pattern from the experimental data (see Table A1 in Appendix). B defines the slope of the sigmoid curve, which represents the slope of the GRF curve during DS: an increase in B makes the curve steeper;

- C (1×1) is the horizontal phase shift of the sigmoid function, which is proportional to the value of t. C was obtained from the experimental data (see equation A2 and Table A1 in Appendix), corresponding the point where DS phase starts.

Modeling 3D-GRF of the first landed foot using sine-sigmoid function

Due to the pattern of the 3D-GRF during DS, we added the sine function to the first part of the sigmoid function to create a model more representative of the 3D-GRF. Typical sigmoid functions start with a horizontal asymptote, however, the GRF starts with an arch pattern (Figure 6a-c). Therefore the first part (t0 to Q) was a combination of sine and sigmoid functions and the rest (Q to t) was fitted by only the sigmoid function. Equations 9-11 were used to compute the parameters. All of the equations were defined based on the kinetic parameters during the singles stance phases before and after double stance provided by force platforms (see Figure 5, Table 2 and Table A1 in Appendix):

Figure 6

a. Antero-posterior (AP) GRF during DS, b.Medio-lateral (ML) GRF during DS, c.Vertical (V) GRF during DS. Recorded trial by force platforms (black), spline interpolation (magenta), parametric curve modeling using 3rd order polynomial function (red), parametric curve modeling using sine function (green), parametric curve modeling using sigmoid function (orange), parametric curve modeling using sine-sigmoid function (blue). SS: The single stance phase, DS: Double stance phase, LFSS: Value of GRF in the last SS instant, t0: Last instant of first single stance (SS) before double stance (DS), tL: The first instant of second SS after DS where the first foot leaves the floor.

where

- S1 (t) is a sigmoid function;

- S2 (t) is a sine function;

Step 3: Computation of GRF under the landing foot during DS

The 3D-GRF of the second foot that lands obtained by subtracting the 3D-GRF of the first foot from the total 3D-GRF and is obtained as follows:

where,

- F [N] is a vector (i×1) of total recorded 3D-GRF by force platforms;

- F [N] is a vector (i×1) of 3D-GRF of second foot that lands;

- F [N] is a vector (i×1) of computed 3D-GRF of the first foot that lands.

It is worth noting that the proposed 3D-GRF decomposition model (step 2 and 3) works in real-time using MATLAB® (MATHWORKS®, USA), with a GRF decomposition time of 0.4, 0.8, and 0.9 ms for V, AP, and ML during one double stance phase. This enables us to implement the developed solution on force platform paths.

Validation of the proposed method

The validation procedure consisted in the analysis of error between the recorded 3D-GRF by force platforms and modeled 3D-GRF. During the DS, the mean absolute error e (1×1) between computed and recorded GRF was defined as

where

- F [N] is a vector (i×1) of GRF computed by the proposed method;

- F [N] is a vector (i×1) of individually measured GRF on each platform during DS.

The best function fitted for each axis was chosen based on the minimal value of mean relative error (E% (1×1)), defined as follow:

where,

- f (N) (1×1) is maximal measured 3D-GRF during DS.

Complementary, the mean global of absolute relative error between measured and computed 3D-GRF were calculated, as bellow:

where,

- e (N) (1×1) is mean of global absolute error;

- f [N] is a vector (i×1) of independent recorded GRF along X, Y and Z axis respectively by force platforms;

- f [N] is a vector (i×1) of computed GRF along X, Y and Z axis respectively.

- E (%) (1×1) is the mean of global relative error;

- f (N) (1×1) is the maximal 3D-GRF during DS phase.

Results

Figures 6a-c illustrate the results of computed 3D-GRF for one subject, as a sample, using the 4 parametric curve investigated in this paper drawn against the Spline interpolation proposed by Ballaz et al. (2013) as well as the recorded 3D-GRF by the force platforms.

Table 3 presents the mean relative errors and the mean absolute errors normalized by body mass between the computed and measured GRF. The mean relative error obtained by using the different functions does not vary greatly, which confirms that these functions were indeed selected properly.

Table 3

Statistical results of the 30 participants. E, e, e and E, respectively, mean relative error between computed and recorded GRF, mean absolute error between computed and recorded GRF normalized to body mass (BM), maximal and mean relative global error . DS: Double stance; ML: Medio-lateral GRF (Fx); AP: Antero-posterior GRF (Fy); V: Vertical GRF (Fz).

The minimal mean relative error calculated for each force component is shown in underlined bold (Table 3). A low average relative error of 2.9% (STD: 1.6%) for the vertical GRF using the sine function. The sine-sigmoid function yielded the lowest error values of 6.3% (STD: 4.3%) and 9.5% (STD: 3.6%) for antero-posterior and medio-lateral GRF respectively. A global 3D-GRF error equal to 3.3% (STD: 1.6%) was obtained by the best function for each axis on global GRF.

As presented in Table 3, the 3rd order polynomial function also provides a mean relative error of 2.9% (STD: ±2.0) for modeling the V GRF, but the sine function shows a 20% smaller standard deviation among 30 subjects (1.6 vs. 2.0). Therefore, we conclude that the sine function is the more app ropriate one to decompose V GRF.

Discussion

The current study proposed a new method to improve the accuracy of decomposition of GRF during double stance phase (DS) of gait using different parametric curves and validated on 30 healthy adults. Comparing our results with the ones reported in the literature, there is a decrease up to 47%[15] and 38%[16] of mean relative errors in decomposing ML and V GRF respectively. The maximum GRF difference normalized to body mass is lower than 1 N/kg (0.1 N/kg, 0.1 N/kg, 0.9 N/kg, for AP, ML and V respectively) which is comparable to the results presented in Villeger et al.[15].

The results of one subject were chosen as a sample to be presented in Figure 5a-c. Observing the behavior of 3D-GRF during DS (Figure 6a-c), shows that the vertical GRF, are less sensitive to the choice of modeling function. It can explain the relative success of different studies in decomposing V GRF solely[9,10,19] in comparison to studies decomposing AP and ML GRF[15,16].

The mean relative error along the Z axis is similar using 3rd order polynomial and sine functions, however, the maximum value of absolute error normalized to body mass and standard deviation are higher for the 3rd order polynomial. The mean relative global error (E) of 3D-GRF using the best fitted functions i.e. sine-sigmoid function for AP and ML GRF and sine function for V-GRF, is 3.3% (STD: 1.6%), which is, to our knowledge, the best in the literature.

Conclusion

The objective of this study was to propose for the first time an automatic method to decompose 3D-GRF along X, Y and Z axes during gait DS. The main results show a mean relative global error on the GRF equal to 3.3±1.6%, which was obtained using the sine-sigmoid function for AP and ML GRF, and the sine function for V GRF. This error is the best one in the literature. Further, the 3D-GRF decomposition process works in real-time, enabling to be implemented on different types of force platforms without the need for any motion capture system. Consequently, the method is recommended to be integrated into devices such as sensory floors[16,28] and treadmills[29] or conventional juxtaposed force platforms to obtain the 3D-GRF on several gait cycles.

The perspectives are to validate the method for various populations with musculo-skeletal disorders, extend it to various movement analyses such as gait with a cane or a walker, and compute the corresponding spatio-temporal parameters without the use of any motion capture system.