1.1 Participants
292 young male participants were recruited (age 33.76 ± 4.59 years, height 176.21 ± 5.33cm, weight 68.44 ± 5.23kg). Participants were required to be physically healthy, without lower limb surgery in the past 3 years, and have normal foot types (arch index between 0.21-0.2613). Prior to testing, participants were informed of the study details and signed an informed consent form. This study was approved by the Ethics Committee of Hebei Normal University (No. 2022LLSC026), and all procedures of this experiment were in accordance with the ethical standards laid out in the 1964 Declaration of Helsinki.
1.2 Data Collection
After familiarizing with the protocol and warming up, participants had their foot pressure data collected during walking using a high-frequency foot pressure plate (RSscan International, Belgium, sampling frequency 126Hz, minimum resolution 0.25N, measuring range 1-60N/cm2, plate length 2m with 1.5m extended runways at start and end). Three trials of foot pressure data were recorded for each participant at their self-selected walking speed. Data from trials with walking speeds between 1.1-1.3m/s and containing complete bilateral footprints were included for processing.
1.3 Data Processing
1.3.1 Data Preprocessing
Based on the collected foot pressure data, the foot regions were adjusted in the instrument's built-in foot pressure analysis software (see Fig. 1). The peak pressures and impulses of different regions were recorded for both feet of each participant. The foot with higher total impulse across 10 regions during the stance phase was identified as the primary force-producing foot. A total of 20 parameters including peak pressures and impulses from the 10 regions of the primary force-producing foot during stance were compiled for all participants. Z-scores were calculated to exclude outliers, with any foot pressure parameter having a Z-score greater than 3 considered an outlier14,15. Data without outliers were retained for further analysis, resulting in a final sample size of 243 participants.
1.3.2 Principal Component Extraction
The preprocessed 20 parameters underwent the Kaiser-Meyer-Olkin (KMO) and Bartlett's tests of sphericity to assess correlations and multicollinearity among the parameters. Principal components with eigenvalues greater than 1 were extracted for dimensionality reduction. To enhance interpretability, the principal components were rotated using the maximum variance method. The scores of the extracted principal components were retained.
1.3.3 Optimized Clustering Algorithm
The K-means + + clustering algorithm was applied to the scores of the extracted principal components. The number of clusters ranged from 2 to n (where n is the number of extracted principal components). Compared to traditional K-means, K-means + + optimizes the selection of initial cluster centers, thereby improving clustering quality16. In this study, the maximum number of iterations was set to 10,000. The silhouette coefficient was introduced to optimize the iteration quality (Eq. 1). The silhouette coefficient ranges from − 1 to 1, where values closer to 1 indicate better clustering and negative values suggest misclassification17,18. The number of misclassified samples (silhouette coefficient < 0) and the average silhouette coefficient (Eq. 2) in each cluster were calculated to optimize each iteration.
The Elbow Method is a technique to determine the optimal number of clusters by finding the "elbow point" of the sum of squared errors (SSE) curve. SSE represents the sum of squared distances between each sample and its nearest cluster centroid (Eq. 3). As the number of clusters increases, SSE gradually decreases, but after the "elbow point," the reduction in SSE becomes insignificant with additional clusters. This elbow point indicates the optimal number of clusters19. In this study, the optimal number of clusters was determined by combining the Elbow Method and silhouette coefficients. The process of the optimized clustering algorithm is shown in Fig. 2.
$${s}_{i}=\frac{({b}_{i}-{a}_{i})}{\text{m}\text{a}\text{x}({a}_{i},{b}_{i})}$$
1
(\({s}_{i}\) is the silhouette coefficient of sample i, \({b}_{i}\) is the distance from i to its cluster centroid, \({a}_{i}\) is the distance from i to the nearest neighboring cluster centroid)
$${s}_{k}=\sum _{i=1}^{n}{s}_{i}$$
2
(Where n is the number of samples, k is the number of clusters. \({s}_{k}\) is the average silhouette coefficient for that number of clusters)
$$SSE=\sum _{k=1}^{n}\sum _{i\in {C}_{k}}{\left|i-{m}_{k}\right|}^{2}$$
3
(Where n is the number of clusters, \({C}_{k}\) is the k-th cluster, i is a sample point in cluster \({C}_{k}\), \({m}_{k}\) is the centroid of the kth cluster)
1.3.4 Feature Index Selection
In this study, feature indexes are selected by combining correlation coefficients and the Fisher criterion20,21. For plantar pressure indices, Pearson's correlation analysis is used to categorize the similarity of indicators, as shown in Eq. (4). If r ≥ 0.60, it was considered that there was redundancy between the indicators, and the highly correlated indices are grouped accordingly.
$$r=\frac{\sum _{i=1}^{n}{(x}_{i}-\stackrel{-}{x}\left)\right({y}_{i}-\stackrel{-}{y})}{\sqrt{\sum _{i=1}^{n}{{(x}_{i}-\stackrel{-}{x})}^{2}}\sqrt{\sum _{i=1}^{n}{{(y}_{i}-\stackrel{-}{y})}^{2}}}$$
4
(where \(\stackrel{-}{x}\) and \(\stackrel{-}{y}\) represent the mean values of the samples within the two feature sets, and the value of r ranges from [-1, 1]).
For indicators with high correlation, the Fisher criterion is used for feature selection in conjunction with the classification results mentioned above. The Fisher criterion extracts the optimal feature indicators for classification by calculating the ratio of between-class variance to within-class variance, achieving the best class discrimination20,22. The calculation method for the Fisher criterion is as follows:
If n samples belong to c categories, the calculation of the between-class distance \({f}_{a}\) is as shown in Eq. (5). The calculation of the between-class distance \({f}_{b}\) for different indicators is as shown in Eq. (6). Finally, the Fisher value \(F\) for each indicator is obtained, as shown in Eq. (7).
$${f}_{a}=\frac{1}{n}\sum _{i=1}^{c}{\sum }_{x\in {n}_{i}}{({x}_{i}^{\left(k\right)}-{m}_{i}^{\left(k\right)})}^{2}$$
5
(where \({x}_{i}^{\left(k\right)}\) denotes the value of the k-th feature for the i-th sample, and \({m}_{i}^{\left(k\right)}\) is the mean value of the k-th feature for the i-th sample).
$${f}_{b}=\sum _{i=1}^{c}\frac{{n}_{i}}{n}{({m}_{i}^{\left(k\right)}-{m}^{\left(k\right)})}^{2}$$
6
(where \({n}_{i}\) is the number of samples in the i-th class, n is the total number of samples, and \({m}^{\left(k\right)}\) is the mean value of the k-th feature across all samples).
$$F=\frac{{f}_{b}}{{f}_{a}}$$
7
Higher F values indicate better separability of features and higher classification performance for the system. Among highly similar indicators, the one with the highest F value, representing the optimal classification performance, is extracted while the remaining indicators are excluded. Under the premise of relative independence among indicators, this completes the extraction of optimal foot pressure feature indicators.
1.3.5 Linear Discriminant Analysis
The principle of linear discriminant analysis is based on multivariate linear regression and discriminant analysis. In this study, stepwise selection was used to exclude redundant independent variables and construct an optimal predictive model22,23. At each step, the contribution of an indicator to classification performance was evaluated based on the change in its F value. If p < 0.05, indicating a significant F statistic, the indicator was included in the discriminant function.
Based on the foot pressure classification results, the extracted foot pressure feature indicators were used as independent variables for linear discriminant analysis to obtain the classification function coefficients and constants for each included indicator. Discriminant functions were constructed for each class. The class corresponding to the maximum value calculated from the discriminant functions determined the class membership of a given sample. This process established the discriminant functions for classifying walking foot pressure and identified the final set of feature indicators.
1.4 Statistical Methods
The optimized K-means + + clustering algorithm was implemented using a custom script in MATLAB R2022b. Principal component analysis and linear discriminant analysis were performed using SPSS 25.0. For different foot pressure classification results, independent samples t-tests or Wilcoxon rank-sum tests were used for between-group comparisons if there were two classes, according to the data distribution characteristics. If there were more than two classes, one-way analysis of variance or Kruskal-Wallis H tests were used for between-group comparisons. The significance level was set at 0.05.